Equilibrium Spacetime Correlations of the Toda Lattice on the Hydrodynamic Scale

We report on molecular dynamics simulations of spacetime correlations of the Toda lattice in thermal equilibrium. The correlations of stretch, momentum, and energy are computed numerically over a wide range of pressure and temperature. Our numerical results are compared with the predictions from linearized generalized hydrodynamics on the Euler scale. The system size is N=3000,4000\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=3000,4000$$\end{document} and time t=600\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t =600$$\end{document}, at which ballistic scaling is well confirmed. With no adjustable parameters, the numerically obtained scaling functions agree with the theory within a precision of less than 3.5%.


Introduction
A central goal of Statistical Mechanics is to explore the structure of equilibrium correlations for observables of physical interest.These could be static correlations, but more ambitiously also correlations in spacetime.An interesting, but very fine-tuned class of hamiltonians are integrable many-body systems, either classical or quantum.This choice restricts us to systems in one dimension.Then, generically, static correlations have exponential decay whether the model is integrable or not.However, the dynamics of correlations is entirely different.In nonintegrable chains correlations propagate as a few narrow peaks at constant speed which then show characteristic sub-ballistic broadening.On the other hand for integrable models correlations still spread ballistically but now with a broad spectrum of velocities.Such behaviour was confirmed through a molecular dynamics (MD) simulation of the Ablowitz-Ladik model [32], an integrable discretization of the nonlinear Schrödinger equation.A further confirmation came from the simulation of the Toda chain [22].On the theoretical side, the 2016 construction of generalized hydrodynamics (GHD) was an important breakthrough [3,6].This theory provides a powerful tool through which, at least in principle, the precise form of the spectrum of correlations can be predicted.With such a development MD simulations can also be viewed as probing the validity of GHD.
From the side of condensed matter physics, integrable quantum models have received considerable attention.Because of size limitations, the simulation of macroscopic profiles are preferred.But time correlations have also been studied through DMRG simulations [4,5,8,34].In recent years, attention has been given to the spacetime spin-spin correlation of the XXZ model at half-filling and at the isotropic point [10,20,25].The same quantity has also been investigated for a discrete classical chain with 3-spins of unit length and interactions such that the model is integrable [7].A comparable situation occurs for the classical sinh-Gordon equation, which is integrable as a nonlinear continuum wave equation and possesses an integrable discretization, see [2] for MD simulations for equilibrium time correlations of the discrete model.
In our contribution we study the correlations of the Toda chain in thermal equilibrium through MD simulations and compare with predictions from GHD.We will comment on the connection to [22] in the last section.To make our article reasonably self-contained we first discuss the Landau-Lifshitz theory for nonintegrable chains.This theory provides the connection between spacetime correlations and linearized hydrodynamics.For the Toda chain, the theory has to be extended so as to accommodate an infinite number of conserved fields.We report on MD simulations of the Toda chain and compare with linearized GHD.

