Ballistic space-time correlators of the classical Toda lattice

The Toda lattice is an integrable system and its natural space-time stationary states are the generalized Gibbs ensembles (GGE). Of particular physical interest are then the space-time correlations of the conserved fields. To leading order they scale ballistically. We report on the exact solution of the respective generalized hydrodynamic equations linearized around a GGE as background state. Thereby we obtain a concise formula for the family of scaling functions.


Introduction
With the discovery of the Toda lattice as an integrable classical field theory an immediate hope was to be able to compute time-correlations in thermal equilibrium. During the 1980ies considerable efforts were invested [1][2][3][4][5][6][7][8][9][10]. But mostly a better understanding of static properties were accomplished [4][5][6][7]. A few at the time pioneering molecular dynamics simulations became available [1,3]. Theoretically low order continued fraction expansions were tried, also expansions close to the two solvable limits, namely harmonic chain and hard rods [2,9]. More recently, with improved computational power, molecular dynamics simulations with approximately 1000 particles have been performed [11]. Timedisplaced correlations as stretch-stretch (the same as free volume) and energy-energy are measured with high precision, which requires not only numerically solving a system of Newton's equations of motion but also generating roughly 10 7 samples so to keep the noise level small. As one main conclusion from these studies, the mentioned correlations scale ballistically over the accessible time span, except for very early times. Thus the challenging goal is to predict, resp. to compute, the scaling function for a specified correlation. Another important advance happened on the theoretical side. Starting with integrable quantum field theories the formalism of generalized hydrodynamics (GHD) has been developed [12,13], which in fact also applies to classical counterparts, in particular to the Toda lattice [14,15]. In the context of time correlations only the equations linearized around thermal equilibrium are of relevance. Some preliminary results are reported in [14,16]. In this contribution we provide the exact solution of linearized hydrodynamics. For quantitative results one still has to numerically determine the density of states (DOS) of the Lax matrix and solve the linear integral equation linked to the dressing transformation. To fully exploit the linear structure all conserved quantities have to be treated on equal footing. Also thermal equilibrium is just one particular case of a generalized Gibbs ensemble (GGE).
The connection between linear response and microscopic correlation functions is standard wisdom. In [17] this theme is studied in detail for a general class of integrable systems governed by GHD. For the one-dimensional sinh-Gordon model molecular dynamics is compared with the predictions from linearized GHD [18]. While our arguments follow a somewhat different route, our final result on time-correlations of the conserved fields is in complete agreement with the general theory developed in [17].
For the Toda lattice, thermal, more generally GGE, average of the conserved fields is most concisely encoded through the DOS of the Lax matrix. But on top there is the conserved stretch which does not seem to be connected to the Lax matrix and plays a special role. The stretch current equals momentum, hence also conserved, which is not the case for all the other conserved fields. As a consequence generalized hydrodynamics consists of two equations, one for the stretch and one for the Lax DOS. The linearized equations then inherit a 2 × 2 matrix structure (with operator entries) and one difficulty in the analysis is to properly handle this feature. Such property appears also, for example, in the XXZ chain in which case the exceptional role is played by the magnetization. The insight from the Toda lattice might thus be helpful for other integrable models.
A basic claim of generalized hydrodynamics is to be valid over the entire parameter range, no exceptions. This claim is particularly intriguing for the Toda chain. In thermal equilibrium, for large N, the length of the system behaves as νN with some coefficient ν(P ), which depends on the applied pressure P . ν can have either sign. With high probability, for small P particles are ordered, while for negative P they are anti-ordered. But this means that there is a specific pressure P c at which ν(P c ) = 0 and the typical distance between particles is of order 1/ √ N only. So the picture that the interactions can be broken up into small groups of particles, which are isolated from the rest of the system and move from an incoming configuration to an outgoing one, becomes questionable. On the Euler level the propagation speed of a perturbation is proportional to ν −1 , thus formally diverges. It is an open problem to understand the dynamical behavior close to ν = 0.
Our predictions are valid for the ballistic scale and one may wonder how much microscopic information is still visible. We recall that for a simple fluid the Euler equations have always the same structure. The particular fluid under consideration is encoded by the thermodynamic pressure as a function of density and internal energy. For an integrable many-body system, apparently the system specific character comes from the two-particle phase shift, whereas the mathematical structure persists over a large variety of models both quantum and classical. For the Toda lattice the phase shift is given by 2 log |v − v ′ | in dependence on the incoming velocities v, v ′ . The entire set of equations of generalized hydrodynamics could be written down by substituting some other function of the momentum transfer. Of course, there is then no reason to have an underlying microscopic particle model, the only confirmed cases being the Toda chain and hard rods, for which the phase shift equals the hard rod length a independently of v, v ′ .
In Sect. 2 we discuss some background material on the Toda chain. A summary on the generalized free energy is shifted to Appendix A. In Sect. 3 we work out a convenient form of the static GGE charge-charge and charge-current correlations. This allows us to state then the exact solution of linearized hydrodynamics with GGE random initial conditions. The required intermediate steps are reported in Sect. 4. There is an alternative way to obtain the static charge-charge correlator, which is related to Dyson's Brownian motion, see Sect. 5. Hard rods turn out to be a useful guiding example, well covered in the literature [19][20][21][22]. To be closer to the Toda lattice, in Sect. 6 we study hard rods as coming from an anharmonic chain with an infinite hard core potential of core size a. Thereby the structure of the corresponding generalized hydrodynamic equations resembles more closely the one of Toda.

