Transition between chaotic and stochastic universality classes of kinetic roughening

The dynamics of non-equilibrium spatially extended systems are often dominated by fluctuations, due to e.g.\ deterministic chaos or to intrinsic stochasticity. This reflects into generic scale invariant or kinetic roughening behavior that can be classified into universality classes defined by critical exponent values and by the probability distribution function (PDF) of field fluctuations. Suitable geometrical constraints are known to change secondary features of the PDF while keeping the values of the exponents unchanged, inducing universality subclasses. Working on the Kuramoto-Sivashinsky equation as a paradigm of spatiotemporal chaos, we show that the physical nature of the prevailing fluctuations (chaotic or stochastic) can also change the universality class while respecting the exponent values, as the PDF is substantially altered. This transition takes place at a non-zero value of the stochastic noise amplitude and may be suitable for experimental verification.

The dynamics of non-equilibrium spatially extended systems are often dominated by fluctuations, due to e.g. deterministic chaos or to intrinsic stochasticity. This reflects into generic scale invariant or kinetic roughening behavior that can be classified into universality classes defined by critical exponent values and by the probability distribution function (PDF) of field fluctuations. Suitable geometrical constraints are known to change secondary features of the PDF while keeping the values of the exponents unchanged, inducing universality subclasses. Working on the Kuramoto-Sivashinsky equation as a paradigm of spatiotemporal chaos, we show that the physical nature of the prevailing fluctuations (chaotic or stochastic) can also change the universality class while respecting the exponent values, as the PDF is substantially altered. This transition takes place at a non-zero value of the stochastic noise amplitude and may be suitable for experimental verification.
Generic scale invariance (GSI) describes the behavior of many spatially extended non-equilibrium systems in which driving and dissipation act at comparable rates, such that strong correlations build up whose space-time behavior lacks characteristic scales [1]. Hence, they are analogous to equilibrium critical systems, a crucial difference being that parameter tuning is not required for criticality in GSI [2]. Thus, GSI is one of the forms in which critical dynamics [3] is being recently generalized to novel non-equilibrium contexts as assessed in physical and non-physical systems, from quantum matter [4,5], to living [6], or social [7] systems.
A key player in recent advances on the understanding of driven systems displaying GSI [8][9][10][11][12] is the celebrated Kardar-Parisi-Zhang (KPZ) equation for the height h(x, t) of an interface at a substrate position x ∈ R d and time t [13], subject to fluctuations, namely (we henceforth set d = 1), where ν, D > 0, and λ are parameters, and η(x, t) is zeroaverage, Gaussian noise. Within the physical image of evolving interfaces, the GSI behavior displayed by the KPZ equation and related stochastic systems is termed kinetic roughening [8,9]. In analogy with equilibrium critical dynamics, it can be classified into universality classes, which are determined by the values of critical exponents and by the probability distribution function (PDF) of, say, height fluctuations [10][11][12]. Indeed, the 1D KPZ universality class is recently being identified in the fluctuation dynamics of a wide range of low-dimensional, strongly-correlated systems, from thin-film growth, random polymers, or randomly-stirred fluids [8,9], to active matter [14], quantum entanglement [15], or spatiotemporal chaos [16], to cite a few. Based on exact solutions of the one-dimensional (1D) KPZ equation, Eq. (1), and related discrete models (see [10][11][12] for reviews), a particularly rich structure is being elucidated * Electronic address: enrodrig@math.uc3m.es † Electronic address: cuerno@math.uc3m.es for GSI universality classes. To begin with, critical exponent values are known not to unambiguously identify a given class. Indeed, explicit examples have been reported of linear height equations with time- [17] or space- [18] correlated noise, which feature the 1D KPZ scaling exponents, but which cannot be in this universality class, their field PDF being Gaussian, while for the nonlinear 1D KPZ equation it is a member of the Tracy-Widom (TW) PDF family [10][11][12].
Moreover, while keeping the same values of the scaling exponents, universality subclasses exist of the 1D KPZ class, which differ by the flavor of the precise TW PDF which occurs: e.g. for globally flat (curved) interfaces growing from a straight line (point), it corresponds to the largest-eigenvalue TW distribution of random matrices in the Gaussian orthogonal (unitary) ensemble [GOE (GUE)] [10][11][12]. Equivalently, finite systems whose size decreases (increases) linearly with time display TW-GOE (GUE) statistics [19][20][21], while analogous transitions have been assessed for changes in the background topology [22] or in the rate of system-size change [23]. Furthermore, the existence of universality subclasses induced by similar changes in geometrical constraints carries over to the main linear [24] and nonlinear [25] universality classes of kinetic roughening other than KPZ, making this a robust trait of (this type of) criticality far from equilibrium.
However, in all these cases the universality subclasses sharing scaling exponent values differ by some global geometrical or topological condition on the system size or background metric, the existence of alternative mechanisms which likewise control the field PDF remaining uncertain. In this Letter, we demonstrate the physical nature of the prevailing system fluctuations as one such mechanism. Specifically, by considering the 1D Kuramoto-Sivashinsky (KS) equation [26,27] for a scalar field u(x, t), which reads where ν 0 , κ 0 , D 0 ,D 0 > 0, and λ 0 are parameters, and η is as in Eq.
(2), we show (for D 0 = 0,D 0 = 0) that the kinetic roughening behavior which it displays is characterized by critical exponent values which are noise-independent, while the field PDF is non-Gaussian (Gaussian) for low (large) noise values, corresponding to dynamics dominated by chaotic (stochastic) fluctuations. This transition in the universality class occurs at a non-zero noise amplitudeD 0 and might be observable in suitable experimental contexts. The deterministic (D 0 =D 0 = 0) KS equation is a paradigm of spatiotemporal chaos (SC) [28,29], where it has become a benchmark to assess novel concepts and tools, like e.g. reservoir computing [30,31]. Either Eq. (3) proper or the equally ubiquitous version of the deterministic KS equation [26,27], satisfied by h(x, t) = x 0 u(y, t)dy [see Eq. (5) below], both provide physical models in many different contexts, from liquid flow down inclines [32,33] to solidification [34]. The stochastic equation for h has been derived in e.g. epitaxial growth [35,36], ion-beam sputtering [37,38], or diffusionlimited growth [39]; in Appendix A we derive Eq. (3) for a falling liquid film under thermal fluctuations [40,41].
The large-scale behavior of the deterministic KS equation, Eq. (3), is known to remarkably coincide [42,43] with that of the stochastic Burgers equation, whose GSI exponents [44] and field PDF [18,45] are known. Likewise, the deterministic [16,46], as well as the stochastic KS equations for h [47,48] are both in the KPZ universality class. However, how and if the nature of the fluctuations, whether deterministic chaos or stochastic noise, reflects into the GSI behavior, has remained overlooked thus far.
We begin by investigating in full detail the universality class of Eq. (3) for the deterministic case, as well as for the stochastic cases with conserved (D 0 = 0,D 0 = 0) and nonconserved (D 0 = 0,D 0 = 0) noise. Note that, as seminally argued for by Yakhot [42], Eq. (3) is expected to renormalize at large scales into an effective stochastic Burgers equation, where notably ν > 0, rendering asymptotically irrelevant the biharmonic term in Eq. We have performed numerical simulations of Eq. (3) using the pseudospectral method in [49] and periodic boundary conditions for system size L = 2048. Initial conditions are random (with 10 −5 amplitude) for the deterministic case and zero otherwise. Parameters are fixed to ν 0 = κ 0 = 1, λ 0 = 10, and space-time discretization steps δx = 1, and δt ∈ [0.01, 0.05]. Under GSI conditions [8,9], the scaling exponents characterizing the universality class can be readily identified in the evolution of the field roughness W (rms deviation of the field fluctuations), which increases with time as W ∼ t β , reaching a saturation value W sat ∼ L α at time t sat ∼ L z , where z = αβ is the dynamic exponent [3] and α is the roughness exponent, related with the fractal dimension of the u(x) profile [8]. The same exponents occur in the two-point statistics, e.g. in the structure factor S(k, t) = |û(k, t)| 2 , where hat denotes space Fourier transform and k is wavenumber (twopoint correlations in real space are similarly assessed in Appendix A). Indeed [8,9], S(k, t) ∼ 1/k 2α+d for t ≫ L z , while S(k, t) ∼ t (2α+d)/z for k ≪ 1 (with d = 1 here). The numerical time evolution of S(k, t) for the various noise conditions is shown in Fig. 1, where the scaling exponents pre- dicted by Yakhot's argument [42] are indeed obtained: both in the deterministic and in the conserved-noise cases, Eq. (3) renormalizes into Burgers equation, Eq. (4), with conserved noise D = 0,D = 0, for which α cn = −1/2 and z cn = 3/2 [18]. While S(k, t) approaches the (k-independent) whitenoise behavior for increasing t in these two cases, the detailed form of the structure factor curves differs noticeably. In the deterministic case, the scaling exponents are obtained via the integrated h field, trivially expected to scale as h ∼ x α+1 [13,18]. In contrast, for non-conserved noise, Eq. (3) now renormalizes into a noisy Burgers equation, Eq. (4), with nonconserved noise D = 0,D = 0, for which α nc = 0 and z nc = 1 [45].
The scaling exponent values already imply that the KS equation with non-conserved noise belongs to a different universality class than the deterministic and conserved-noise equations, which share the same values of α and z. Hence, for the remainder of this work we set D 0 = 0 and focus on the latter two cases. Still, as noted above, the values of the two independent scaling exponent are currently known [17,18] not to necessarily fix the GSI universality class unambiguously. Moreover, the 1D KPZ universality class also illustrates the fact [10][11][12] that the field PDF can differ in the nonlinear regime prior to saturation to a steady state, as compared with the PDF of fluctuations around such a steady state after it has been reached. Hence, we next investigate the fluctuation statistics for the u field in Eq. (3) within the nonlin- ear regime prior to saturation, setting D 0 = 0 and considering different values of the conserved-noise amplitudeD 0 ; results are shown in Fig. 2. In the deterministicD 0 = 0 case, the rescaled fluctuations of u around its space averageū, defined as χ = (u −ū)/std(u), exhibit a symmetric probability density function (PDF) whose tails decay much faster than those of a Gaussian distribution. Hence, the kurtosis is much smaller than 3, see While the PDFs of the deterministic and the (large) conserved-noise cases of Eq. (3) are both even in χ, they are obviously different, specially with respect to the occurrence of "typical" fluctuation values. One could speak of two different subclasses of a single universality class which additionally features α cn = −1/2 (white noise) and z cn = 3/2 (superdiffusive spread of correlations, as in the 1D KPZ equa- tion). However, further dynamical properties suggest that we rather speak of a change in the universality class, with a transition at a well-defined value ofD 0 separating the predominance of chaotic or of stochastic fluctuations. This interpretation is supported by a study of the behavior of the finitesize Lyapunov exponents of the system [50,51] as a function of the conserved-noise amplitude. Specifically, we measure the so-called scale-dependent Lyapunov exponent Λ(ǫ) [51,52], where ǫ is a distance between trajectories extracted from the time series u(x 0 , t) where x 0 is any fixed value of x and ∆t is a sampling time. These trajectories are efficient tools to reconstruct the attractor of the dynamics [50]. Once all (i, j) pairs of close-by trajectories are selected for a Results are averaged for different values of x 0 and, after binning, an ǫ-dependent Lyapunov exponent Λ(ǫ) is obtained. For purely stochastic systems, Λ(ǫ) has been shown to display monotonic power-law decay with ǫ [52]; in contrast, chaotic systems are characterized by a well-defined plateau where Λ(ǫ) shows ǫindependent behavior, for not too large ǫ values [52].
The plateau width decreases for increasingD 0 . In contrast, for "large"D 0 = 1, Λ(ǫ) decays monotonically with ǫ. We consider that a transition takes place when the plateau first vanishes, forD 0 ≃D 0,c = 1.6·10 −3 . Actually, starting at this value the kurtosis K(D 0 ) departs from its deterministic value, approaching Gaussian behavior. Moreover, the L-dependence of the threshold valueD 0,c seems weak, see Fig. 3.
We have additionally assessed this transition between chaotic and stochastic GSI behavior in the KS equation satisfied by the space integral of u defined above, namely, where we only consider the non-conserved noise case, relevant within Yakhot's argument [16,42,46]. In Fig. 4 we confirm that, for an increasing noise amplitude D 0 , Λ(ǫ) behaves similarly to Eq. (3): the (D 0 = 0) deterministic system shows a plateau of ǫ-independent behavior which disappears for D 0 0.01, beyond which stochastic behavior ensues. However, now the transition does not reflect into a change between two different fluctuation PDFs. Indeed, in the chaotic case (D 0 = 0) fluctuations of Eq. (5) are known to be Tracy-Widom distributed in the nonlinear regime prior to saturation [16], while Fig. 4 indicates that so do fluctuations for a "large" value of stochastic noise. Scaling exponent values are known to be 1D KPZ, independently of the value of D 0 [16,37].
In summary, we have seen that Yakhot's classic argument on the asymptotic equivalence between SC and GSI holds only partly for Eqs. (3) and (5); specifically, for the conserved KS equation, Eq. (3), different universality classes occur with different field PDF, albeit with the same scaling exponent values. This fact is not captured by Yakhot's argument which, while incorporating the basic system symmetries (conserved vs non-conserved dynamics, etc.), misses key differences between PDFs which are otherwise consistent with the former. The dynamical role of symmetries can be subtle indeed in the present class of nonequilibrium critical systems [4,45,53].
Crucially, each one of the GSI universality classes occurring in Eq. (3) correlates with the nature (chaotic or stochastic) of the mechanism controlling fluctuations in the system, the transition between them nontrivially occurring at a nonzero stochastic noise amplitude. This transition could be experimentally verified, for instance in epitaxial growth of vicinal surfaces [36], where the KS equation describes the dynamics of atomic steps separating terraces under non-negligible adatom desorption [35]. In such a case the equation for the step slope is Eq. (3), whereD 0 scales as an inverse power of the characteristic desorption time.
For the non-conserved KS equation, Eq. (5), although an analogous transition takes place in the dominance of chaotic or stochastic fluctuations, on both sides of the transition the field PDF and the scaling exponent are those of the 1D KPZ universality class. Indeed, the TW distribution is not only relevant to the stochastic 1D KPZ class, but also describes the fluctuations of deterministic chaotic systems [54]. This coincidence might well be accidental and limited to 1D systems. Its exploration in 2D might provide some clue on the relation between the deterministic KS and the stochastic KPZ equa- tions in higher dimensions [55], an open challenge in the fundamental understanding of spatiotemporal chaos [56].
With respect to the specifics of the present transition, it would be interesting to obtain analytical estimates on the threshold noise amplitude and to assess nontrivial consequences on physical quantities beyond the field PDF. The behavior discussed above for S(k, t) already indicates differences in the equal-time two-point statistics, but two-time statistics may introduce additional novelties. In this process, it would be interesting to find analogous transitions, but in which the chaotic and stochastic "phases" might also differ by the values of the scaling exponents.
Finally, the correlation between the field PDF and the nature of the fluctuations underscores the importance of assessing the PDF explicitly, to correctly identify the GSI universality class. This is particularly critical in view of the plethora of experimental complex systems that can be described by a paradigmatic model like the KS equation, in its different (conserved or non-conserved; deterministic or stochastic) forms, especially when both, chaotic and stochastic fluctuations may be operative at comparable space-time scales.