Landau-Lifshitz theory
The dynamics of the Toda chain is governed by the Hamiltonian H = j∈Z 1 2 p 2 j + exp(−(q j+1 − q j )) , where (q j , p j ) ∈ R 2 are position and momentum of the j-th particle [43,44].Introducing the j-th stretch (free volume) through r j = q j+1 − q j , the equations of motion read By tradition, one introduces coefficients for the range and strength of the interaction potential through (g/γ) exp(−γ(q j+1 − q j )).However, by a suitable change of spacetime scales, the form (2) can be regained, see the discussion in Section 5.The Toda hamiltonian has no free parameters.Since the equilibrium measure for (1) is of product form, static correlations are easily accessible.Time correlations are more challenging, see [36,37] for early attempts.A novel approach has been developed, now known as GHD.The guiding idea is to first identify the hydrodynamic equations for the Toda chain, which by necessity are a set of nonlinear coupled hyperbolic conservation laws.
Given such an input one can construct the corresponding Landau-Lifshitz theory [13,24], as based on linearized GHD.Before entering into details, it will be useful to first recall the Landau-Lifshitz theory for a chain with a generic interaction potential, denoted by V (for the Toda lattice V (x) = e −x ), see [38] and references listed therein.Thus in (1) the interaction term reads V (q j+1 − q j ) and the equations of motion become To define spacetime correlations we first have to specify the random initial data modelling thermal equilibrium.By Galileian invariance one restricts to the case of zero average momentum.Then the Gibbs states are characterized by the inverse temperature β > 0 and a parameter P such that the physical pressure equals P/β.For simplicity, we will also refer to P as pressure.The allowed range of P depends on V .If V diverges faster than |x| for |x| → ∞, then P ∈ R. For the Toda lattice P > 0 because of the one-sided divergence of the exponential.In thermal equilibrium {(r j , p j ), j ∈ Z} are a collection of i.i.d.random variables with single site probability density Here Z 0 (P, β) is the normalizing partition function.Note that, with our convention, P and β appear linearly in the exponent.Expectations with respect to such i.i.d.random variables are denoted by • P,β .We also shorten the notation for the covariance through X 1 X 2 c P,β = X 1 X 2 P,β − X 1 P,β X 2 P,β , where the particular random variables X 1 , X 2 will be obvious from the context.
For general V , the conserved fields are stretch, momentum, and energy with densities using as shorthand V j = V (r j ).Q is a three-vector with components labeled by n = 0, 1, 2. The static space correlator is defined through and the static susceptibility by summing over space, Since the underlying measure is product, only the j = 0 term is nonvanishing and the zero entries resulting from p 0 P,β = 0, p 3 0 P,β = 0, and r 0 , p 0 being independent random variables.Later on we will need the statistics of the conserved fields on the hydrodynamic scale.More precisely, for smooth test functions f , we consider the random field Then, by the central limit theorem for independent random variables, where the limit field u(x) is a Gaussian random field on R with mean zero, E( u(x)) = 0, and covariance in other words, u(x) is Gaussian white noise with correlated components.Microscopically, spacetime correlations are defined by evolving one of the observables to time t which yields S m,n (j, t) = Q m (j, t)Q n (0, 0) c P,β .
Note that the Gibbs measure is spacetime stationary and thus without loss of generality both arguments in Q n in (7) can be taken as (0, 0).To understand the structure of S m,n one has to rely on approximations.For the long time ballistic regime a standard scheme is the Landau-Lifshitz theory, which views Q n (0, 0) as a small perturbation of the initial Gibbs measure at the origin.This perturbation will propagate and is then probed by the average of Q m at the spacetime point (j, t).For large (j, t) the microscopic dynamics is approximated by the Euler equations, but only in their linearized version since the perturbation is small.More concretely, the approximate theory will be a continuum field u(x, t) over R × R, which is governed by with random initial conditions as specified in (6).The 3 × 3 matrix A is constant, i.e. independent of (x, t).To explain the structure of A requires some further efforts.We refer to [38] for more details and proofs of the key identities.
From the equations of motion one infers that to each density Q n (j, t) there is a current density Explicitly, the current densities are where we adopted the convention that omission of time argument t means time 0 fields.One then defines the static current-conserved field correlator and the corresponding susceptibility Despite its asymmetric looking definition, As a general property, Euler equations are built on thermally averaged currents.Linearizing them with respect to the average fields yields Here B appears when differentiating the average currents with respect to the chemical potentials and C −1 when switching from intensive to extensive variables.By construction C = C T and C > 0, in addition B = B T according to (11).Hence which ensures that A has real eigenvalues and a complete set of left-right eigenvectors.Anharmonic lattices are symmetric under time reversal, which implies the eigenvalues c = (−c, 0, c), with c > 0 the isentropic speed of sound.We denote the right, resp.left eigenvectors of A by |ψ α and ψα |, α = 0, 1, 2. With this input the solution to (8) with initial conditions (6) reads with m, n = 0, 1, 2. There are three δ-peaks, the heat peak standing still and two sound peaks propagating in opposite directions with speed c.Specifying m, n, each peak has a signed weight which depends on C and the left-right eigenvectors of A.
The Landau-Lifshitz theory asserts that the microscopic correlator for j = xt , • denoting integer part, with t sufficiently large.The reader might be disappointed by the conclusion.But with such basic information the fine-structure of the peaks can be investigated, in particular their specific sub-ballistic broadening and corresponding scaling functions [31,38,39].When turning to the Toda lattice, the conservation laws are now labeled by n = 0, 1, ... and thus A, B, C become infinite dimensional matrices.The corresponding Landau-Lifshitz theory has been worked out in [40].As to be discussed in the following section, with appropriate adjustments Eq. ( 12) is still valid.

