Equilibrium Spacetime Correlations of the Toda Lattice on the Hydrodynamic Scale

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 [5, 4, 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 [20, 25, 10]. 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.

The dynamics of the Toda chain is governed by the Hamiltonian

$$H=\sum_{j\in\mathbb{Z}}\big{(}\tfrac{1}{2}p^{2}_{j}+\exp(-(q_{j+1}-q_{j}))\big{)},$$

(1)

where $(q_{j},p_{j})\in\mathbb{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

$$\frac{d}{dt}r_{j}=p_{j+1}-p_{j}\,,\qquad\frac{d}{dt}p_{j}=-\mathrm{e}^{-r_{j}}+\mathrm{e}^{-r_{j-1}},\qquad j\in\mathbb{Z}.$$

(2)

By tradition, one introduces coefficients for the range and strength of the interaction potential through $(g/\gamma)\exp(-\gamma(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 [37, 36] 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 [24, 13], 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)=\mathrm{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

$$\frac{d}{dt}r_{j}=p_{j+1}-p_{j}\,,\qquad\frac{d}{dt}p_{j}=V^{\prime}(r_{j})-V^{\prime}(r_{j-1}).$$

(3)

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 $\beta>0$ and a parameter $P$ such that the physical pressure equals $P/\beta$. 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|\to\infty$, then $P\in\mathbb{R}$. For the Toda lattice $P>0$ because of the one-sided divergence of the exponential. In thermal equilibrium $\{(r_{j},p_{j}),j\in\mathbb{Z}\}$ are a collection of i.i.d. random variables with single site probability density

$$Z_{0}(P,\beta)^{-1}\exp\!\big{(}-\beta\big{(}\tfrac{1}{2}p_{0}^{2}+V(r_{0})\big{)}-Pr_{0}\big{)}.$$

(4)

Here $Z_{0}(P,\beta)$ is the normalizing partition function. Note that, with our convention, $P$ and $\beta$ appear linearly in the exponent. Expectations with respect to such i.i.d. random variables are denoted by $\langle\cdot\rangle_{P,\beta}$. We also shorten the notation for the covariance through $\langle X_{1}X_{2}\rangle_{P,\beta}^{\mathrm{c}}=\langle X_{1}X_{2}\rangle_{P,\beta}-\langle X_{1}\rangle_{P,\beta}\langle X_{2}\rangle_{P,\beta}$, 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

$$\vec{Q}(j)=\big{(}r_{j},p_{j},e_{j}\big{)},\qquad e_{j}=\tfrac{1}{2}p^{2}_{j}+V_{j},$$

(5)

using as shorthand $V_{j}=V(r_{j})$. $\vec{Q}$ is a three-vector with components labeled by $n=0,1,2$. The static space correlator is defined through

$$\mathsfit{C}_{m,n}(j)=\langle Q_{m}(j)Q_{n}(0)\rangle_{P,\beta}^{\mathrm{c}}$$

(6)

and the static susceptibility by summing over space,

$$\mathsfit{C}_{m,n}=\sum_{j\in\mathbb{Z}}\langle Q_{m}(j)Q_{n}(0)\rangle_{P,\beta}^{\mathrm{c}},$$

(7)

$m,n=0,1,2$. Since the underlying measure is product, only the $j=0$ term is nonvanishing and

$$\mathsfit{C}=\begin{pmatrix}\langle r_{0}r_{0}\rangle_{P,\beta}^{\mathrm{c}}&0&\langle r_{0}e_{0}\rangle_{P,\beta}^{\mathrm{c}}\\[2.15277pt] 0&\langle p_{0}p_{0}\rangle_{P,\beta}^{\mathrm{c}}&0\\[2.15277pt] \langle r_{0}e_{0}\rangle_{P,\beta}^{\mathrm{c}}&0&\langle e_{0}e_{0}\rangle_{P,\beta}^{\mathrm{c}}\end{pmatrix},$$

(8)

the zero entries resulting from $\langle p_{0}\rangle_{P,\beta}=0$, $\langle p_{0}^{3}\rangle_{P,\beta}=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

$$\vec{\xi}_{\epsilon}(f)=\sqrt{\epsilon}\sum_{j\in\mathbb{Z}}f(\epsilon j)\big{(}\vec{Q}(j)-\langle\vec{Q}(0)\rangle_{P,\beta}\big{)}.$$

(9)

Then, by the central limit theorem for independent random variables,

$$\lim_{\epsilon\to 0}\vec{\xi}_{\epsilon}(f)=\int_{\mathbb{R}}\mathrm{d}xf(x)\vec{u}(x),$$

(10)

where the limit field $\vec{u}(x)$ is a Gaussian random field on $\mathbb{R}$ with mean zero, $\mathbb{E}(\vec{u}(x))=0$, and covariance

$$\mathbb{E}(u_{m}(x)u_{n}(x^{\prime}))=\mathsfit{C}_{m,n}\delta(x-x^{\prime}),$$

(11)

in other words, $\vec{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

$$\mathsfit{S}_{m,n}(j,t)=\langle Q_{m}(j,t)Q_{n}(0,0)\rangle_{P,\beta}^{\mathrm{c}}.$$

(12)

Note that the Gibbs measure is spacetime stationary and thus without loss of generality both arguments in $Q_{n}$ in (12) can be taken as $(0,0)$. To understand the structure of $\mathsfit{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 $\vec{u}(x,t)$ over $\mathbb{R}\times\mathbb{R}$, which is governed by

$$\partial_{t}\vec{u}(x,t)+\mathsfit{A}\partial_{x}\vec{u}(x,t)=0\,,$$

(13)

with random initial conditions as specified in (11). The $3\times 3$ matrix $\mathsfit{A}$ is constant, i.e. independent of $(x,t)$. To explain the structure of $\mathsfit{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 $J_{n}(j,t)$ such that

$$\frac{d}{dt}Q_{n}(j,t)+J_{n}(j+1,t)-J_{n}(j,t)=0.$$

(14)

Explicitly, the current densities are

$$\vec{J}(j)=-(p_{j},V^{\prime}_{j-1},p_{j}V^{\prime}_{j-1}),$$

(15)

where we adopted the convention that omission of time argument $t$ means time $0$ fields. One then defines the static current-conserved field correlator

$$\mathsfit{B}_{m,n}(j)=\langle J_{m}(j)Q_{n}(0)\rangle_{P,\beta}^{\mathrm{c}},$$

(16)

and the corresponding susceptibility

$$\mathsfit{B}_{m,n}=\sum_{j\in\mathbb{Z}}\langle J_{m}(j)Q_{n}(0)\rangle_{P,\beta}^{\mathrm{c}}.$$

(17)

Despite its asymmetric looking definition,

$$\mathsfit{B}_{m,n}=\mathsfit{B}_{n,m}.$$

(18)

As a general property, Euler equations are built on thermally averaged currents. Linearizing them with respect to the average fields yields

$$\mathsfit{A}=\mathsfit{B}\mathsfit{C}^{-1}.$$

(19)

Here $\mathsfit{B}$ appears when differentiating the average currents with respect to the chemical potentials and $\mathsfit{C}^{-1}$ when switching from intensive to extensive variables. By construction $\mathsfit{C}=\mathsfit{C}^{\mathrm{T}}$ and $\mathsfit{C}>0$, in addition $\mathsfit{B}=\mathsfit{B}^{\mathrm{T}}$ according to (18). Hence

$$\mathsfit{A}=\mathsfit{C}^{1/2}\mathsfit{C}^{-1/2}\mathsfit{B}\mathsfit{C}^{-1/2}\mathsfit{C}^{-1/2},$$

(20)

which ensures that $\mathsfit{A}$ has real eigenvalues and a complete set of left-right eigenvectors. Anharmonic lattices are symmetric under time reversal, which implies the eigenvalues $\vec{c}=(-c,0,c)$, with $c>0$ the isentropic speed of sound. We denote the right, resp. left eigenvectors of $\mathsfit{A}$ by $|\psi_{\alpha}\rangle$ and $\langle\tilde{\psi}_{\alpha}|$, $\alpha=0,1,2$. With this input the solution to (13) with initial conditions (11) reads

	$\displaystyle\mathsfit{S}^{\mathrm{LL}}_{m,n}(x,t)=\mathbb{E}\big{(}u_{m}(x,t)u_{n}(0,0)\big{)}$
	$\displaystyle\hskip 31.0pt=(\delta(x-\mathsfit{A}t)\mathsfit{C})_{m,n}=\sum_{\alpha=0}^{2}\delta(x-c_{\alpha}t)(|\psi_{\alpha}\rangle\langle\tilde{\psi}_{\alpha}|\mathsfit{C})_{m,n}$		(21)

with $m,n=0,1,2$. There are three $\delta$-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 $\mathsfit{C}$ and the left-right eigenvectors of $\mathsfit{A}$.

The Landau-Lifshitz theory asserts that the microscopic correlator

$$\mathsfit{S}_{m,n}(j,t)\simeq\mathsfit{S}^{\mathrm{LL}}_{m,n}(x,t)$$

(22)

for $j=\lfloor xt\rfloor$, $\lfloor\cdot\rfloor$ 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 [38, 31, 39].

When turning to the Toda lattice, the conservation laws are now labeled by $n=0,1,...$ and thus $\mathsfit{A},\mathsfit{B},\mathsfit{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. (2) is still valid.

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

$$a_{j}=\mathrm{e}^{-r_{j}/2}.$$

(23)

Then the equations of motion become

$$\frac{d}{dt}a_{j}=\tfrac{1}{2}a_{j}(p_{j}-p_{j+1}),\quad\frac{d}{dt}p_{j}=a_{j-1}^{2}-a_{j}^{2}.$$

(24)

The Lax matrix, $L$, is defined by

$$L_{j,j}=p_{j},\qquad L_{j,j+1}=L_{j+1,j}=a_{j},$$

(25)

$j\in\mathbb{Z}$, and $L_{i,j}=0$ otherwise. Clearly $L=L^{\mathrm{T}}$. The conserved fields are labelled by nonnegative integers and have densities given by

$$Q_{0}(j)=r_{j},\qquad Q_{n}(j)=(L^{n})_{j,j}\,,$$

(26)

with $n\geq 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

$$J_{0}(j)=-p_{j},\qquad J_{n}(j)=(L^{n}L^{\scriptscriptstyle\downarrow})_{j,j},\quad n=1,2,...\,,$$

(27)

where $L^{\scriptscriptstyle\downarrow}$ is the lower triangular part of $L$. Then under the Toda dynamics

$$\frac{d}{dt}Q_{n}(j,t)+J_{n}(j+1,t)-J_{n}(j,t)=0,$$

(28)

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 (5), (15). However, for $n=2$ one obtains $(L^{2})_{0,0}=p_{0}^{2}+a_{-1}^{2}+a_{0}^{2}$ and $(L^{2}L^{\scriptscriptstyle\downarrow})_{0,0}=a_{-1}^{2}(p_{-1}+p_{0})$, which differs from (5), (15) 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 (26) and (27), 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 (6) and the current-field correlator as in (16). In particularly, $B=B^{\mathrm{T}}$. Of course, $C,B$ are now matrices in the Hilbert space of sequences indexed by $\mathbb{N}_{0}$, i.e. the space $\ell_{2}(\mathbb{N}_{0})$. To distinguish $3\times 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

$$S_{m,n}(j,t)=\langle Q_{m}(j,t)Q_{n}(0,0)\rangle_{P,\beta}^{\mathrm{c}}.$$

(29)

and we plan to construct its Landau-Litshitz approximation. In essence this amounts to an analysis of

$$\big{(}\mathrm{e}^{At}C\big{)}_{m,n},\qquad A=BC^{-1}.$$

(30)

While we are mainly interested in the physical fields corresponding to the indices $m,n=0,1,2$, for the operator in (30) an understanding of the infinite dimensional matrices is required.

Starting from the basics, the free energy of the Toda lattice is given by

$$F_{\mathrm{eq}}(P,\beta)=\log\sqrt{\beta/2\pi}+P\log\beta-\log\Gamma(P).$$

(31)

In particular, the average stretch, $\nu$, is determined through

$$\nu(P,\beta)=\partial_{P}F_{\mathrm{eq}}(P,\beta)=\langle Q_{0}(0)\rangle_{P,\beta}=\log\beta-\psi(P),$$

(32)

with $\psi$ the digamma function. Expectations of higher order fields can be written as moments of a probability measure denoted by $\nu\rho_{\mathsf{p}}$,

$$\kappa_{n}=\langle Q_{n}(0)\rangle_{P,\beta}=\int_{\mathbb{R}}\mathrm{d}w\nu\rho_{\mathsf{p}}(w)w^{n},$$

(33)

$n\geq 1$. $\rho_{\mathsf{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

$$Tf(w)=2\int_{\mathbb{R}}\mathrm{d}w^{\prime}\log|w-w^{\prime}|f(w^{\prime}),$$

(34)

$w\in\mathbb{R}$, considered as an operator on $L^{2}(\mathbb{R},\mathrm{d}w)$ and define the number density

$$\rho_{\mathsf{n}}(w)=\mathrm{e}^{-\varepsilon(w)},$$

(35)

with quasi-energies $\varepsilon$. The quasi-energies satisfy the TBA equation

$$\varepsilon(w)=\tfrac{1}{2}\beta w^{2}-\mu-(T\mathrm{e}^{-\varepsilon})(w),$$

(36)

where the chemical potential $\mu$ has to be adjusted such that

$$\int_{\mathbb{R}}\mathrm{d}w\rho_{\mathsf{n}}(w)=P.$$

(37)

Thereby the number density depends on the parameters $P$ and $\beta$.

The TBA equation is closely connected to the $\beta$-ensemble of random matrix theory. We rewrite (36) as

$$-\log\rho_{\mathsf{n}}(w)=\tfrac{1}{2}\alpha w^{2}-\mu-\alpha P(T\rho_{\mathsf{n}})(w).$$

(38)

As $\alpha\to\infty$, 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\sqrt{P},2\sqrt{P}]$. The Lax DOS is the $P$-derivative of $\rho_{\mathsf{n}}$, which diverges as $(w\pm 2\sqrt{P})^{-1/2}$ at the two borders. As $\alpha$ 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 [35, 1, 12]. Denoting the solution of (36) for $\beta=1$ and the constraint (37) by $\rho_{\mathsf{n}}^{*}$ one has

$$\rho_{\mathsf{n}}^{*}(w)=\frac{\mathrm{e}^{-w^{2}/2}}{\sqrt{2\pi}|\hat{f}_{P}(w)|^{2}},\quad\hat{f}_{P}(w)=\int_{0}^{\infty}\mathrm{d}tf_{P}(t)\mathrm{e}^{\mathrm{i}wt},\quad f_{P}(t)=\sqrt{2}\pi^{-1}\Gamma(P)^{-1/2}t^{P-1}\mathrm{e}^{-\frac{1}{2}t^{2}}.$$

(39)

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, 30, 29, 27].

Next we introduce the dressing transformation of some function $f$ by

$$f^{\mathrm{dr}}=\big{(}1-T\rho_{\mathsf{n}}\big{)}^{-1}f$$

(40)

with $\rho_{\mathsf{n}}$ regarded as a multiplication operator. Then number and particle density are related as

$$\rho_{\mathsf{n}}(w)=\frac{\rho_{\mathsf{p}}(w)}{1+T\rho_{\mathsf{p}}(w)}$$

(41)

with inverse

$$\rho_{\mathsf{p}}=(1-\rho_{\mathsf{n}}T)^{-1}\rho_{\mathsf{n}}=\rho_{\mathsf{n}}\varsigma_{0}^{\mathrm{dr}},$$

(42)

using the convention $\varsigma_{n}(w)=w^{n}$.

For the average currents similar identities are available. The central novel quantity is the effective velocity

$$v^{\mathrm{eff}}=\frac{\varsigma_{1}^{\mathrm{dr}}}{\varsigma_{0}^{\mathrm{dr}}},$$

(43)

see [3, 6, 41, 45]. Then

$$\langle J_{0}(0)\rangle_{P,\beta}=-\kappa_{1},$$

(44)

and, for $n\geq 1$,

$$\langle J_{n}(0)\rangle_{P,\beta}=\int_{\mathbb{R}}\mathrm{d}w\rho_{\mathsf{p}}(w)(v^{\mathrm{eff}}(w)-\kappa_{1})w^{n}.$$

(45)

In thermal equilibrium we have $\kappa_{1}=0$.

Since in the following there will be many integrals over $\mathbb{R}$, let us first introduce the abbreviation

$$\langle f\rangle=\int_{\mathbb{R}}\mathrm{d}wf(w).$$

(46)

With this notation the $C$ matrix turns out to be of the form

	$\displaystyle C_{0,0}=\nu^{3}\langle\rho_{\mathsf{p}}\varsigma_{0}^{\mathrm{dr}}\varsigma_{0}^{\mathrm{dr}}\rangle,$
	$\displaystyle C_{0,n}=C_{n,0}=-\nu^{2}\langle\rho_{\mathsf{p}}\varsigma_{0}^{\mathrm{dr}}(\varsigma_{n}-\kappa_{n}\varsigma_{0})^{\mathrm{dr}}\rangle,$
	$\displaystyle C_{m,n}=\nu\langle\rho_{\mathsf{p}}(\varsigma_{m}-\kappa_{m}\varsigma_{0})^{\mathrm{dr}}(\varsigma_{n}-\kappa_{n}\varsigma_{0})^{\mathrm{dr}}\rangle,$		(47)

$m,n\geq 1$. Note that the matrix $C$ has the block structure

$\displaystyle C=\begin{pmatrix}C_{0,0}&C_{0,n}\\ C_{m,0}&C_{m,n}\end{pmatrix},$

(48)

in the sense that $C_{m,n}$ for $m,n\geq 1$ follows a simple pattern. This structure will reappear for $B$ and $\mathrm{e}^{At}C$.

The field-current correlator $B$ can be computed in a similar fashion with the result

	$\displaystyle B_{0,0}=\nu^{2}\langle\rho_{\mathsf{p}}(v^{\mathrm{eff}}-\kappa_{1})\varsigma_{0}^{\mathrm{dr}}\varsigma_{0}^{\mathrm{dr}}\rangle,$
	$\displaystyle B_{0,n}=B_{n,0}=-\nu\langle\rho_{\mathsf{p}}(v^{\mathrm{eff}}-\kappa_{1})\varsigma_{0}^{\mathrm{dr}}(\varsigma_{n}-\kappa_{n}\varsigma_{0})^{\mathrm{dr}}\rangle,$
	$\displaystyle B_{m,n}=\langle\rho_{\mathsf{p}}(v^{\mathrm{eff}}-\kappa_{1})(\varsigma_{m}-\kappa_{m}\varsigma_{0})^{\mathrm{dr}}(\varsigma_{n}-\kappa_{n}\varsigma_{0})^{\mathrm{dr}}\rangle.$		(49)

As in (2), we want to determine the propagator of the Landau-Lifshitz theory, denoted by $S^{\mathrm{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

$$\nu^{-1}(v^{\mathrm{eff}}(w)-\kappa_{1})$$

(50)

in $L^{2}(\mathbb{R},\mathrm{d}w)$. Here $v^{\mathrm{eff}}$ is the effective velocity defined in (43). Using the block convention as in (48), the spacetime correlator in the Landau-Lifshitz approximation is given by

	$\displaystyle S^{\mathrm{LL}}(x,t)=\int_{\mathbb{R}}\mathrm{d}w\delta\big{(}x-t\nu^{-1}(v^{\mathrm{eff}}(w)-\kappa_{1})\big{)}\nu\rho_{\mathsf{p}}(w)$		(51)
	$\displaystyle\times\begin{pmatrix}\nu^{2}\varsigma_{0}^{\mathrm{dr}}(w)^{2}&\nu\varsigma_{0}^{\mathrm{dr}}(w)(\varsigma_{n}-\kappa_{n}\varsigma_{0})^{\mathrm{dr}}(w)\\[4.30554pt] \nu\varsigma_{0}^{\mathrm{dr}}(w)(\varsigma_{m}-\kappa_{m}\varsigma_{0})^{\mathrm{dr}}(w)&(\varsigma_{m}-\kappa_{m}\varsigma_{0})^{\mathrm{dr}}(w)(\varsigma_{n}-\kappa_{n}\varsigma_{0})^{\mathrm{dr}}(w)\end{pmatrix}.$

Note that $S^{\mathrm{LL}}(x,0)=\delta(x)C$. As a property of the Euler equations, the expression (51) possesses exact ballistic scaling,

$$S^{\mathrm{LL}}_{m,n}(x,t)=\frac{1}{t}S^{\mathrm{LL}}_{m,n}(x/t,1).$$

(52)

The correlator $S_{m,n}(j,t)$ is computed in our MD simulations which will then be compared with $S^{\mathrm{LL}}_{m,n}(x,t)$.

For a molecular dynamics simulation one has to first specify a finite ring $[1,\dots,N]$ with suitable boundary conditions. For the dynamics of positions $q_{j}$ and momenta $p_{j}$ one imposes

$$q_{N+1}=q_{1}+\nu N.$$

(53)

The parameter $\nu$ 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},\dots,r_{N},p_{N}\}$ are i.i.d. random variables with a single site distribution as specified in (4). Then the deterministic time evolution is governed by (24) with boundary conditions

$$r_{0}=r_{N},\qquad p_{N+1}=p_{1}.$$

(54)

In fact, the boundary condition in (53) amounts to the micro-canonical constraint

$$\sum_{j=1}^{N}r_{j}=\nu N.$$

(55)

If one sets $\nu=\langle Q_{0}(0)\rangle_{P,\beta}$, 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 $\Delta j/\Delta t=2$ for the entire range of $(P,\beta,m,n)$ displayed in these figures. Taking into account that the initial correlations are proportional to $\delta_{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.

As can be seen from Table 1, we picked the intermediate pressure such that $\nu\simeq 0$. In the particle picture $\nu=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 $\sqrt{N}$. Thus the free volume is of order $1/\sqrt{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. $\nu=0$ poses no particular difficulties for MD simulations. In GHD the factor $1/\nu$ appears in the expression for $v^{\mathrm{eff}}$, see Eq. (51). This makes the numerical scheme slow and only values close to $\nu=0$ are accessible. However the correlator $S$ changes smoothly through $\nu=0$. GHD also covers this seemingly singular point.

Simultaneously A. Kundu [21] posted a somewhat puzzling note. He considers the parameter values $\beta=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 $\delta$-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

$$V_{\mathrm{kd}}(x)=\frac{g}{\gamma}\mathrm{e}^{-\gamma x}$$

(59)

with arbitrary $\gamma,g>0$. To distinguish their parameters from ours, the variables in [22] are here denoted by $\bar{t},\bar{r},\bar{P},\bar{\beta}$. $\bar{P}$ is the physical pressure and, comparing the Gibbs weights, one obtains the relations

$$\beta=\frac{g}{\gamma}\bar{\beta},\qquad P=\frac{1}{\gamma}\bar{P}\bar{\beta}.$$

(60)

From the equations of motions one deduces

$$\bar{t}=\frac{1}{\sqrt{\gamma g}}t,\quad r(t)=\gamma\bar{r}(\bar{t}),\quad p(t)=\frac{g}{\gamma}\bar{p}(\bar{t}).$$

(61)

Thus, translating to our units, the MD simulations reported in [22] are (i) $P=0.01$, $\beta=0.01$, $N=1024$, $t=400$, (ii) $P=1$, $\beta=1$, $N=1024$, $t=200,300$, and (iii) $P=400$, $\beta=400$, $N=256$, $t=80$. In fact, in all three cases the time scales are identical, $t=\bar{t}$. Since GHD was not available yet, no comparison could have been attempted.

Case (i) is a very dilute chain. In this limit $\nu\rho_{\mathsf{p}}$ is a unit Gaussian. The dressed functions become polynomials as $\varsigma_{0}^{\mathrm{dr}}(w)=a_{0}$, $\varsigma_{1}^{\mathrm{dr}}(w)=a_{1}w$, and $\varsigma_{2}^{\mathrm{dr}}(w)=a_{2}w^{2}+a_{3}$ with coefficients $a_{0},...,a_{3}$ depending on $(P,\beta)$. 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 $\beta=1,P=0.40$ and $\beta=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 $\beta\to\infty$ will freeze any motion. But the simultaneous limit $\beta\to\infty$ with $P=\bar{P}\beta$ at fixed physical pressure $\bar{P}$ is meaningful, at least statistically. In this limit $\nu>0$ always. Also the density of states converges to the arcsine distribution,

$$\lim_{\beta\to\infty}\nu\rho_{\mathsf{p}}(w)=\frac{1}{\pi\sqrt{4\bar{P}-w^{2}}},\quad|w|\leq 2\sqrt{\bar{P}}.$$

(62)

To understand the dynamical behavior, the effective potential is expanded as

$$\mathrm{e}^{-r}+\bar{P}r\simeq\tfrac{1}{2}\bar{P}(r-r_{0})^{2}+c_{0}$$

(63)

at its minimum $r_{0}$. Since $\beta$ is large, the initial fluctuations are of order $1/\sqrt{\beta}$. Therefore the dynamics can be approximated by a harmonic chain with $\omega^{2}=\bar{P}$. The equilibrium time correlations of the harmonic chain have intricate oscillatory behavior [14], which in the large $\beta$ 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.

This material is based upon work supported by the National Science Foundation under Grant No. 1440140, while five of the authors were in residence at the Mathematical Sciences Research Institute in Berkeley, California, during the fall semester of 2021.

The authors would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme “Dispersive hydrodynamics: mathematics, simulation and experiments, with applications in nonlinear waves” where some work on this paper was undertaken. This work was supported by EPSRC grant no EP/R014604/1. TG acknowledges the support of the European Union’s H2020 Marie Skłodowska–Curie grant No. 778010 IPaDEGAN, of INdAM/GNFM and of the research project Mathematical Methods in NonLinear Physics (MMNLP), Gruppo 4-Fisica Teorica of INFN. GM is financed by the KAM grant number 2018.0344. KTRM was supported by a Visiting Wolfson research fellowship from the Royal Society.