---
title: 'Leibenson Process: Nonlinear Diffusion Analysis'
url: https://www.emergentmind.com/topics/leibenson-process
type: topic
---

# Leibenson Process: Nonlinear Diffusion Analysis

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 [2601.20640] [2508.12979]. 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 [2505.08304].

## 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 [2601.20640].

Several classical reductions are immediate. When $p=2$, one obtains
$$
\partial_t u=\Delta(u^q),
$$
which is the porous medium equation for $q>1$ and the fast diffusion equation for $0<q<1$. When $q=1$, one recovers the $p$-Laplacian evolution
$$
\partial_t u=\Delta_p u.
$$
In the mixed reaction–diffusion setting
$$
u_t=\Delta_p(u^m)+u^q,
$$
the diffusion acts on the composite variable $u^m$, not on $u$; this distinguishes it both from $\Delta_p u$ and from $\Delta(u^m)$. That diffusion term is identified as the Leibenson equation modeling filtration of turbulent compressible fluids in a porous medium, while the term $u^q$ is a source term [2505.08304].

The terminology “Leibenson process” therefore has two related uses. In deterministic nonlinear diffusion it denotes the evolution generated by $\partial_t u=\Delta_p(u^q)$; in the McKean–Vlasov framework it denotes the nonlinear Markov process canonically associated with the Barenblatt solution flow of the same equation [2508.12979].

## 2. Parameter regimes and qualitative behavior

The principal qualitative discriminator is
$$
D:=1-q(p-1).
$$
When $q(p-1)>1$, equivalently $D<0$, 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
$$
R(t)=\Big(\frac{C}{\kappa}\Big)^{\frac{p-1}{p}} t^{1/\beta},
\qquad
\beta=p+d\bigl(q(p-1)-1\bigr),
$$
with
$$
\gamma=\frac{p-1}{q(p-1)-1},
\qquad
\kappa=\frac{q(p-1)-1}{pq}\,\beta^{-\frac{1}{p-1}}.
$$
When $q(p-1)<1$, equivalently $D>0$, the equation is in the fast-diffusion regime; on manifolds this regime exhibits infinite speed of propagation and algebraic spatial decay. The borderline case $q(p-1)=1$ separates these behaviors [2508.12979] [2603.27791].

The same slow/fast dichotomy appears in the generalized diffusion $\partial_t u=\Delta_p(u^m)$, where the sign of $m(p-1)-1$ distinguishes slow diffusion, fast diffusion, and the pseudo-linear threshold $m(p-1)=1$. 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 [2505.08304].

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(x,t)\le C\,t^{-1/\beta}\Bigl(1+\frac{|x|}{t^{1/\beta}}\Bigr)^{-p/D},
\qquad
\beta=(p-1)(1+q),
$$
where $|x|=d(x,\operatorname{supp}u_0)$. This matches the spatial and temporal decay of Barenblatt-type self-similar solutions in the Euclidean case and on certain model manifolds [2603.27791].

## 3. Weak solutions, comparison principles, and existence theory

On an arbitrary geodesically complete Riemannian manifold, a central existence theorem states that if $p>1$, $q>0$, and $pq\ge1$, then for every nonnegative initial datum
$$
u_0\in L^1(M)\cap L^\infty(M),
$$
the Cauchy problem for
$$
\partial_t u=\Delta_p(u^q)
$$
has a nonnegative bounded weak solution satisfying
$$
u\in C([0,T];L^{q+1}(M)),
\qquad
u^q\in L^p((0,T);W^{1,p}(M))
\quad\forall T<\infty.
$$
In distributional form, the weak solution identity is
$$
-\int_0^\infty\!\!\int_M u\,\partial_t\varphi\,d\mu\,dt
+\int_0^\infty\!\!\int_M |\nabla(u^q)|^{p-2}\nabla(u^q)\cdot\nabla\varphi\,d\mu\,dt
=\int_M u_0(x)\varphi(x,0)\,d\mu
$$
for all $\varphi\in C_c^\infty(M\times(0,\infty))$ [2601.20640].

