Further evidence: The KPZ Equation of Kinetic Interface Roughening: A Variational Perspective

1. Introduction

Originally formulated in 1986 to describe kinetic interface roughening, the nonlinear stochastic partial differential equation known as the Kardar–Parisi–Zhang (KPZ) equation readily became a paradigm of non-equilibrium dynamic critical scaling systems, defining as its canonical member the rather comprehensive KPZ universality class. This foundational work [1] established a profound connection between the dynamics of interface growth and the broader field of non-equilibrium physics. The mathematical formulation of the KPZ equation initially posed significant challenges, as its strong nonlinearity exceeded the scope of existing methods for studying stochastic partial differential equations. A major breakthrough in addressing the “well-definedness” of this equation occurred in 2013 [2]. The resolution of this fundamental mathematical problem was essential for advancing the analytical understanding of complex stochastic systems and extracting deeper physical insights.

Universality is a cornerstone concept in statistical physics, asserting that diverse physical systems and mathematical models can exhibit identical statistical behaviors in their long-time and large-scale limits, irrespective of their specific microscopic details. The KPZ equation describes and connects a broad spectrum of significant stochastic growth-like processes across physics, chemistry, and biology, spanning both classical and quantum systems. This wide applicability means that the KPZ universality class offers a powerful coarse-grained predictive framework, enabling the inference of macroscopic behaviors across diverse systems without requiring detailed knowledge of their microscopic interactions. Researchers can infer macroscopic behaviors and statistical properties across vastly different physical, chemical, and biological systems without needing to delve into their specific microscopic interactions. This “surprising and broad universality” of the one-dimensional KPZ equation (even appearing in deterministic Hamiltonian systems) further elevates its status from a mere model to a fundamental principle that unifies the study of non-equilibrium phenomena, suggesting that common mathematical structures dictate emergent behavior irrespective of the underlying physical substrate.

The study of the KPZ equation has revealed deep and unexpected mathematical connections that have been pivotal in its analytical understanding and in broadening its perceived universality, due to its relation with both the Burgers equation of turbulence and the restricted partition function of directed polymers in random media (DPRM). In particular, the KPZ equation and its associated universality class have remarkably extended their reach into the domain of biological growth phenomena, providing a framework for understanding complex pattern formation and spatial dynamics in living systems. The KPZ equation is being increasingly applied to movement ecology due to its ability to model complex stochastic processes that resemble ecological patterns. Since the scaling behavior of the surface of Escherichia coli and Bacillus subtilis colonies growing on agar plates was established to be self-affine and consistent with models in the so-called “KPZ universality class” [3], there have been numerous applications of the KPZ equation in movement ecology [4,5], among which we might cite the following:

• Modeling nonlinearities in ecological systems. The KPZ equation’s inclusion of nonlinear growth terms is particularly useful for modeling ecological processes where interactions between individuals or between species and their environments lead to complex, emergent movement patterns. For example, the propagation of invasive species, herd dynamics, or predator–prey interactions often exhibit such nonlinearities, and KPZ can capture these effects efficiently. A foundational exploration of movement patterns, emphasizing the need for nonlinear models like the KPZ to understand large-scale spatial patterns, can be found in [6].

• Spatial pattern formation in ecosystems. KPZ-like models can explain spatial clustering and aggregation, phenomena common in ecosystems due to habitat fragmentation or social behaviors in species. In this sense, the KPZ equation aligns well with empirical observations in movement ecology, where spatially heterogeneous environments cause uneven dispersal patterns. Book [7] elaborates on spatial pattern formation and how stochastic models, including those similar to KPZ, apply to understanding the distribution and aggregation of organisms.

• Capturing anomalous diffusion in animal movement. Traditional diffusion models often fail to capture the complexity of animal movement, which frequently exhibits subdiffusive or superdiffusive behavior due to factors like resource hotspots or behavioral adaptations. The KPZ equation, with its growth and roughening terms, provides a more flexible framework for modeling these deviations from classical diffusion. The article [8] reviews the inadequacy of simple diffusion models and how alternative stochastic approaches like KPZ are better suited to describe animal search behaviors that diverge from Brownian assumptions.

• Modeling stochastic processes in animal movement. The KPZ equation is an excellent tool for capturing random yet structured movement patterns, such as those seen in animal migrations or dispersal behaviors. Unlike simple random walks, the KPZ framework can describe anisotropic movement, a key aspect of ecological processes where environmental and physiological conditions lead to non-uniform patterns. The article [9] discusses the role of stochastic models, like KPZ, in modeling animal movement patterns beyond simple Brownian motion.