Conserved charges, currents, and normal form
The Toda lattice is an anharmonic chain, for which the interaction potential is specified to be exponential. The corresponding hamiltonian is written as Here q j , p j is the position and momentum of the j-th particle. The positional increments, r j , are called stretches. We consider the infinitely extended lattice and indicate suitable finite number approximations whenever required. Then the equations of motion read which is viewed as a discrete nonlinear wave equation. Let us introduce the Flaschka variables [23], In fact the original Flaschka variables carry both a prefactor 1 2 and, in principle, any prefactor could be used. However, only with our choice the generalized free energy has a particularly simple form, see Appendix A. The Lax matrix is the tridiagonal real symmetric matrix with matrix elements (2.4) For the finite lattice [1, ...., N] the N eigenvalues of L are conserved. As functions on phase space they are non-local, their local version being tr[L n ], n = 1, 2, ... [24]. In generalized hydrodynamics such conserved fields are usually called charges or conserved charges and it is convenient to follow this practice. The locally conserved charges of the Toda lattice have a strictly local density given by with j ∈ Z. In addition, the stretch is conserved with density The respective current densities are of the form where L off denotes the off-diagonal part of L [16].
For the Lax matrix we adopted the convention L = 2L Flaschka . This may look like an arbitrary choice, but is so to speak dictated by the structure of the Toda generalized free energy. This point was slightly overlooked in [16], see also the arXiv version 4, and for convenience we briefly repeat in Appendix A.
As a next step we have to introduce a few standard items from generalized hydrodynamics. A generalized Gibbs state (GGE) of the Toda lattice is characterized by the pressure P > 0 and the chemical potential V (w), see Appendix A. P is the thermodynamic dual to the stretch and V (w) the dual to the collection of conserved charges, where instead of the discrete label n a continuous label w ∈ R will be more convenient. To compute averages of the charges and their currents one starts from the free energy functional (2.9) Introducing the Lagrange multiplier, µ, the Euler-Lagrange equation for the minimizer, Setting ρ µ = e −ε with quasi-energies ε, (2.10) turns into the classical version of the TBA equation, Let us define the integral operator Then the TBA equation is rewritten as and one introduces the dressing of a function ψ by On the right hand side ρ µ is regarded as multiplication operator, (ρ µ ψ)(w) = ρ µ (w)ψ(w). The DOS (density of states) of the Lax matrix under GGE is given by νρ p , such that where r 0 P,V stands for the average with respect to the GGE state labeled by the parameters P, V , while f is merely a short hand for R dwf (w). Differentiating TBA with respect to µ we conclude Here [1] stands for the constant function, ψ(w) = 1, and similarly [w n ] for the n-th power, ψ(w) = w n . We also use [w n ] dr2 = ([w n ] dr ) 2 . A tracer pulse with bare velocity v in fact travels with an effective velocity, which is determined through the integral equation More concisely the effective velocity can be written as With this input the GGE averaged charge densities are 0 P,V = ν ρ p w n = q n , n = 1, 2, ... , (2.19) with the corresponding average current densities The status of both identities is distinctly different. Under suitable conditions on the confining potential, (2.19) is established in [14,16], while (2.20) should be viewed as a conjecture, however with strongly supporting numerical evidence [25].
Remark on notations. In [14] B. Doyon investigates the generalized hydrodynamics for the Toda lattice, both in the fluid picture and the chain picture, the latter being adopted here. Our notations differ only minimally. More precisely the definitions of ν, T , v eff , and ρ p are identical. n in [14] becomes ρ µ here, because n is used already otherwise. ✷ Note that through (2.17) v eff becomes a functional of ρ p . Thus the Euler type equations for the Toda chain read We claim that the nonlinear transformation ρ p → ρ µ , defined by results in the normal form the continuity equation yields Using this identity together with (2.22), we obtain The effective velocity is solution of the integral equation Hence Here t refers to the Toda time evolution with the convention Q j (t = 0). Through a Galilei transformation one can always achieve q 1 = 0, which we adopt from now on. An exact solution is of out reach and we focus on the hydrodynamic prediction for this correlator. On the Euler scale the initial randomness is merely propagated. The noise resulting from the dynamics becomes visible only on the diffusive scale. Therefore the linearized equation has to be solved with random initial data as deduced from the GGE. On the hydrodynamic scale they are delta-correlated in space with non-trivial charge correlations as determined through the static charge-charge covariance Here we use the shorthand As to be explained in Sect. 4, for all m, n ≥ 1, For later computations the basis consisting of moments is not so convenient. In addition the index 0 plays a special role. For these reasons we introduce the two vector (r, φ) with r ∈ R and φ a real-valued function, formally φ(w) = ∞ n=1 a n w n with some real coefficients a n . The matrix then inherits a two-block structure. We define and where T stands for transpose and we freely use the Dirac notation for row and column vectors. F is a linear operator acting on the functions to the right. For example F T hF φ means first to compute F φ, then multiply the result by the function h, and finally act with F T . Note that F 1 = 0. We also introduce the static charge-current correlator On abstract grounds, see Sect. 4, the matrix B is symmetric. In addition the column m = 0 is already determined by C, since J [0] = −Q [0] . As to be discussed, one arrives at (3.9) Using [w] dr = v eff [1] dr , one obtains the more concise operator form

