Research Article

Fibonacci family of dynamical universality classes

Vladislav Popkov, Andreas Schadschneider, Johannes Schmidt, and Gunter M. Schütz
PNAS October 13, 2015 112 (41) 12645-12650; first published September 30, 2015; https://doi.org/10.1073/pnas.1512261112

  1. Edited by Giorgio Parisi, University of Rome, Rome, Italy, and approved August 31, 2015 (received for review June 23, 2015)

Significance

Universality is a well-established central concept of equilibrium physics. It asserts that, especially near phase transitions, the properties of a physical system do not depend on its details such as the precise form of interactions. Far from equilibrium, such universality has also been observed, but, in contrast to equilibrium, a deeper understanding of its underlying principles is still lacking. We show that the two best-known examples of nonequilibrium universality classes, the diffusive and Kardar−Parisi−Zhang classes, are only part of an infinite discrete family. The members of this family can be identified by their dynamical exponent, which, surprisingly, can be expressed by a Kepler ratio of Fibonacci numbers. This strongly indicates the existence of a simpler underlying mechanism that determines the different classes.

Abstract

Universality is a well-established central concept of equilibrium physics. However, in systems far away from equilibrium, a deeper understanding of its underlying principles is still lacking. Up to now, a few classes have been identified. Besides the diffusive universality class with dynamical exponent z=2, another prominent example is the superdiffusive Kardar−Parisi−Zhang (KPZ) class with z=3/2. It appears, e.g., in low-dimensional dynamical phenomena far from thermal equilibrium that exhibit some conservation law. Here we show that both classes are only part of an infinite discrete family of nonequilibrium universality classes. Remarkably, their dynamical exponents zα are given by ratios of neighboring Fibonacci numbers, starting with either z1=3/2 (if a KPZ mode exist) or z1=2 (if a diffusive mode is present). If neither a diffusive nor a KPZ mode is present, all dynamical modes have the Golden Mean z=(1+5)/2 as dynamical exponent. The universal scaling functions of these Fibonacci modes are asymmetric Lévy distributions that are completely fixed by the macroscopic current density relation and compressibility matrix of the system and hence accessible to experimental measurement.

The Golden Mean, φ=1/2+5/2≈1.61803..., also called Divine Proportion, has been an inspiring number for many centuries. It is widespread in nature; i.e., arrangements of petals of the flowers and seeds in the sunflower follow the golden rule (1). Being considered an ideal proportion, the Golden Mean appears in famous architectural ensembles such as the Parthenon in Greece, the Giza Great Pyramids in Egypt, or Notre Dame de Paris in France. Ideal proportions of the human body follow the Golden Rule.

Mathematically, the beauty of the Golden Mean number is expressed in its continued fraction representation: All of the coefficients in the representation are equal to unity,φ=1+11+ 11+ 11+⋱.[1]Systematic truncation of the above continued fraction gives the so-called Kepler ratios, 1/1,2/1,3/2,5/3,8/5,..., which approximate the Golden Mean. Subsets of denominators (or numerators) of the Kepler ratios form the celebrated Fibonacci numbers, Fi=1,1,2,3,5,8,.., such that Kepler ratios are ratios of two neighboring Fibonacci numbers. As well as the Golden Mean, Fibonacci ratios and Fibonacci numbers are widespread in nature (1).

The occurrence of the Golden Mean is not only interesting for aesthetic reasons but often indicates the existence of some fundamental underlying structure or symmetry. Here we demonstrate that the Divine Proportion as well as all of the truncations (Kepler ratios) of the continued fraction (Eq. 1) appear as universal numbers, namely, the dynamical exponents, in low-dimensional dynamical phenomena far from thermal equilibrium. The two well-known paradigmatic universality classes, Gaussian diffusion with dynamical exponent z=2 (2, 3) and the Kardar−Parisi−Zhang (KPZ) universality class with z=3/2 (4), enter the Kepler ratios hierarchy as the first two members of the family.