Toda lattice, linearized generalized hydrodynamics
The conservation laws of the Toda lattice are obtained from a Lax matrix [11,26].For this purpose, we first introduce the Flaschka variables

Then the equations of motion become
The Lax matrix, L, is defined by L j,j = p j , L j,j+1 = L j+1,j = a j , j ∈ Z, and L i,j = 0 otherwise.Clearly L = L T .The conserved fields are labelled by nonnegative integers and have densities given by with n ≥ 1.Note that Q n (j) is local in the sense that it depends only on the variables with indices in the interval [j − n, j + n].An explicit expression for these quantities is given in [15].For the current densities one obtains where L ↓ is the lower triangular part of L. Then under the Toda dynamics which is the n-th conservation law in local form.
The first items in the list are stretch and momentum for which our current definitions agree with those in (4), (9).However, for n = 2 one obtains ( 4), ( 9) on two accounts.First there is the trivial factor of 2. In our numerical plots we use the physical energy density e j .The second point is more subtle.Densities are not uniquely defined, since one can add a difference of some local function and its shift by one.When summing a particular choice for the density over some spatial interval, the result differs from another choice of the density by a boundary term only.Thus the bulk term will have a correction of order 1/(length of interval), which does not affect the hydrodynamic equations.For the currents the difference can be written as a total time derivative which is again a boundary term when integrating over some time interval.In this section we adopt the conventions ( 14) and (15), since the analysis heavily relies on the Lax matrix.Beyond n = 2, while the fields no longer have a name, they still have to be taken into account in a hydrodynamic theory.
The infinite volume static field-field correlator is defined as in ( 5) and the current-field correlator as in (10).In particularly, B = B T .Of course, C, B are now matrices in the Hilbert space of sequences indexed by N 0 , i.e. the space 2 (N 0 ).To distinguish 3 × 3 matrices from their infinite dimensional counterparts, for the latter we use standard italic symbols.The spacetime correlator of the Toda lattice is defined by and we plan to construct its Landau-Litshitz approximation.In essence this amounts to an analysis of e At C m,n , While we are mainly interested in the physical fields corresponding to the indices m, n = 0, 1, 2, for the operator in (17) an understanding of the infinite dimensional matrices is required.
Starting from the basics, the free energy of the Toda lattice is given by F eq (P, β) = log β/2π + P log β − log Γ(P ).
In particular, the average stretch, ν, is determined through with ψ the digamma function.Expectations of higher order fields can be written as moments of a probability measure denoted by νρ p , n ≥ 1. ρ p is called particle density.To determine this density one first has to solve the thermodynamic Bethe equations (TBA).For this purpose we introduce the integral operator w ∈ R, considered as an operator on L 2 (R, dw) and define the number density with quasi-energies ε.The quasi-energies satisfy the TBA equation where the chemical potential µ has to be adjusted such that Thereby the number density depends on the parameters P and β.
The TBA equation is closely connected to the β-ensemble of random matrix theory.We rewrite (21) as As α → ∞, the entropy term on the lefthand side can be neglected and one recognizes the defining equation for the Wigner semi-cirle law on the interval [−2 √ P , 2 √ P ].The Lax DOS is the Pderivative of ρ n , which diverges as (w ± 2 √ P ) −1/2 at the two borders.As α is lowered the borders become smeared to eventually cross over to a Gaussian.
In practice, the TBA equation has to be solved numerically.But for thermal equilibrium an exact solution is available [1,12,35].Denoting the solution of (21) for β = 1 and the constraint (22) by ρ * n one has In our numerical simulations it is of advantage to use the exact solution.
The TBA equation is a standard tool from GHD as one way to write the Euler-Lagrange equations for the variational principle associated with the generalized free energy.For the Toda lattice such a variational formula was obtained in [9,42].Proofs using methods from the theory of large deviations and transfer operator method have also become available [16,27,29,30].
Next we introduce the dressing transformation of some function f by with ρ n regarded as a multiplication operator.Then number and particle density are related as with inverse using the convention ς n (w) = w n .For the average currents similar identities are available.The central novel quantity is the effective velocity see [3,6,41,45].Then J 0 (0) P,β = −κ 1 , and, for n ≥ 1, In thermal equilibrium we have κ 1 = 0. Since in the following there will be many integrals over R, let us first introduce the abbreviation With this notation the C matrix turns out to be of the form in the sense that C m,n for m, n ≥ 1 follows a simple pattern.This structure will reappear for B and e At C. The field-current correlator B can be computed in a similar fashion with the result As in (12), we want to determine the propagator of the Landau-Lifshitz theory, denoted by S LL m,n (x, t).In principle, all pieces have been assembled.However to figure out the exponential of A requires its diagonalization.Details can be found in [40] and we only mention that one constructs a linear similarity transformation, R, such that R −1 AR is multiplication by in L 2 (R, dw).Here v eff is the effective velocity defined in (26).Using the block convention as in (28), the spacetime correlator in the Landau-Lifshitz approximation is given by .
Note that S LL (x, 0) = δ(x)C.As a property of the Euler equations, the expression (31) possesses exact ballistic scaling, The correlator S m,n (j, t) is computed in our MD simulations which will then be compared with S LL m,n (x, t).