Main result
The particular structure of the matrices C, B allows for an educated guess of the full space-time correlators in the hydrodynamic approximation. As can be seen from the following section, it will require actually some efforts to confirm this guess. Let us start with a generic example of κ conservation laws linearized around a spatially homogeneous equilibrium state, which are governed by A n,m ∂ x u m (x, t) = 0, n = 1, ..., κ. (3.11) Here x ∈ R is the continuum approximation for the lattice Z. The κ × κ linearization matrix A is x-independent, generically not symmetric, but ensured to have a right and left system of eigenvectors, |ψ α , ψ α |, α = 1, ..., κ, and has real eigenvalues c α , α = 1, ..., κ.
(3.11) has to be solved with random initial conditions with mean zero and covariance where to distinguish from other averages we use E(·) for expectation over the initial noise. As before, C denotes the static correlator and δ(x) reflects the exponential decay of spatial correlations for the underlying microscopic model. Since the spatial part in (3.11) is merely a translation proportional to t, the space-time correlator is given by which has a simple interpretation: The α-th peak travels with velocity c α and has a weight depending on A and C. The static field-current correlator equals (3.14) Returning to the Toda chain, the charge-charge correlator C has been computed in (3.7) and charge-current B in (3.10). To be in agreement with (3.14), the natural guess for the full solution reads This expression is our main result, which is exact on the ballistic Euler scale. The right hand side for t = 0 equals δ(x)C and differentiating with respect to t at t = 0 yields the required identity B = AC. Still to be shown is the exponential form, e At C, of the correlator, see Section 4. Manifestly, the right hand side of Eq. (3.15) scales ballistically. The left hand side is the true correlator of the Toda chain which exhibits ballistic scaling only for sufficiently large j, t both of the same order. This is the meaning of ≃, where one should identify x with j up to a suitable scale factor.
(3.15) is an operator identity. For specific correlations one has to compute the respective matrix elements. As examples, we list the time-dependent stretch-stretch and momentum-momentum correlation, where in the second line q 1 = 0 has been used.