The universal dynamical exponents in the present context characterize the self-similar space−time fluctuations of locally conserved quantities, characterizing, e.g., mass, momentum, or thermal transport in one-dimensional systems far from thermal equilibrium (5). The theory of nonlinear fluctuating hydrodynamics (NLFH) has recently emerged as a powerful and versatile tool to study space−time fluctuations, and specifically the dynamical structure function that describes the behavior of the slow relaxation modes, and from which the dynamical exponents can be extracted (6).

The KPZ universality class has been shown to explain the dynamical exponent observed in interface growth processes as diverse as the propagation of flame fronts (7, 8), the growth of bacterial colonies (9), or the time evolution of droplet shapes such as coffee stains (10) where the Gaussian theory fails. For a nice introduction into the KPZ class and its relevance, we refer to ref. 11. Recent reviews (12, 13) provide a more detailed account of theoretical and experimental work on the KPZ class. The dynamical structure function originating from the one-dimensional KPZ equation has a nontrivial scaling function obtained exactly by Prähofer and Spohn from the totally asymmetric simple exclusion process (TASEP) and the polynuclear growth model (14, 15) and was beautifully observed in experiments on turbulent liquid crystals (16, 17). The theoretical treatment, both numerical and analytical, of generic model systems with Hamiltonian dynamics (18), anharmonic chains (19, 20), and lattice models for driven diffusive systems (21, 22) have demonstrated an extraordinary robust universality of fluctuations of the conserved slow modes in one-dimensional systems.

Despite this apparent ubiquity, dynamical exponents different from z=2 or z=3/2 were observed frequently. Usually, it is not clear whether this corresponds to genuinely different dynamical critical behavior or is just a consequence of imperfections in the experimental setting. Moreover, recently, a new universality class with dynamical exponent z=5/3 for the heat mode in Hamiltonian dynamics (18) was discovered, followed by the discovery of some more universality classes in anharmonic chains (19, 20) and lattice models for driven diffusive systems (21, 22). What is lacking, even in the conceptually simplest case of the effectively one-dimensional systems that we are considering, is the understanding of the plethora of dynamical nonequilibrium universality classes within a larger framework. Such a framework exists, e.g., for 2D critical phenomena in equilibrium systems where the spatial symmetry of conformal invariance together with internal symmetries give rise to discrete families of universality classes in which all critical exponents are simple rational numbers.

It is the aim of this article to demonstrate that discrete families of universality classes with fractional critical exponents appear also far from thermal equilibrium. This turns out to be a hidden feature of the NLFH equations that we extract using mode coupling theory. It is remarkable that one finds dynamical exponents zα, which are ratios of neighboring Fibonacci numbers {1,1,2,3,5,8,…} defined recursively as Fn=Fn−1+Fn−2. The first two members of this family are diffusion (z=2=F3/F2) and KPZ (z=3/2=F4/F3). The corresponding universal scaling functions are computed and shown to be (in general asymmetric) zα-stable Lévy distributions with parameters that can be computed from the macroscopic current density relation and compressibility matrix of the corresponding physical system and which thus can be obtained from experiments without detailed knowledge of the microscopic properties of the system. The theoretical predictions, obtained by mode coupling theory, are confirmed by Monte Carlo simulations of a three-lane asymmetric simple exclusion process, which is a model of driven diffusive transport of three conserved particle species.

Nonlinear Fluctuating Hydrodynamics

We consider a rather general interacting nonequilibrium system of length L described macroscopically by n conserved order parameters ρλ(x,t) with stationary values ρλ and associated macroscopic stationary currents jλ(ρ1,…,ρn) and compressibility matrix K with matrix elements Kλμ=(1/L)〈(Nλ−ρλL)(Nμ−ρμL)〉 where Nλ=∫0Ldxρλ(x,t) are the time-independent conserved quantities.

