Fluctuating hydrodynamics approach to equilibrium time correlations for anharmonic chains
FFluctuating hydrodynamics approach to equilibriumtime correlations for anharmonic chains
Herbert Spohn
Zentrum Mathematik and Physik Department, Technische Universit¨at M¨unchen,Boltzmannstraße 3, 85747 Garching, Germany. [email protected]
Abstract.
Linear fluctuating hydrodynamics is a useful and versatile tool for describingfluids, as well as other systems with conserved fields, on a mesoscopic scale. In onespatial dimension, however, transport is anomalous, which requires to develop a nonlinearextension of fluctuating hydrodynamics. The relevant nonlinearity turns out to be thequadratic part of the Euler currents when expanding relative to a uniform background.We outline the theory and compare with recent molecular dynamics simulations.1 a r X i v : . [ c ond - m a t . s t a t - m ec h ] J a n able of contents
1. Introduction, long time tails for simple fluids2. Anharmonic chains2.1 Conservation laws and equilibrium time correlations2.2 Linearized hydrodynamics3. Nonlinear fluctuating hydrodynamics3.1 Euler currents to second order3.2 Stationary measure for the physical fields3.3 Transformation to normal modes4. Mode-coupling theory4.1 Decoupling hypothesis4.2 One-loop, diagonal, and small overlap approximations4.3 Numerical simulations of the mode-coupling equations4.4 Asymptotic self-similarity4.5 No signal beyond the sound cone4.6 Dynamical phase diagram5. Molecular dynamics simulations6. Total current correlations7. Other Hamiltonian chains7.1 Coupled rotators7.2 Lattice nonlinear Schr¨odinger equation2
Introduction, long time tails for simple fluids
In the mid 1950ies, Green [1] and Kubo [2, 3] discovered that transport coefficients forsimple fluids can be obtained through a time-integral over the respective total currentcorrelation function. For tracer diffusion such a connection is more immediate and wasunderstood much earlier. But the then novel insight was that collective transport co-efficients, such as viscosity and thermal conductivity, follow the same pattern. Thus itbecame a central issue to determine the time decay of such current correlations. Withessentially no tools available this amounted to an impossible task. The static equilibriumcorrelations were known to decay exponentially fast, as confirmed by a convergent seriesexpansion. But for the dynamics one would have to deal with a huge set of coupleddifferential equations. At the time the only theoretical tool available was the Boltzmannequation valid at low density. Kinetic theory predicts an exponential decay for the currenttime correlations and it was tacitly assumed that such behavior would extend to moderatedensities. Alder and Wainwright [4] tried to check the situation in a pioneering moleculardynamics simulation of 500 hard disks, resp. hard spheres, at periodized volume 2, 3,and 5 times larger than close packing. For tracer diffusion they convincingly observed apower law decay as t − d/ , dimension d = 2 ,
3, which was baptized “long time tail”. Theyalso argued that the same behavior should hold for collective transport. Theory quicklyjumped in and predicted a decay as t − d/ for viscosity and thermal conductivity [5, 6].There are several theoretical schemes and they all arrive at the same prediction, whichof course increases their confidence level. As recognized for some time, the most directapproach is linear fluctuating hydrodynamics plus small nonlinear perturbations. Werefer to the recent monograph [7] for a comprehensive discussion. Here I provide only arough sketch of the method with the purpose to explain why one dimension is so special.In the physical dimension a simple fluid has five conservation laws and correspondinglyfluctuating hydrodynamics has to deal with a five component field, where the momentumcomponents are odd, density and energy are even under time reversal. As well known[8, 7], the full structure has to be used in order to arrive at quantitative predictions. Butthe argument becomes even more direct for the (unphysical) case of a single conservationlaw.Let us thus consider the scalar field ρ ( x , t ), which for concreteness is called density.Space is x ∈ R d and time is t ∈ R . ρ ( x , t ) is a fluctuating field. On the macroscopic scalefluctuations are not visible and ρ satisfies the hyperbolic conservation law ∂ t ρ + ∇ · (cid:126)j ( ρ ) = 0 , (1.1)where (cid:126)j is the density current. To also include mesoscopic details, in particular to in-corporate fluctuations, one argues that the current has, in addition to the deterministicpart, also a random contribution which is essentially uncorrelated in space-time. Sincefluctuations are always associated with dissipation, on a more refined scale Eq. (1.1)becomes ∂ t ρ + ∇ · (cid:0) (cid:126)j ( ρ ) − D ∇ ρ + σ(cid:126)ξ (cid:1) = 0 . (1.2)3ere (cid:126)ξ ( x , t ) is Gaussian white noise with mean zero and covariance (cid:104) ξ α ( x , t ) ξ α (cid:48) ( x (cid:48) , t (cid:48) ) (cid:105) = δ αα (cid:48) δ ( x − x (cid:48) ) δ ( t − t (cid:48) ) , α, α (cid:48) = 1 , ..., d . (1.3) σ is the noise strength and D is the diffusion constant. They are both treated as numbers.Physically they will depend on the density. But this would be higher order effects, whichare ignored in our discussion.The goal is to compute the density time correlations in the stationary regime, whichis no easy task, since (1.2) is nonlinear. But correlations can be thought of as imposingat t = 0 a small density perturbation at the origin and then record how the perturbationpropagates in space-time. For this purpose one might hope to get away with linearizing(1.2) at the uniform background density ρ as ρ + (cid:37) ( x , t ), which yields ∂ t (cid:37) + ∇ · (cid:0) (cid:126)j (cid:48) ( ρ ) (cid:37) − D ∇ (cid:37) + σ(cid:126)ξ (cid:1) = 0 , (1.4)where (cid:48) refers to differentiation w.r.t ρ . (1.4) is a linear Langevin equation, hence solvedeasily. Since (cid:37) is a deviation from the uniform background, we are interested in the space-time stationary process with zero mean. First note that (1.4) has a unique time-stationaryzero mean measure, which is Gaussian white noise in the spatial variable, (cid:104) (cid:37) ( x ) (cid:37) ( x (cid:48) ) (cid:105) = ( σ / D ) δ ( x − x (cid:48) ) . (1.5)For the stationary space-time covariance one obtains (cid:104) (cid:37) ( x , t ) (cid:37) (0 , (cid:105) = ( σ / D ) p ( x − (cid:126)c t, t ) , (1.6)where p is the Gaussian transition kernel, p ( x , t ) = (4 πD | t | ) − d/ exp (cid:0) − x / D | t | (cid:1) , (1.7)and (cid:126)c = (cid:126)j (cid:48) ( ρ ). By (1.4) the fluctuating current is given by (cid:126) J = (cid:126)c(cid:37) − D ∇ (cid:37) + σ(cid:126)ξ . (1.8)For the stationary correlations of the total current one arrives at (cid:90) R d dx (cid:104) (cid:126) J ( x , t ) · (cid:126) J (0 , (cid:105) = dσ δ ( t ) . (1.9)No surprise. A density fluctuation propagates with velocity (cid:126)c and spreads diffusively.The currents are delta-correlated, which should translate into exponential decay for theunderlying microscopic system. But before jumping at such conclusions one has to studythe stability of (1.9) against including higher orders in the expansion. By power countingthe next to leading term is the nonlinear current at second order, which amounts to ∂ t (cid:37) + ∇ · (cid:0) (cid:126)c(cid:37) + (cid:126)G(cid:37) − D ∇ (cid:37) + σ(cid:126)ξ (cid:1) = 0 , (cid:126)G = (cid:126)j (cid:48)(cid:48) ( ρ ) . (1.10)The task is to compute the current correlation (1.9) perturbatively in (cid:126)G .4e first remove (cid:126)c by switching to a moving frame of reference. Secondly we notethat, quite surprisingly, white noise is still time-stationary under the nonlinear Langevinequation (1.10). Since, as argued above, spatial white noise is time-stationary underthe linear part, one only has to check its invariance under the deterministic evolution ∂ t (cid:37) = − (cid:126)G · ∇ (cid:37) . Formally its vector field is divergence free, since one can choose adivergence free discretization [9], compare with the discussion in Section 3 below. Hencethe “volume measure” is preserved and one has to check only the time change of thelogarithm of the stationary density, ddt (cid:90) R d dx (cid:37) ( x ) = − (cid:90) R d dx(cid:37) ( x ) (cid:126)G · ∇ (cid:37) ( x ) = (cid:90) R d dx (cid:126)G · ∇ (cid:37) ( x ) = 0 . (1.11)For an expansion in (cid:126)G · ∇ (cid:37) it is most convenient to use the Fokker-Planck generator,denoted by L = L + L , where L corresponds to the Gaussian part and L to thenonlinear flow. We define S ( x , t ) = (cid:104) (cid:37) ( x , t ) (cid:37) (0 , (cid:105) = (cid:104) (cid:37) ( x )e Lt (cid:37) (0) (cid:105) , (1.12)average with respect to the time-stationary Gaussian measure, and plan to use the generalsecond moment sum rule d dt (cid:90) R d dx x S ( x , t ) = (cid:90) R d dx (cid:104) (cid:126) J ( x , t ) · (cid:126) J (0 , (cid:105) , (1.13)which follows directly from the conservation law. By second order time-dependent per-turbation theory, S ( x , t ) = (cid:104) (cid:37) ( x ) e L t (cid:37) (0) (cid:105) + (cid:90) t dt (cid:104) (cid:37) ( x )e L ( t − t ) L e L t (cid:37) (0) (cid:105) + (cid:90) t dt (cid:90) t dt (cid:104) (cid:37) ( x )e L ( t − t ) L e L ( t − t ) L e L t (cid:37) (0) (cid:105) + O ( | (cid:126)G | ) . (1.14) L (cid:37) is linear, while L (cid:37) is quadratic in (cid:37) . Thus the second term on the right is odd in (cid:37) and vanishes. For any functional, F , (cid:104) (cid:37) ( x ) L F ( (cid:37) ) (cid:105) = −(cid:104) ( L (cid:37) ( x )) F ( (cid:37) ) (cid:105) , (1.15)see (4.8), and the left L is swapped over to act on e L ( t − t ) (cid:37) ( x ). By translation invariance, (cid:126)G · ∇ can be pulled in front as acting on x . Combining the terms one arrives at S ( x , t ) = (cid:104) (cid:37) ( x )e L t (cid:37) (0) (cid:105) + (cid:90) t dt (cid:90) t dt ( (cid:126)G · ∇ ) × (cid:90) R d dx (cid:90) R d dx p ( x − x , t − t ) p ( x , t ) (cid:104) (cid:37) ( x ) e L ( t − t ) (cid:37) ( x ) (cid:105) . (1.16)The Gaussian average is computed as (cid:104) (cid:37) ( x ) e L t (cid:37) ( x ) (cid:105) = 2 (cid:0) ( σ / D ) p ( x − x , t ) (cid:1) + s-c . (1.17)5he self-contraction does not depend on x , hence vanishes when applying (cid:126)G · ∇ . Weinsert in (1.13). Working out the integrals yields, including second order, (cid:90) R d dx (cid:104) (cid:126) J ( x , t ) · (cid:126) J (0 , (cid:105) = dσ δ ( t ) + 4( σ / D ) | (cid:126)G | p (0 , t ) , (1.18)which is the claimed long time tail for a scalar conserved field in d dimensions.Eq. (1.10) is singular at short distances, the worse the higher the dimension. Inphysical systems there is a natural cut-off at the microscopic scale which is simply ignoredin (1.10). There are many possibilities to improve, but the serious constraint is to obtaina still manageable nonlinear stochastic equation. The noise should remain δ -correlatedin time so to preserve the Markov property of the time evolution. One could smoothenin space, but thereby may loose the information on the time-stationary measure. To myexperience the best compromise is to adopt the obvious spatial discretization. Then onehas a set of coupled stochastic differential equations. For them one can check rigorouslythe time-stationary measure and identities as e.g. (1.15). On the perturbative levelthis then leads to the continuum expressions as (1.16). For a simple fluid the currentcorrelations are bounded and will not diverge near t = 0. So in (1.18) only the long timeprediction can be trusted.What can be learned from the long time tails? In dimension d ≥ d = 2 is marginal. The diffusivity is weaklydivergent. In principle, say, a system of hard disks has infinite viscosity and thermalconductivity. But since the divergence is only logarithmic one could convert the resultinto a weak system size dependence. In one dimension the conductivity is truly infinite.Obviously, not even the power law as based on the perturbative expansion (1.16) can betrusted and one needs to develop non-perturbative techniques.Eq. (1.10) for d = 1 is known as stochastic Burgers equation, which we record forlater reference as ∂ t (cid:37) − ∂ x (cid:0) λ(cid:37) + ν∂ x (cid:37) + σξ (cid:1) = 0 , x ∈ R . (1.19)One can introduce a potential through (cid:37) = ∂ x h . Then h satisfies the one-dimensionalversion of the Kadar-Parisi-Zhang (KPZ) equation [10], ∂ t h = λ ( ∂ x h ) + ν∂ x h + σξ , (1.20)for the moment using the more conventional symbols for the coefficients. Over the lastfifteen years many properties, including rigorous results, of the KPZ equation have beenobtained. While this is not the place to dive into details, I note that the solution iscontinuous in x , but so singular that ( ∂ x h ) is ill-defined. However it is proved that theultraviolet divergence is very mild and can be tamed by what would be an infinite energyrenormalization in a 1+1 dimensional quantum field theory, compare with the discussionin [11]. More precisely one chooses a regularized version of (1.20) by replacing ξ ( x, t )through the spatially smoothed version ξ ϕ ( x, t ) = (cid:82) ϕ ( x − x (cid:48) ) ξ ( x (cid:48) , t ) dx (cid:48) with ϕ ≥
0, even,smooth, of rapid decay at infinity, and normalized to 1. It can be proved that then the6olution to (1.20) is well-defined. One introduces an ultraviolet cut-off of spatial size (cid:15) bychoosing the δ -function sequence ϕ (cid:15) ( x ) = (cid:15) − ϕ ( (cid:15) − x ) and substituting the noise ξ by itssmoothed version ξ (cid:15) = ξ ϕ (cid:15) . Let us denote the corresponding solution of the KPZ equationby h (cid:15) ( x, t ). Then h (cid:15) ( x, t ) − v (cid:15) t , i.e. h (cid:15) ( x, t ) viewed in the frame moving with the divergingvelocity v (cid:15) = (cid:15) − (cid:82) ϕ ( x ) dx , has a non-degenerate limit as (cid:15) →
0, which coincides withthe Cole-Hopf solution of the KPZ equation [12, 13].In this context Hairer recently developed a solution theory, for which he was awardedthe 2014 Fields Medal. His theory works for the KPZ equation, as well as a large classof other singular stochastic partial differential equations, and for general cutoffs [14, 15].Of course, the solution theory studies the small scale structure of solutions, and not thelarge scale, where the universal behaviors of interest are observed.There is an exact formula for S ( x, t ), which involves Fredholm determinants. In thelong time limit S ( x, t ) (cid:39) ( σ / ν )( √ λ | t | ) − / f KPZ (cid:0) √ λ | t | ) − / x (cid:1) . (1.21)The scaling function f KPZ will reappear below, where more details are given. The Burgerscurrent reads J ( x, t ) = − ( λ(cid:37) + ν∂ x (cid:37) + σξ (cid:1) (1.22)and for its total correlation function one obtains, using the sum rule (1.13), (cid:90) R dx (cid:104)J ( x, t ) J (0 , (cid:105) (cid:39) (cid:16) ( σ / ν )( √ λ ) (cid:90) R dx x f KPZ ( x ) (cid:17) ( √ λt ) − / (1.23)valid for large t . Note that the true decay turns out to be t − / , to be contrasted withthe perturbative result t − / . In fact the power 2 / f KPZ and the pure number √ cknowledgements . First of all I owe my thanks to Christian Mendl. Only throughhis constant interest, his help in checking details, and his superb numerical efforts myproject could advance to its current level. Various collaborations developed and I ammost grateful to my coauthors A. Dhar, P. Ferrari, D. Huse, M. Kulkarni, T. Sasamoto,and G. Stoltz. At various stages of the project, I greatly benefited from discussions withH. van Beijeren, C. Bernardin, J. Krug, J.L. Lebowitz, S. Lepri, R. Livi, S. Olla, H. Posch,and G. Sch¨utz. Conservation laws and equilibrium time correlations . We consider a classical fluidconsisting of particles with positions q j and momenta p j , j = 1 , ..., N , q j , p j ∈ R , possibleboundary conditions to be discussed later on. We use units such that the mass of theparticles equals 1. Then the Hamiltonian is of the standard form, H fl N = N (cid:88) j =1 12 p j + N (cid:88) i (cid:54) = j =1 V ( q i − q j ) , (2.1)with pair potential V ( x ) = V ( − x ). The potential may have a hard core and otherwiseis assumed to be short ranged. The dynamics for long range potentials is of independentinterest [24], but not discussed here. Furthermore the potential is assumed to be thermo-dynamically stable, meaning the validity of a bound as (cid:80) Ni (cid:54) = j =1 V ( q i − q j ) ≥ A − BN forsome constants A and B >
0. A substantial simplification is achieved by assuming a hardcore of diameter a , i.e. V ( x ) = ∞ for | x | < a , and restricting the range of the smoothpart of the potential to at most 2 a . Then the particles maintain their order, q j ≤ q j +1 ,and in addition only nearest neighbor particles interact. Hence H fl N simplifies to H N = N (cid:88) j =1 12 p j + N − (cid:88) j =1 V ( q j +1 − q j ) . (2.2)As a, at first sight very different, physical realization, we could interpret H N as describingparticles in one dimension coupled through anharmonic springs which is then usuallyreferred to as anharmonic chain.In the second interpretation the spring potential can be more general than anticipatedso far. No ordering constraint is required and the potential does not have to be even.To have well defined thermodynamics the chain is pinned at both ends as q = 0 and q N +1 = (cid:96)N . It is convenient to introduce the stretch r j = q j +1 − q j . Then the boundarycondition corresponds to the microcanonical constraint N (cid:88) j =1 r j = (cid:96)N . (2.3)8witching to canonical equilibrium according to the standard rules, one the arrives at theobvious condition of a finite partition function Z ( P, β ) = (cid:90) dx e − β ( V ( x )+ P x ) < ∞ , (2.4)using the standard convention that the integral is over the entire real line. Here β > P is the thermodynamically conjugate variable to the stretch.By partial integration P = − Z ( P, β ) − (cid:90) dxV (cid:48) ( x ) e − β ( V ( x )+ P x ) , (2.5)implying that P is the average force in the spring between two adjacent particles, henceidentified as thermodynamic pressure. To have a finite partition function, a naturalcondition on the potential is to be bounded from below and to have a one-sided linearbound as V ( x ) ≥ a + b | x | for either x > x < b >
0. Then there is anon-empty interval I ( β ) such that Z ( P, β ) < ∞ for P ∈ I ( β ). For the particular case ofa hard-core fluid one imposes P > Note:
The sign of P is chosen such that for a gas of hard-point particles one has thefamiliar ideal gas law P = 1 /β(cid:96) . The chain tension is − P .Famous examples are the harmonic chain, V ha ( x ) = x , the Fermi-Pasta-Ulam (FPU)chain, V FPU ( x ) = x + αx + βx , in the historical notation [25], and the Toda chain[26], V ( x ) = e − x , in which case P > V hc ( x ) = ∞ for | x | < a and V hc ( x ) = 0 for | x | ≥ a , are infact integrable systems which have a very different correlation structure and will not bediscussed here. Except for the harmonic chain, one simple way to break integrability isto assume alternating masses, say m j = m for even j and m j = m for odd j .We will mostly deal with anharmonic chains described by the Hamiltonian (2.2), in-cluding one-dimensional hard-core fluids with a sufficiently small potential range. Thereare several good reasons. Firstly on the level of fluctuating hydrodynamics a genericone-dimensional fluid cannot be distinguished from an anharmonic chain. Thus withthe proper translation of the various terms we would also predict the large scale corre-lation structure of one-dimensional fluids. The second reason is that in the large bodyof molecular dynamics simulations there is not a single one which deals with an “hon-est” one-dimensional fluid. To be able to reach large system sizes all simulations areperformed for anharmonic chains. In addition, from a theoretical perspective, the equilib-rium measures of anharmonic chains are particularly simple in being of product form inmomentum and stretch variables. Thus material parameters, as compressibility and soundspeed, can be expressed in terms of one-dimensional integrals involving the Boltzmannfactor e − β ( V ( x )+ P x ) , V ( x ), and x .Anharmonic chains should be thought of as a particular class of 1+1 dimensionalfield theories. Thus q j is viewed as the displacement variable at lattice site j and notnecessarily the physical position of the j -th particle on the real line. There is a simpletranslation between both pictures, but we will stick to the field theory point of view. For9 fluid with unlabeled particles and, say, a bounded potential, the equivalence is lost andonly the fluid picture can be used.This being said, we follow the standard rules. We write down the dynamics of stretchesand momenta and identify the conserved fields. From there we infer the microcanonicaland canonical equilibrium measures. For slowly varying equilibrium parameters we de-duce the Euler equations. In particular, their version linearized at uniform equilibriumwill constitute the backbone in understanding the equilibrium time correlations of theconserved fields.The dynamics of the anharmonic chain is governed by ddt q j = p j , ddt p j = V (cid:48) ( q j +1 − q j ) − V (cid:48) ( q j − q j − ) . (2.6)For the initial conditions we choose a lattice cell of length N and require q j + N = q j + (cid:96)N , p j + N = p j (2.7)for all j ∈ Z . This property is preserved under the dynamics and thus properly mimics asystem of finite length N . The stretches are then N -periodic, r j + N = r j , and the singlecell dynamics is given by ddt r j = p j +1 − p j , ddt p j = V (cid:48) ( r j ) − V (cid:48) ( r j − ) , (2.8) j = 1 , . . . , N , together with the periodic boundary conditions p N = p , r = r N andthe constraint (2.3). Through the stretch there is a coupling to the right neighbor andthrough the momentum a coupling to the left neighbor. The potential is defined only upto translations, since the dynamics does not change under a simultaneous shift of V ( x )to V ( x − a ) and r j to r j + a , in other words, the potential can be shifted by shiftingthe initial r -field. Note that our periodic boundary conditions are not identical to fluidparticles moving on a ring, but they may become so for large system size when lengthfluctuations become negligible. Both equations are already of conservation type and weconclude that ddt N (cid:88) j =1 r j = 0 , ddt N (cid:88) j =1 p j = 0 . (2.9)We define the local energy by e j = p j + V ( r j ) . (2.10)Then its local conservation law reads ddt e j = p j +1 V (cid:48) ( r j ) − p j V (cid:48) ( r j − ) , (2.11)implying that ddt N (cid:88) j =1 e j = 0 . (2.12)10t this point we assume that there are no further local conservation laws. Unfortu-nately our assumption, while reasonable, is extremely difficult to check. It certainly rulesout the integrable chains, which have N conservation laws. There is no natural exampleknown with, say, seven conservation laws. But the mere fact that there are exceptionsimplies that close to integrability the predictions from fluctuating hydrodynamics couldbe on time scales which are not accessible. The parameters entering in fluctuating hy-drodynamics depend smoothly on the potential. Thus as one approaches, for example,the Toda potential no abrupt changes will be detected. In this sense, the theory cannotdistinguish the Toda chain from a FPU chain both at moderate temperatures. There isanother limitation which can be addressed more quantitatively. If V ( x ) + P x has a uniqueminimum, then at very low temperatures the potential is close to a harmonic one. Thisfeature will be properly reflected by fluctuating hydrodynamics through the temperaturedependence of the coupling coefficients and the noise strength.The microcanonical equilibrium state is defined by the Lebesgue measure constrainedto a particular value of the conserved fields as N (cid:88) j =1 r j = (cid:96)N , N (cid:88) j =1 p j = u N , N (cid:88) j =1 (cid:0) p j + V ( r j ) (cid:1) = e N (2.13)with (cid:96) the stretch, u the momentum, and e the total energy per particle. In our contextthe equivalence of ensembles holds and computationally it is of advantage to switch to thecanonical ensemble with respect to all three constraints. Then the dual variable for thestretch (cid:96) is the pressure P , for the momentum the average momentum, again denoted by u , and for the total energy e the inverse temperature β . For the limit of infinite volumethe symmetric choice j ∈ [ − N, ..., N ] is more convenient. In the limit N → ∞ eitherunder the canonical equilibrium state, trivially, or under the microcanonical ensemble, bythe equivalence of ensembles, the collection ( r j , p j ) j ∈ Z are independent random variables.Their single site probability density is given by Z ( P, β ) − e − β ( V ( r j )+ P r j ) (2 π/β ) − / e − β ( p j − u ) . (2.14)Averages with respect to (2.14) are denoted by (cid:104)·(cid:105) P,β, u . The dependence on the averagemomentum can be removed by a Galilei transformation. Hence we mostly work with u = 0, in which case we merely drop the index u . We also introduce the internal energy, e , through e = u + e , which agrees with the total energy at u = 0. The canonical freeenergy, at u = 0, is defined by G ( P, β ) = − β − (cid:0) − log β + log Z ( P, β ) (cid:1) . (2.15)Then (cid:96) = (cid:104) r (cid:105) P,β , e = ∂ β (cid:0) βG ( P, β ) (cid:1) − P (cid:96) = 12 β + (cid:104) V ( r ) (cid:105) P,β . (2.16)The relation (2.16) defines ( P, β ) (cid:55)→ ( (cid:96) ( P, β ) , e ( P, β )), thereby the inverse map ( (cid:96), e ) (cid:55)→ ( P ( (cid:96), e ) , β ( (cid:96), e )), and thus accomplishes the switch between the microcanonical thermo-dynamic variables (cid:96), e and the canonical thermodynamic variables P, β .11t is convenient to collect the conserved fields as the 3-vector (cid:126)g = ( g , g , g ), (cid:126)g ( j, t ) = (cid:0) r j ( t ) , p j ( t ) , e j ( t ) (cid:1) , (2.17) (cid:126)g ( j,
0) = (cid:126)g ( j ). Then the conservation laws are combined as ddt(cid:126)g ( j, t ) + (cid:126) J ( j + 1 , t ) − (cid:126) J ( j, t ) = 0 (2.18)with the local current functions (cid:126) J ( j ) = (cid:0) − p j , − V (cid:48) ( r j − ) , − p j V (cid:48) ( r j − ) (cid:1) . (2.19)Our prime interest are the equilibrium time correlations of the conserved fields, which aredefined by S αα (cid:48) ( j, t ) = (cid:104) g α ( j, t ) g α (cid:48) (0 , (cid:105) P,β − (cid:104) g α (0 , (cid:105) P,β (cid:104) g α (cid:48) (0 , (cid:105) P,β , (2.20) α, α (cid:48) = 1 , ,
3. The infinite volume limit has been taken already and the average is withrespect to thermal equilibrium at u = 0. It is known that such a limit exists [28]. Also thedecay in j is exponentially fast, but with a correlation length increasing in time. Often itis convenient to regard S ( j, t ), no indices, as a 3 × S ( j, t ) has certainsymmetries, the first set resulting from space-time stationarity and the second set fromtime reversal, even for α = 1 ,
3, odd for α = 2, S αα (cid:48) ( j, t ) = S α (cid:48) α ( − j, − t ) , S αα (cid:48) ( j, t ) = ( − α + α (cid:48) S αα (cid:48) ( j, − t ) . (2.21)At t = 0 the average (2.20) reduces to a static average, which is easily computed. Thefields are uncorrelated in j , i.e. S ( j,
0) = δ j C (2.22)with the static susceptibility matrix C = (cid:104) r ; r (cid:105) P,β (cid:104) r ; V (cid:105) P,β β − (cid:104) r ; V (cid:105) P,β β − + (cid:104) V ; V (cid:105) P,β . (2.23)Here, for X , Y arbitrary random variables, (cid:104) X ; Y (cid:105) = (cid:104) XY (cid:105) − (cid:104) X (cid:105)(cid:104) Y (cid:105) denotes the secondcumulant and V = V ( r ), following the same notational convention as for e . Theconservation law implies the zeroth moment sum rule (cid:88) j ∈ Z S αα (cid:48) ( j, t ) = (cid:88) j ∈ Z S αα (cid:48) ( j,
0) = C αα (cid:48) . (2.24)An explicit computation of S ( j, t ) is utterly out of reach. But with the current com-puter power MD simulations have become an essential source of information. A broadercoverage will be provided in Section 5. Just to have a first impression I show in Fig. 1a recent MD simulation of the correlator for a FPU chain. One notes the central peak,12alled heat peak, which is standing still, and two symmetrically located peaks, calledsound peaks which move outwards with the speed of sound. The peaks broaden with acertain power law which will have to be discussed. One expects, better hopes for, self-similar shape functions, at least for sufficiently long times. We do not know yet theirform. But the central peak seems to have fat tails while the sound peaks fall off morerapidly, at least towards the outside of the sound cone. The area under each peak ispreserved in time and normalized to 1 in our plot. If the chain is initially perturbed near0 and the response in one of the conserved fields is observed at ( j, t ), then one records asignal, which is a linear combination of the peaks in Fig. 1, the coefficients depending onthe initial perturbation. It might happen that one peak is missing. If the perturbationis orthogonal to all three physical fields, then there is no peak at all, only low amplitudenoise. -4000 -2000 0 2000 4000x00.0010.0020.0030.004 S ( x , t ) t=800t=1300t=2700 Figure 1: Heat peak and sound peaks, area normalized to 1, at times t = 800 , , α = 2, β = 1, pressure P = 1, and temperature β − = 0 . S αα (cid:48) ( j, t ) encoding thepropagation of local perturbations of the equilibrium state. On the crudest level, theyshould be captured by linearized hydrodynamics, to which we turn next. Linearized hydrodynamics . We start the dynamics from a product measure of the form(2.14), but replace the uniform P, u , β by slowly varying spatial fields P ( (cid:15)j ) , u ( (cid:15)j ) , β ( (cid:15)j )with (cid:15) (cid:28) (cid:15) − is the macroscopic scale measured in lattice units. Equivalently we mayregard (cid:15) as the lattice spacing. The relation (2.16) then defines also the slowly varyingfields (cid:96) ( (cid:15)j ) , u ( (cid:15)j ) , e ( (cid:15)j ). Because of the conservation laws the time change of such a state isslow and varies only over microscopic times of order (cid:15) − t with macroscopic t of order 1. Weaverage Eq. (2.18) over the slowly varying initial state, which gives then the time changeof the average locally conserved fields. The expectation of the currents is more difficult.13he difference in j becomes (cid:15)∂ x on the macroscopic scale. Since the currents are functionsdepending only on at most two neighboring lattice sites, to lowest approximation theiraverage can be computed in the equilibrium state with the corresponding local values ofthe fields, i.e. (cid:96) ( (cid:15)j, (cid:15) − t ) , u ( (cid:15)j, (cid:15) − t ) , e ( (cid:15)j, (cid:15) − t ). Therefore we define the hydrodynamicEuler currents by the equilibrium averages (cid:104) (cid:126) J ( j ) (cid:105) (cid:96), u , e = (cid:0) − u , P ( (cid:96), e − u ) , u P ( (cid:96), e − u ) (cid:1) = (cid:126) j ( (cid:96), u , e ) (2.25)with P ( (cid:96), e ) defined implicitly through (2.16). Our argument then leads to the macroscopicEuler equations ∂ t (cid:96) − ∂ x u = 0 , ∂ t u + ∂ x P ( (cid:96), e − u ) = 0 , ∂ t e + ∂ x (cid:0) u P ( (cid:96), e − u ) (cid:1) = 0 . (2.26)We refer to a forthcoming monograph [28], where the validity of the Euler equations isproved up to the first shock. Since, as emphasized already, it is difficult to deal withdeterministic chaos, the authors add random velocity exchanges between neighboringparticles which ensure that the dynamics locally enforces the microcanonical state.We are interested here only in small deviations from equilibrium and therefore linearizethe Euler equations as (cid:96) + u ( x ), 0 + u ( x ), e + u ( x ) to linear order in the deviations (cid:126)u ( x ). This leads to the linear equation ∂ t (cid:126)u ( x, t ) + ∂ x A(cid:126)u ( x, t ) = 0 (2.27)with A = − ∂ (cid:96) P ∂ e P P . (2.28)Here, and in the following, the dependence of A , C and similar quantities on the back-ground values (cid:96), u = 0 , e , hence on P, β , is suppressed from the notation. Beyond (2.24)there is the first moment sum rule which states that (cid:88) j ∈ Z jS ( j, t ) = AC t . (2.29)A proof, which in essence uses only the conservation laws and space-time stationarity ofthe correlations, is given in [29], see also see [30, 31]. Microscopic properties enter onlyminimally. However, since C = C T and S ( j, t ) T = S ( − j, − t ), Eq. (2.29) implies theimportant relation AC = ( AC ) T = CA T , (2.30)with T denoting transpose. Of course, (2.30) can be checked also directly from the defini-tions. Since C > A is guaranteed to have real eigenvalues and a nondegenerate systemof right and left eigenvectors. For A one obtains the three eigenvalues 0 , ± c with c = − ∂ (cid:96) P + P ∂ e P > . (2.31)14hus the solution to the linearized equation has three modes, one standing still, one rightmoving with velocity c and one left moving with velocity − c . Hence we have identifiedthe adiabatic sound speed as being equal to c .(2.27) is a deterministic equation. But the initial data are random such that withinour approximation (cid:104) u α ( x, u α (cid:48) ( x (cid:48) , (cid:105) = C αα (cid:48) δ ( x − x (cid:48) ) . (2.32)To determine the correlator S ( x, t ) with such initial conditions is most easily achieved byintroducing the linear transformation R satisfying RAR − = diag( − c, , c ) , RCR T = 1 . (2.33)Up to trivial phase factors, R is uniquely determined by these conditions. Explicit for-mulas are found in [29]. Setting (cid:126)φ = A(cid:126)u , one concludes ∂φ α + c α ∂ x φ α = 0 , α = − , , , (2.34)with (cid:126)c = ( − c, , c ). By construction, the random initial data have the correlator (cid:104) φ α ( x, φ α (cid:48) ( x (cid:48) , (cid:105) = δ αα (cid:48) δ ( x − x (cid:48) ) . (2.35)Hence (cid:104) φ α ( x, t ) φ α (cid:48) (0 , (cid:105) = δ αα (cid:48) δ ( x − c α t ) . (2.36)We transform back to the physical fields. Then in the continuum approximation, at thelinearized level, S ( x, t ) = R − diag (cid:0) δ ( x + ct ) , δ ( x ) , δ ( x − ct ) (cid:1) R − T (2.37)with R − T = ( R − ) T .Rather easily we have gained a crucial insight. S ( j, t ) has three peaks which separatelinearly in time. For example, S ( j, t ) has three sharp peaks moving with velocities ± c, R − , usually called Landau-Plazcek ratios. A Landau-Placzek ratio could vanish,either accidentally or by a particular symmetry. An example is the momentum correlation S ( j, t ). Since ( R − ) = 0 always, its central peak is absent.For integrable chains each conservation law generates a peak. Thus, e.g. , S ( j, t ) ofthe Toda chain is expected to have a broad spectrum expanding ballistically, rather thanconsisting of three sharp peaks. Euler currents to second order . The broadening of the peaks results from randomfluctuations in the currents, which tend to be uncorrelated in space-time. Therefore the15rudest model would be to assume that the current statistics is space-time Gaussianwhite noise. In principle, the noise components could be correlated. But since the stretchcurrent is itself conserved, its fluctuations will be taken care of by the momentum equation.Momentum and energy currents have different signature under time reversal, hence theircross correlation vanishes. As a result, there is a fluctuating momentum current of strength σ u and an independent energy current of strength σ e . According to Onsager, noise is linkedto dissipation as modeled by a diffusive term. Thus the linearized equations (2.27) areextended to ∂ t (cid:126) u ( x, t ) + ∂ x (cid:0) A(cid:126) u ( x, t ) − ∂ x D(cid:126) u ( x, t ) + B(cid:126)ξ ( x, t ) (cid:1) = 0 . (3.1)Here (cid:126)ξ ( x, t ) is standard white noise with covariance (cid:104) ξ α ( x, t ) ξ α (cid:48) ( x (cid:48) , t (cid:48) ) (cid:105) = δ αα (cid:48) δ ( x − x (cid:48) ) δ ( t − t (cid:48) ) (3.2)and, as argued, the noise strength matrix is diagonal as B = diag(0 , σ u , σ e ) . (3.3)To distinguish the linearized Euler equations (2.27) from the Langevin equations (3.1),we use (cid:126) u = ( u , u , u ) for the fluctuating fields.From the introduction, we know already that a Gaussian fluctuation theory will fail.But still, it is useful to first explore the structure of the Langevin equation (3.1). Thestationary measures for (3.1) are spatial white noise with arbitrary mean. Since smalldeviations from uniformity are considered, we always impose mean zero. Then the com-ponents are correlated as (cid:104) u α ( x ) u α (cid:48) ( x (cid:48) ) (cid:105) = C αα (cid:48) δ ( x − x (cid:48) ) . (3.4)Stationarity relates the linear drift and the noise strength through the steady state co-variance as − ( AC − CA T ) ∂ x + ( DC + CD T ) ∂ x = BB T ∂ x . (3.5)The first term vanishes by (2.30) and the diffusion matrix is uniquely determined as D = D u D e D e . (3.6)with ˜ D e = −(cid:104) r ; V (cid:105) P,β (cid:104) r ; r (cid:105) − P,β D e . Here D u > D e > σ u = (cid:104) p ; p (cid:105) P,β D u , σ e = (cid:104) e ; e (cid:105) P,β D e . (3.7)We still have to establish that the stationary measure (3.4) is unique and is approachedin the limit t → ∞ . For this purpose it suffices that the 3 × πkA − (2 πk ) D hasits eigenvalues in the left hand complex plane, where for convenience we have switched16o Fourier space with respect to x . If one drops D , then i2 πkA has the eigenval-ues i2 πkc ( − , , − (2 πk ) D , which is given by (cid:104) ˜ ψ α , Dψ α (cid:105) , where Aψ α = c α ψ α and A T ˜ ψ α = c α ˜ ψ α are theright and left eigenvectors of A , as listed in [29]. One simply has to follow the definitionsand express every term through the cumulants of r and V . As to be expected, the ma-trix elements from above have a definite value and the eigenvalues are shifted into the lefthand complex plane. Similarly, − D has eigenvalues 0 , − D u , − D e and the zero eigenvalueis shifted to the left by second order perturbation in i(2 πk ) − A . The only condition isthe strict positivity of D u , D e .Based on (3.1) one computes the stationary space-time covariance, which most easilyis written in Fourier space, S αα (cid:48) ( x, t ) = (cid:104) u α ( x, t ) u α (cid:48) (0 , (cid:105) = (cid:90) dk e i2 πkx (cid:0) e − i t πkA −| t | (2 πk ) D C (cid:1) αα (cid:48) . (3.8)To extract the long time behavior it is convenient to transform to normal modes. Butbefore, we have to introduce a more systematic notation. We will use the superscript (cid:93) for a normal mode quantity. Thus for the anharmonic chain S (cid:93) ( j, t ) = RS ( j, t ) R T , S (cid:93)αα (cid:48) ( j, t ) = (cid:104) ( R(cid:126)g ) α ( j, t )( R(cid:126)g ) α (cid:48) (0 , (cid:105) P,β . (3.9)The hydrodynamic fluctuation fields are defined on the continuum, thus functions of x, t ,and we write S αα (cid:48) ( x, t ) = (cid:104) u α ( x, t ) u α (cid:48) (0 , (cid:105) , S (cid:93) ( x, t ) = RS ( x, t ) R T . (3.10)Correspondingly A (cid:93) = RAR − = diag( − c, , c ), D (cid:93) = RDR − , B (cid:93) = RB . Note that (cid:126) u ( x, t ) will change its meaning when switching from linear to nonlinear fluctuating hydro-dynamics.In normal mode representation Eq. (3.8) becomes S (cid:93) ( x, t ) = (cid:90) dk e i2 πkx e − i t πkA (cid:93) −| t | (2 πk ) D (cid:93) . (3.11)The leading term, i t πkA (cid:93) , is diagonal, while the diffusion matrix D (cid:93) couples the compo-nents. But for large t the peaks are far apart and the cross terms become small. Moreformally we split D (cid:93) = D dia + D off and regard the off-diagonal part D off as perturbation.When expanding, one notes that the off-diagonal terms carry an oscillating factor withfrequency c α − c α (cid:48) , α (cid:54) = α (cid:48) . Hence these terms decay quickly and S (cid:93)αα (cid:48) ( x, t ) (cid:39) δ αα (cid:48) (cid:90) dk e i2 πkx e − i t πkc α −| t | (2 πk ) D (cid:93)αα (3.12)for large t . Each peak has a Gaussian shape function which broadens as ( D (cid:93)αα | t | ) / .Besides the peak structure, we have gained a second important insight. Since thepeaks travel with distinct velocities, on the linearized level the three-component systemdecouples into three scalar equations, provided it is written in normal modes.17he linear fluctuation theory should be tested against adding nonlinear terms. Thecomputation of the Introduction indicates that in one dimension the quadratic part ofthe Euler current is a relevant perturbation. This can be seen even more directly byrescaling (cid:126) u ( x, t ) to large space-time scales as (cid:126) u (cid:15) ( x, t ) = (cid:15) − b (cid:126) u (cid:15) ( (cid:15) − x, (cid:15) − z t ) and counting thepowers of the nonlinear terms. To have the correct t = 0 covariance requires b = 1 / z = 3 /
2, and the scalingexponents of cubic terms are only marginally relevant, while a possible dependence of
D, B on (cid:126) u ( x, t ) can be ignored. Thus, we retain the linear part (3.1) but expand the Eulercurrents including second order in (cid:126)u , which turns (3.1) into the equations of nonlinearfluctuating hydrodynamics, ∂ t u − ∂ x u = 0 ,∂ t u + ∂ x (cid:0) ( ∂ (cid:96) P ) u + ( ∂ e P ) u + ( ∂ (cid:96) P ) u − ( ∂ e P ) u + ( ∂ e P ) u + ( ∂ (cid:96) ∂ e P ) u u − D u ∂ x u + σ u ξ (cid:1) = 0 ,∂ t u + ∂ x (cid:0) P u + ( ∂ (cid:96) P ) u u + ( ∂ e P ) u u − ˜ D e ∂ x u − D e ∂ x u + σ e ξ (cid:1) = 0 . (3.13)To explore their consequences is a more demanding task than solving the linear Langevinequation and the results of the analysis will be more fragmentary. Stationary measure for the physical fields . Adding quadratic terms could changedrastically the character of the solution. To find out we first attempt to investigate thetime-stationary measure. The vector field of the nonlinear part of (3.13) reads (cid:126)F = − ∂ x (cid:0) , ( ∂ (cid:96) P ) u − ( ∂ e P ) u + ( ∂ e P ) u + ( ∂ (cid:96) ∂ e P ) u u , ( ∂ (cid:96) P ) u u + ( ∂ e P ) u u (cid:1) . (3.14)Formally the drift is divergence free, since (cid:88) α =1 (cid:90) dx δF α ( x ) δ u α ( x ) = ∂ e P (cid:90) dx ( − ∂ x u + ∂ x u ) = 0 (3.15)and the infinite dimensional Lebesgue measure is invariant under the flow generated by (cid:126)F . Since the equilibrium measure is a product, a natural ansatz for the invariant measureis Gaussian white noise retaining the physical susceptibility,exp (cid:16) − (cid:88) α,α (cid:48) =1 12 ( C − ) αα (cid:48) (cid:90) dx u α ( x ) u α (cid:48) ( x ) (cid:17) (cid:89) α,x d u α ( x ) . (3.16)As established before, for the linear Langevin equation this measure is stationary. Thus,to find out whether it is also stationary for (3.13), we only have to check the invarianceunder the nonlinear flow0 = ddt (cid:88) α,α (cid:48) =1 ( C − ) αα (cid:48) (cid:90) dx u α ( x ) u α (cid:48) ( x ) = 2 (cid:88) α =1 (cid:90) dx ( C − (cid:126) u ) α ( x ) ∂ t u α ( x )= (cid:90) dxa u ∂ x (cid:0) ( ∂ (cid:96) P ) u − ( ∂ e P ) u + ( ∂ e P ) u + 2( ∂ (cid:96) ∂ e P ) u u (cid:1) +2 (cid:90) dx ( a u + a u ) ∂ x (cid:0) ( ∂ (cid:96) P ) u u + ( ∂ e P ) u u (cid:1) , (3.17)18here a = ( C − ) , a = ( C − ) , a = ( C − ) . The term cubic in u vanishes by thesame argument as for the one-component case. All other terms are linear in u , thus theirsum has to vanish point-wise,0 = ( a ∂ (cid:96) P − a ∂ (cid:96) P ) ∂ x u + ( a ∂ e P − a ∂ e P ) ∂ x u +2 (cid:0) a ( ∂ (cid:96) ∂ e P ) ∂ x ( u u ) − a ( ∂ e P ) u ∂ x u − a ( ∂ (cid:96) P ) u ∂ x u (cid:1) . (3.18)This leads to the constraints on the coefficients as a ∂ (cid:96) P = a ∂ (cid:96) P , a ∂ e P = a ∂ e P , a ∂ (cid:96) ∂ e P = a ∂ (cid:96) P , a ∂ (cid:96) ∂ e P = a ∂ e P . (3.19)There are four constraints and five partial derivatives which may be regarded as inde-pendent parameters. Thus one would expect that there is a sub-manifold in
V, P, β spacefor which (3.19) can be satisfied. We will come back to another representation of theseconstraints below. Away from the special subset of invariant Gaussian measures, we haveno tools, but one would hope that the invariant measure still has a finite correlation lengthand exponential mixing. Based on the mechanical model, it is suggestive to assume u to be independent of u , u and to have white noise statistics. But this forces again theconstraints (3.19) and results in the same Gaussian measure as before.A basic property of the mechanical model is invariance under time reversal. On thelevel of fluctuating hydrodynamics this translates to the following property: We fix a timewindow [0 , T ]. Then, in the stationary process, the trajectories( u ( t ) , u ( t ) , u ( t )) and ( u ( T − t ) , − u ( T − t ) , u ( T − t )) (3.20)with 0 ≤ t ≤ T have the same probability. To check (3.20) requires some informationon the invariant measure. But under the Gaussian measure (3.16), hence assuming thevalidity of the constraints, it can be shown that time reversal invariance indeed holds. Transformation to normal modes . To proceed further, it is convenient to write (3.13)in vector form, ∂ t (cid:126) u ( x, t ) + ∂ x (cid:0) A(cid:126) u ( x, t ) + (cid:104) (cid:126) u , (cid:126)H(cid:126) u (cid:105) − ∂ x ˜ D(cid:126) u ( x, t ) + ˜ B(cid:126)ξ ( x, t ) (cid:1) = 0 , (3.21)where (cid:126)H is the vector consisting of the Hessians of the currents with derivatives evaluatedat the background values ( (cid:96), , e ), H αγγ (cid:48) = ∂ u γ ∂ u γ (cid:48) j α , (cid:104) (cid:126) u , (cid:126)H(cid:126) u (cid:105) = (cid:88) γ,γ (cid:48) =1 (cid:126)H γγ (cid:48) u γ u γ (cid:48) . (3.22)As for the linear Langevin equation we transform to normal modes through (cid:126)φ = R(cid:126) u . (3.23)Then ∂ t φ α + ∂ x (cid:0) c α φ α + (cid:104) (cid:126)φ, G α (cid:126)φ (cid:105) − ∂ x ( D (cid:93) (cid:126)φ ) α + ( B (cid:93) (cid:126)ξ ) α (cid:1) = 0 . (3.24)19y construction B (cid:93) B (cid:93) T = 2 D (cid:93) . The nonlinear coupling constants, denoted by (cid:126)G , aredefined by G α =
12 3 (cid:88) α (cid:48) =1 R αα (cid:48) R − T H α (cid:48) R − (3.25)with the notation R − T = ( R − ) T .Since derived from a chain, the couplings are not completely arbitrary, but satisfy thesymmetries G αβγ = G αγβ , G σαβ = − G − σ − α − β , G σ − = G σ ,G σσ = − G − σ − σ , G αβ = 0 otherwise . (3.26)In particular note that G = 0 , (3.27)always, while G = − G − − − are generically different from 0. This property signals thatthe heat peak will behave differently from the sound peaks. The (cid:126)G -couplings are listedin [29] and as a function of P, β expressed in cumulants up to third order in r , V . Thealgebra is somewhat messy. But there is a short MATHEMATICA program available [32]which, for given P, β, V , computes all coupling constants, including the matrices
C, A, R .We return to the issue of Gaussian time-stationary measures, where we regard thecoefficients (cid:126)G, D (cid:93) , B (cid:93) as arbitrary, up to 2 D (cid:93) = B (cid:93) B (cid:93) T , ignoring for a while their particularorigin. The Langevin equation (3.24) is slightly formal. To have a well-defined evolution,we discretize space by a lattice of N sites. The field (cid:126)φ ( x, t ) then becomes (cid:126)φ j ( t ) withcomponents φ j,α ( t ), j = 1 , . . . , N , α = 1 , ,
3. The spatial finite difference operator isdenoted by ∂ j , ∂ j f j = f j +1 − f j , with transpose ∂ T j f j = f j − − f j . Then the discretizedequations of fluctuating hydrodynamics read ∂ t φ j,α + ∂ j (cid:0) c α φ j,α + N j,α + ∂ T j D (cid:93) φ j,α + B (cid:93) ξ j,α (cid:1) = 0 (3.28)with (cid:126)φ j = (cid:126)φ N + j , (cid:126)ξ = (cid:126)ξ N , where ξ j,α are independent Gaussian white noises with covari-ance (cid:104) ξ j,α ( t ) ξ j (cid:48) ,α (cid:48) ( t (cid:48) ) (cid:105) = δ jj (cid:48) δ αα (cid:48) δ ( t − t (cid:48) ) . (3.29)The diffusion matrix D (cid:93) and noise strength B (cid:93) act on components, while the differenceoperator ∂ j acts on the lattice site index j . N j,α is quadratic in φ . But let us first consider the case N j,α = 0. Then φ j,α ( t ) is aGaussian process. The noise strength has been chosen such that one invariant measure isthe Gaussian N (cid:89) j =1 3 (cid:89) α =1 exp[ − φ j,α ](2 π ) − / dφ j,α = ρ G ( φ ) N (cid:89) j =1 3 (cid:89) α =1 dφ j,α . (3.30)Because of the conservation laws, the hyperplanes N (cid:88) j =1 φ j,α = N ρ α , (3.31)20re invariant and on each hyperplane there is a Gaussian process with a unique invariantmeasure given by (3.30) conditioned on that hyperplane. For large N it would becomeindependent Gaussians with mean ρ α , our interest being the case of zero mean, ρ α = 0.The generator of the diffusion process (3.28) with N j,α = 0 is given by L = N (cid:88) j =1 (cid:16) − (cid:88) α =1 ∂ j (cid:0) c α φ j,α + ∂ T j D (cid:93) φ j,α (cid:1) ∂ φ j,α + (cid:88) α,α (cid:48) =1 ( B (cid:93) B (cid:93) T ) αα (cid:48) ∂ j ∂ φ j,α ∂ j ∂ φ j,α (cid:48) (cid:17) . (3.32)The invariance of ρ G ( φ ) can be checked through L ∗ ρ G ( φ ) = 0 , (3.33)where ∗ is the adjoint with respect to the flat volume measure. Furthermore linear func-tions evolve to linear functions according toe L t φ j,α = N (cid:88) j (cid:48) =1 3 (cid:88) α (cid:48) =1 (e A t ) jα,j (cid:48) α (cid:48) φ j (cid:48) ,α (cid:48) , (3.34)where the matrix A = − ∂ j ⊗ diag( c , c , c ) − ∂ j ∂ T j ⊗ D (cid:93) , the first factor acting on j andthe second on α .We now add the nonlinearity N j,α . In general, this will modify the time-stationarymeasure and we have little control how. Therefore we propose to choose N j,α such thatthe corresponding vector field ∂ j N j,α is divergence free [9]. If N j,α depends only on thefield at sites j and j + 1, then the unique solution reads N j,α =
13 3 (cid:88) γ,γ (cid:48) =1 G αγγ (cid:48) (cid:0) φ j,γ φ j,γ (cid:48) + φ j,γ φ j +1 ,γ (cid:48) + φ j +1 ,γ φ j +1 ,γ (cid:48) (cid:1) . (3.35)For ρ G to be left invariant under the deterministic flow generated by the vector field − ∂ j N requires L ρ G = 0 , L = − N (cid:88) j =1 3 (cid:88) α =1 ∂ j N j,α ∂ φ j,α , (3.36)which implies N (cid:88) j =1 3 (cid:88) α =1 φ j,α ∂ j N j,α = 0 (3.37)and thus the constraints G αβγ = G βαγ (cid:0) = G αγβ (cid:1) (3.38)for all α, β, γ = 1 , ,
3, where in brackets we added the symmetry which holds by definition.Denoting the generator of the Langevin equation (3.28) by L = L + L , (3.39)21ne concludes L ∗ ρ G = 0, i.e. the time-invariance of ρ G .In the continuum limit the condition (3.37) reads (cid:88) α,β,γ =1 G αβγ (cid:90) dxφ α ( x ) ∂ x (cid:0) φ β ( x ) φ γ ( x ) (cid:1) = 0 , (3.40)where G αβγ = G αγβ . By partial integration2 (cid:88) α,β,γ =1 G αβγ (cid:90) dxφ α ( x ) φ β ( x ) ∂ x φ γ ( x ) = − (cid:88) α,β,γ =1 G αβγ (cid:90) dxφ β ( x ) φ γ ( x ) ∂ x φ α ( x ) . (3.41)Hence (3.40) is satisfied only if G γβα = G αβγ , which is the condition (3.38) claimed for thediscrete setting.(3.38) is the generalization of (3.19), which is specific for the anharmonic chain. Infact, while abstractly true, it is not so easy to verify directly. But now we can argue moreconvincingly why one should be allowed to continue with assuming the validity of theconstraints (3.38). As we will discuss in the next section, the leading coupling constantsare of the form G ααα , while the sub-leading couplings have equal lower indices, G αγγ , γ (cid:54) = α .The off-diagonal matrix elements are irrelevant for the large scale behavior. When onedoes the counting, all leading and sub-leading couplings can be chosen freely and theirrelevant couplings can be adjusted so that the constraint (3.38) is satisfied. Appealingto universality, the large space-time behavior should not depend on that particular choice.We expect that for general (cid:126)G the true time-stationary measure will have short rangecorrelations and nonlinear fluctuating hydrodynamics remains a valid approximation tothe dynamics of the anharmonic chain.In related problem settings, a different point of view has been suggested [33, 34].Firstly one notes that the Gaussian stationary measure (3.16), hence also (3.30), is simplyinherited from the canonical equilibrium measure. In this respect there is no choice.Also the nonlinear Euler currents are on the safe side. But the remaining terms arephenomenological to some extent. D (cid:93) , B (cid:93) could depend on (cid:126)φ itself. One could also includehigher derivative terms. In fact, one could try to choose the nonlinearities precisely insuch a way that the dynamics is invariant under time-reversal and leaves the Gaussianmeasure invariant. The program as such may be easily endorsed, but so far I have notseen a convincing handling of the details. Decoupling hypothesis . For the linear equations the normal modes decouple for longtimes. The hypothesis claims that such property persists when adding the quadraticnonlinearities. For the precise phrasing, we have to be somewhat careful. We consider afixed component, α , in normal mode representation. It travels with velocity c α , which isassumed to be distinct from all other mode velocities. If G ααα (cid:54) = 0, then for the purpose ofcomputing correlations of mode α at large scales, one can use the scalar conservation law ∂ t φ α + ∂ x (cid:0) c α φ α + G ααα φ α − D (cid:93)αα ∂ x φ α + B (cid:93)αα ξ α (cid:1) = 0 , (4.1)22hich coincides with the stochastic Burgers equation (1.19). If decoupling holds, onehas the exact asymptotics as stated in (1.21) with λ = 2 √ | G ααα | . The universal scalingfunction f KPZ is tabulated in [35], denoted there by f . f KPZ ≥ (cid:82) dxf KPZ ( x ) = 1, f KPZ ( x ) = f KPZ ( − x ), (cid:82) dxf KPZ ( x ) x (cid:39) . f KPZ looks like a Gaussian with a large | x | decay as exp[ − . | x | ]. Plots are provided in [17, 35].For an anharmonic chain, G = 0 always and the decoupling hypothesis applies onlyto the sound peaks, provided G = − G − − − (cid:54) = 0 which generically is the case. If G (cid:54) = 0,then the exact scaling form is S (cid:93)σσ ( x, t ) ∼ = ( λ s t ) − / f KPZ (cid:0) ( λ s t ) − / ( x − σct ) (cid:1) , λ s = 2 √ | G σσσ | , (4.2) σ = ±
1. To find out about the scaling behavior of the heat mode other methods have tobe developed.For one-dimensional fluids, van Beijeren [36] follows the scheme developed in [6] andarrived first at the prediction (4.2) together with the L´evy 5/3 heat peak to be discussedbelow. In [36] no Langevin equations appear. I regard them as a useful intermediatestep valid on a mesoscopic scale. In the Langevin form the theory can be applied to alarge class of one-dimensional systems. As a tool, fluctuating hydrodynamics has beenproposed considerably earlier [37] and used to predict the t − / decay of the total energycurrent correlation. One-loop, diagonal, and small overlap approximations . We return to the Langevinequation (3.28) and consider the mean zero, stationary φ j,α ( t ) process with ρ G as t = 0measure. The stationary covariance reads S (cid:93)αα (cid:48) ( j, t ) = (cid:104) φ j,α ( t ) φ ,α (cid:48) (0) (cid:105) = (cid:104) φ ,α (cid:48) e Lt φ j,α (cid:105) eq , t ≥ . (4.3)On the left, (cid:104)·(cid:105) denotes the average with respect to the stationary φ j,α ( t ) process and onthe right (cid:104)·(cid:105) eq refers to the average with respect to ρ G . By construction S (cid:93)αα (cid:48) ( j,
0) = δ αα (cid:48) δ j . (4.4)The time derivative reads ddt S (cid:93)αα (cid:48) ( j, t ) = (cid:104) φ ,α (cid:48) (e Lt L φ j,α ) (cid:105) eq + (cid:104) φ ,α (cid:48) (e Lt L φ j,α ) (cid:105) eq . (4.5)We insert e Lt = e L t + (cid:90) t ds e L ( t − s ) L e Ls (4.6)in the second summand of (4.5). The term containing only e L t is cubic in the time zerofields and hence its average vanishes. Therefore one arrives at ddt S (cid:93)αα (cid:48) ( j, t ) = A S αα (cid:48) ( j, t ) + (cid:90) t ds (cid:104) φ ,α (cid:48) e L ( t − s ) L (e Ls L φ j,α ) (cid:105) eq . (4.7)For the adjoint of e L ( t − s ) we use (3.34) and for the adjoint of L we use (cid:104) φ j,α L F ( φ ) (cid:105) eq = −(cid:104) ( L φ j,α ) F ( φ ) (cid:105) eq , (4.8)23hich both rely on (cid:104)·(cid:105) eq being the average with respect to ρ G . Furthermore L φ j,α = − ∂ j N j,α . (4.9)Inserting in (4.7) one arrives at the identity ddt S (cid:93)αα (cid:48) ( j, t ) = A S (cid:93)αα (cid:48) ( j, t ) − (cid:90) t ds (cid:104) (e A T ( t − s ) ∂ j N ,α (cid:48) )(e Ls ∂ j N j,α ) (cid:105) eq . (4.10)To obtain a closed equation for S (cid:93) we note that the average (cid:104) ∂ j (cid:48) N j (cid:48) ,α (cid:48) e L s ∂ j N j,α (cid:105) eq = (cid:104) ∂ j N j,α ( s ) ∂ j (cid:48) N j (cid:48) ,α (cid:48) (0) (cid:105) (4.11)is a four-point correlation. We invoke the Gaussian factorization as (cid:104) φ ( s ) φ ( s ) φ (0) φ (0) (cid:105) ∼ = (cid:104) φ ( s ) φ ( s ) (cid:105)(cid:104) φ (0) φ (0) (cid:105) + 2 (cid:104) φ ( s ) φ (0) (cid:105)(cid:104) φ ( s ) φ (0) (cid:105) . (4.12)The first summand vanishes because of the difference operator ∂ j . Secondly we replacethe bare propagator e A ( t − s ) by the interacting propagator S (cid:93) ( t − s ), which correspondsto a partial resummation of the perturbation series in (cid:126)G . Finally we take a limit ofzero lattice spacing. This step could be avoided, and is done so in our numerical schemefor the mode-coupling equations. We could also maintain the ring geometry which, forexample, would allow to investigate collisions between the moving peaks. Universality isonly expected for large j, t , hence in the limit of zero lattice spacing. The continuum limitof S (cid:93) ( j, t ) is denoted by S (cid:93) ( x, t ), x ∈ R . With these steps we arrive at the mode-couplingequation ∂ t S (cid:93)αβ ( x, t ) = (cid:88) α (cid:48) =1 (cid:16)(cid:0) − c α δ αα (cid:48) ∂ x + D αα (cid:48) ∂ x (cid:1) S (cid:93)α (cid:48) β ( x, t )+ (cid:90) t ds (cid:90) R dyM αα (cid:48) ( y, s ) ∂ x S (cid:93)α (cid:48) β ( x − y, t − s ) (cid:17) (4.13)with the memory kernel M αα (cid:48) ( x, t ) = 2 (cid:88) β (cid:48) ,β (cid:48)(cid:48) ,γ (cid:48) ,γ (cid:48)(cid:48) =1 G αβ (cid:48) γ (cid:48) G α (cid:48) β (cid:48)(cid:48) γ (cid:48)(cid:48) S (cid:93)β (cid:48) β (cid:48)(cid:48) ( x, t ) S (cid:93)γ (cid:48) γ (cid:48)(cid:48) ( x, t ) . (4.14)In numerical simulations of both, the mechanical model of anharmonic chains and themode-coupling equations, it is consistently observed that S (cid:93)αα (cid:48) ( j, t ) becomes approximatelydiagonal fairly rapidly. To analyse the long time asymptotics on the basis of (4.13) wetherefore rely on the diagonal approximation S (cid:93)αα (cid:48) ( x, t ) (cid:39) δ αα (cid:48) f α ( x, t ) . (4.15)Then f α ( x,
0) = δ ( x ) and the f α ’s satisfy ∂ t f α ( x, t ) = ( − c α ∂ x + D (cid:93)αα ∂ x ) f α ( x, t ) + (cid:90) t ds (cid:90) R dy∂ x f α ( x − y, t − s ) M αα ( y, s ) , (4.16)24 = − , ,
1, with memory kernel M αα ( x, t ) = 2 (cid:88) γ,γ (cid:48) =0 , ± ( G αγγ (cid:48) ) f γ ( x, t ) f γ (cid:48) ( x, t ) . (4.17)The solution to (4.16) has two sound peaks centered at ± ct and the heat peak sittingat 0. All three peaks have a width much less than ct . But then, in case γ (cid:54) = γ (cid:48) , theproduct f γ ( x, t ) f γ (cid:48) ( x, t ) (cid:39) t . Hence for the memory kernel (4.17) we invoke asmall overlap approximation as M αα ( x, t ) (cid:39) M dg α ( x, t ) = 2 (cid:88) γ =0 , ± ( G αγγ ) f γ ( x, t ) , (4.18)which is to be inserted in Eq. (4.16). Numerical simulations of the mode-coupling equations . When starting this projecttogether with Christian Mendl, in the summer of 2012 we spent many days in numericallysimulating the mode coupling equations with initial conditions S (cid:93)αα (cid:48) ( j,
0) = δ αα (cid:48) δ j . Onlya few plots are in print [38], simply because there is such a large parameter space and itis not clear where to start and where to end. Still, for our own understanding this periodwas extremely helpful. Mostly we simulated in Fourier space. System size was up to 400.Speeds were of order 1, thereby limiting the simulation time to about 200, the time of thefirst peak collision. For such sizes the simulations are fast and many variations could beexplored. We started from the scalar equation, to be discussed below, moved up to twomodes, and eventually to three modes with parameters taken from an actual anharmonicchain. | (cid:126)G αα (cid:48) α (cid:48) | was either 0 or somewhere in the range 0 . . D (cid:93) is a free parameterwhich was varied from 0 to | (cid:126)G αα (cid:48) α (cid:48) | /
2. We always simulated the complete matrix-valuedmode-coupling equations (4.13). Our main findings can be summarized as: (i)
For a large range of parameters, the diagonal approximation in generally failed forshort times, but was quickly restored with the off-diagonal elements being at most 10%of the diagonal ones. (ii)
The results were fairly insensitive to the choice of D (cid:93) . In fact, D (cid:93) = 0 works also well.Apparently the memory term generates already enough dissipation. (iii) We varied the overlap coefficients G αγγ (cid:48) with γ (cid:54) = γ (cid:48) . Over the time scale of thesimulation no substantial changes were observed.All these findings confirm the approximations proposed.As our biggest surprise, except for trivial cases we never reached the asymptotic regime.The peak structure develops fairly rapidly. The peak shape then changes slowly, roughlyconsistent with the predicted scaling exponents, but it does not reach a self-similar form.For example in the case of the sound peak, on the scale t / , rather than being symmetric,as claimed by (4.2), it is still badly distorted, tilted away from the central peak with rapiddecay outside the sound cone but rather slow power law type of decay towards the centralpeak. To improve one would have to simulate larger system sizes and longer times. Butthen numerical simulations become heavy and the fun evaporates. More attention can beachieved by molecular dynamics (MD) simulation of the mechanical chain.25or given parameters V, P, β one easily computes all the required coefficients. So onegoal was to run a MD and put the results on top of the ones from a simulation of the modecoupling equations. For this to be a reasonable program, one would have to simulate themode-coupling equations for sizes of N = 4000 and more, which we never attempted.For the scalar case the situation is much simpler. The mode-coupling equation takesthe form ∂ t f ( x, t ) = D∂ x f ( x, t ) + 2 G (cid:90) t ds (cid:90) R dy∂ x f ( x − y, t − s ) f ( y, s ) , (4.19)which is the one-loop approximation for the stochastic Burgers equation [39] . For large x, t , its solution with initial condition f ( x,
0) = δ ( x ) takes the scaling form f ( x, t ) ∼ = ( λ s t ) − / f mc (cid:0) ( λ s t ) − / x (cid:1) . (4.20)Inserting in (4.19), one first finds the non-universal scaling coefficient λ s = 2 √ | G σσσ | . (4.21)Secondly ˆ f mc , the Fourier transform of f mc , is defined as solution of the fixed point equation ˆ f (cid:48) mc ( w ) = − π w (cid:90) ds ˆ f mc ((1 − s ) / w ) (cid:90) R dq ˆ f mc ( s / ( w − q )) ˆ f mc ( s / q ) (4.22)with w ≥ f mc (0) = 1, ˆ f (cid:48) mc (0) = 0.Eq. (4.22) is based on the closure assumption (4.12) and there is no reason to inferthat it is exact. However from our numerical simulations we conclude that f KPZ differsfrom f mc by a few percent only. We regard the scalar case as a strong support for theentire approach. But for several components the large finite size effects prohibit one toarrive at a similarly simple claim. Asymptotic self-similarity . Within mode-coupling the asymptotic shape function forthe sound peaks is given by f σ ( x, t ) ∼ = ( λ s t ) − / f mc (cid:0) ( λ s t ) − / ( x − σct ) (cid:1) , (4.23) σ = ±
1. For the heat peak we employ (4.16) together with (4.18), using as an input thatthe asymptotic form of f σ is known already. In fact the scaling exponent for f σ is crucial,but the precise shape of f σ enters only mildly. Hence, again switching to Fourier space,one has to solve ∂ t ˆ f ( k, t ) = − D (cid:93) (2 πk ) ˆ f ( k, t ) − (cid:88) σ = ± ( G σσ ) (2 πk ) (cid:90) t ds ˆ f ( k, t − s ) (cid:90) R dq ˆ f σ ( k − q, s ) ˆ f σ ( q, s ) , (4.24)ˆ f ( k,
0) = 1. For ˆ f σ one inserts the asymptotic result (4.20). (4.24) is a linear equationwhich is solved through Laplace transform with the resultˆ f ( k, t ) ∼ = e −| k | / λ h t , (4.25)26here λ h = λ − / ( G σσ ) (4 π ) (cid:90) ∞ dtt − / cos(2 πct ) (cid:90) R dxf mc ( x ) = λ − / ( G σσ ) (4 π ) (2 πc ) − / π ) 1cos( π ) (cid:90) R dxf mc ( x ) (4.26)and we used the symmetry G σσ = − G − σ − σ , see [29] for details. (4.25) is the Fouriertransform of the symmetric α -stable distribution with exponent α = 5 /
3, also known asL´evy distribution. In real space the asymptotics reads, for | x | ≥ ( λ h t ) / , f ( x, t ) (cid:39) π − λ h t | x | − / . (4.27) f mc is a smooth function with rapid decay. On the other hand, f has fat tails andits variance is divergent. According to (4.27), at x = σct the heat peak f ( σct, t ) ∼ = π − λ h c − / t − / . This explains why there is still coupling between f and f σ , despite thelarge spatial separation. In fact, numerically one observes that beyond the sound cone, x = ± ct , the solution decays exponentially fast. As t becomes large the tails of f arebuild up between the two sound peaks, so to speak they unveil the L´evy distribution.In (4.24) we could also insert for f σ the exact scaling function f KPZ , which wouldslightly modify λ h . f L´evy is an approximation, just as f mc . But the MD simulationsdisplay so convincingly the L´evy distribution that one might be willing to regard it asexact. If so, the exact λ h must be based on f KPZ . To obtain the correlations of thephysical fields, one has to use S ( j, t ) = R − S (cid:93) ( j, t ) R − T . (4.28)In particular the correlations of the physical fields are given through S αα ( j, t ) = (cid:88) σ =0 , ± | ( R − ) ασ | f σ ( j, t ) , (4.29)where for f σ the asymptotic scaling form is inserted. Then asymptotically the (cid:96) - (cid:96) andthe e - e correlations show generically all three peaks. However, for the u - u correlations thecentral peak is missing asymptotically, since ( R − ) = 0 independently of the interactionpotential V .We note that the coefficient D (cid:93) does not appear in the asymptotic scaling form, ofcourse neither B (cid:93) B (cid:93) T = 2 D (cid:93) . This result is consistent with the picture that noise and dis-sipation are required to maintain the correct local stationary measure with susceptibility C . The long time asymptotics is however governed by the nonlinearities. No signal beyond the sound cone . Physically the sound speed is an upper limitfor the propagation of small disturbances. Since the initial state has a finite correlationlength, one would expect that towards the outside of the sound cone correlations decayexponentially, while inside the sound cone there seems to be no particular restriction. Asa consistency check, one would hope that such a general feature is properly reproduced27y mode-coupling. Their numerical solutions conform with this expectation, at least forthe small system sizes explored. But for the scaling limit one has to let t → ∞ and thedecay information seems to be lost. However there is still a somewhat subtle trace.To explain, we have to first recall some properties of L´evy stable distributions. Exceptfor trivial rescalings, they are characterized by two parameters, traditionally called α, β ,where for simplicity we momentarily stick to this convention, without too much risk ofconfusion. The probability density has a simple form in Fourier space,ˆ f L´evy ,α,β ( k ) = exp (cid:0) − | k | α (cid:2) − i β tan( πα )sgn( k ) (cid:3)(cid:1) . (4.30)The parameter α controls the steepness, 0 < α <
2, while β controls the asymmetry, | β | ≤
1. For | β | > α = 2 only β = 0 is admitted and the probability density is a Gaussian. If | β | <
1, the asymptotic decay of f L´evy ,α,β ( x ) is determined by α and is given by | x | − α − for | x | → ∞ . At | β | = 1 the two tails show different decay. The functions correspondingto β = 1 and β = − β = 1 the slow decay being for x → −∞ and still as | x | − α − . For 0 < α ≤ f L´evy ,α, ( x ) = 0 for x >
0, while for 1 < α < − c x α/ (1 − α ) ) with known constant c . Werefer to [40] for more details.For the heat peak we obtained the symmetric L´evy distribution because the soundpeaks are reflection symmetric, implying c = − c − and ( G ) = ( G − − ) . If hypo-thetically we would choose distinct couplings, or c (cid:54) = − c − , then this imbalance wouldproduce a β (cid:54) = 0. If one of the sound peaks would be completely missing, as in the casefor a system with only two conserved fields, then necessarily | β | = 1. In accordance withthe physical principle, the sign of β is such that the fast decay of f L´evy ,α, ± ( x ) is towardsthe outside of the sound cone, while the slow decay is towards the single sound peak. Forfinite t , this slow decay will be cut by the sound peak. Thus the scaling solution of themode-coupling equations reproduces the rapid decay towards the exterior of the soundcone. This is a completely general fact, any number of components and any (cid:126)G [41]. Dynamical phase diagram . As already indicated through the particular case G = 0,the large scale structure of the solution depends on whether G αγγ = 0 or not. One extremecase would be G ααα (cid:54) = 0 for all α , implying that the three peaks have KPZ scaling be-havior. The other extreme is G αγγ = 0 for all α, γ , resulting in all peaks to have diffusivebroadening. For the case of only two modes, the full phase diagram has seven distinctphases, with unexpected details worked out in [42]. For the general case of n componentsthe long time asymptotics is completely classified in [43, 41]. Anharmonic chains havespecial symmetries and not all possible couplings (cid:126)G can be realized. Given that G = 0and because the sound peaks are symmetric, to have a distinct scaling requires G = 0 , (4.31)which can be realized. The behavior is then determined by the value of the remainingdiagonal matrix elements. For the central peak one finds σG σσ > , (4.32)28hile for the sound peak diagonals, G − − = − G − , G = − G − , there seems to be noparticular restriction. In principle, there could be sort of accidental zeros of G αγγ whichare then difficult to locate. A more direct approach starts from the observation that the (cid:126)G coefficients are expressed through cumulants in r , V . If the integrands are antisymmetricunder reflection, many terms vanish. The precise condition on the potential is to havesome a , P such that V ( x − a ) + P x = V ( − x − a ) − P x (4.33)for all x . Then for arbitrary β and P = P , one finds G = 0 , G − − = − G − = 0 , G = − G − = 0 , (4.34)while σG σ σ (cid:48) >
0. The standard examples for (4.33) to hold are the FPU chain with nocubic interaction term, the β -chain, and the square well potential with alternating masses,both at zero pressure.Under (4.34) the heat mode is coupled to the sound mode, but there is no back reactionfrom the sound mode. Hence the sound peak is diffusive with scaling function f σ ( x, t ) = 1 √ πD s t e − ( x − σct ) / D s t . (4.35) D s is a transport coefficient. It can defined through a Green-Kubo formula, which alsomeans that no reasonably explicit answer can be expected. The feed back of the soundpeak to the central peak follows by the same computation as before, with the resultˆ f ( k, t ) = e −| k | / λ h t , (4.36)where λ h = ( D σ ) − / ( G σσ ) (4 π ) (2 πc ) − / (cid:90) ∞ dt t − / cos( t )(2 √ π ) − . (4.37)Since 3 / < /
3, the density f ( x, t ) turns out to be broader than the L´evy 5 / G (cid:54) = 0.In testing nonlinear fluctuating hydrodynamics almost subconsciously one tries toconfirm (or not) the scaling exponents, resp. functions. This can be difficult becauseof limited size. The dynamical phase diagram offers a different option. For exceptionalpoints in the phase diagram, without too precise a verification of the scaling, one shouldfind that the standard scaling exponent does not properly fit the data. Such qualitativeproperty is possibly more easy to access. In Fig. 2 we display heat and sound peak for aFPU chain with G = 0. In 1953 Fermi, Pasta and Ulam, technically supported by Tsingou, simulated 32 particlesinteracting through a quartic potential at the extremely low energy of e = 5 × − per29
10 -5 0 5 10 x/ (λ t) ( λ e t ) / S ( x , t ) t=800t=1200t=1600Levy 3/2 -4 -2 0 2 4 (x+ct)/( λ t) ( λ s t ) / S - - ( x , t ) t=800t=1200t=1600Gaussian Figure 2: Scaling plot of heat and sound peak for a FPU chain with N = 8192, potentialparameters α = 0, β = 1, pressure P = 1, and temperature β − = 1.particle (above the ground state energy) [25]. They injected energy in the highest Fouriermode and were looking for equipartition of the modes at long times. However they foundquasi-periodic motion with time averages settling to some definite value different fromequipartition. The observed quasi-periodicity triggered the connection to KAM tori, thediscovery of integrable systems with many degrees of freedom, and to the developmentof the theory of solitons and breathers [44, 45]. For sure, a rich harvest, see [46] for ahistorical perspective. Later on Izrailev and Chirikov [47] repeated the simulation at thehigher energy e = 5 × − and observed equipartition.Anomalous transport surfaced much later [48], see the reviews [49, 50]. One connectsthe two ends of the chain to thermal reservoirs. To explain, the end particles are tieddown as q = 0, q N +1 = 0, and to the equations of motion for the boundary particles oneadds Langevin terms as¨ q = − γp + (cid:112) γT − ξ − , ¨ q N = − γp N + (cid:112) γT + ξ + , (5.1)where γ is a friction constant, T ± are the boundary temperatures, and ξ ± ( t ) are indepen-dent standard Gaussian white noises. For T − = T + the system settles in the canonicalequilibrium state. But for T − (cid:54) = T + there is a non-trivial steady state with a non-zeroenergy flux j e ( N ) depending on the length, N , of the chain. For regular heat transport,Fourier’s law implies j e ( N ) (cid:39) c N − . However for FPU chains one finds an enhancedtransport as j e ( N ) (cid:39) c N − α ( T − − T + ) (5.2)with an exponent α characterizing the anomaly. [A further α , but better to stick tostandard conventions]. Over the last two decades many MD simulations have been imple-mented for a wide variety of one-dimensional systems. Early results indicated α = 2 / α = 1 / (cid:96) (0 , t ) = 0 = (cid:96) ( L, t ), u (0 , t ) = 0 = u ( L, t ), but e (0 , t ) = e − and e ( L, t ) = e + and tries to investigate the steadystate. Unfortunately, at least for the moment, we have no powerful techniques to dealwith this problem. On the other hand it is argued [49] that the energy flux is related toa Green-Kubo formula by j e ( N ) ∼ (cid:90) N/c dt (cid:0) (cid:104)J ( t ); J (0) (cid:105) − (cid:104)J ( ∞ ); J (0) (cid:105) (cid:1) . (5.3)Under the time integral appears the total energy current correlation in thermal equilibriumwith its possibly non-zero value at t = ∞ subtracted, see Section 6 for more explanations.The decay of such correlation can be predicted by mode-coupling. In fact, we will confirmthe value α = 1 /
3. But the argument is subtle because it is only indirectly related to thespreading of the heat peak which is on scale t / .Because of the relation (5.3), in many MD simulations the total energy current corre-lation (cid:104)J ( t ); J (0) (cid:105) is measured as an addition to steady state transport [51, 52]. Thereare also MD simulations exclusively focussed on momentum and energy current correla-tions [53]. However simulations of correlations of the conserved fields have been fairlyscarce until recently. The peak structure was noted already in [54], see also [55, 56]. Butsurprisingly enough, even such a basic issue as the quantitative comparison between themeasured speed of sound and formula (2.31) is apparently not a routine check. So farthere have been three independent sets of simulations with the specific aim to check thepredictions from mode-coupling. For the details the reader is encouraged to look at theoriginal papers. I will try to compare so to reach some sort of conclusion.The first and second set are FPU chains with either symmetric or asymmetric poten-tial, both the αβ and the pure β chain [57, 58]. In this case one has to integrate numer-ically the differential equations governing the evolution, for which both a velocity-Verletalgorithm and a fourth order symplectic Runge-Kutta algorithm are used. The third setconsists of chains with a piecewise constant potential [59]. We call them hard-collision,since the force is zero except for δ -spikes and the dynamics proceeds from collision to colli-sion. Now one has to develop an efficient algorithm by which one finds the time-wise nextcollision. Except for rounding, there is no discretization of time. The simplest exampleis the hard-point potential, V hc ( x ) = ∞ for x < V hc ( x ) = 0 for x ≥
0. A variantis the infinite square well potential, V sw ( x ) = ∞ for x < x > a and V sw ( x ) = 0 for0 ≤ x ≤ a [60, 59]. In this case two neighboring particles at the maximal distance a arereflected inwards as if connected by a massless string of finite length a . For both modelsthe dynamics remembers the initial velocities. To have only the standard conservationlaws, one imposes alternating masses, say m j = m for even j and m j = m for odd j . Inboth models the unit cell then contains two particles and the scheme explained before hasto be extended. But at the very end the difference is minimal. A further variant is theshoulder potential V sh ( x ), for which V sh ( x ) = ∞ for | x | < , V sh ( x ) = (cid:15) for ≤ | x | ≤ V sh ( x ) = 0 for | x | >
1. The potential is either repulsive, (cid:15) >
0, or attractive (cid:15) < P = 0 constitute a distinct dynamical phase.Such phase is absent for the hard-point and the square shoulder potentials. But the squarewell at P = 0 has the same properties as can be seen from taking a = a/ N = 2 to 2 , even size 2 = 65 ,
536 has been attempted[61]. Periodic boundary conditions are imposed. Time is restricted to t ≤ t max = N/ c ,the time of the first collision between sound peaks. For given potential, one has to decideon the thermodynamic parameters, P, β . One constraint is to have c approximately in therange 1 ... (cid:126)G , which should notbe too small, at least for the relevant couplings. For the models from above they aretabulated. The relevant couplings show quite some variability taking values in the range0 . ... .
4. At the very end, one has to make a physically reasonable choice, perhaps usethe same parameters as previous MD simulations so to have the possibility to compare. Asystematic study of the dependence on
P, β seems to be too costly and most likely not sointeresting. But it does make sense to probe values at the border. For some of the FPUsimulations the energy per particle is chosen as e = 0 .
1, at which point nonlinearities aresmall [57]. There are also very extensive simulations at the even lower energy e = 5 × − ,with the goal to explore the route to equipartition, which is a somewhat distinct issue[62].Once all parameters are fixed, there are several options to run the simulation. Sincein canonical equilibrium { r j , p j , j = 1 , ..., N } are independent random variables, one cansample the initial conditions through a random number generator. For the hard-collisionpotentials the geometric constraints are still simple enough for allowing one to generatethe microcanonical ensemble by Monte Carlo methods. In our simulations the correlator S hardly depends on the choice of the ensemble. With such generated random initialdata the equations of motion are simulated up to t max . A single run is noisy and onehas to repeat many times, order 10 . The much more common choice is to start from areasonable nonequilibrium configuration and to equilibrate before measuring correlations.Usually one then simulates very long trajectories, up to times of order 2 , over which thetime lag g α ( j, t + τ ) g α (cid:48) (0 , t ) is sampled, τ ≤ t max . In addition one averages over a smallnumber of runs, of order 10 . The total number of samples is roughly the same in bothmethods. The random number generator method produces the thermal average with ahigher reliability.The sampled S αα (cid:48) ( j, t ) can be Fourier transformed in the spatial variable and/or in thetime variable. One can also transform to normal coordinates. These are linear operationswhich can be done for each sample or only after averaging. Should one keep the fullresolution or only some data points? Of course it depends somewhat on the goals. In[58, 59], the full 3 × R matrix to obtain S (cid:93) ( j, t ). This approach allows to test diagonality. Becausethe peak structure is most easily seen in the space-like j coordinate, maximal resolution32 x/( λ e t) ( λ e t ) / S ( x , t ) t=800t=1300t=2700Levy (a) -4 -2 0 2 4 (x+ct)/( λ s t) ( λ s t ) / S - - ( x , t ) t=500t=800t=1300t=2700KPZ(b) Figure 3: Scaling plot of heat and sound peak for a FPU chain with N = 8192, potentialparameters α = 2, β = 1, pressure P = 1, and temperature β − = 0 . j is retained and only three times ( t = 250 , , t . In [57] the lowest Fourier modes are plotted as a function of the frequency ω . Aseparate issue are the much simulated total current correlations. The total currents aresampled directly, in the most complete version momentum, energy, and cross correlations,and then plotted as a function of t or ω .We reproduce only a few figures. Many more details can be found in the originalpapers. In Fig. 3 we display the data for the FPU chain [58] with V FPU ( x ) = x + x + x , β = 2, P = 1, N = 8192, and t max = 2700. The sound speed is c = 1 .
45. Note thatthe sound peak is somewhat distorted, not symmetric relative to ct , but has a rapid fall-offaway from the sound cone. As only fit parameter one uses λ h , resp. λ s . The optimal fitat the longest available time is denoted by λ emph , resp. λ emps , standing for empirical value.In most cases there is also a theoretical value based on decoupling and/or mode-coupling,which is indicated in square brackets. For the FPU simulations the results are for theheat peak λ emph = 13 . . λ emps = 2 .
05 [0 . N = 4096 with t max = 1024. The following parameters have been chosen, shoulder: (cid:15) = 1, P = 1 . β = 2, c = 1 .
74, hard-point: m /m = 3, P = 2, β = 0 . c = 1 . m /m = 3, a = 1, P = 0, β = 2, c = 1 .
73. The fit to the predicted scalingfunction has an error less than 5%. For the hard-collision models the simulation resultsare, shoulder: λ emps = 1 .
62 [1 . λ emph = 1 .
44 [1 . λ emps = 1 .
42 [2 . λ emph = 1 .
04 [0 . λ emph = 0 .
95 [1 . -
10 0 10 20 x10 - - - - H l t L (cid:144) S HH l t L (cid:144) x,t L (a) shoulder, heat - -
10 0 10 20 30 x10 - - - H l t L (cid:144) S HH l t L (cid:144) x,t L (b) hard-point gas, heat - - - - - - H l t L (cid:144) S HH l t L (cid:144) x,t L (c) square-well, heat - - - - H l t L (cid:144) S HH l t L (cid:144) x + ct,t L (d) shoulder, sound - - - H l t L (cid:144) S HH l t L (cid:144) x + ct,t L (e) hard-point gas, sound - - - - - - H l t L (cid:144) S HH l t L (cid:144) x + ct,t L (f) square-well, sound Figure 4: Heat and sound peak for shoulder, hard-point, and square well potential inlogarithmic scale.at P = 0 is in a distinct dynamical universality class, hence the different scaling expo-nents. For this model λ s is related to a diffusion constant, which can be obtained onlynumerically.The MD simulation [57] is at low temperatures and considers the positional correla-tions, which in Fourier space differ from the stretch correlations only by a k -dependentprefactor. Thus sampled is the correlator(1 − cos(2 πk )) − (cid:90) dt e i ωt ˆ S ( k, t ) = (1 − cos(2 πk )) − ˆ S ( k, ω ) . (5.4)For low temperatures the area under the central peak is a factor 10 smaller than the oneunder the sound peak. Hence only the sound peak is explored. Its asymptotic scalingform is ˆ S ( k, t ) = cos(2 π i kct ) C ˆ f KPZ ( k ( λ h | t | ) / ) . (5.5)Considering only the right moving sound peak by setting ω max = 2 πkc ,ˆ S ( k, ω + ω max ) = (cid:90) dt e i ωt C ˆ f KPZ ( k ( λ h | t | ) / )= (cid:90) dt e i( ω/λ h | k | / ) t ( λ h | k | / ) − C ˆ f KPZ ( | t | / ) . (5.6)Thus defining h KPZ ( ω ) = (cid:90) dt e i ωt ˆ f KPZ ( | t | / ) , (5.7)34ne arrives at ˆ S ( k, ω + ω max ) = C ( λ | k | / ) − h KPZ ( ω/λ | k | / ) . (5.8)If one normalizes the maximum of ˆ S to 1, then h KPZ ( ω ) is replaced by h KPZ ( ω ) /h KPZ (0),which amounts to setting the prefactor in (5.8) equal to 1. In Fig. 5 the spectrum at twodifferent choices of the asymmetry parameter is displayed.In the simulation the potential is chosen as V ( x ) = x + α x + x , where theasymmetry varies from 0 to 2. Increasing α , the pressure P increases from 0 to 0 . .
55 to 9 .
75. The sound speed c (cid:39) .
1. In frequency spacethe peak moves linearly in k at around ω = 0 .
01. For the shape function one uses ˆ S ( k, ω )for k = 1 , , , ,
16 to generate a scaling plot. Over the whole range of α ’s the fit withthe scaling function (5.7) is fairly convincing. The optimal fit parameter starts from λ emps = 0 .
02 [0 .
04] at α = 0 . λ emps = 0 .
08 [0 .
37] at α = 1 . λ emps = 0 .
07 [0 . α = 2 . -40 -20 0 20 40( ω −ω max ) / λ num q S ( q , ω ) / S m a x q=0.01227q=0.02454q=0.04908h( ω)/ h(0) KPZ scaling - FPU αβ e=0.5 α=0.1 β=1 λ num = 0.0055 Figure 5: Scaling plot of the sound peak in ( k = q, ω ) variables of a FPU chain with N = 1024, potential parameters α = 0 . β = 1, pressure P = − .
04, and temperature β − = 0 . c = 1 . (i) The separation into three peaks is a fast process. In normal mode representationthe off-diagonal matrix elements are indeed small. There are no correlations beyond thesound cone. (ii)
The distinction between G (cid:54) = 0 and G = 0 is seen very convincingly.35 iii) The central peak is on the t / scale and adjusts well to the predicted L´evy dis-tribution with one caveat. As seen from the center, the shape function fairly rapidlyswitches into the power law decay. However at the location of the sound peaks there aresome wiggles and beyond only small amplitude noise is observed. In this sense, the L´evydistribution gets uncovered as time progresses. The L´evy distribution is a fairly directconsequence of G = 0 and thus one of the strongest supports for the theory. (iv) The sound peaks are mirror images of each other. One plots them on the t / scale,but then there is still a slow change in time. For the hard-collision models the shape isalmost perfect, but the λ s parameter is dropping in time. For the FPU models the peakis distorted and still away from the symmetric shape predicted by the theory. To theoutside of the sound cone there is the rapid fall-off in accordance with the KPZ scalingfunction. But towards the heat peak there is slow decay. The sound peak is tilted awayfrom the central peak. Apparently there is still a strong interaction between the peaks. (v) In the majority of the simulations the peaks vary slowly on the scale t / for sound,resp. t / for the heat peak. Thus it becomes meaningful to use as a fit the theoreticalscaling function with λ s , resp. λ h , as only free parameter. This optimal choice has beendenoted above as empirical value λ emps / h . The error between the measured and theoreticalshape function is less than 5%. However one observes that λ emps and λ emph are still changingin time signaling that the simulation has not yet reached the truly asymptotic regime. Insome simulations λ emps / h drops monotonically in time and differs not too strongly from λ s / h .One is then willing to believe that for even longer simulation times the asymptotic valueis reached. But in other simulations there is a much stronger discrepancy, which asks formore explanations.Mode-coupling is not specific to anharmonic chains. In principle any one-dimensionalsystem with conserved fields can be handled by the same scheme. This offers the possibilityto test the theory through other models, possibly finding systems with less strong finitetime effects. One obvious choice are stochastic lattice gases with several type of particleslike several lane TASEP [64, 43] and the AHR model [65, 66]. For them the couplings canbe more easily adjusted than for anharmonic chains, which offers the possibility to test thedynamical phase diagram. Also anharmonic chains with a stochastic collision mechanism,respecting the conservation laws, have been studied in considerable detail [67, 68, 42]. The total current is a fluctuation observable, in contrast to S αα (cid:48) ( j, t ) which refers to theaverage of the product of two local observables. Thus we need some additional consider-ations to establish the link to nonlinear fluctuating hydrodynamics. For a ring of size N ,Λ N = [1 , ..., N ], the total currents are defined by (cid:126) J tot , Λ N ( t ) = 1 √ N N (cid:88) j =1 (cid:126) J ( j, t ) (6.1)36nd the total current covariance readsΓ Λ N ,αα (cid:48) ( t ) = (cid:104)J tot , Λ N ,α ( t ); J tot , Λ N ,α (cid:48) (0) (cid:105) P,β, Λ N = N (cid:88) j =1 (cid:104)J α ( j, t ); J α (cid:48) (0 , (cid:105) P,β, Λ N . (6.2)The cumulant (cid:104)· ; ·(cid:105) means that the static average is subtracted and system size is indicatedexplicitly. In the limit N → ∞ Γ αα (cid:48) ( t ) = (cid:88) j ∈ Z (cid:104)J α ( j, t ); J α (cid:48) (0 , (cid:105) P,β . (6.3)For fixed t the integrand decays exponentially in j , but with a correlation length increasingin time. Γ ( t ) (a) momentum current correlations
10 100 1000 t10 - - - - - Γ ( t )- m c / β (b) logarithmic plot of Γ ∆ ( t ) Γ ( t ) (c) energy current correlations Γ ( t )- P /( m β ) (d) logarithmic plot of Γ ∆ ( t ) Figure 6: Total momentum and energy correlations for a hard-point particles with alter-nating massesBefore we emphasized the transformation to normal coordinates. For the currents westick to the physical fields. The current J tot , Λ N , ( t ) is itself conserved. Thus only the (2 , αα ( t ) = Γ αα ( − t ), while the off-diagonal elements are odd and satisfyΓ ( t ) = − Γ ( − t ) = Γ ( − t ) = − Γ ( t ) . (6.4)37ode-coupling predicts that this matrix element vanishes. In our simulations we observean exponential decay with a decay time of order 20 to 30. Thus the correlations of realinterest are Γ ( t ) and Γ ( t ). For chains the latter is the most frequently simulatedequilibrium time correlation. The momentum current correlations have been measured in[53, 69] and the momentum-energy cross correlations only recently in [70].Γ Λ N ,αα (cid:48) ( t ) is a fluctuation observable, for which the equivalence of ensembles does nothold. But the time-independent difference between microcanonical and canonical averagecan be computed explicitly. For the microcanonical ensemble, lim t →∞ Γ Λ N , micro ( t ) = 0.On the other hand, there is no reason for Γ Λ N ( t ) to vanish asymptotically if the canonicalaverage, as in (6.2), is used. In fact, for infinite volume,lim t →∞ Γ ( t ) = β − c , lim t →∞ Γ ( t ) = β − P . (6.5)The asymptotic values in (6.5) are called Drude weight, which has received a lot ofattention in the context of current correlations for integrable quantum chains. A non-zero Drude weight indicates that the correlator of the corresponding conserved field has aballistically moving component, but it cannot resolve the structure of this component. Fornon-integrable anharmonic chains, the ballistic pieces are just the sharply concentratedsound peaks, while in the integrable case one expects to have a broad spectrum whichexpands ballistically.The link to mode-coupling is achieved through the general observation that the currentcorrelations are proportional to the memory kernel. We use its diagonal approximationand insert the asymptotic form of f α . The memory kernel is in normal mode representa-tion. Thus we still have to transform back to the physical fields through the R matrix.The computation can be found in [70] with the resultΓ ∆ ( t ) = Γ ( t ) − β − c (cid:39) ( λ h t ) − / (cid:104) ψ , H u ψ (cid:105) (cid:90) dx f L´evy , / ( x ) + ( λ s t ) − / (cid:104) ψ , H u ψ (cid:105) (cid:90) dx f KPZ ( x ) Γ ∆ ( t ) = Γ ( t ) − β − P (cid:39) c β − ( λ s t ) − / (cid:90) dx f KPZ ( x ) , (6.6)where { ψ α } are the eigenvectors of A , Aψ = 0, Aψ = cψ , see [29]. If (4.33) is satisfied, e.g. an even potential at P = 0, then the sound peak is diffusive, see Eq. (4.35), and thecentral peak is L´evy 3 /
2. Furthermore (cid:104) ψ , H u ψ (cid:105) = 0 = (cid:104) ψ , H u ψ (cid:105) . Hence Γ ( t ) − β − c is expected to decay integrably, whileΓ ( t ) − β − P (cid:39) c β − (8 πD s t ) − / . (6.7)The energy current correlation is predicted to decay as t − / which has been reportedin MD simulations already more than 15 years ago [51, 52]. The true mechanism behindthe decay is actually somewhat subtle. From the conservation law it follows that thesecond moment of the heat peak is related to the second time derivate of the currentcorrelation. Using the asymptotic form (4.27), including the cut-off at the sound peak,one arrives at d dt (cid:90) ct − ct dx x f ( x, t ) (cid:39) π ( λ h ) / c / ( λ h t ) − / . (6.8)38his argument overlooks that the scaling form is for the normal mode representation, whileΓ ( t ) refers to the physical energy current. The complete computation leads however tothe same power law except for a different prefactor [70]. Γ ( t ) (a) momentum current correlations - - - - Γ ( t )- c / β (b) logarithmic plot of Γ ∆ ( t ) Γ ( t ) (c) energy current correlations × - × - Γ ( t )- P / β (d) logarithmic plot of Γ ∆ ( t ) Figure 7: Total momentum and energy current correlations for a hard collision modelwith shoulder potential.The momentum current correlation should decay as t − / , which is a recent finding.However, its prefactor could vanish, in principle. Mode-coupling with the currently avail-able precision would not provide an answer, then. In contrast the prefactor for Γ ( t ) isstrictly positive.Our simulation results [70] are shown in Fig. 6,7. The parameters are as before, N = 4096, t max = 1024, shoulder: (cid:15) = 1, P = 1 . β = 2, c = 1 .
74, and hard-point: m /m = 3, P = 2, β = 0 . c = 1 .
73. The red lines indicate the predictions based onmode-coupling. It is interesting to note that for the shoulder potential the evidence fora t − / decay is not so overwhelming as one might have anticipated and, by looking ata different time window, one could as well fit to a slightly different exponent. On theother hand the hard-point potential with alternating masses shows a very clean powerlaw decay. Through MD with shoulder potential the predicted decay of Γ ( t ) is wellconfirmed. However, for the hard-point potential it so happens that both prefactors, (cid:104) ψ , H u ψ (cid:105) and (cid:104) ψ , H u ψ (cid:105) , vanish. Numerically we estimate a decay as t − . Again thisis a strong qualitative support of nonlinear fluctuating hydrodynamics. One might have39hought that all hard-collision potentials have the same asymptotic power law for themomentum current correlation. This expectation is born out under the proviso that therespective prefactors do not vanish. Matrix elements, as (cid:104) ψ , H u ψ (cid:105) , must come from amicroscopic computation and cannot be deduced by a mere inspection of the potential.An additional confirmation of (6.5), (6.7) has been accomplished recently [71]. Forthe potential V ( x ) = ax + x , β = 1 and at P = 0 .
59 with a = −
2, resp. at P = − . a = 1 .
89, up to very small errors the signature of the (cid:126)G matrices is identical to an evenpotential at P = 0. In the MD simulation the energy current is found to decay as t − / and the momentum current seems to be integrable. Keeping all parameters fixed andshifting slightly to a = − .
7, resp. a = 2, the decay as stated in (6.6) is restored. For the anharmonic chains studied so far, the potential depends only on q j +1 − q j and henceremains without change under spatial translations. Physically this property is obviousand seems hard to avoid. However there could be a substrate potential which forcesthe particles preferentially to particular locations. One could consider a two-componentsystem, which then has an acoustic and an optical mode. The latter would be comparableto a one-component system with an on-site potential. Such considerations lead to the moregeneral class of Hamiltonians H os = N (cid:88) j =1 (cid:0) p j + V os ( q j ) (cid:1) + N − (cid:88) j =1 V ( q j +1 − q j ) (7.1)with some confining on-site potential V os . The only conserved field is the energy. Thereis a unique equilibrium measure. The Euler currents vanish. From the perspective offluctuating hydrodynamics all evidence points towards diffusive energy transport. For thecase of a quadratic V and for V os = V FPU very detailed MD simulation confirm diffusivetransport [72].More interesting are models with two conservation laws. We discuss separately coupledrotators, which can be thought of as a classical limit of a quantum Heisenberg chain, andthe discrete nonlinear Schr¨odinger equation on a lattice, which is the classical field theoryfor lattice bosons.
Coupled rotators . The Hamiltonian of the rotator chain reads H CR = N (cid:88) j =1 (cid:0) p j + V ( ϕ j +1 − ϕ j ) (cid:1) (7.2)with periodic boundary conditions, ϕ N +1 = ϕ . At first glance we have only rebaptized q j as ϕ j . But the ϕ j ’s are angles and the p j ’s angular momenta. Hence the phase space is( S × R ) N with S denoting the unit circle. The standard choice for V is V ( ϑ ) = − cos ϑ ,but in our context any 2 π -periodic potential is admitted. The equations of motion are ddt ϕ j = p j , ddt p j = V (cid:48) ( ϕ j +1 − ϕ j ) − V (cid:48) ( ϕ j − ϕ j − ) . (7.3)40bviously angular momentum is locally conserved with the angular momentum current J ( j ) = − V (cid:48) ( ϕ j − ϕ j − ) . (7.4)As local energy we define e j = p j + V ( ϕ j +1 − ϕ j ). Then e j is locally conserved, since ddt e j = p j +1 V (cid:48) ( ϕ j +1 − ϕ j ) − p j V (cid:48) ( ϕ j − ϕ j − ) , (7.5)from which one reads off the energy current J ( j ) = − p j V (cid:48) ( ϕ j − ϕ j − ) . (7.6)For the angles, in analogy to the stretch, one defines the phase difference ˜ r j = Θ( ϕ j +1 − ϕ j ),where Θ is 2 π -periodic and Θ( x ) = x for | x | ≤ π . Because of the jump discontinuity thestretch is not conserved. A rotator chain has only two conserved fields.To apply nonlinear fluctuating hydrodynamics one has to compute the Euler currentsin local equilibrium. Since there are two conserved fields the canonical equilibrium statereads 1 Z N N (cid:89) j =1 exp (cid:2) − β (cid:0) ( p j − u ) + V ( ϕ j +1 − ϕ j ) (cid:1)(cid:3) dϕ j dp j (7.7)with u the average angular momentum. Now (cid:104)J ( j ) (cid:105) N = −(cid:104) V (cid:48) ( ϕ j − ϕ j − ) (cid:105) N , (cid:104)J ( j ) (cid:105) N = − u (cid:104) V (cid:48) ( ϕ j − ϕ j − ) (cid:105) N , (7.8)average with respect to the canonical ensemble (7.7). We claim thatlim N →∞ (cid:104) V (cid:48) ( ϕ j − ϕ j − ) (cid:105) N = 0 . (7.9)For this purpose we expand in Fourier series as e − V ( ϑ ) = (cid:88) m ∈ Z a ( m ) e − i mϑ , f ( ϑ ) = (cid:88) m ∈ Z ˆ f ( m ) e − i mϑ . (7.10)Then, working out all Kronecker deltas from the integration over the ϕ j ’s, one arrives at (cid:104) f ( ϕ j +1 − ϕ j ) (cid:105) N = (cid:16) (cid:88) m ∈ Z a ( m ) N (cid:17) − (cid:88) m ∈ Z a ( m ) N (cid:16) a ( m ) − (cid:88) (cid:96) ∈ Z ˆ f ( (cid:96) − m ) a ( (cid:96) ) (cid:17) . (7.11)Since a (0) > | a ( m ) | for all m (cid:54) = 0,lim N →∞ (cid:104) f ( ϕ j +1 − ϕ j ) (cid:105) N = a (0) − (cid:88) (cid:96) ∈ Z ˆ f ( (cid:96) ) a ( (cid:96) ) = 1 Z (cid:90) π − π dϑf ( ϑ ) e − βV ( ϑ ) . (7.12)For f ( ϑ ) = V (cid:48) ( ϑ ), the latter integral vanishes because of periodic boundary conditions in ϑ . We conclude that both currents vanish on average.41s before we consider the infinite lattice with thermal expectation (cid:104)·(cid:105) u,β and form theequilibrium time correlations as S αα (cid:48) ( j, t ) = (cid:104) g α ( j, t ) g α (cid:48) (0 , (cid:105) u,β − (cid:104) g α (0) (cid:105) u,β (cid:104) g (cid:48) α (0) (cid:105) u,β , (7.13) (cid:126)g ( j, t ) = (cid:0) p j ( t ) , e j ( t ) (cid:1) . p j is odd and e j is even under time reversal. Hence in the Green-Kubo formula the cross term vanishes. Thus, for large j, t , fluctuating hydrodynamicspredicts S αα (cid:48) ( j, t ) = δ αα (cid:48) (4 πD α t ) − / f G ((4 πD α t ) − / j ) , (7.14)where f G is the unit Gaussian. D α is the diffusion coefficient of mode α . Of course, it canbe written as a time-integral over the corresponding total current-current correlation, butits precise value has to be determined numerically. This has been done for the standardchoice V ( x ) = − cos x , to which we specialize now. For β = 1, with a lattice size N = 500,the diffusive peaks are well established at t = 2000 [73, 74]. Energy diffusion has beenconfirmed much earlier [75, 76].At low temperatures one finds a different, perhaps more interesting scenario. At zerotemperature, there is the one-parameter family of ground states with ϕ j = ¯ ϕ , p j = 0.When heating up, under the canonical equilibrium measure, the phase ϕ j jumps to ϕ j +1 with a jump size O (1 / √ β ). Next we have to understand how the conservation of ˜ r j field isbroken. In a pictorial language, the event that | ϕ j +1 ( t ) − ϕ j ( t ) | = π is called an umklappfor phase difference ˜ r j or an umklapp process to emphasize its dynamical character. Atlow temperatures a jump of size π has a small probability of order e − β ∆ V with ∆ V = 2the height of the potential barrier. Hence ˜ r j is locally conserved up to umklapp processesoccurring at a very small frequency only. This can be measured more quantitatively byconsidering the average Γ uk ( t ) = (cid:88) j ∈ Z (cid:0) (cid:104) ˜ r j ( t )˜ r (0) (cid:105) u,β − (cid:104) ˜ r (cid:105) u,β (cid:1) . (7.15)At β = 1, Γ uk ( t ) decays exponentially due to umklapp. But at β = 5 the decay rate isalready very much suppressed, see [74].In the low temperature regime it is tempting to use an approximation, where thepotential V ( x ) = − cos x is Taylor expanded at the minimum x = 0. But such procedurewould underestimate the regime of low temperatures, as can be seen from the example ofa potential, still with ∆ V = 2, but several shallow minima. The proper small parameteris β − such that β ∆ V >
1. To arrive at an optimal low temperature hamiltonian, wefirst parametrize the angles ϕ , . . . , ϕ N through r j = ϕ j +1 − ϕ j with r j ∈ [ − π, π ]. Todistinguish, we denote the angles in this particular parametrization by φ j . The dynamicsgoverned by H CR corresponds to periodic boundary conditions at r j = ± π . For a lowtemperature description we impose instead specular reflection, i.e., if r j = ± π , then p j , p j +1 are scattered to p (cid:48) j = p j +1 , p (cid:48) j +1 = p j . By fiat all umklapp processes are nowsuppressed, while between two umklapp events the CR dynamics and the low temperaturedynamics are identical. The corresponding hamiltonian reads H CR , lt = N (cid:88) j =1 (cid:16) p j + ˜ V ( φ j +1 − φ j ) (cid:17) (7.16)42ith ˜ V ( x ) = − cos x for | x | ≤ π , ˜ V ( x ) = ∞ for | x | > π , (7.17)periodic boundary conditions φ N +1 = φ being understood. The pair ( φ j , p j ) are canon-ically conjugate variables. Note that as weights exp[ − βH CR ] = exp[ − βH CR , lt ]. Thus allequilibrium properties of the coupled rotators remain untouched.The hamiltonian H CR , lt is a variant of the hard collision model with square well poten-tial as discussed before, see [59]. The dynamics governed by H CR , lt has three conservedfields, the stretch r j = φ j +1 − φ j , the momentum p j , and the energy e j = p j + ˜ V ( r j ).Because of φ = φ N +1 , one has (cid:80) Nj =1 r j = 0. The model is in the dynamical phasecharacterized by an even potential at zero pressure.We claim that, for β ∆ V >
1, the CR equilibrium time correlations are well approx-imated by those of H CR , lt , provided the time of comparison is not too long. The lattercorrelations can be obtained within the framework of nonlinear fluctuating hydrodynam-ics. Thereby one arrives at fairly explicit dynamical predictions for the low temperatureregime of the CR model.One physically interesting information concerns the Landau-Placzek ratio at low tem-peratures. We use (4.29) and expand R − in 1 /β . To lowest order it suffices to use theharmonic approximation, ˜ V ( x ) = ax , a = 1 for the cosine potential. Using the formulasin Appendix A of [29] one arrives at the following value for the Landau-Placzek ratios, r - r : (2 aβ ) − (1 , , , p - p : (2 β ) − (1 , , , e - e : (2 β ) − ( a(cid:96) , β − , a(cid:96) ) . (7.18)The correlations are small, order β − . For the stretch correlations there is no centralpeak, to this order, and for the energy correlations the central peak is down by a factor β − relative to the sound peaks.To have a unified picture we add ˜ r j to the list of fields of physical interest. At hightemperatures ˜ r j is not conserved and one has diffusive spreading of the conserved fields.At low temperatures ˜ r j is conserved up to small errors and the conventional three-peakstructure, including universal shape functions, results. For extremely long times umklappprocesses will happen and one expects that they force a cross over to the Gaussian scaling(7.14). The precise dynamical structure of such cross over still needs to be investigated. Nonlinear Schr¨odinger equation on a lattice . A further example with a dynam-ically distinct low temperature phase is the nonlinear Sch¨odinger equation on the one-dimensional lattice. In this case the lattice field is ψ j ∈ C , for which real and imaginarypart are the canonically conjugate fields. The Hamiltonian reads H = N (cid:88) j =1 (cid:0) | ψ j +1 − ψ j | + g | ψ j | (cid:1) (7.19)with periodic boundary conditions and coupling g >
0. The sign of the hopping term playsno role, since it can be switched through the gauge transformation ψ j (cid:59) e i πj ψ j . The chainis non-integrable and the locally conserved fields are the number density ρ j = | ψ j | and43he local energy e j = | ψ j +1 − ψ j | + g | ψ j | . Hence the canonical equilibrium state isgiven by Z − e − β ( H − µ N ) N (cid:89) j =1 dψ j dψ ∗ j , N = N (cid:88) j =1 | ψ j | , (7.20)with the chemical potential µ . We assume β >
0. But also negative temperature states, inthe microcanonical ensemble, have been studied [77, 78]. Then the dynamics is dominatedby a coarsening process mediated through breathers. In equilibrium, the ψ -field hashigh spikes at random locations embedded in a low noise background, which is verydifferent from the positive temperature states considered here. For them the densityand energy currents are symbolically of the form i( z − z ∗ ), hence their thermal averagevanishes. Both fields are expected to have diffusive transport. In fact, this is confirmedby MD simulations [79]. They also show Gaussian cross-correlations, which is possiblesince density and energy are both even under time reversal. In the previous studies [80]transport coefficients have been measured in the steady state set-up.To understand the low temperature phase, it is convenient to transform to the newcanonical pairs ρ j , ϕ j through ψ j = √ ρ j e i ϕ j . (7.21)In these variables the Hamiltonian becomes H = N (cid:88) j =1 (cid:0) − √ ρ j +1 ρ j cos( ϕ j +1 − ϕ j ) + ρ j + gρ j (cid:1) . (7.22)The equations of motion read then ∂ t ϕ j = − ∂ ρ j H , ∂ t ρ j = ∂ ϕ j H . (7.23) ϕ j takes values on the circle S and ρ j ≥
0. From the continuity of ψ j ( t ) when movingthrough the origin, one concludes that at ρ j ( t ) = 0 the phase jumps from ϕ j ( t ) to ϕ j ( t )+ π .One recognizes the similarity to the coupled rotators (7.2). But now the equilibriummeasure carries a nearest neighbor coupling. For µ >
0, in the limit β → ∞ the canonicalmeasure converges to the one-parameter family of ground states with ρ j = ¯ ϕ = µ/g , ϕ j = ¯ ϕ with ¯ ϕ uniformly distributed on S . At low temperatures the field of phase dif-ferences ˜ r j is approximately conserved. The low temperature hamiltonian is constructedin such a way that the equilibrium ensemble remains unchanged while all umklapp pro-cesses are suppressed. To achieve our goal we follow verbatim the CR blueprint. Thephases are parametrized such that ϕ j +1 − ϕ j lies in the interval [ − π, π ] and this particularparametrization denoted by φ j . Umklapp is a point at the boundary of this interval. Now( φ j , ρ j ) are a pair of canonically conjugate variables, only ρ j ≥ p j ∈ R . Thusthe proper low temperature hamiltonian reads H lt = N − (cid:88) j =0 (cid:16) √ ρ j +1 ρ j U ( φ j +1 − φ j ) + V ( ρ j ) (cid:17) , (7.24)44here U ( x ) = − cos( x ) for | x | ≤ π , U ( x ) = ∞ for | x | > π , (7.25)and V ( x ) = x + g x for x ≥ , V ( x ) = ∞ for x < . (7.26)The low temperature Hamiltonian has a nearest neighbor coupling, which complicatesthe scheme through which the (cid:126)G matrices are determined [79]. Progress is achievedthrough the miraculous identity (cid:126) j = ( µ, P, µP ) for the Euler currents. The (cid:126)G coefficientsare evaluated at P = 0. As for a generic anharmonic chain, one finds that G = 0 and G (cid:54) = 0. Thus the heat peak is predicted to be L´evy 5/3 and the sound peaks to beKPZ. The sound peak for the density-density correlation was first observed in [81] using k, ω space, see also [82]. In [79] we use normal mode representation, as explained in thisarticle. The sound peaks fit nicely with KPZ, but the normalized heat peak is very broadand noisy, still with a shape not unlikely L´evy 5/3. References α - β model. Master’sthesis, University of Florence[58] Das, S.G., Dhar, A., Saito, K., Mendl, C.B., Spohn, H. (2014) Numerical test ofhydrodynamic fluctuation theory in the Fermi-Pasta-Ulam chain. Phys. Rev. E 90:pp. 012124 4859] Mendl, C.B., Spohn, H. (2014) Equilibrium time-correlation functions for one-dimensional hard-point systems. Phys. Rev. E 90: pp. 012147[60] Delfini, L., Denisov, S., Lepri, S., Livi, R., Mohanty, P.K. (2007) Energy diffusion inhard-point systems. Eur. Phys. J. 146: pp. 21–35[61] Das, S.G., Dhar, A., Narayan, O. (2013) Heat conduction in the α - β -Fermi-Pasta-Ulam chain. J. Stat. Phys. 154: pp. 204–213[62] Benettin, G., Ponno, A., Christodoulidi, H. (2013) The Fermi–Pasta–Ulam problemand its underlying integrable dynamics. Journ. Stat. Phys. 152: pp. 195–212[63] van Beijeren, H., Posch, H., private communication, TU M¨unchen, June 2013[64] Popkov, V., Schmidt, J., Sch¨utz, G. (2014) Non-KPZ modes in two-species drivendiffusive systems. Phys. Rev. Lett. 112: pp. 200602[65] Arndt, P.F., Heinzel, T., Rittenberg, V. (1999) Spontaneous breaking of translationalinvariance and spatial condensation in stationary states on a ring. I. The neutralsystem. J. Stat. Phys. 97: pp. 1–65[66] Ferrari, P.L., Sasamoto, T., Spohn, H. (2013) Coupled Kardar-Parisi-Zhang equa-tions in one dimension. J. Stat. Phys. 153: pp. 377–399[67] Bernardin, C., Gon¸calves, P., Jara, M., 3 / φ4