Towards the exact solution 4.1 The Cand B-matrix
We first collect a few identities, which will be used to compute the derivatives of the average charges and currents with respect to the pressure P and the chemical potential V (w). The most basic identity is which follows from expanding ρ p = ∂ µ ρ µ = (1−ρ µ T ) −1 ρ µ in a power series. Hence, setting Furthermore, for a general function f , The variational derivative with respect to V (w) is slightly more complicated, since the constraint (2.9) has to be respected. The corresponding directional derivative, oriented along g(w), is denoted by D g = dwg(w)δ/δV (w). Then, for general functions g, f,f, one obtains The average charges are 0 P,V = ν ρ p w n . (4.5) By suitable linear combinations the powers w n are replaced by a general function f (w).
The covariance matrix C is obtained by taking derivatives of (4.5) with respect to P and V .
and for C 0,n , n ≥ 1 we use (4.3) to arrive at as claimed in (3.7). For m, n ≥ 1 we use (4.4) for the special choicef = 1. Then (4.8) in agreement with (3.7). The B-matrix is symmetric. Since this property is a useful control check, for conciseness, we repeat the argument in our specific context, see also [20]. By stationarity S m,n (j, t) = S n,m (−j, −t) and hence j∈Z jS m,n (j, t) = − j∈Z jS n,m (j, −t). (4.9) Because of the conservation law, j (t)); Q Upon inserting in (4.9) one arrives at B m,n t = B n,m t.
Returning to Toda, the average currents are The second version will be more convenient for us. Note that while one can shift to q 1 = 0, when taking derivatives this term has to be kept. As before the powers w n are replaced by a general function f (w). Since J [0] = −Q [1] , the border matrix elements are already determined. As control we still repeat with the result The two middle terms vanish, since q 1 = 0, and the remainder agrees with (3.9). Finally, using (4.4) withf (w) = w, we turn to (4.14) Since q 1 = 0, we have obtained the (2, 2) matrix element of (3.10).

Linearized transformation to normal modes
In (2.22) we introduced a nonlinear map which transforms the system of conservation laws into its quasi-linear version in such a way that the linearized operator A is manifestly diagonal. One would expect this property to persist under linearization. In the ν, h variables the map is We linearize on both sides as ρ µ + ǫg, ν + ǫr, h + ǫφ, φ = 0. To first order in ǫ this yields the linear map R : g → (r, φ) given by where Note that indeed F T ψ = 0. (4.15) can be inverted as thereby deriving, by a similar argument as before, Indeed one checks that the first 1 standing for the identity operator as a 2 × 2 block matrix and the second 1 for the identity operator in the space of scalar functions. The next step is to write the C, B matrices in the new basis. Using one arrives at There are two cancellations which yield, as operators, Correspondingly for the B-matrix, As before there are two cancellations yielding In the new basis, as anticipated, R −1 AR is simply multiplication by ν −1 v eff . We conclude that e At C = Re (v eff /ν)t R −1 C. (4.26) Working out the algebra, one arrives at As explained in the beginning of Sect. 3.2, with this input one can write the solution to the hydrodynamic equations (2.21) linearized relative to some precsribed GGE and with random initial conditions having covariance matrix δ(x)C. The result is the right hand side of (3.15).
While not used, just for completeness the matrix A is recorded as (4.28)