The starting point for investigating density fluctuations uλ(x,t):=ρλ(x,t)−ρλ in the nonequilibrium steady state are the NLFH equations (5)∂tu→=−∂x(Ju→+12〈u|H→|u〉−∂xDu→+Bξ→)[2]where J is the current Jacobian with matrix elements Jλμ=∂jλ/∂ρμ, H→ is a column vector whose entries (H→)λ=Hλ are the Hessians with matrix elements Hμνλ=∂2jλ/(∂ρμ∂ρν), and the bra-ket notation represents the inner product in component space 〈u|(H→)λ|u〉=u→THλu→=∑μνuμuγHμνλ with 〈u|=u→T and |u〉=u→. The diffusion matrix D is a phenomenological quantity. The noise term Bξ→ does not appear explicitly below, but plays an indirect role in the mode coupling analysis. The product JK of the Jacobian with the compressibility matrix K is symmetric (23), which guarantees a hyperbolic system of conservation laws (24). We ignore possible logarithmic corrections arising from cubic contributions (25).

This system of coupled noisy Burgers equations is conveniently treated in terms of normal modes ϕ→=Ru→, where RJR−1=diag(vα) and the transformation matrix R is normalized such that RKRT=1. The eigenvalues vα of J are the characteristic velocities of the system. From Eq. 2, one thus arrives at ∂tϕα=−∂x(vαϕα+〈ϕ|Gα|ϕ〉−∂x(D˜ϕ→)α+(B˜ξ→)α), with D˜=RDR−1, B˜=RB, and the mode coupling matricesGα=12∑βRαβ(R−1)THβR−1,[3]whose matrix elements we denote by Gβγα.

Computation of the Dynamical Structure Function

The dynamical structure function describes the stationary fluctuations of the conserved slow modes and is thus a key ingredient for understanding the interplay of noise and nonlinearity and their role for transport far from equilibrium. We focus on the case of strict hyperbolicity where all vα are pairwise different and study the large-scale behavior of the dynamical structure function Sαβ(x,t)=〈ϕα(x,t)ϕβ(0,0)〉. Because all modes have different velocities, only the diagonal elements Sα(x,t):=Sαα(x,t) are nonzero for large times. Mode coupling theory yields (5)∂tSα(x,t)=(−vα∂x+Dα∂x2)Sα(x,t)+∫0tds∫ℝdySα(x−y,t−s)∂y2Mαα(y,s)[4]with the diagonal element Dα:=D˜αα of the phenomenological diffusion matrix for the eigenmodes and the memory kernel Mαα(y,s)=∑β,γ(Gβγα)2Sβ(y,s)Sγ(y,s). The task therefore is to extract for arbitrary n the large-time and large-distance behavior from this nonlinear integro-differential equation.

Remarkably, these equations can be solved exactly in the long-wavelength limit and for t→∞ by Fourier and Laplace transformation (see Materials and Methods and SI Text). Using a suitable scaling ansatz for the transformed structure function then allows analysis of the small-p behavior from which the dynamical exponents can be determined. We find that different conditions arise depending on which diagonal elements of the mode coupling matrices vanish.

Fibonacci Family of Dynamical Universality Classes

Fibonacci Case.

First, we consider the case where the self-coupling Gααα is nonzero for one mode only, e.g., G111≠0. For all other modes α > 1, we assume a single nonzero coupling to the previous mode, so Gα−1,α−1α≠0, and Gβ,βα=0 for β≠α−1. Then, as follows from our analysis (see Materials and Methods), we find the following recursion for the dynamical exponents:zα=1+1zα−1[5]with z1=3/2.