The same theory establishes comparison and uniqueness through an $L^2$-contraction of the positive part:
$$
\int_{\Omega}(v-u)_+^2(\cdot,t)\,d\mu
\le
\int_{\Omega}(v_0-u_0)_+^2\,d\mu,
$$
for subsolutions $v$ and supersolutions $u$ 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 $\ell\ge1$ the map
$$
t\mapsto \|u(\cdot,t)\|_{L^\ell(\Omega)}
$$
is monotone decreasing; the manifold theory explicitly notes that there is no mass conservation in general [2601.20640].

The construction proceeds by solving truncated auxiliary Dirichlet problems on precompact balls, using a finite-dimensional Galerkin scheme for the flux
$$
A(u,\nabla u)
=
q^{p-1}(u^N)^{(q-1)(p-1)}|\nabla u|^{p-2}\nabla u,
$$
establishing uniform bounds, identifying the nonlinear limit via monotonicity of the $p$-Laplacian, and then exhausting the manifold by balls [2601.20640].

## 4. Manifold geometry and reaction-driven global existence

For the reaction–diffusion equation
$$
u_t=\Delta_p(u^m)+u^q
$$
on a complete, non-compact Riemannian manifold $(M,g)$ of infinite volume and dimension $N>1$, the geometry enters through global functional inequalities. The Sobolev inequality is
$$
\|v\|_{L^{p^*}(M)}\le \frac1{\mathcal S_p}\|\nabla v\|_{L^p(M)},
\qquad
p^*=\frac{pN}{N-p},
$$
and the Poincaré inequality is
$$
\|v\|_{L^p(M)}\le \frac1{\mathcal C_p}\|\nabla v\|_{L^p(M)}.
$$
These hold, for example, on Cartan–Hadamard manifolds, with the Poincaré inequality available under curvature bounded above by a negative constant [2505.08304].

Under the baseline assumptions
$$
1<p<N,\qquad m>1,\qquad m(p-1)\ge1,\qquad q>m(p-1),
$$
global existence for small data depends sharply on the geometry. Under Sobolev alone, the threshold is
$$
q>m(p-1)+\frac pN,
$$
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
$$
q>m(p-1).
$$
In the Sobolev-only case, if the data are sufficiently small in $L^1(M)$ and $L^s(M)$ for some $s>q_0$, where
$$
q_0:=\frac{[q-m(p-1)]N}{p},
$$
then the solution exists globally and satisfies the explicit smoothing estimate
$$
\|u(t)\|_{L^\infty(M)}
\le
c\,t^{-N/(N[m(p-1)-1]+p)}\,
\|u_0\|_{L^1(M)}^{\,p/(N[m(p-1)-1]+p)}.
$$
With Sobolev and Poincaré, global existence continues to hold for sufficiently small data in the larger range $q>m(p-1)$, together with local $L^\infty$ control on balls [2505.08304].

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 $\mathcal C_p>0$, enhances dissipative control and can suppress Fujita-type blow-up mechanisms that dominate bounded-domain theory [2505.08304].

## 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 $\delta:=q(p-1)-1>0$, the pressure is
$$
v=\frac{q(p-1)}{\delta}\,u^{\delta/(p-1)},
\qquad
F=\frac{|\nabla v|^p}{v}-\frac{\partial_t v}{v}.
$$
In the fast-diffusion case $D:=1-q(p-1)>0$ with the additional condition $p-nD>0$, the pressure is
$$
v=\frac{q(p-1)}{D}\,u^{-D/(p-1)},
\qquad
F=\frac{|\nabla v|^p}{v}+\frac{\partial_t v}{v}.
$$
Under the structural bounds
$$
\Lambda_{\min}\le |\nabla v|^{p-2}v\le \Lambda_{\max},
$$
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 [2506.07221].

