The stochastic porous medium equation in one dimension

Introduction. The porous medium equation (PME) $\partial_{t}h=\nabla D(h)\nabla h$ describes the evolution of a conserved scalar field $h$ in a $d$ dimensional non-linear medium Vazquez (2006). As such it appears ubiquitously in physics, e.g. to model heat transfer in presence of non linear thermal conductivity, or density flows of an ideal gas in porous media with non-linear equations of state Barenblatt (1952); Zel’dovich and Kompaneets (1950). The PME equation was proved to arise as an hydrodynamic limit of interacting particle systems Varadhan (1991); Oelschläger (1990), in kinetically constrained lattice gases Gonçalves et al. (2007); Cardoso et al. (2023). It also arises in the conserved dynamics of magnetic interfaces Spohn (1993a, b) and in sandpile models Bántay and Jánosi (1992). More recently it was found to describe superdiffusion of energy in non integrable Luttinger liquids Bulchandani et al. (2020). The PME is often studied for $D(h)\propto|h|^{s-1}$ a power law, the case $s\to 0$ being the logarithmic PME Pedron et al. (2005). For general $s$ it admits special front solutions, known as Barenblatt’s solutions Barenblatt (1952); Zel’dovich and Kompaneets (1950); Vazquez (2006), spreading as $x\sim t^{1/(2+d(s-1))}$ subdiffusively for $s<1$, and superdiffusively for $s>1$.

In this Letter we study a stochastic version of the porous medium equation (SPME) in presence of additive non-conservative noise. We will consider it as a model for the growth of an interface of height $h(x,t)$ with evolution equation

where $\eta(x,t)$ is a Gaussian white noise uncorrelated in time and space, $V(h)$ is an increasing function of $h$, and the Laplacian term can also be written as $\nabla D(h)\nabla h$ with $D(h)=V^{\prime}(h)$, i.e. $V(h)\propto h|h|^{s-1}$ in the power law case. For $s=1$ it reduces to the celebrated Edwards-Wilkinson (EW) equation Edwards and Wilkinson (1982), i.e. the simplest linear equation describing the equilibrium growth of an interface. In that case the height fluctuations scale as $h\sim x^{\alpha}$ with roughness exponent $\alpha=1/2$, while the space-time scaling $t\sim x^{z}$ corresponds to dynamical exponent $z=2$. For $s\neq 1$ the non-linearities make the problem non-equilibrium, as for the Kardar-Parisi-Zhang equation Kardar et al. (1986), although here the universality class is different since the equation is conservative.

An anisotropic version of the SPME in space dimension $d=2$ has been studied before in the context of erosion of tilted landscapes Pastor-Satorras and Rothman (1998a, b); Antonov and Kakin (2017); Birnir and Rowlett (2013); Pelletier (1999) and of continuous models for self-organized criticality Hwa and Kardar (1989). More recently the stochastic PME has also been studied in the context of the Wilson-Cowan equation in critical neural networks models of the brain dynamics Tiberi et al. (2022) as well as in active matter spe . Furthermore the stochastic PME has also been studied in mathematics Barbu et al. (2016) where existence of solutions was obtained both for the additive noise Da Prato and Röckner (2004); Kim (2006) and multiplicative noise Barbu et al. (2009); Dareiotis et al. (2021).

Despite these works, there is no quantitative theory for the SPME in $d=1$, unlike the KPZ equation. The aim of the present work is thus to determine the scaling exponents $\alpha$ and $z$, and to numerically check whether standard scaling applies to this highly non-linear stochastic growth equation nag . Surprisingly, we find that the scaling for $s\neq 1$ is anomalous, suggesting the presence of a second roughness exponent $\alpha_{\text{loc}}$. Although anomalous scaling has been observed in other systems López (1999); Krug (1994); Bernard et al. (2024); Lopez and Schmittbuhl (1998), its origin is fundamentally different in the context of SPME.

In particular, we consider the following discrete version of the SPME in $d=1$ with additive white noise:

Note that in the above equation, the time dependence of $h_{n}$ is kept implicit. This equation represents the time evolution of a spring chain characterized by non-homogeneous and $h_{n}$-dependent elastic constants $D(h_{n})$. Different behavior of $D(h_{n})$ can be considered with the only constraint that $D(h_{n})>0$. In this paper, we study the power-law case, $D(h_{n})\sim|h_{n}|^{s-1}$ for large $h_{n}$ with $s>0$ and $D(h_{n})$ even. At the origin, the behavior is regularized with $D(0)=\nu$ (See SM for the precise form of $D(h)$ used in numerical simulation). In Fig.1, we show two interfaces in the stationary regime. The interface with $s>1$ is flatter compared to the case $s<1$. Indeed, for large $h$, since $D(h)\sim|h|^{s-1}$ the springs become stiff when $s>1$ and soft when $s<1$. One of the main result of this paper is the following prediction for the roughness and dynamical exponents