Numerical simulations
For a molecular dynamics simulation one has to first specify a finite ring [1, . . ., N ] with suitable boundary conditions.For the dynamics of positions q j and momenta p j one imposes The parameter ν fixes the free volume per particle and can have either sign.In our simulation, we actually allow for a fluctuating free volume by choosing random initial conditions such that {r 1 , p 1 , . . ., r N , p N } are i.i.d.random variables with a single site distribution as specified in (3).
Then the deterministic time evolution is governed by (13) with boundary conditions In fact, the boundary condition in (33) amounts to the micro-canonical constraint If one sets ν = Q 0 (0) P,β , then for large N , by the equivalence of ensembles, the two schemes for sampling the correlator S m,n (j, t) should differ by the size of statistical fluctuations.For a few representative examples we checked that indeed the equivalence of ensembles holds for the particular observables under study.
Returning to the choice of system size there is an important physical constraint.In all simulations one observes a sharp right and left front, which travel with constant speed and beyond which spatial correlations are exponentially small.On a ring necessarily the two fronts will collide after some time.Such an encounter has a noticeable effect on the molecular dynamics which is not captured by the linearized GHD analysis.Therefore the simulation time is limited by the time of first collision.Indeed, we note in Figures 1-3 that both linearized GHD and MD clearly display maximal speeds of at most ∆j/∆t = 2 for the entire range of (P, β, m, n) displayed in these figures.Taking into account that the initial correlations are proportional to δ 0j , we conclude that for a ring of size N = 3000 there will be no collision of the two fronts up to time t = 750 which is larger than t = 600 used in our simulations.
Before displaying and discussing our results, we provide more details on numerically solving the TBA equations and on the actual scheme used for MD.

Details of the numerical implementation 4.1.1 Solving linearized GHD
To numerically solve the linearized GHD equations, we use a numerical method similar to the one from [33].First, Eq. ( 23) can be expressed in terms of the parabolic cylinder function D ν (z), which is readily available in Mathematica.This provides the solution to the TBA equations ( 21), (22).
Then, we use a simple finite element discretization of the w-dependent functions by hat functions, resulting in piecewise linear functions on a uniform grid.After precomputing the integral operator T in (20) for such hat functions, the dressing transformation (24) becomes a linear system of equations, which can be solved numerically.This procedure yields ς dr n , and subsequently ρ p via (25) and v eff via (26).The moments can be computed from κ n = R dwνρ n (w)ς dr n (w), or (equivalently) Eq. (19).
To evaluate the correlator in (31), we note that the delta-function in the integrand results in a parametrized curve, with the first coordinate (corresponding to x/t) equal to ṽeff from (30), and the second coordinate equal to the remaining terms in the integrand divided by the Jacobi factor | d dw ṽeff (w)| resulting from the delta-function.

