Leibenson Process: Nonlinear Diffusion Analysis

• The Leibenson process is a doubly nonlinear diffusion equation where the p-Laplacian acts on u^q, generalizing porous medium and p-Laplace models.
• Parameter regimes dictate slow diffusion with finite propagation versus fast diffusion with infinite speed, influencing self-similar behavior.
• The equation supports rigorous weak solutions on Riemannian manifolds and connects to nonlinear Markov processes via a McKean–Vlasov framework.

The Leibenson process is the doubly nonlinear diffusion governed by

$$\partial_t u=\Delta_p(u^q), \qquad \Delta_p f=\operatorname{div}\!\bigl(|\nabla f|^{p-2}\nabla f\bigr),$$

posed on Euclidean space or, more generally, on a Riemannian manifold. In the PDE literature this evolution is also called the Leibenson equation or a doubly nonlinear evolution equation; in recent probabilistic work, the same term is used more specifically for the nonlinear Markov process whose one-dimensional marginals are Barenblatt solutions of the equation (Sürig, 28 Jan 2026, Barbu et al., 18 Aug 2025). A closely related reaction–diffusion variant,

$u_t=\Delta_p(u^m)+u^q,$

is studied as a Leibenson-type equation on non-compact Riemannian manifolds and preserves the same doubly nonlinear diffusion mechanism while adding a source term (Meglioli et al., 13 May 2025).

1. Definition, nomenclature, and model class

The defining feature of the Leibenson equation is that the $p$-Laplacian acts on a nonlinear constitutive variable $u^q$, rather than on $u$ itself. On a Riemannian manifold $(M,g)$, the gradient and divergence are taken with respect to the metric $g$, and integrals are taken with respect to the Riemannian measure. In this form, the equation simultaneously generalizes the porous medium equation and the evolutionary $p$-Laplace equation (Sürig, 28 Jan 2026).

Several classical reductions are immediate. When $p=2$, one obtains

$\partial_t u=\Delta(u^q),$

which is the porous medium equation for $u_t=\Delta_p(u^m)+u^q,$0 and the fast diffusion equation for $u_t=\Delta_p(u^m)+u^q,$1. When $u_t=\Delta_p(u^m)+u^q,$2, one recovers the $u_t=\Delta_p(u^m)+u^q,$3-Laplacian evolution

$u_t=\Delta_p(u^m)+u^q,$4

In the mixed reaction–diffusion setting

$u_t=\Delta_p(u^m)+u^q,$5

the diffusion acts on the composite variable $u_t=\Delta_p(u^m)+u^q,$6, not on $u_t=\Delta_p(u^m)+u^q,$7; this distinguishes it both from $u_t=\Delta_p(u^m)+u^q,$8 and from $u_t=\Delta_p(u^m)+u^q,$9. That diffusion term is identified as the Leibenson equation modeling filtration of turbulent compressible fluids in a porous medium, while the term $p$0 is a source term (Meglioli et al., 13 May 2025).

The terminology “Leibenson process” therefore has two related uses. In deterministic nonlinear diffusion it denotes the evolution generated by $p$1; in the McKean–Vlasov framework it denotes the nonlinear Markov process canonically associated with the Barenblatt solution flow of the same equation (Barbu et al., 18 Aug 2025).

2. Parameter regimes and qualitative behavior

The principal qualitative discriminator is

$p$2

When $p$3, equivalently $p$4, the equation lies in the slow-diffusion regime. In this range Barenblatt solutions have finite propagation speed and compact support, and the support radius has the form

$p$5

with

$p$6