The dynamical structure function in momentum space is proportional to the zα-stable Lévy distribution with maximal asymmetry σα=±1; see ref. 26 and Eq. 13 below. The sign of the asymmetry depends on whether the mode (α−1) has bigger or smaller velocity than the mode α, σα=−sgn(vα−vα−1). The dynamical exponents (Eq. 5) form a sequence of rational numbers,zα=Fα+3Fα+2,[6]which are consecutive ratios of neighboring Fibonacci numbers Fα, defined by Fα=Fα−1+Fα−2 with initial values F0=0,F1=1, which converge exponentially to the Golden Mean φ:=12(1+5)≈1.618, as first observed by Kepler in 1611 in a treatise on snowflakes. In a model with n conservation laws, one has the Fibonacci modes with dynamical exponents {3/2,5/3,8/5, …, zn}.

Finally, we remark that if mode 1 is diffusive rather than KPZ, then we find the same sequence (Eq. 6) of exponents, except that it starts with z1=F2=2.

In Fig. 1, we show some representative examples of the scaling functions that are quite different in shape. Furthermore, the relation between the exponents zα, determined by Eq. 11 in Materials and Methods, and the mode coupling matrices Gα is illustrated for the case n=2.

Fig. 1.
Fig. 1.

The scaling functions (Bottom) and dynamical exponents are related to the structure of the mode coupling matrices Gα (Top). The table shows the dynamical exponents zα in the case n=2; see Eq. 11. The symbols ∗ and ⋆ denote nonzero elements. Red symbols correspond to self-coupling; black symbols correspond to couplings to other modes. Matrix elements not indicated can take any value. The colors in the table correspond to the colors of the graphs of the scaling functions.

Golden Mean Case.

As a second representative example, we consider the case where all self-coupling coefficients vanish, Gααα≡0 for all α, while each mode has at least one nonzero coupling to another mode, Gββα≠0 for some β≠α. Then, Eq. 5 reduces to zα=1+1/zβ for all modes α, β. The unique solution of this equation is the Golden Mean zα=φ=(1+5)/2 for all α. The scaling functions (see SI Text) are proportional to φ-stable Lévy distributions with parameters fixed by the collective velocities and the mode coupling coefficients. The asymmetry of the fastest right-moving (left-moving) mode is predicted to be β=−1 (β=1).

Simulation Results

To check the theoretical predictions for the two cases, we simulate mass transport with three conservation laws, i.e., three distinct species of particles. To maintain a far-from-equilibrium situation, a driving force is applied that leads to a constant drift superimposed on undirected diffusive motion. This is a natural setting for transport of charged particles in nanotubes (see Fig. 2 for an illustration), where a direct measurement of the stationary particle currents is experimentally possible (27). However, due to the universal applicability of NLFH, the actual details of the interaction of the particles with their environment and the driving field are irrelevant for the theoretical description of the large-scale dynamics. Hence, for good statistics, we simulate a lattice model for transport that represents a minimal realization of the essential ingredients, namely, a nonlinear current density relation for all three conserved masses.

Fig. 2.
Fig. 2.

Schematic drawing of three particle species drifting inside a nanotube. Due to the interaction between the particles and with the walls, one expects a nonlinear current density relation.

Our model is the three-species version of the multilane TASEP (28). Particles hop randomly in field direction on three lanes to their neighboring sites on a periodic lattice of 3×L sites with rates that depend on the nearest-neighbor sites. Lane changes are not allowed so that the total number of particles on each lane is conserved. Due to excluded volume interaction, each lattice site can be occupied by at most one particle. Thus, the occupation numbers nk(λ) of site k on lane λ take only values 0 or 1. The hopping rate rk(λ) from site k on lane λ to site k+1 on the same lane is given byrk(λ)=bλ+12∑{μ:μ≠λ}γλμ(nk(μ)+nk+1(μ))[7]with a species-dependent drift parameter bλ and symmetric interaction constants γλμ=γμλ. Hopping attempts onto occupied sites are rejected. The conserved quantities are the three total numbers of particles Nλ on each lane with corresponding densities ρλ=Nλ/L.