• Interfacing with statistical mechanics for predictive ecology. KPZ’s roots in statistical mechanics make it a natural fit for ecological models that require predictions based on system dynamics under uncertainty. By using KPZ, researchers can leverage known statistical properties to predict critical transitions, tipping points, or responses to environmental changes. The article [10] discusses the interface between movement ecology and statistical mechanics, highlighting how frameworks like KPZ facilitate a predictive approach to ecological modeling.

• In sum, the KPZ equation offers a robust framework for addressing the complex, nonlinear, and stochastic characteristics of movement patterns in ecology. Its applicability extends beyond traditional diffusion models, making it valuable for studying the spatial and temporal heterogeneity observed in animal populations. By leveraging the KPZ equation, ecologists can create models that more accurately reflect real-world movement dynamics, supporting better predictive and management tools in ecological research [11,12,13,14]. The KPZ universality class is characterized by definite values of either the roughness exponent α or the dynamical exponent z, and by the relation $\alpha +z=2$. If L denotes the substrate’s length along any dimension, then the global width

  $$W(L,t):={\langle \left\{{[h(x,t)-\left\{h\right\}\left(t\right)]}^2\right\}\rangle }^{1/2}$$

  scales as $L^{\alpha }$. Here h is the local height (flat substrates) or radius (curved substrates), and $\{\dots \}$ denote spatial average and $\langle \dots \rangle$ the ensemble one. Moreover, the correlation length $\xi \left(t\right)$ on the interface scales as $t^{1/z}$. For finite L, the scaling regime ends at saturation, when $\xi \sim L$. As $L\to \infty$, it is a dynamical critical phenomenon. Putting it all together, it turns out that in the scaling regime $W(L,t)\sim t^{\beta }$, where $\beta :=\alpha /z$ is called the growth exponent.

Of course, ${\nabla }^{4}h$ becomes the most relevant operator for $\nu =0$. Golubović and Bruinsma [16] studied the effect of this on the lower critical dimension and the exponent $\alpha$, and concluded that it leads to a different universality class than the KPZ one, showing that the ${(\nabla h)}^2$ term does not suffice (either by itself or through its interplay with the noise) to determine the universality class, but it is their interplay with interface relaxation that does.

In Section 2, we review a variational formulation of the KPZ equation and related ones and point out some issues like the necessity of Galilean invariance in discrete versions, the Markov character of process $h\left(x\right)$, and the existence of steady-state distribution (before reaching saturation!) in $d=1$ substrates. In Section 3, as a prolegomenon to Section 4, we recall previous numerical results on the time evolution of the NEP in $d=1$, 2, and 3, as well as its scaling behavior with the nonlinearity parameter $\lambda$. In Section 4 (always from a variational perspective), we show the intimate connection between the NEP and entropy production (both from a path-integral formulation or starting directly from the KPZ equation) and address the memory of process $h\left(x\right)$ with respect to initial conditions, both in KPZ and GB. In Section 5, we highlight some open questions that can be approached from a variational perspective.

2. The Variational Formulation

2.1. Model A—Like Variational Formulation

Although describing a non-equilibrium process, Equation (2) originated in purely kinetic models. For some time, it was believed that its drift $\nu {\nabla }^{2}h(\mathbf{x},t)+{\displaystyle \frac{\lambda }{2}}{[\nabla h(\mathbf{x},t)]}^2$ cannot be derived from an effective free energy [17].

Using the Hopf–Cole transformation

$$\Psi (\mathbf{x},t):=exp\left[{\displaystyle \frac{\lambda }{2\nu }}h(\mathbf{x},t)\right]$$

(mapping h onto a field $\Psi$ obeying a linear equation, but with multiplicative noise), it was proved in [18] that for any d, Equation (2) can be cast as model A (namely, a relaxational system without conservation laws, in the nomenclature of [19])

$${\displaystyle \frac{\partial h(\mathbf{x},t)}{\partial t}}=-\Gamma \left[h\right]{\displaystyle \frac{\delta \mathcal{F}\left[h\right]}{\delta h(\mathbf{x},t)}}+\xi (\mathbf{x},t),$$(4)

with functional $\mathcal{F}\left[h\right]$ given by Equation (5), and functional $\Gamma \left[h\right]$ (known as the motility) by