This prediction verifies the exact exponent relationWolf and Villain (1990)

valid for conservative dynamics with non conserved noise ¹¹1Indeed, integrating the SPME equation (1) over the system volume leads to $h\sim t^{1/2}/L^{d/2}$, and thus to the exponent relation.. First, we perform a functional RG calculation and find stable fixed points that preserve the same power law behavior $D(h)\sim|h|^{s-1}$ at large $h$. This leads to the prediction Eq. (3) for $d\leq 2$, obtained to first order in a $\epsilon=2-d$ expansion, but argued to be exact for $s>0$. However, we also find that local observables exhibit anomalous scaling. For instance, for $\ell\ll L,$ and $n\sim L$, the $q$ moments of the stationary equal time height-to-height correlation read:

where $\langle\dots\rangle$ indicates the average over the noise. In Eq.(5), $\alpha_{\rm loc}(q)$ is commonly referred to as the local roughness exponent López (1999); Krug (1994); Bernard et al. (2024); Lopez and Schmittbuhl (1998) and can depend on the value of $q$. The anomalous scaling arises due to the presence of inhomogeneous increments, $h_{n+1}-h_{n}$, along the interface.

Third, in the final part of the Letter, we show with heuristic arguments and numerical simulations, that the stationary interface of the SPME, provided that $L$ is large, has the same statistical properties as the following random walk:

where $\chi_{n}$ is a centered unit Gaussian random variable, independent of $\tilde{h}_{n}$. Using this mapping, we can recover the roughness exponent from Eq.(3) as well as the anomalous scaling. In particular, we predict that $\alpha_{\text{loc}}(q)=1/2$ for $0<s<1$ while for $s>1$, one has

with $q<q_{c}=2+2/(s-1)$.

Renormalization group. To perform the RG analysis of the SPME (1) we use the Martin-Siggia-Rose dynamical field theory in space dimension $d$, near the upper critical dimension $d=2$, in a perturbative expansion in $\epsilon=2-d$. The RG analysis of this model was first performed in Antonov and Vasil’Ev (1995) within a polynomial expansion of the function $V(h)$ keeping only first few terms. Here we analyze it in the full functional form, that is necessary to identify all possible fixed points which are determined by the large $h$ behavior of $V(h)$. We introduce an infrared renormalization scale $\mu$ (an inverse length scale which can be taken as $L^{-1}$). We denote the bare potential in the dynamical action by $V_{0}(h)$, where $V^{\prime}_{0}(0)=\nu$ is the bare diffusivity, and the $\mu$-dependent renormalized potential at scale $\mu$ by $V(h)$, defined from the effective action. The functional RG flow to one loop (i.e. to first order in $\epsilon>0$) then reads SM

We introduce $u=\frac{1}{4\pi}\int_{\mu}^{\mu_{0}}d\mu^{\prime}/(\mu^{\prime})^{1+\epsilon}$ so that $u\sim\mu^{-\epsilon}/\epsilon$ for $\epsilon>0$ and $u=\frac{1}{4\pi}\log(\mu_{0}/\mu)$ in $d=2$. Let us denote the running $V^{\prime}(h)$ at the scale $\mu$ by $D_{u}(h)$. The RG equation can be put in the form

which, interestingly, is itself a logarithmic PME. We can look for scale invariant fixed point solutions of the RG flow in the form

Once we find such a fixed point solution, since $u\sim\mu^{-\epsilon}\sim L^{\epsilon}$, one can identify $\alpha=a\epsilon$ where $\alpha$ is the anomalous scaling dimension of $h\sim L^{\alpha}$. The corresponding time scale is $t\sim L^{2}/\nu(L)$ where $\nu(L)=D_{u}(0)\sim L^{\epsilon(1-2a)}$ is the renormalized diffusivity, leading to $t\sim L^{z}$ with $z=2-\epsilon(1-2a)$. Hence the exact relation $z=2\alpha+d$ follows and there is a single exponent (as the noise is not corrected under RG) see SM for details. Inserting (10) into (9) we obtain the following fixed point equation for $F(H)$ with $H=h/u^{a}$

Since here we consider only the case where $V^{\prime}(h)$ is finite at $h=0$ and even, we must have $F^{\prime}(0)=0$. Moreover, if $F(H)$ is solution of (11), so is $b^{2}F(Hb)$ for any $b$. We therefore choose $F(0)=1$. We have found through analytical and numerical analysis that there are four types of such solutions to the equation (11):

(i) a line of fixed points with $0<a<1$ and $F(H)\sim_{H\gg 1}H^{-2+1/a}$.

(ii) a special solution $F(H)=1/(1+H^{2}/2)$ associated to $a=1$.