The stationary distribution of our model factorizes (28) and thus allows for the exact computation of the macroscopic current density relations jλ(ρ1,ρ2,ρ3) and the compressibility matrix K(ρ1,ρ2,ρ3). Furthermore, because there is no particle exchange between lanes, the compressibility matrix is diagonal, with elements denoted by κλ. One hasjλ=ρλ(1−ρλ)(bλ+∑{μ:μ≠λ}γλμρμ)[8]κλ=ρλ(1−ρλ).[9]The diagonalization matrix R and the mode coupling matrices Gα are fully determined by these quantities.

According to mode coupling theory, three different Fibonacci modes with z1=3/2,z2=5/3,z3=8/5 occur, e.g., when G111≠0, G112≠0, G223≠0, and G222=G332=G333=0. For our simulation, we compute numerically densities, bare hopping rates, and interaction parameters to satisfy these properties as described in Materials and Methods. For this choice, the velocities of the normal modes are v1=0.592315, v2=0.0281578, and v3=1.58226, which ensures a good spatial separation after quite small times. The propagation of the three normal modes (Fig. 3) with the predicted velocities is observed with an error of less than 10−3. Moreover, the numerically obtained dynamical structure function for mode 3 shows a startling agreement with the theoretically predicted Lévy scaling function with z=8/5 and maximal asymmetry (see Fig. 4). It takes longer for the other two modes (KPZ mode and Lévy stable 5/3 mode) to reach their asymptotic form, which we argue is due to the much smaller respective couplings, (G111/G223)2≪1, (G112/G223)2≪1.

Fig. 3.
Fig. 3.

Space−time propagation of three normal modes in the three-lane model. The modes (from left to right) are the Fibonacci mode with z=8/5 (mode 3), the KPZ mode with z=3/2 (mode 1), and the Fibonacci mode with z=5/3 (mode 2). The physical and simulation parameters are given in Materials and Methods.

Fig. 4.
Fig. 4.

(Left) Vertical least squares fit of the numerically obtained dynamical structure function for the Fibonacci 8/5 mode (points), at time t=1,000 with an 8/5-stable Lévy distribution, maximal asymmetry −1, and theoretical center of mass (line), predicted by the mode coupling theory. The only fit parameter is the scale parameter of the Lévy stable distribution. The simulation results agree very well with the asymptotic theoretical result already for moderate times. (Right) Close-ups of the peak region and tail regions, according to a color code. Every tenth data point is plotted, to improve the visibility of the data. The statistical error ϵ99% with 99% confidence bound is for every data point smaller than 1.6299⋅10−5.

To observe three Golden Mean modes, it is sufficient to require that each mode has zero self-coupling and at least one nonzero coupling to other modes. This can be achieved with the set of parameters given in Materials and Methods, which lead to the velocities v1=1.83149, v2=0.762688, and v3=0.326778 of the normal modes. The propagation of the three normal modes with the predicted velocities is observed, approaching, for large times, a very small relative error of about 10−4. The structure function for the fastest mode 1 converges to its asymptotic shape faster than for the other modes, due to the large coupling coefficient G331. In Fig. 5, we show a scaling plot of the measured structure function for mode 1 with dynamical exponent z≡φ=(1+5)/2 together with a fitted to a φ-stable Lévy function (Eq. 13) with maximal asymmetry β=−1 as predicted by the theory. The data collapse shows a striking agreement between the measured and theoretical scaling function. Alternatively, the dynamical exponent zα can be derived from the maximum of the structure function, which scales as max(S1(x,t))=const⋅t−1/z. We obtain z≈1.63, which differs from the predicted value z≡φ=(1+5)/2 by less than 0.8%.

Fig. 5.
Fig. 5.

Scaling plot of the measured structure function of mode 1 with dynamical exponent z≡φ=(1+5)/2 for the Golden Mean case, fitted to a φ-stable Lévy distribution with maximal asymmetry −1 (see Eq. 13). The scale parameter E1 for the Lévy-stable distribution and the center of mass velocity v1 are obtained by a vertical least square fit. Fitted parameters are v1,fit=1.83107±0.00009 and E1,fit=1.1±0.01. The fitted velocity v1,fit differs by 0.02% from the theoretical velocity.