The C-matrix from Dyson's Brownian motion
The variational problem (2.8) is closely linked to Dyson's Brownian motion with confining potential V . To make this section self-contained, some background material is required. The precise connection to Toda will be at the end of this section. We consider the stochastic particle system on R governed by dt + √ 2db j (t), j = 1, ..., N, α ≥ 0, (5.1) with {b j (t), j = 1, ..., N} a collection of independent standard Brownian motions. This is Dyson's Brownian motion in an external potential V . The interaction has strength α/N, which corresponds to a standard mean-field limit. (As to be discussed, the proper identification will α = P ). Let us introduce the empirical density If at the initial time t = 0, ρ N (x, 0) converges in the limit N → ∞ to a non-random density ρ(x, 0), then such a limit will hold for any t > 0 and the limit density satisfies the nonlinear Fokker-Planck equation where the effective potential is defined by with T is the linear operator defined in (2.12). For this convergence result even a proof is available [27]. If V is suitably confining, then Eq. (5.3) has a unique stationary solution ρ s , which can also be characterized as the minimizer of Let us now consider the stationary dynamics with ρ(x, 0) = ρ s (x) and study the fluctuations of the density. It is convenient to integrate against some smooth test function f . Then the scaled density fluctuations are For them a standard central limit theorem holds, compare with [28], where φ(x, t) is governed the linear Langevin equation Here ξ(x, t) is normalized space-time white noise and ∂ x D the linearized evolution operator with V ′ eff is defined as in (5.4) with ρ(x, t) substituted by ρ s (x). For the following argument V ′ eff and ρ s act as multiplication operators, while ∂ = ∂ x is the differentiation operator, which commutes with T , T ∂ = ∂T .
The stationary solution to Eq. (5.8) is a Gaussian measure, its covariance denoted by C ♯ , which is determined by where ·, · denotes the standard L 2 inner product, f, g = dxf (x)g(x). We claim that as an operator The second term ensures that there are no fluctuations in the number of particles, C ♯ [1] = 0.

(5.13)
Since ρ s is stationary, we have which when inserted in (5.13) leads to the condition It follows from using T ∂ = ∂T , thereby confirming our claim. ✷ The see the connection to the (2, 2) matrix element of (3.7), one recalls that because of linear ramp argument the covariance for the conserved charges is given by ∂ α (αC ♯ ), an expression depending only on αρ s . Thus differentiating with respect α becomes identical to differentiating with respect to P . Using this observation in (5.11), one arrives at the four terms defining the (2, 2) matrix element, thus providing an alternative derivation for this particular contribution. The stretch does not seem to be explicitly linked to the Lax matrix and for the (1, 2) matrix element of C one has to rely on the thermodynamic reasoning.

Hard rods viewed as lattice model
Hard rods are more naturally viewed as a one-dimensional fluid. To make the connection with the Toda chain, we discuss here the chain point of view. The rod length is denoted by a. The r j 's are changing linearly in time at constant p j 's. At the instant when r j = a with incoming momenta p j , p j+1 the momenta are simply interchanged, resulting in d dt r j > 0. The Euler equations could be easily guessed. But in our context it is more instructive to follow the systematic route outlined in Sect. 2. ν denotes the mean distance between hard rods, ν > a, and h(v) is the normalized velocity distribution. Thus ρ p = ν −1 h and u = hv . The phase shift equals −a and hence Since ρ p = [1] dr ρ µ , one arrives at and, using (2.18), v eff = (ν − a) −1 (νv − au). Thus the equations of GHD read with the normal mode transform The factor (ν − a) −1 is the equilibrium density at contact. The static covariance equals and the charge-current correlator, setting u = 0, The linearized operator A is explicit, The limit t → 0 and its first derivative yield C and B, respectively.

A Free energy
We briefly explain how the Toda generalized free energy is linked to the variational problem (2.8), more details being provided in [16]. For free boundary conditions, (L N ) 1,N = 0, and in the variables of (2.3), the Toda partition function reads where 1 2 w 2 +Ṽ (w) = V (w). Here P is the pressure and V the chemical potential. More precisely, the n-th charge is controlled by the chemical potential µ n . The grand-canonical weight is therefore the exponential of ∞ n=1 µ n Q [n] = ∞ n=1 µ n tr[L n ]. Introducing V (w) = − ∞ n=1 µ n w n the exponent can be written more concisely as −tr[V (L N )]. We compare Z toda with the Dumitriu-Edelman partition function,  with ρ * the minimizer of the mean field free energy functional. It is of advantage to absorb P into ̺ through P F MF P (P −1 ̺) = F (̺) − P log P . Then with F as in (2.8). Further properties are discussed in [16]. In particular F toda (P ) = µ(P ). (A.8)