(iii) $F(H)\sim_{H\gg 1}e^{-C_{1}H}/H^{2}$ associated to $a>1$.

(iv) solutions diverging at finite $H=H_{0}$ as $\sim 1/(H_{0}-H)$ for $a<0$ and $\sim 1/(H_{0}-H)^{2}$ for $a=0$ .

We show by linear stability analysis, see SM , that the solutions (i) which are relevant for this paper correspond to attractive fixed points of the RG flow Eq.(9). In order to study in more details the basin of attraction of these fixed points, we solve Eq.(9) with several initial conditions (See Fig.2). For $s>0$, the large $|h|$ behavior of $D_{u}(h)\sim|h|^{s-1}$ is conserved by the RG flow, which converges towards the FP (i) with $a=1/(1+s)$. Finally, one can argue that the higher order loop corrections to the RG equations (8) and (11) cannot change the large $h$ behavior of $V(h)$ for $s>0$, leading to our conjecture that the prediction (3) for the exponents is exact.

Numerical results. We use an Euler scheme to integrate Eq.(2), from a flat initial condition, see SM for details and for our choice of function $D(h)$. We fix $h_{0}=0$ and let the right boundary free ($D(h_{L})=0$). We first compute the time evolution of the width $W(L,t)=\sqrt{\frac{1}{L}\langle\sum_{n=1}^{L}(h_{n}(t)-\bar{h}(t))^{2}\rangle}$, where $\bar{h}(t)=\frac{1}{L}\sum_{n=1}^{L}h_{n}(t)$ denotes the center of mass. Fig.3 shows a very good agreement with Family-Vicsek scaling:

Here the values of the exponents $\alpha$ and $z$ are the values predicted by the RG (Eq.(3)). The scaling function $f_{w}(u)$ depends on the boundary conditions, with $f_{w}(u)\to{\rm cst}$ for $u\gg 1$ and $f_{w}(u)\sim u^{\alpha/z}$ for $u\ll 1$. Indeed, the interface is rough up to a length $\xi(t)\sim t^{1/z}$, and is flat above this scale. As a consequence, at short times ($u\ll 1$), the width $W(L,t)$, as well as the typical height $h_{n}(t)$, grow as $\sim\xi(t)^{\alpha}\sim t^{\beta}=t^{\alpha/z}$.

The width is a global observable which probes the behavior of the entire interface. We have also computed local observables, particularly the $q$-th moments of the stationary height-to-height equal-time correlation. In the inset of Fig. 3, we show, for $s=1/2$, that it exhibits anomalous scaling behavior for $\ell\ll L$, as displayed in Eq.(5), with the exponent $\alpha_{\rm loc}(q)=1/2$. Here, the average is taken over both the noise but also over $n\lesssim L$. Although Eq.(5) holds for the stationary state, the anomalous scaling can be extended to finite times by replacing $L$ by $\xi(t)$. We will argue that $\alpha_{\rm loc}(q)=1/2$ for $s<1$, while for $s>1$ it depends explicitly on $q$ (See Eq.(7)). Finally, we have also studied the structure factor of the height field which also displays the anomalous scaling SM .

Connection to the random walk. To clarify the relation between the random walk (RW) model defined in Eq.(6) and the stationary interfaces of the SPME, we recall that if the spring constants were independent of $h$, i.e. for the EW model, the stationary limit would be described by the equipartition theorem. This theorem predicts Gaussian increments with zero mean and variance $1/D$. This reasoning can be extended to cases where $D(h_{n})$ varies slowly along the $x$-axis, namely $D(h_{n})\simeq D(h_{n+1})$. Using a Taylor expansion, we can rewrite the latter condition as a condition for the increments:

where the asymptotic behavior $\sim|h_{n}|$ holds only for large $h_{n}$. For small increments, the distribution of $h_{n+1}$ conditioned to $h_{n}$ becomes Gaussian of mean $h_{n}$ and variance $1/D(h_{n})$. For large increments, deviations from the Gaussian behavior are observed. However, as can be seen in the inset of Fig.4, these deviations disappear when $|h_{n}|\gg 1$, in which case the condition (13) is satisfied ²²2The condition Eq.(13) for the increments can be rewritten as a condition for $h_{n}$. Indeed, for small increments, we expect $|h_{n+1}-h_{n}|\sim 1/\sqrt{D(h_{n})}$, hence from Eq.(13), the condition becomes $|h_{n+1}-h_{n}|\sim 1/\sqrt{D(h_{n})}\sim|h_{n}|^{(1-s)/2}\ll D(h_{n})/|D^{\prime}(h_{n})|\sim|h_{n}|$. This condition is satisfied for $|h_{n}|\gg 1$.. This is also confirmed in the inset of Fig.1, where we show that $\langle(h_{n+1}-h_{n})^{2}\rangle=1/D(h_{n})$ for large $|h_{n}|$. In the stationary regime, for $n=O(L)$, the height of the interface scale as $h_{n}\sim L^{\alpha}\gg 1$, hence, in the limit of large system size, we may expect that the statistical properties of the SPME are described by the random walk model of Eq.(6) at all scales. This is indeed what we observe, as detailed in SM and illustrated in Fig. 4.