When $p$7, equivalently $p$8, the equation is in the fast-diffusion regime; on manifolds this regime exhibits infinite speed of propagation and algebraic spatial decay. The borderline case $p$9 separates these behaviors (Barbu et al., 18 Aug 2025, Grigor'yan et al., 29 Mar 2026).

The same slow/fast dichotomy appears in the generalized diffusion $u^q$0, where the sign of $u^q$1 distinguishes slow diffusion, fast diffusion, and the pseudo-linear threshold $u^q$2. In Euclidean space, the slow-diffusion regime is associated with finite speed of propagation and Barenblatt self-similarity, whereas the fast-diffusion regime has infinite speed (Meglioli et al., 13 May 2025).

For long-distance behavior in the fast regime on complete manifolds satisfying a relative Faber–Krahn inequality and a lower volume bound, bounded non-negative weak subsolutions satisfy the sharp estimate

$u^q$3

where $u^q$4. This matches the spatial and temporal decay of Barenblatt-type self-similar solutions in the Euclidean case and on certain model manifolds (Grigor'yan et al., 29 Mar 2026).

3. Weak solutions, comparison principles, and existence theory

On an arbitrary geodesically complete Riemannian manifold, a central existence theorem states that if $u^q$5, $u^q$6, and $u^q$7, then for every nonnegative initial datum

$u^q$8

the Cauchy problem for

$u^q$9

has a nonnegative bounded weak solution satisfying

$u$0

In distributional form, the weak solution identity is

$u$1

for all $u$2 (Sürig, 28 Jan 2026).

The same theory establishes comparison and uniqueness through an $u$3-contraction of the positive part:

$u$4

for subsolutions $u$5 and supersolutions $u$6 of the auxiliary problems. In particular, two solutions with the same initial data coincide almost everywhere. Nonnegativity is preserved, the solutions are globally bounded in time, and for every $u$7 the map

$u$8

is monotone decreasing; the manifold theory explicitly notes that there is no mass conservation in general (Sürig, 28 Jan 2026).

The construction proceeds by solving truncated auxiliary Dirichlet problems on precompact balls, using a finite-dimensional Galerkin scheme for the flux

$u$9

establishing uniform bounds, identifying the nonlinear limit via monotonicity of the $(M,g)$0-Laplacian, and then exhausting the manifold by balls (Sürig, 28 Jan 2026).

4. Manifold geometry and reaction-driven global existence

For the reaction–diffusion equation

$(M,g)$1

on a complete, non-compact Riemannian manifold $(M,g)$2 of infinite volume and dimension $(M,g)$3, the geometry enters through global functional inequalities. The Sobolev inequality is

$(M,g)$4

and the Poincaré inequality is

$(M,g)$5

These hold, for example, on Cartan–Hadamard manifolds, with the Poincaré inequality available under curvature bounded above by a negative constant (Meglioli et al., 13 May 2025).

Under the baseline assumptions

$(M,g)$6

global existence for small data depends sharply on the geometry. Under Sobolev alone, the threshold is

$(M,g)$7

and this threshold is stated to be sharp in the class of manifolds supporting only Sobolev. Under the additional Poincaré inequality, the global result extends to the full supercritical interval

$(M,g)$8

In the Sobolev-only case, if the data are sufficiently small in $(M,g)$9 and $g$0 for some $g$1, where

$g$2

then the solution exists globally and satisfies the explicit smoothing estimate

$g$3

With Sobolev and Poincaré, global existence continues to hold for sufficiently small data in the larger range $g$4, together with local $g$5 control on balls (Meglioli et al., 13 May 2025).

The geometric significance is explicit: the paper states that this larger global-existence interval has no Euclidean counterpart on bounded domains. The argument is that non-compactness, infinite measure, and global geometric functional inequalities alter the balance between diffusion and reaction. Negative curvature, encoded analytically through $g$6, enhances dissipative control and can suppress Fujita-type blow-up mechanisms that dominate bounded-domain theory (Meglioli et al., 13 May 2025).

5. Gradient bounds, extinction, and propagation on manifolds

For positive solutions of the Leibenson equation on a geodesically complete manifold with Ricci curvature bounded from below by a non-positive constant, Li–Yau type gradient estimates are formulated through pressure variables. In the slow-diffusion case $g$7, the pressure is

$g$8

In the fast-diffusion case $g$9 with the additional condition $p$0, the pressure is

$p$1

Under the structural bounds

$p$2

the theory yields local and global estimates for these functionals, with curvature contributions controlled by time-dependent barrier functions and local Sobolev inequalities on geodesic balls (Sürig, 8 Jun 2025).

In the fast-diffusion regime for the weighted equation

$p$3

finite extinction occurs under a weighted Sobolev inequality and the integrability condition

$p$4

If $p$5 and

$p$6

with

$p$7

then every bounded nonnegative weak subsolution extinguishes in finite time. Writing

$p$8

the proof derives

$p$9

and hence

$p=2$0

The geometric input is entirely through the weighted Sobolev inequality (Sürig, 2024).

Propagation properties also depend sharply on the regime. On arbitrary geodesically complete manifolds, if $p=2$1 and the initial datum vanishes in a geodesic ball $p=2$2, then the solution remains zero in $p=2$3, where

$p=2$4

Thus the front propagates with finite speed, and the rate depends on the relative Faber–Krahn geometry of the ball (Sürig, 28 Jan 2026).

6. Nonlinear Fokker–Planck structure and the McKean Leibenson process

A distinctive recent development identifies the Leibenson equation as a nonlinear Fokker–Planck equation. The exact reformulation is

$p=2$5

with

$p=2$6

and

$p=2$7

The associated McKean–Vlasov SDE is

$p=2$8

with marginal constraint

$p=2$9

The coefficients therefore depend pointwise on the current marginal density and its first derivatives, while the drift involves second derivatives through the spatial derivative of the nonlinear diffusion coefficient (Barbu et al., 18 Aug 2025).

Within the Barenblatt class, this yields a nonlinear Markov process in the sense of McKean. Under

$\partial_t u=\Delta(u^q),$0

and, if $\partial_t u=\Delta(u^q),$1, the additional condition

$\partial_t u=\Delta(u^q),$2

the family of path laws with one-dimensional marginals given by Barenblatt densities forms a unique nonlinear Markov process. The time-homogeneous family $\partial_t u=\Delta(u^q),$3 corresponding to $\partial_t u=\Delta(u^q),$4 is called the Leibenson process (Barbu et al., 18 Aug 2025).

A further theorem proves strong well-posedness from every strictly positive time. If

$\partial_t u=\Delta(u^q),$5

then the weak solution with marginals $\partial_t u=\Delta(u^q),$6 is a probabilistically strong solution, and pathwise uniqueness holds on every finite interval within the class of solutions with those fixed marginals. The proof combines a restricted Yamada–Watanabe theorem, Crippa–De Lellis type maximal-function estimates for BV vector fields, and Muckenhoupt $\partial_t u=\Delta(u^q),$7 weighted estimates. The same work states that this is the first McKean–Vlasov model whose coefficients depend pointwise on time marginals and their first and second spatial derivatives, and it leaves extension from explicit Barenblatt flows to general initial data as an open problem (Barbu et al., 18 Aug 2025).