Acknowledgments
This work has been supported by Ministerio de Economía y Competitividad, Agencia Estatal de Investigación, and Fondo Europeo de Desarrollo Regional (Spain and European Union) through grant No. PGC2018-094763-B-I00. E. R.-F. also acknowledges financial support by Ministerio de Educación, Cultura y Deporte (Spain) through Formación del Profesorado Universitario scolarship No. FPU16/06304. While the Kuramoto-Sivashinsy (KS) equation has been derived as a physical model in many different contexts, either in the deterministic case or subject to non-conserved noise, there does not seem to be so many analogous explicit derivations in which this equation comes out perturbed by conserved noise. In what follows, we provide one such derivation in the context of falling liquid films, not far from e.g. a similar one reported in [33] for the case of non-conserved noise.
Consider a liquid film which is falling down an inclined plane (see Fig. 5) and is so thin that thermal fluctuations can no longer be neglected [40,41]. Assuming fluid incompressibility, the evolution equation for the film thickness, h(x, t), can be obtained from the mass balance h t +  5 for coordinate conventions. Actually, the full velocity field (u, v) can be obtained from the balance of linear momentum [33], namely, x y h(x,t) where ρ is the liquid density, µ is the liquid viscosity (both assumed constant), p denotes hydrostatic pressure, g is the acceleration of gravity, and S ij are the components of a symmetric, zero-mean, delta-correlated fluctuation tensor, S. implementing thermal fluctuations in the stress tensor as in classical stochastic hydrodynamics [41], and subindices denote partial derivatives. We consider non-slip, no-penetration boundary conditions u = v = 0 at the planar rigid substrate (y = 0) and a simple stress balance at the free surface of the film (y = h(x, t)), namely, || n T n|| = γ C and || n T t|| = 0, where C is the curvature of the free surface, γ is surface tension, assumed isotropic [33], the unit normal and tangential vector are the stress tensor for this Newtonian fluid, including thermal fluctuations, reads and Π is the disjoining pressure [40,41]. Now, we consider the average thickness, h 0 , of the liquid layer as a typical length scale, w 0 = ρgh 2 0 sin θ/2µ as a velocity scale, w 0 /h 0 as a time scale, and µw 0 /h 0 as a representative scale for pressure and stress, and use all these to rewrite the previous equations in dimensionless units. The resulting momentum balance equations become where Re = ρw 0 h 0 /µ is the Reynolds number. The stress balance at the free surface (y = h) yields Now, we introduce a small parameter ǫ and the new variables x ′ = ǫx, t ′ = ǫt, and v ′ = v/ǫ, adapted to a lubrication approximation [33] within which the cross-stream dimension of the film will be considered much smaller than its streamwise extent. We consider the capillary number, Ca = µw 0 /γ, to be order ǫ 2 and define Ca ′ = Ca/ǫ 2 . We expand u = u 0 + ǫu 1 + O(ǫ 2 ), v ′ = v ′ 0 + ǫv ′ 1 + O(ǫ 2 ) and p = p 0 + ǫp 1 + O(ǫ 2 ), and consider S xx , S yy to be O(ǫ −2 ) and S xy , S yx to be O(ǫ −1 ) [41]. Last, by defining S xx′ = S xx /ǫ 2 , S yy ′ = S yy /ǫ 2 , S xy′ = S xy /ǫ, and S yx′ = S yx /ǫ the momentum balance equations and the surface boundary conditions become, respectively, We can now compute the velocity profile u = u 0 + ǫu 1 + O(ǫ 2 ). At O(1), Eq. (A7) becomes u 0yy = −2. As u 0y = 0 at the fluid surface y = h [leading order of Eq. (A10)] and u 0 = 0 at the substrate y = 0, we obtain u 0 = 2 hy − y 2 /2 . Considering the fluid film to be ultrathin, Re ≪ 1 can be neglected, and Eq. (A8) at O(ǫ) becomes u 1yy = p 0x ′ − S xy′ y . Here we have u 1y = −S yx [Eq. (A10) at O(ǫ)] and u 1 = 0 as boundary conditions at the fluid surface and the substrate, respectively, which allow us to obtain the profile for u 1 = −p 0x ′ hy − y 2 /2 − y 0 S yx′ dy. The p 0 contribution can be obtained from Eq. (A7) at O(1), p 0y = −2 cot θ, with p 0 = −Π − h x ′ x ′ /Ca ′ as boundary condition at the fluid surface, Finally, mass balance reads h t ′ + h 0 u 0 + ǫ u 1 dy yields the evolution equation Taking y 0 where η is zero average, Gaussian white noise [41] and Π = −φ ′ , where φ is the interface potential [41], we finally obtain Finally, a weakly-nonlinear expansion allows us to get the KS equation with conserved noise from Eq. (A13). Considering very small fluctuations around the flat film solution, h = 1 + ǫh, Eq. (A13) becomes where By defining κ 0 = 1/(3Ca ′ ), ν 0 = φ ′′ (1) − 2 cot θ/3, andη = η/ǫ,  As can be seen, the data collapse in both cases to the expected scaling form, C(x, t) = t 2β c(x/t 1/z ), with c(y) ∼ cst. − y 2α for y ≪ 1 and 0 for y ≫ 1 [8,9], using the corresponding values of the scaling exponents as determined in the MT. For the h field, collapse is to the exact covariance of the Airy 1 process, as expected in the growth regime for 1D KPZ scaling with periodic boundary conditions, see references, e.g. in [12]. In general, we can observe that the fluctuation distributions remain largely unchanged over time. Panels (a-c) actually show that, in the corresponding systems (and at variance with the behavior of the 1D KPZ equation), the PDF within the nonlinear time evolution actually coincides with the corresponding PDF at saturation. Such a saturation PDF is reported in [43] for the deterministic KS equation [panel (a)], while in the cases of the stochastic KS equation with conserved [panel (b)] and nonconserved [panel (c)] noise results are only available for the corresponding stochastic Burgers equation with which they share asymptotic scaling behavior (in terms of scaling exponents and, presumably, of fluctuation PDF, given that these are stochasticsdominated systems with noise amplitudes above threshold, see Figs. 2 and 3), reported in [18] and [45], respectively. On the other hand, while the full PDF corresponding to panels (a,b,d) is provided in the MT, this is not the case for the stochastic Kuramoto-Sivashinsky equation with non-conserved noise [panel (c)], whose PDF in the nonlinear regime prior to saturation is presented in Fig. 8.
More specifically, the time evolution of the skewness and excess kurtosis shown in Fig. 7 implies a PDF which exhibits a symmetric, notably platykurtic, non Gaussian behavior for the deterministic conserved KS equation [panel (a)], and almost Gaussian behavior in the conserved KS equation with conserved [panel (b)] or non-conserved noise [panel (c)]. Finite-size deviations of the excess kurtosis from their exact zero values seen in panels (b,c) are comparable to similar deviations in the corresponding cases of the stochastic Burgers equation with conserved or non-conserved noise, respectively [18,45]. Finally, the PDF reproduces closely the expected GOE-TW behavior for the non-conserved KS equation with non-conserved noise, Eq. (5), see Fig. 7(d). However, this case is well-known to feature quite different PDF behavior at steady state, as is the case in the 1D KPZ universality class [12].  c. Transition between chaotic and stochastic PDF Figure 9 depicts a more detailed view than that provided by Fig. 2, on the transition between the deterministic (chaotic) and the stochastic PDF, which takes place in the KS equation, Eq. (3), for increasing values of the conserved-noise amplitude. In the MT, we have identified the threshold valueD 0,c ≃ 0.0016 as that value ofD 0 above which the kurtosis departs from its deterministic (D 0 = 0) value, see Fig. 3. Moreover, forD 0 >D 0,c the scale-dependent Lyapunov exponent Λ(ǫ) changes qualitative behavior from chaos-to stochastic-dominated fluctuations, see Fig. 4. In contrast, Fig. 9 illustrates how this change is more difficult to see by naked-eye inspection of the form of the field PDF. Note that for all cases considered in this figure, D 0 >D 0,c . The two leftmost panels display PDF with inflection points which are akin those characteristic of the PDF for purely chaotic fluctuations (D 0 = 0). Nevertheless, such inflection points disappear once the stochastic-noise amplitude increases even further aboveD 0,c , beyond which the fluctuation PDF eventually reaches fully-Gaussian form, see Fig. 2.