Molecular dynamics simulations
We approximate the expectation value that is contained in the MD-definition of the correlations S m,n in equation ( 16) by the following numerical scheme, whose implementation program is written in Python, and can be found at [28].First, we generate the random initial conditions distributed according to the Gibbs measure, as given by (3) for the i.i.d.random variables (r j , p j ) 1≤j≤N .Specifically, the variables p j are distributed according to a standard normal random variable, that we generate with Numpy v1.23's native function random.defaultrng().normal[18], times 1/ √ β.It takes a brief calculation to see that r j can be chosen to be − ln(X/(2β)) where X is chi-square distributed with shape parameter 2P .We obtain the random variable X using Numpy v1.23's native function random.defaultrng().chisquare.Having chosen the initial conditions in such a manner, we solve equation (2).
For the evolution, we adapt the classical Störmer-Verlet algorithm [17] of order 2 to work with the variables (p, r).Specifically, we used a time step equal to δ = 0.05, and, given the solution (r(t), p(t)) at time t, we approximate the solution at time t + δ through the following scheme, for all j = 1, . . ., N .In this part of the implementation, we extensively used the library Numba [23] to speed up the computations.Our approximation for the expectation S m,n is then extracted from 3 × 10 6 trials with independent initial conditions.Here we take the empirical mean of all trials where for each trial we also take the mean of the N = 3000 sets of data that are generated by choosing each site on the ring for j = 0.
To evaluate the quality of our numerical simulations, we have repeated the numerical experiments up to five times including variations for the length of the ring and evaluating the solutions at more intermediate time steps than displayed in the figures below.Furthermore, we have compared the results with the corresponding outcomes obtained by a MATLAB program that has been developed independently from the Python program, and that follows a different numerical scheme.It uses MATLAB's random number generators randn for initial momenta and rand combined with the rejection method to produce initial stretches.The dynamics is then evaluated by the solver ode45, which exploits the Runge-Kutta method to numerically solve the Hamiltonian system associated with (1) on the ring.We found that the deviations between different experiments are comparable to the size of the amplitudes of the high frequency oscillations that are present in figures 4-5.These oscillations are due to the random fluctuations of the empirical means around their expectation values S m,n .Agreement of different experiments up to the order of these oscillations therefore shows the consistency of the corresponding numerical results.
We also want to mention that all the pictures that appeared in this paper are made using the library matplotlib [19].

Comparison of linearized GHD with MD at time t = 600
We compare the GHD predictions with MD simulations for three different temperatures that correspond to β = 0.5 (Fig. 1), β = 1 (Fig. 2), and β = 2 (Fig. 3).For each β we choose three different values for the pressure parameter P in such a way that the corresponding mean stretches, given by ( 18), are positive (≈ 2.57) for low pressure, negative (≈ −0.42) for high pressure and approximately zero for medium pressure.We summarize their values in Table 1.In each of the nine cases we have evaluated the Landau-Lifshitz approximations S LL m,n (•, 1), see (31), of the correlators for all 0 ≤ n ≤ m ≤ 2 using the numerical scheme described in Section 4.1.1.Their graphs are displayed in Figures 1-3 as dashed lines.Note that the speeds of the sound peaks depend significantly on both pressure and temperature.Moreover, the predicted fine-structure of both the heat and the sound peaks are quite different for low pressure when compared to medium and high pressure.
The colored lines in Figures 1-3 show our numerical results for the corresponding molecular dynamics.According to the predicted ballistic scaling (32) we plot tS m,n (j, t) as a function of j/t for t = 600.Here the values of S m,n (j, t) are approximated using the numerics explained in Section 4. 1.2.The agreement between linearized GHD and MD is striking, in particular since there are no adjustable parameters.In all of the 54 comparisons shown in Figures 1-3 the GHD predictions for the fine-structure of heat and sound peaks are in excellent agreement with the ones observed from molecular dynamics at time t = 600.As we show in more detail in the next subsection the largest deviations occur mostly near the sound peaks and do not exceed 3.5% of the peaks' maximal values.