This RW model allows for several predictions that we now detail. Since the $\xi_{n}$ are i.i.d one obtains the second moment of the height difference (removing the tilde symbol)

At large $n$, we assume the scaling $h_{n}\sim n^{\alpha}{\sf h}$, where $\alpha$ is the Hurst exponent and ${\sf h}$ is the rescaled height, a random variable independent of $n$. Using $D(h)\simeq c|h|^{s-1}$ we obtain for $n\gg 1$ (see SM for details)

To determine $\alpha$ we consider $\ell\sim n$, in which case the r.h.s behaves as $n^{\alpha(1-s)+1}$ while the l.h.s behaves as $n^{2\alpha}$. Equating the two gives $\alpha=1/(1+s)$ for the RW model, which coincides with the prediction (3) for the SPME in $d=1$. Furthermore, one can study the local behavior for $\ell\ll n$. From (15) we obtain

which shows that the RW model exhibits anomalous scaling with $\alpha_{\text{loc}}(q=2)=1/2$.

To obtain all the moments of the height difference $\Delta=h_{n+\ell}-h_{n}$ we need the distribution of the rescaled height at large $n$. This distribution can be obtained using the continuum version of the RW model, analyzed in SM . There we show that this continuum model is related to a Bessel process in the variable $Y=|h|^{(s+1)/2}$. It allows to compute the distribution of the height difference, for given $n$ and $\ell\ll n$ SM . Its expression depends on the value of $s$. For $s>1$, it behaves as

with a typical value $\Delta_{\text{typ}}=n^{\alpha-1/2}\ell^{1/2}$ and an upper cutoff $\Delta_{c}=\ell^{\alpha}$. From (17), we can calculate the $q$-th moment, which shows a multiscaling behavior

with $q_{c}=2+\frac{2}{s-1}$. Indeed, for $q<q_{c}$, the average is dominated by the typical values, namely $\sim\Delta_{\rm typ}^{q}$, while for $q>q_{c}$, it is dominated by the upper cutoff $\sim\Delta_{c}^{q+1}P_{n,\ell}(\Delta_{c})$. For $s<1$, the power law decay of $P_{n,\ell}(\Delta)$ for $\Delta_{\rm typ}\ll\Delta\ll\Delta_{c}$ is replaced by a stretched exponential decay. As a consequence the moments do not display multiscaling and behave as

Averaging over $1\leq n\lesssim L$, we obtain $\alpha_{\rm loc}(q)=1/2$ for $s<1$, while, for $s>1$, we recover the multifractal $\alpha_{\rm loc}(q)$ given in (7). The continuum model enables further predictions, including the determination of the one-point distribution of the rescaled height $\sf h$ (see also Fa and Lenzi (2003); giu )

For $0<s<1$, this distribution exhibits an integrable divergence at the origin. Conversely, for $s>1$, the distribution vanishes at the origin and develops a bimodal shape. Furthermore, the two-point height probability, derived in Eq. (S76), is depicted in Fig. 4 and shows excellent agreement with the SPME. This mapping thus provides a detailed characterization of the statistical properties of the stationary interface.

Conclusion. In summary, we have studied the scaling behavior of the stochastic porous medium equation with a diffusion coefficient $D(h)\sim|h|^{s-1}$ for $s>0$, and obtained the critical exponents. We have predicted and verified numerically that it obeys anomalous scaling with non-trivial multiscaling properties for $s>1$. We have found that the stationary interfaces are remarkably well described by a random walk model, which enables many analytical predictions. Its large scale continuum limit is related to a Bessel process, allowing for calculation of, e.g., the stationary height distribution as well as the multifractal exponents. Using functional RG method we have shown that for arbitrary even functions $D(h)$ there exists only four families of fixed points finite at $h=0$. While only one of them is relevant for the scaling behavior of the SPME, the others could also have applications. In particular, our functional RG equation can be used to study the models of brain dynamics, such as considered in Tiberi et al. (2022).

Acknowledgments: We thank Bjorn Birnir for discussions. AAF acknowledges support from the ANR Grant No. ANR-18-CE40-0033 (DIMERS), AR from the ANR Grant No. ANR-23-CE30-0031-04 and PLD from the ANR Grant ANR-23-CE30-0020-01 (EDIPS). This work was performed using HPC resources from GENCI-IDRIS (Grant 2024-AD011015488).

.

In this Supplemental Material, we give the principal details of the calculations described in the main text of the Letter.