$$\Gamma \left[h\right]:={\left({\displaystyle \frac{2\nu }{\lambda }}\right)}^{2}exp\left[-{\displaystyle \frac{\lambda }{\nu }}h(\mathbf{x},t)\right].$$

$$\mathcal{F}\left[h\right]:=\int d\mathbf{x}exp\left[{\displaystyle \frac{\lambda }{\nu }}h(\mathbf{x},t)\right]{\displaystyle \frac{\lambda }{2\nu }}\left\{-Fh(\mathbf{x},t)+{\displaystyle \frac{\lambda }{4}}{\left[\nabla h(\mathbf{x},t)\right]}^2\right\}$$(5)

The Lyapunov property $\dot{\mathcal{F}}\left[h\right]\le 0$ is easy to check, and the minimum of $\mathcal{F}\left[h\right]$ is reached by constant functions. In addition, the effective potential $f\left(h\right)$—the restriction of $\mathcal{F}$ to $h(\mathbf{x},t)=\mathrm{const}.$—is minimum at $h=0$ for $F=0$ (Figure 1). Hence, $\mathcal{F}\left[h\right]$ is a Lyapunov functional for the KPZ problem, and the deterministic KPZ equation displays simple dynamics: the asymptotic stability of constant solutions indicates that arbitrary initial conditions approach constant profiles for long times. While the deterministic case is simple, the stochastic dynamics remain highly non-trivial, producing self-affine fractal profiles. Notably, the Lyapunov functional offers little intuition about these stochastic features.

As shown in [18], $\mathcal{F}\left[h\right]$ and $\Gamma \left[h\right]$ (together with the density $\tilde{\mathcal{F}}[h,\nabla h]$, defined by $\mathcal{F}\left[h\right]:=\int \tilde{\mathcal{F}}[h,\nabla h]d\mathbf{x}$) fulfill the relations postulated in [20] on the grounds of global shift invariance.

It is worth noting that several studies [21,22,23,24] have shown that in oscillator systems subjected to noise, synchronization is achieved when the coupling is strong enough. This dynamic process can be related to a roughness process and exhibit, in one dimension, scale invariance within the KPZ universality class. Since Kuramoto-type models can be obtained from energy functionals, this provides further support for the idea of the existence of such a functional for KPZ.

2.2. Discrete Version

Discretization is a necessity if the substrate dimensionality is $d>1$, where the DRGFP is not reachable perturbatively and one must resort to numerical integration to find the critical exponents. In the past, a criterion to judge the quality of a discretization of Equation (2) was that it should preserve the two main symmetries of Equation (2): Galilean invariance (GI) and a fluctuation–dissipation relation peculiar of 1d (1dFDR). Note that although there is an explicit FDR for diffusion phenomena, there is in principle no such thing for growth processes. Recently, however, a kind of FDR—similar to the diffusive one—was found for growth phenomena [25,26]. In turn, Galilean invariance implies that the physics of interface growth remains the same when viewed from a reference frame moving laterally at constant speed, provided the surface height is appropriately tilted. This symmetry protects parameter $\lambda$ from being renormalized and imposes the constraint $\alpha +z=2$.

Consistency in employing the Hopf–Cole transformation constraints the discrete version of ${(\nabla h)}^2$ to have the form ${\displaystyle \frac{1}{2}}[{(h_{i+1}-h_i)}^2+{(h_{i-1}-h_i)}^2]$. In [27], scaling was obtained to a very good approximation (Figure 2) using this discretization, which breaks GI and 1dFDR explicitly.

2.3. Memory of Process $h(\mathbf{x},t)$

From the Family–Vicsék scaling of roughening processes [29], it turns out that

$$h(\mathbf{x},t)\sim v_{\infty }t+A{t}^{\beta }\chi \left({\mathbf{x}}^{\prime }\right),{\mathbf{x}}^{\prime }=B{\displaystyle \frac{\mathbf{x}}{t^{1/z}}}.$$(6)

Over the last 15 years or so, numerous works by mathematicians and physicists [2,30,31,32,33,34,35] have shown that the asymptotic (in the scaling regime) probability distribution $P\left(\chi \right)$ (for any x) not only has skewness, but shows a remarkable dependence on the initial conditions. In other words, process $h\left(x\right)$ has memory.