Deviation of linearized GHD from MD at times t = 150 and t = 600
The purpose of this subsection is twofold.On the one hand we have a look at the small differences between GHD predictions and molecular dynamics simulations that can hardly be detected in left) and S 1,0 (right) for medium pressure and increasing temperatures (top to bottom).For each value of β and P the top panels show GHD prediction vs. numerical simulations as in Figures 1-3 but with the the molecular dynamics evaluated at two times t = 150 and t = 600.The bottom panels display the differences between the GHD prediction and numerical simulations at time t = 150 (red) and at time t = 600 (green).Number of trials: 3 × 10 6 .For each value of β and P the top panels show GHD prediction vs. numerical simulations as in Figure 2 but with the the molecular dynamics evaluated at two times t = 150 and t = 600.The bottom panels display the differences between the GHD prediction and numerical simulations at time t = 150 (red) and at time t = 600 (green).Number of trials: 3 × 10 6 .
Figures 1-3.On the other hand we indicate how these differences evolve in time by including time t = 150 for the molecular dynamics.Recall that the GHD predictions are time-invariant in the scaling y → tS m,n (yt, t) we have chosen, see (32).
From the 54 comparisons that are displayed in Figures 1-3 we select 12 cases that are representative and show all the phenomena that we have observed.In Figure 4 we consider correlations S 1,1 and S 1,0 at medium pressure (cf.Table 1) for all three values of β.The small scale fluctuations displayed in the bottom panels are due to the approximation of expectation values by empirical averages.Their amplitudes become smaller if one increases the number of trials.Note that the difference in amplitudes of these fluctions between t = 150 and t = 600 is mostly due to the scaling y → tS m,n (yt, t) that we use.This implies that the values of the correlations are multiplied by a factor that is 4 times larger at the later time.The same holds for the plots in Figure 5 where the correlations S 0,0 and S 2,0 are shown for fixed β = 1 and our three different choices for pressure.We now summarize our main findings: 1.The deviations occur mostly near the sound peaks and amount to 1.5%-3.5% of the peaks' maximal values at time t = 600.
2. There appear to be small but systematic deviations concerning the shape of the sound peak in all cases.One would need to conduct experiments with a higher resolution, i.e. more sites and consequently larger times and more trials, to determine whether there is indeed such a systematic deviation.With the resolution present in our experiments the question of a systematic deviation with respect to the shape of the peak cannot be decided.
3. In some of the experiments the maximal deviations would be significantly smaller if a constant only depending on the values of β, P , m, n is added to all values of S (j, t), see e.g.correlations S 0,0 and S 2,0 for β = 1, P = 0.4 in Figure 5.This seems to be related to the approximation errors for the means r , p , and e , that appear to be less pronounced in the case of momentum p.We have observed that these deviations decrease as the number of trials is increased and we do not expect a systematic deviation between GHD and MD in this respect.
4. For (β; P ) ∈ {(0.5; 0.95), (0.5; 1.21)} we observe that the size of the deviations is essentially the same for times t = 150 and t = 600 whereas for (β; P ) ∈ {(0.5; 0.32), (1; 0.4), (2; 0.52), (2; 2.55), (2; 3.53)} these deviations are significantly larger at the smaller time.The remaining two cases (β; P ) ∈ {(1; 1.5), (1; 2)} are somewhat in between, also depending on the correlation function that is considered, see Figure 5.This is an indication that the speed of convergence of tS m,n (yt, t) to the GHD prediction S LL m,n (y, 1) as t → ∞ depends on the values of β and P .As a rule we have observed that both increasing temperature or increasing pressure leads to a faster speed of convergence.