Discussion

Our work demonstrates that nonequilibrium phenomena are much richer than just the diffusive and KPZ universality suggest. We have established that in nonequilibrium phenomena governed by NLFH with n conservation laws, mode coupling theory predicts a family of dynamical universality classes with dynamical exponents given by the sequence of consecutive Kepler ratios (Eq. 6) of Fibonacci numbers. With slightly modified initial conditions on G111, this result is easily generalized for the case when the first mode α=1 is diffusive. Then the sequence of dynamical exponents becomes shifted by one unit with respect to Eq. 6. On the other hand, if all self-couplings vanish, but at least one other diagonal element Gββα of the mode coupling matrix is non-zero, one has, as a unique solution for all modes α, the fixed point value zα=z∞=φ, which is the Golden Mean.

For general mode coupling matrices, all critical exponents can be computed (from Eq. 11 in Materials and Methods). The scaling functions of the nondiffusive and non-KPZ modes are asymmetric Lévy distributions whose parameters are completely determined by the macroscopic current density relation and compressibility matrix of the system.

For 1+1 dimensional systems out of equilibrium, this is the first time, to our knowledge, that an infinite family of discrete universality classes is found. Recalling that 1 + 1 dimensional nonequilibrium systems with short-range interactions can be mapped onto 2D equilibrium systems (with the time evolution operator playing the role of the transfer matrix), one is reminded of the discrete families of conformally invariant critical equilibrium systems in two space dimensions (29, 30). We do not know whether there is any mathematical link, but the analogy is suggestive in so far as conformal invariance is a local symmetry of spatially isotropic systems with z=1 (which happens to be the lowest-order Kepler ratio) whereas z>1 corresponds to strongly anisotropic systems for which local symmetry groups are also known to exist (31).

Because an infinite number of lanes of coupled one-dimensional systems correspond to a 2D system, it is intriguing to observe that the Golden Mean is close to the numerical value z=1.612−1.618 of the dynamical exponent of the 2+1-dimensional KPZ equation (32, 33). The scaling function of the 2+1-dimensional KPZ equation, however, is not Lévy (34).

To observe and distinguish between the different new classes, highly precise experimental data will be required. For example, in the Fibonacci case, the dynamical exponents converge quickly to the Golden Mean. A feature that might be easier to observe experimentally is the scaling function itself, which, for higher Fibonacci ratios 5/3,8/5, …, usually has a strong asymmetry (see Figs. 1, 4, and 5), whereas KPZ and Gauss scaling functions are symmetric. Growth processes that can be mapped on exclusion processes with several conservation laws might be potentially suitable candidates for an experimental verification; see, e.g., refs. 16 and 17 for an example of a system with one conservation law.

Materials and Methods

Computation of the Dynamical Structure Function.

The mode coupling Eq. 4 can be solved in the scaling limit by applying a Fourier transform (FT) f(x)→f˜(p) and a Laplace transform (LT) f(t)→f˜(ω). For more details, we refer to ref. 22, where the case n=2 of two conservation laws has been treated. After making the scaling ansatzS˜α(p,ω˜α)=p−zαgα(ζα)[10]