In the fast-diffusion regime for the weighted equation
$$
\rho\,\partial_t u=\Delta_p(u^q),
$$
finite extinction occurs under a weighted Sobolev inequality and the integrability condition
$$
\Bigl\|\frac{\omega}{\rho}\Bigr\|_{L^\kappa(M,d\mu)}<\infty.
$$
If $D=1-q(p-1)>0$ and
$$
u_0\in L^\sigma(M,d\mu)\cap L^\infty(M),
$$
with
$$
\sigma\ge \max\left\{p,\ pq,\ -D,\ \frac{D}{\kappa-1}\right\},
$$
then every bounded nonnegative weak subsolution extinguishes in finite time. Writing
$$
\Phi(t)=\int_M u^{\sigma+D}(\cdot,t)\,d\mu,
$$
the proof derives
$$
\frac{d}{dt}\Phi(t)\le -c\,\Phi(t)^\beta,
\qquad
\beta>1,
$$
and hence
$$
T_*\le \frac{\sigma+D}{cD}
\left(\int_M u_0^{\sigma+D}\,d\mu\right)^{D/(\sigma+D)}.
$$
The geometric input is entirely through the weighted Sobolev inequality [2412.06496].

Propagation properties also depend sharply on the regime. On arbitrary geodesically complete manifolds, if $q(p-1)>1$ and the initial datum vanishes in a geodesic ball $B_0=B(x_0,R)$, then the solution remains zero in $\tfrac12B_0\times[0,t_0]$, where
$$
t_0
=
c_*\,
\bigl(\mathcal F(B_0)\bigr)^{\frac{q(p-1)-1}{\sigma}}
R^p
\|u_0\|_{L^\sigma(M)}^{-[q(p-1)-1]}.
$$
Thus the front propagates with finite speed, and the rate depends on the relative Faber–Krahn geometry of the ball [2601.20640].

## 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
$$
\partial_t u
=
\partial_{ij}\bigl(a_{ij}(u,x)\,u\bigr)
-
\partial_i\bigl(b_i(u,x)\,u\bigr),
$$
with
$$
a_{ij}(u,x)
=
q^{p-1}\delta_{ij}\,|\nabla u(x)|^{p-2}\,u(x)^{(p-1)(q-1)},
$$
and
$$
b_i(u,x)
=
q^{p-1}\,\partial_i\!\bigl(|\nabla u(x)|^{p-2}\,u(x)^{(p-1)(q-1)}\bigr).
$$
The associated McKean–Vlasov SDE is
$$
\begin{aligned}
dX(t)
&=
q^{p-1}\nabla\!\Bigl(|\nabla u(t,X(t))|^{p-2}u(t,X(t))^{(p-1)(q-1)}\Bigr)\,dt \\
&\quad+
\sqrt{2q^{p-1}|\nabla u(t,X(t))|^{p-2}u(t,X(t))^{(p-1)(q-1)}}\,dW(t),
\end{aligned}
$$
with marginal constraint
$$
\mathcal L_{X(t)}(dx)=u(t,x)\,dx.
$$
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 [2508.12979].

Within the Barenblatt class, this yields a nonlinear Markov process in the sense of McKean. Under
$$
d\ge2,\qquad p>\frac{d}{d-1},\qquad q>\frac1{p-1},
$$
and, if $p<2$, the additional condition
$$
q>\frac{2-p+d}{d(p-1)},
$$
the family of path laws with one-dimensional marginals given by Barenblatt densities forms a unique nonlinear Markov process. The time-homogeneous family $(P_y)_{y\in\mathbb R^d}$ corresponding to $\zeta=\delta_y$ is called the Leibenson process [2508.12979].

A further theorem proves strong well-posedness from every strictly positive time. If
$$
d\ge2,\qquad p>\frac{d-1}{d},\qquad q>\frac{|p-2|+d}{d(p-1)},
$$
then the weak solution with marginals $w^y(t+\delta,\cdot)\,dx$ 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 $A_2$ 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 [2508.12979].

Source: https://www.emergentmind.com/topics/leibenson-process