Conclusions and outlook
As can be seen from Table 1, we picked the intermediate pressure such that ν 0. In the particle picture ν = 0 corresponds to the boundary condition q 1 = q N .In thermal equilibrium the positions then perform an unbiased random walk with typical excursions of order √ N .Thus the free volume is of order 1/ √ N .The particles are extremely dense and the picture of successive pair collisions breaks down completely.So one might wonder whether GHD is still valid under such extreme conditions.ν = 0 poses no particular difficulties for MD simulations.In GHD the factor 1/ν appears in the expression for v eff , see Eq. ( 31).This makes the numerical scheme slow and only values close to ν = 0 are accessible.However the correlator S changes smoothly through ν = 0. GHD also covers this seemingly singular point.
Simultaneously A. Kundu [21] posted a somewhat puzzling note.He considers the parameter values β = 1, P = 1.When cutting the matrices C m,n and A m,n at low orders, the resulting S m,n consists of a few δ-peaks which move at constant velocity.After ballistic scaling, with high precission they turn out to lie on the curve obtained from GHD.A theoretical explanations seems to be missing.
In [22] the molecular dynamics of Toda lattice correlations are simulated for the potential with arbitrary γ, g > 0. To distinguish their parameters from ours, the variables in [22] are here denoted by t, r, P , β. P is the physical pressure and, comparing the Gibbs weights, one obtains the relations β = g γ From the equations of motions one deduces t = 1 √ γg t, r(t) = γr( t), p(t) = g γ p( t).
Thus, translating to our units, the MD simulations reported in [22] are (i) P = 0.01, β = 0.01, N = 1024, t = 400, (ii) P = 1, β = 1, N = 1024, t = 200, 300, and (iii) P = 400, β = 400, N = 256, t = 80.In fact, in all three cases the time scales are identical, t = t.Since GHD was not available yet, no comparison could have been attempted.Case (i) is a very dilute chain.In this limit νρ p is a unit Gaussian.The dressed functions become polynomials as ς dr 0 (w) = a 0 , ς dr 1 (w) = a 1 w, and ς dr 2 (w) = a 2 w 2 + a 3 with coefficients a 0 , ..., a 3 depending on (P, β).Note that for a noninteracting fluid a 3 would vanish.As a result S 0,0 is Gaussian, S 1,1 has two peaks, and S 2,2 has either two or three peaks.This is in good agreement with [22] and explains our motivation not to venture into the low density regime.Case (ii) interpolates between our β = 1, P = 0.40 and β = 1, P = 1.5.Note that now S 0,0 has a local minimum at w = 0, which is very different from the structure in the dilute regime.On the other hand, S 2,2 has a local maximum at w = 0, as is the case for low density/high temperature.
The most interesting parameter value is (iii), which deserves more detailed studies.The issue is the behavior of the Toda chain at very low temperatures.Simply letting β → ∞ will freeze any motion.But the simultaneous limit β → ∞ with P = P β at fixed physical pressure P is meaningful, at least statistically.In this limit ν > 0 always.Also the density of states converges to the arcsine distribution, To understand the dynamical behavior, the effective potential is expanded as e −r + P r 1 2 P (r − r 0 ) 2 + c 0 at its minimum r 0 .Since β is large, the initial fluctuations are of order 1/ √ β.Therefore the dynamics can be approximated by a harmonic chain with ω 2 = P .The equilibrium time correlations of the harmonic chain have intricate oscillatory behavior [14], which in the large β limit should match with the Toda lattice, as partially evidenced through case (iii).Clearly, GHD cannot reproduce such fine details.Still, when averaged on suitable scales, the gross behavior of the harmonic chain oscillations might be visible.

Figure 2 :
Figure 2: Toda correlation functions: GHD predictions y → S LL m,n (y, 1) vs. numerical simulations of the molecular dynamics y → tS m,n (yt, t) at t = 600 for β = 1.0 with low pressure (top), medium pressure (middle) and high pressure (bottom).Numerical simulations are colored according to the legend, the corresponding GHD predictions are displayed by dashed lines.Number of trials: 3 × 10 6 .

Figure 3 :
Figure 3: Toda correlation functions: GHD predictions y → S LL m,n (y, 1) vs. numerical simulations of the molecular dynamics y → tS m,n (yt, t) at t = 600 for β = 2.0 with low pressure (top), medium pressure (middle) and high pressure (bottom).Numerical simulations are colored according to the legend, the corresponding GHD predictions are displayed by dashed lines.Number of trials: 3 × 10 6 .

Figure 4 :
Figure4: Toda correlation functions S 1,1 (left) and S 1,0 (right) for medium pressure and increasing temperatures (top to bottom).For each value of β and P the top panels show GHD prediction vs. numerical simulations as in Figures1-3but with the the molecular dynamics evaluated at two times t = 150 and t = 600.The bottom panels display the differences between the GHD prediction and numerical simulations at time t = 150 (red) and at time t = 600 (green).Number of trials: 3 × 10 6 .

Figure 5 :
Figure5: Toda correlation functions S 0,0 (left) and S 2,0 (right) for β = 1 and increasing pressure (top to bottom).For each value of β and P the top panels show GHD prediction vs. numerical simulations as in Figure2but with the the molecular dynamics evaluated at two times t = 150 and t = 600.The bottom panels display the differences between the GHD prediction and numerical simulations at time t = 150 (red) and at time t = 600 (green).Number of trials: 3 × 10 6 .

Table 1 :
Values for β and P and the corresponding mean stretches used in experiments