for the transformed dynamical structure function where S^α(p,0)=1/2π and ζα=ω˜α|p|−zα, we are in a position to analyze the small-p behavior. One has to search for dynamical exponents for which the limit p→0 is nontrivial, which requires a self-consistent treatment of all modes. We find that different conditions arise depending on which diagonal elements of the mode coupling matrices vanish. To characterize the possible scenarios, we define the set Iα:={β:Gββα≠0} of nonzero diagonal mode coupling coefficients. Through power counting one obtainszα={2ifIα=Ø3/2ifα∈Iαminβ∈Iα[(1+1zβ)]else[11]

and the domain1<zα≤2 ∀α.[12]

for the possible dynamical exponents.

In the Fibonacci case, the dynamical structure function of mode α in momentum space has the scaling formS^α(p,t)=12πe−ivαpt−Eα|p|zαt(1−iσpαβtan(πzα2))[13]

with inverse time scales Eα. The dynamical exponents then satisfy the recursion (Eq. 5). Up to the normalization 1/2π the scaling form (Eq. 13) is an α-stable Lévy distribution (26).

Simulation Algorithm.

For the Monte Carlo simulation of the model, we choose a large system size L≥5⋅105, which avoids finite-size effects. At time t=0, Nλ particles are placed on each lane according to the desired initial state. One Monte Carlo time unit consists of 3⋅L⋅r* random sequential update steps where r*=max{rk(λ)}: In each update step, a bond (k(λ),k(λ)+1) is chosen randomly with uniform distribution. If nk(λ)(1−nk+1(λ))=1, then the particle at site k is moved to k+1 with probability rk(λ)/r* where r* is the maximal value that the rk(λ) can take among all possible particle configurations on the neighboring lanes. If nk(λ)(1−nk+1(λ))=0, the particle configuration remains unchanged.

Simulation of the Dynamical Structure Function.

To determine the dynamical structure function, we initialize the system by placing Nλ particles uniformly on each lane λ. This yields a random initial distribution drawn from the stationary distribution of the process. No relaxation is required.

Then we use translation invariance and compute the space and time averageσL,kλμ(M,τ,t)=1M∑j=1M1L∑l=1Lnl+k(λ)(jτ+t)nl(μ)(jτ)−ρλρμ.[14]

To avoid noisy data of σL,kλμ, we take in Eq. 14 the system size L and the time average parameter M sufficiently large. To obtain Skλμ(t), we average over P independently generated and propagated initial configurations of σL,kλμ. The error estimates for Skλμ(t) are calculated from the P independent measurements. From Skλμ(t), we compute the structure function of the normal modes by transformation with the diagonalizing matrix R determined by Eqs. 8 and 9.

To obtain model parameters for three different Fibonacci modes with z1=3/2,z2=5/3,z3=8/5, we solve the equations given in the text after Eq. 9 numerically with a C program performing direct minimization of the absolute values of the targeted G elements until the given tolerance value (10−6) is reached. The data shown here for the three-mode case have been obtained from simulations with densities ρ1=0.2, ρ2=0.25, ρ1=0.3, bare hopping rates b1=0.613185, b2=0.425714, b3=0.799831, and interaction parameters γ12=1.36145, γ23=3.69786, γ13=0.143082 for which the needed relations are satisfied. This choice of parameters yields G111=0.322507, G112=−0.15, G223=1.04547, while the absolute values of G222, G332, G333 are smaller than 10−6. Besides these physical parameters, the simulation parameters for the Fibonacci modes (Figs. 3 and 4) are L=5·105, τ=250, M=1,400, P=98.

For the Golden Mean case (Fig. 5), the set of parameters is ρ1=0.2, ρ2=ρ3=0.25, γ12=0.0082334758646, γ23=1.68447706968, γ13=3.72140740146, and b1=0.905073261248, b2=0.86, and b3=1.18875738638. This leads to G221=0.405702, G331=0.929315, G112=−0.104141, G332=−0.208477, G113=−0.182467, G223=0.271246, while the absolute values of G111,G222,G333 are smaller than 10−6. The simulation parameters for the Golden Mean case are L=5·106 τ=750, M=30 and P=303.

SI Text

Remarkably, the mode coupling equations (Eq. 4) can be solved exactly in the scaling limit by FT and LT. To this end, we define the FT asf^(p):=12π∫−∞∞dx e−ipxf(x)[S1]and the LT asf˜(ω):=∫0∞dt e−ωtf(t).[S2]With D^α(p)=ivαp+Dαp2, we obtain from Eq. 4, in momentum−frequency space,S˜α(p,ω)=S^α(p,0)ω+D^α(p)+C˜αα(p,ω)[S3]with memory kernelC˜αα(p,ω)=2∑β,γ(Gβγα)2p2∫0∞ds e−ωs∫ℝdqS^β(q,s)S^γ(p−q,s).[S4]and S^α(p,0)=1/2π.

Next, we introduce ω˜α:=ω+ivαp and make the scaling ansatzS˜α(p,ω˜α)=p−zαgα(ζα)[S5]with ζα=ω˜α|p|−zα. Having in mind systems with short-range interactions, we anticipate that all modes spread subballistically, i.e., zα>1 for all α. Using strict hyperbolicity, one obtains, after some substitutions of variables,gα(ζα)=limp→0[ζα+Dα|p|2−zα+Qααζα2−zα−1zα|p|3−2zα+∑β≠αQαβ(−ivpαβ)1zβ−1|p|1+1zβ−zα]−1.[S6]with vpαβ:=|vα−vβ|sgn[p(vα−vβ)] andQαβ=2(Gββα)2Γ(1−1zβ)Ω[S^β]≥0,[S7]whereΩ[f^]=∫−∞∞dp f^(p)f^(−p).[S8]With σpαβ=sgn[p(vα−vβ)], one has(−ivpαβ)1zβ−1=sin(π2zβ)|vα−vβ|1zβ−1⋅[1−iσpαβtan((1+1zβ)π2)].[S9]

Now we are in a position to analyze the small-p behavior. One has to search for dynamical exponents for which the limit p→0 is nontrivial, which is determined by the smallest power of p in Eq. S6. This has to be done self-consistently for all modes. We find that different conditions arise depending on which diagonal elements of the mode coupling matrices vanish. To characterize the possible scenarios, we define the set Iα:={β:Gββα≠0} of nonzero diagonal mode coupling coefficients. One obtains, from Eq. S6, through power counting the system of equations,zα={2ifIα=Ø3/2ifα∈Iαminβ∈Iα[(1+1zβ)]else[S10]and the domain1<zα≤2 ∀α[S11]for the possible dynamical exponents.

Acknowledgments

We thank Herbert Spohn for helpful comments on a preliminary version of the manuscript. This work was supported by Deutsche Forschungsgemeinschaft under Grant SCHA 636/8-1.

Footnotes

  • ↵1Present address: Helmholtz Institut für Strahlen- und Kernphysik, University of Bonn, 53115 Bonn, Germany.

  • ↵2To whom correspondence should be addressed. Email: as{at}thp.uni-koeln.de.

  • Author contributions: A.S. and G.M.S. designed research; V.P., J.S., and G.M.S. performed research; V.P., A.S., J.S., and G.M.S. analyzed data; and V.P., A.S., and G.M.S. wrote the paper.

  • The authors declare no conflict of interest.

  • This article is a PNAS Direct Submission.

  • This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1512261112/-/DCSupplemental.

View Abstract
Back to top
Article Alerts
Email Article
Citation Tools
Request Permissions
Fibonacci family of dynamical universality classes
Vladislav Popkov, Andreas Schadschneider, Johannes Schmidt, Gunter M. Schütz
Proceedings of the National Academy of Sciences Oct 2015, 112 (41) 12645-12650; DOI: 10.1073/pnas.1512261112
Digg logo Reddit logo Twitter logo Facebook logo Google logo Mendeley logo
Proceedings of the National Academy of Sciences: 112 (41)
Table of Contents

Submit

Article Classifications

Jump to section

You May Also be Interested in

Large piece of gold
News Feature: Tracing gold's cosmic origins
Astronomers thought they’d finally figured out where gold and other heavy elements in the universe came from. In light of recent results, they’re not so sure.
Image credit: Science Source/Tom McHugh.

Similar Articles