Equation (2) is an extended Langevin equation; thus, the joint process $\left[h\right]\left(t\right)$ obeys a Fokker–Planck equation and is Markovian. Now $h\left(x\right)$ obeys a marginal distribution and is generically non-Markovian (unless a conservation law, like, e.g., $\int d\mathbf{x}{\nabla }^{2}h=0$ in the EW model, forces it to be Markovian).

Besides the model A variational formulation of Equations (4) and (5), Equation (2) can also be cast (again, for any d) as [18]

$${\displaystyle \frac{\partial h(\mathbf{x},t)}{\partial t}}=-{\displaystyle \frac{\delta \Phi \left[h\right]}{\delta h(\mathbf{x},t)}}+\xi (\mathbf{x},t),$$(7)

with

$$\Phi \left[h\right]:={\displaystyle \frac{\nu }{2}}\int d\mathbf{x}{\left(\nabla h\right)}^2-F\int d\mathbf{x}h(\mathbf{x},t)-{\displaystyle \frac{\lambda }{2}}\int d\mathbf{x}{\int }_{h_{\mathrm{ref}}}^{h(\mathbf{x},t)}\mathcal{D}\psi {\left(\nabla \psi \right)}^2.$$(8)

Again, $\dot{\Phi }\left[h\right]\le 0$ and its effective potential is minimum at $h=0$ for $F=0$ (Figure 3). But at variance with $\mathcal{F}\left[h\right]$, $\Phi \left[h\right]$ is not bounded from below, so it is not a Lyapunov functional for $h(\mathbf{x},t)$. However, $\Phi \left[h\right]$ does provides intuition on the stochastic dynamics (for instance, it resembles an extended escape problem).

In particular, Equation (8) can be transformed into

$$\Phi \left[h\right]={\displaystyle \frac{\nu }{2}}\int d\mathbf{x}{\left(\nabla h\right)}^2-F\int d\mathbf{x}h(\mathbf{x},t)-{\displaystyle \frac{\lambda }{2}}\int d\mathbf{x}{\int }_0^{t}ds\left[{\partial }_{s}h(\mathbf{x},s)\right]{\left[\nabla h(\mathbf{x},s)\right]}^2.$$(9)

which shows explicitly that $\Phi \left[h\right]$ depends on the whole trajectory $\left[h\right]\left(s\right)$, from 0 to t.

As in the case of reaction–diffusion systems, Equation (2) looks formally as a stochastically forced gradient flow (in the sense that the coupling terms have the same form for all sites), but it is not strictly gradient. Equation (7) is regarded as a gradient flow in the sense that its motility is $\Gamma \left[h\right]=1$.

3. KPZ NEP’s Time Behavior and Its $\lambda$-Dependence

We recall our numerical results on the time behavior of the NEP in dimensions $d=1,2,3$ [37]. These results (i) provide the first systematic characterization of the dimensional dependence of NEP, (ii) quantify memory effects in NEP that had not been reported, and (iii) establish explicit connections with stochastic thermodynamics.

Figure 7 and Figure 8 display the mean $\lambda$-dependence over 100 samples of the NEP’s asymptotic slope and of the time of occurrence of the NEP’s maximum, respectively, for $d=1$, 2, and 3. Whereas the first magnitude shows no dependence whatsoever on d, a mild $d$-dependence is seen in the second one.

We feel that, while RG fails for $d>1$, the variational approach could offer a qualitative or complementary insight for obtaining the critical exponents.

4. Novel Results

4.2. Memory of Initial Conditions in KPZ

The dependence of the asymptotic h-statistics on initial conditions is an established fact [30,31,32,33,44]. A numerical test, undertaken with initial conditions of the form

$$\begin{array}{ccc}\hfill h\left(0\right)& =& {\displaystyle \frac{1}{a}}\mathrm{abs}(x-L/2),\hfill \\ \hfill a& \in & \{0.5,1.5,3,20\}.\hfill \end{array}$$

shows that the NEP also displays memory effects.

Within the studied interval, the slopes show a strong dependence on the initial conditions, and in some cases, the curves also display concavity (Figure 9).

4.3. The Related Golubović–Bruinsma Model

As indicated in the Introduction, Golubović and Bruisnma [16] have analyzed an extension of the KPZ model. The equation for the the growth rate they have found is as follows:

$${\displaystyle \frac{\partial h(\mathbf{x},t)}{\partial t}}=\nu {\nabla }^{2}h(\mathbf{x},t)-{\left(K{\nabla }^2\right)}^{2}h(\mathbf{x},t)+{\displaystyle \frac{\lambda }{2}}{\left(\nabla h(\mathbf{x},t)\right)}^2+\xi (\mathbf{x},t).$$(21)

As in the KPZ case, this equation can be derived from a NEP, given by

$${\Phi }_{GB}\left[h\right]:={\displaystyle \frac{\nu }{2}}\int d\mathbf{x}{\left(\nabla h\right)}^2+{\displaystyle \frac{K}{2}}\int d\mathbf{x}{\left({\nabla }^{2}h\right)}^2-{\displaystyle \frac{\lambda }{2}}\int d\mathbf{x}{\int }_0^{t}ds{\partial }_{s}h{\left(\nabla h\right)}^2,$$(22)

such that

$${\displaystyle \frac{\partial h(\mathbf{x},t)}{\partial t}}=-{\displaystyle \frac{\delta {\Phi }_{GB}\left[h\right]}{\delta h(\mathbf{x},t)}}+\xi (\mathbf{x},t).$$(23)

Equation (23) expresses the Golubovic–Bruinsma model as a stochastic gradient flow driven by the functional derivative of the non-equilibrium potential ${\Phi }_{GB}\left[h\right]$. This formulation highlights the shared variational structure between this model and the standard KPZ equation, allowing for a unified analysis of their dynamical and statistical properties within the same theoretical framework.

From this form of the potential, we can immediately predict that, for plane initial conditions, the long time behavior will be dominated by a linear dependence on t, with a slope depending on ${\lambda }^2$. We can also anticipate the existence of memory effects.

In the following figures, we show the time behavior of the indicated functional as well as, changing initial conditions, we prove the existence of memory effects. In all these figures, the parameters we have used are $\nu =\lambda =K=1$.

From Figure 10, it is clear that the linear contributions have a similar behavior than in the EW case, while the nonlinear contribution, dominating at long time, reaches a linear dependence on t. The slope should scale with ${\lambda }^2$, which is an aspect to be discussed in a forthcoming work.

In Figure 11, the memory effects are clearly shown when changing the initial conditions. A detailed study will be the subject of a forthcoming work.

Considering the associated Fokker–Planck equation, and considering the stationary case we observe that, independently of the bc, the exponential

$$P_{BG}^{\mathrm{stat}}\left[h\right]\propto exp\{-{\Phi }_{KPZ}\left[h\right]/\epsilon \},$$

fullfils the stationary condition

$${\displaystyle \frac{\partial }{\partial t}}{P}_{BG}^{\mathrm{stat}}\left[h\right]=0,$$

with ${\Phi }_{BG}^{\mathrm{stat}}\left[h\right]$ given by Equation (22). Clearly, as in the KPZ case, this is a formal solution since $P_{BG}^{\mathrm{stat}}\left[h\right]$ cannot be normalized.

An interesting aspect was studied in [45] and corresponds to the large deviations of the interface height in this model. Such a study follows [46,47] on similar aspects for the KPZ equation.

5. Conclusions and Outlook

Here we highlighted, from a variational approach, several aspects of the KPZ equation as a theory of rough interface growth. We have addressed several misconceptios about KPZ—such as the alleged absence of a free-energy-type functional (the NEP), the supposed role of Galilean invariance in fixing the universality class, the Markovian character of the process, and the conjecture linking the stationary EW pdf to $d=1$ KPZ with periodic boundaries—and dismantled them.

In addition, we presented numerical results of the NEP (from which one can extract information on the $\lambda$ dependence of its long-time behavior and the crossover time from EW to KPZ regime), and have shown numerically that also the NEP displays memory of initial conditions (a fact established in several theoretical works). We have also referred to our work on the stochastic thermodynamic description of KPZ based on a variational approach, allowing us to numerically obtain several known fluctuation theorems and TUR for entropy production.

We have also shown the existence of a NEP in a related model, the Golubovic–Bruinsma one. We have shown the time dependence of such a functional, as well as the existence of memory effects. We expect to extend this approach to other related systems. We also expect to exploit our functional approach to investigate large deviation functions as discussed in [45,46,47].

We are currently working on the possibility of obtaining information about the critical exponents and the existence of an upper critical dimension, using the time behavior of the NEP. While RG fails for $d>1$, the variational approach offers a qualitative or complementary insight, even if exact exponents are still an open challenge.