The GGE averaged currents of the classical Toda chain

TL;DR

This paper investigates the averaged currents in the classical Toda chain with random initial data, providing numerical evidence supporting the collision-rate assumption and suggesting a universal structure for integrable systems.

Contribution

The study offers the first numerical validation of the collision-rate assumption for the Toda chain's averaged currents, reinforcing the generalized Euler equations framework.

Findings

Numerical evidence supports the collision-rate assumption in the Toda chain.

Averaged currents align with predictions from generalized Euler-type equations.

Supports the universality of integrable systems' hydrodynamic descriptions.

Abstract

The Toda chain with random initial data is studied. Of particular interest are generalized Gibbs ensembles, their averaged conserved fields, and the averages of the corresponding currents. While averaged fields are well-understood, the description of averaged currents has hitherto relied on the collision-rate assumption. For the Toda chain, the rate assumption can be investigated numerically. Here, we provide convincing evidence for the validity of the rate assumption. This lends further support to the idea that generalized Euler-type equations have a structure common to all integrable extensive systems.

Equations

$$\displaystyle v^{\mathrm{eff}}(v)\Delta t=v\Delta t+\Big{(}\int_{v}^{\infty}dw\Phi(v,w)\rho_{\mathrm{p}}(w)(v^{\mathrm{eff}}(w)-v^{\mathrm{eff}}(v))$$

$$\displaystyle\hskip 80.0pt+\int_{-\infty}^{v}dw\Phi(v,w)\rho_{\mathrm{p}}(w)(v^{\mathrm{eff}}(v)-v^{\mathrm{eff}}(w)\Big{)}\Delta t.$$

$$H=\sum _{j=1}^N\frac{1}{2}{p}_j^2+\sum _{j=1}^{N-1}{V}_{\mathrm{ch}}(q_{j+1}-q_j),$$

$${\ddot{q}}_j(t)=V_{\mathrm{ch}}^{\prime }(q_{j+1}(t)-q_j(t))-V_{\mathrm{ch}}^{\prime }(q_j(t)-q_{j-1}(t)).$$

$$r_j=q_{j+1}-q_j,$$

$${\dot{r}}_j=p_{j+1}-p_j,{\dot{p}}_j=V_{\mathrm{ch}}^{\prime }(r_j)-V_{\mathrm{ch}}^{\prime }(r_{j-1}),$$

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

$$\big{(}-p_{j},-V_{\mathrm{ch}}^{\prime}(r_{j-1}),-p_{j}V_{\mathrm{ch}}^{\prime}(r_{j-1})\big{)}.$$

$$Z^{-1}\exp\big{[}-\beta\big{(}\tfrac{1}{2}(p_{0}-u)^{2}+V_{\mathrm{ch}}(q_{0})+\tilde{P}q_{0})\big{)}\big{]}\,.$$

$$\tilde{P}=-⟨V_{\mathrm{ch}}^{\prime }{⟩}_{\tilde{P},\beta }$$

$$⟨(r_j,p_j,e_j){⟩}_{\tilde{P},u,\beta }=(\ell ,u,\mathfrak{e}),\mathfrak{e}=\frac{1}{2}{\beta }^{-1}+⟨V_{\mathrm{ch}}{⟩}_{\tilde{P},\beta },$$

$$⟨(-p_j,-V_{\mathrm{ch}}^{\prime }(r_{j-1}),-p_j{V}_{\mathrm{ch}}^{\prime }(r_{j-1})){⟩}_{\tilde{P},u,\beta }=(-u,\tilde{P},u\tilde{P}).$$

$${\partial }_t\ell +{\partial }_x{\mathsf{j}}_{\ell }=0,{\partial }_t\mathsf{u}+{\partial }_x{\mathsf{j}}_{\mathsf{u}}=0,{\partial }_t\mathfrak{e}+{\partial }_x{\mathsf{j}}_{\mathfrak{e}}=0,$$

$$(\mathsf{j}_{\ell},\mathsf{j}_{\mathsf{u}},\mathsf{j}_{\mathfrak{e}})=\big{(}-\mathsf{u},\tilde{P}(\ell,\mathfrak{e}-\tfrac{1}{2}\mathsf{u}^{2}),\mathsf{u}\tilde{P}(\ell,\mathfrak{e}-\tfrac{1}{2}\mathsf{u}^{2})\big{)}.$$

$$V_{\mathrm{ch}}(r)={\mathrm{e}}^{-r},$$

$$a_j=\frac{1}{2}{\mathrm{e}}^{-r_j/2},b_j=\frac{1}{2}{p}_j.$$

$$(L_N{)}_{j,j}=b_j,(L_N{)}_{j,j+1}=(L_N{)}_{j+1,j}=a_j,$$

$$Q^{[n],N}=\mathrm{tr}\big{[}(L_{N})^{n}\big{]}$$

$$Q_j^{[n]}=L_{j,j}.$$

$$J^{[n],N}=\mathrm{tr}\big{[}(L_{N})^{n}L_{N}^{\mathrm{off}}\big{]}$$

$$J_j^{[n]}=(L^n{L}_{\mathrm{off}}{)}_{j,j}.$$

$$Q_j^{[\mathrm{s}]}=r_j,J_j^{[\mathrm{s}]}=-p_j.$$

$$(Z_{\mathrm{toda},N}{)}^{-1}{\mathrm{e}}^{-\mathrm{tr}[V(L_N)]}\prod _{j=1}^N{\mathrm{e}}^{-P{r}_j}\mathrm{d}{r}_j\mathrm{d}{p}_j,P>0.$$

$$F_{\mathrm{toda}}(P,V)=-\underset{N\to \infty }{\lim }\frac{1}{N}\log Z_{\mathrm{toda},N}.$$

$$L_N{\psi }_{\ell }(j)={\lambda }_{\ell }{\psi }_{\ell }(j),\ell =1,...,N.$$

$$\frac{1}{N}{Q}^{[n],N}={\int }_{\mathbb{R}}dx{x}^n\frac{1}{N}\sum _{\ell =1}^N\delta (x-{\lambda }_{\ell })x^n,$$

$$\frac{1}{N}{J}^{[n],N}={\int }_{\mathbb{R}}dx{x}^n\frac{1}{N}\sum _{\ell =1}^N\sum _{j=1}^{N}2a_j{\psi }_{\ell }(j){\psi }_{\ell }(j+1)\delta (x-{\lambda }_{\ell }).$$

$$\frac{1}{N}\sum _{\ell =1}^N\delta (x-{\lambda }_{\ell })\to {\rho }_Q(x),\frac{1}{N}\sum _{\ell =1}^N\sum _{j=1}^{N}2a_j{\psi }_{\ell }(j){\psi }_{\ell }(j+1)\delta (x-{\lambda }_{\ell })\to {\rho }_J(x).$$

$$\underset{N\to \infty }{\lim }\frac{1}{N}⟨Q^{[n],N}{⟩}_{N,P,V}=⟨(L^n{)}_{0,0}{⟩}_{P,V}={\int }_{\mathbb{R}}dx{\rho }_Q(x)x^n$$

$$\underset{N\to \infty }{\lim }\frac{1}{N}⟨J^{[n],N}{⟩}_{N,P,V}=⟨(L^n{L}^{\mathrm{off}}{)}_{0,0}{⟩}_{P,V}={\int }_{\mathbb{R}}dx{\rho }_J(x)x^n.$$

$${\mathcal{F}}_P^{\mathrm{MF}}(\rho )={\int }_{\mathbb{R}}\mathrm{d}x\rho (x)V(x)-P{\int }_{\mathbb{R}}\mathrm{d}x{\int }_{\mathbb{R}}\mathrm{d}y\log |x-y|\rho (x)\rho (y)+{\int }_{\mathbb{R}}\mathrm{d}x\rho (x)\log \rho (x),$$

Full text

The GGE averaged currents of the classical Toda chain

Xiangyu Cao

Department of Phyics, University of California, Berkeley, 94720 California, United States of America.

Vir B. Bulchandani

Department of Phyics, University of California, Berkeley, 94720 California, United States of America.

Herbert Spohn

Zentrum Mathematik and Physik Department, TUM, Boltzmannstr. 3, 85747 Garching, Germany.

Abstract. The Toda chain with random initial data is studied. Of particular interest are generalized Gibbs ensembles, their averaged conserved fields, and the averages of the corresponding currents. While averaged fields are well-understood, the description of averaged currents has hitherto relied on the collision-rate assumption. For the Toda chain, the rate assumption can be investigated numerically. Here, we provide convincing evidence for the validity of the rate assumption. This lends further support to the idea that generalized Euler-type equations have a structure common to all integrable extensive systems.

1 Introduction, collision rate assumption

Euler equations are based on microscopic conservation laws, with possible additions when symmetries are broken. One might wonder whether such a structure remains valid for one-dimensional integrable many-body systems, for which the number of conservation laws is extensive. There has been a lot of interest in this topic, mostly driven by the investigation of quantum integrable systems. We refer to a few fundamental contributions [1, 2, 3, 4, 5, 6, 7, 8, 9], in which the reader can find a more comprehensive list.

The basic ingredients are easily recalled. For a given model one has to determine all local (perhaps quasi-local) conservation laws. Then one writes down the generalized Gibbs ensemble (GGE). In many cases its free energy can be obtained from a Bethe ansatz. The GGE averaged conserved fields are then first derivatives of this free energy. Determining GGE averaged currents is a more complicated story. In most models it is already difficult to write down an expression for the microscopic current and computing the GGE average current can be accomplished only in very exceptional cases. One thus has to a rely on an educated guess that we call the collision rate assumption, which is first explained in general terms.

Consider a uniform, one-dimensional fluid, consisting of quasi-particles with bare velocity $v\in\mathbb{R}$. Quasi-particles of velocity $v$ have a uniform spatial density $\rho_{\mathrm{p}}(v)$, that does not change under the dynamics. Two colliding quasi-particles retain their velocity but undergo a spatial shift. Quantum mechanically this would correspond to the two-particle phase shift. When viewed over longer time intervals, the intrinsic velocity $v$ of a quasiparticle is modified to an effective velocity $v^{\mathrm{eff}}(v)$ through collisions. To compute this velocity, we add a tracer quasi-particle with velocity $v$ at the origin. Upon colliding with a fluid quasi-particle of velocity $w$, the position of the tracer is shifted by $\Phi(v,w)$. Let us now consider a time span $\Delta t$, which is short on the macro-scale but long on the micro-scale. Then $v^{\mathrm{eff}}(v)\Delta t$ is the true displacement of the tracer particle up to time $\Delta t$. $v\Delta t$ is its bare displacement to which one has to add the contribution from collisions. The collision rate is $\rho_{\mathrm{p}}(w)|v^{\mathrm{eff}}(w)-v^{\mathrm{eff}}(v)|$, which leads to

[TABLE]

As a consequence the macroscopic current, used in a hydrodynamic description, is given by $v^{\mathrm{eff}}(v)\rho_{\mathrm{p}}(v)$.

According to our current understanding the relation (1) holds for all one-dimensional integrable many-particle systems, classical and quantum. There also seems to be no small parameter, except that the conserved fields are assumed to vary slowly on the microscopic scale. Of course, the labelling of the quasi-particles could be more complicated, and the same is true for the bare velocity. But $\Phi$ is always deduced from the microscopic two-particle phase shift. So what is the evidence for the collision rate assumption (1)?

Much work has been done for quantum systems, such as the integrable spin-$1/2$ XXZ chain [3, 4, 5, 7, 8]. One typically starts from a natural initial state, such as a domain wall or a localized cloud in infinite space, and solves numerically the corresponding kinetic equation, which is based on the collision rate assumption. Some specific averages of interest are then compared with DMRG simulations of the microscopic dynamics. While the agreement is often excellent, the test as such involves many intermediate steps, although recent analytical work [8] has demonstrated that the collision rate assumption for averaged spin and energy currents in the XXZ model follows directly from thermodynamic Bethe ansatz (TBA). On the classical side, an early work discusses the KdV equation [10]. One prepares initially a low density of solitons, adds a tracer soliton, and verifies quite accurately the properly adapted Eq. (1) [11]. The regime of high soliton density is addressed in [12], where the collision rate assumption is verified by an analytical argument involving the thermodynamic limit of the Whitham equations (see also  [13]). For a fluid of hard rods with diameter $a$, Eq. (1) can be easily verified, since the tracer particle collisions turn out to be statistically independent and, by definition of the model, the shift $\Phi=\pm a$ is independent of $v,w$ [14]. The interpretation of the quantum rate assumption in terms of classical flea-gas dynamics was developed in [15, 16], while the connection to the Whitham modulation theory of classical soliton lattices was discussed in [5, 17].

There is one further classical system for which (1) has been established recently [18]. It is a cellular automaton known as the particle-box system [19]. One considers the one-dimensional lattice $\mathbb{Z}$ (= boxes) and each site is either occupied by a particle (= ball) or empty. Thus a configuration is an infinite string of $0,1$’s. To define a single time step we assume a finite number of particles and introduce a carrier which starts to the left with load 0 and moves step-by-step to the right. At the first encounter with a particle, say at site $j$, the site $j$ switches to empty and the carrier load increases to $1$. If now site $j+1$ happens to be empty, then the carrier drops its particle, i.e. the occupation at $j+1$ switches to $1$, while the carrier load decreases to [math]. Otherwise, if site $j+1$ happens to be occupied, then the carrier picks up this particle, i.e. the occupation at site $j+1$ switches to [math] and the carrier load increases to $2$, etc., until the whole sweep is completed resulting in the time $1$ configuration of particles and carrier load 0. A $k$-soliton is a configuration of the form $111...000...$ each block of length $k$, $k=1,2,...$ . The bare velocity of the soliton is simply $k$. If a $k$-soliton bumps into a slower $m$ soliton, $m<k$, then the $k$-soliton is shifted by $2m$ and the $m$-soliton by $-2m$. The obvious adaptation of (1) is proved for a large class of space-time stationary probability measures on configurations. The proof is difficult and uses a particular representation of the dynamics, which is hidden behind the bare definition the box-ball system.

In our paper we investigate the collision rate assumption for the classical Toda lattice [20, 21, 22]. This model has several simplifying features. Firstly, it is classical, so that the GGE can be sampled by Monte-Carlo methods. Secondly, there is a simple and, in a sense, very explicit expression for the microscopic currents. Furthermore, the GGE free energy comes from a variational principle, whose minimizer can be obtained from Dyson’s Brownian motion, which simplifies the task of finding solutions to the TBA equation. Our goal is to check the collision rate assumption directly, by computing the GGE average of the microscopic currents.

2 Hydrodynamics for non-integrable chains

We briefly recall the hydrodynamics of a non-integrable classical chain [23]. Considered are $N$ particles in one dimension with positions, $q_{j}$, momenta $p_{j}$, $j=1,\ldots,N$. Their mass is set equal to one. Nearest neighbours in index space are coupled through a chain potential $V_{\mathrm{ch}}$. The Hamiltonian of the chain then reads

[TABLE]

resulting in the equations of motion

[TABLE]

We regard (2.2) as a discrete nonlinear wave equation. One could also think of particles moving on the real line. Their positions are not ordered and particles have a somewhat peculiar interaction, which is not invariant under relabelling. This second picture is referred to as “fluid”. As discussed at length in [21], one can switch back and forth between the fluid and chain pictures. We adopt here the latter picture, which turns out to be somewhat simpler.

It will be convenient to introduce the positional differences

[TABLE]

more physically referred to as the stretch, which can be negative. $r_{j}$ is the free “volume” between particles $j+1$ and $j$. The equations of motion (2.2) then turn into

[TABLE]

$j=1,\ldots,N$, where periodic boundary conditions have been adopted, $p_{N+1}=p_{1}$, $r_{0}=r_{N}$. Viewed as a classical lattice field theory, the underlying lattice is $[1,\ldots,N]$ with periodic boundary conditions and the field variables are $(r_{j},p_{j})$. Dynamically, $r_{j}$ is coupled to its right and $p_{j}$ to its left neighbour. From the dynamics, one can read off the local conservation laws, namely stretch, momentum, and energy,

[TABLE]

and their currents

[TABLE]

The assumption of being non-integrable can be rephrased as having no further local conservation laws than already listed.

Hydrodynamics is based on the assumption that local equilibrium is approximately maintained over the course of time evolution. While plausible, such a dynamical property is difficult to deduce from the equations of motion. Nevertheless, to write down the correct macroscopic equations one merely has to average the local fields and currents in an equilibrium state. With our choice of variables, the canonical equilibrium states are of product form, with one factor given by

[TABLE]

To yield a well-defined partition function, the chain potential should be bounded from below and have an at least one-sided, say, linearly increasing lower bound. The pressure, $\tilde{P}$, which might have to be restricted to a half-line, is dual to the stretch. $\beta$ is the inverse temperature and $u$ the mean velocity. Note that

[TABLE]

is the average force between neighbouring particles. The required thermal averages are easily accomplished. For the fields, one obtains

[TABLE]

and their currents,

[TABLE]

To arrive at an explicitly closed system, one still has to reexpress the average currents in terms of the conserved fields. In other words, the intensive parameters of the canonical ensemble have to be substituted by the extensive ones, which can always be done due to the convexity of the free energy. If we use $x\in\mathbb{R}$ for the corresponding continuum approximation, the conserved fields are the local stretch $\ell(x)$, the local momentum $\mathsf{u}(x)$, and the local total energy $\mathfrak{e}(x)=\frac{1}{2}\mathsf{u}(x)^{2}+\mathsf{e}(x)$ per particle, with internal energy $\mathsf{e}$. The microscopic conservation laws then turn into the Euler equations of the nonlinear chain, as

[TABLE]

where the hydrodynamic currents are given by

[TABLE]

Here $\tilde{P}(\ell,\mathsf{e})$ results from inverting the map $(\tilde{P},\beta)\mapsto(\ell,\mathsf{e})$.

The Toda lattice is a special anharmonic chain with interaction potential

[TABLE]

which requires the restriction $\tilde{P}>0$. With $N$ degrees of freedom, the Toda chain has $N$ independent local conservation laws. To write down its Euler-type equations, we will still use the blue-print provided by the non-integrable case.

3 Lax matrix and GGE averaged fields

The conserved fields and currents of the Toda chain may be obtained concisely from the Lax matrix $L_{N}$ [24]. We define

[TABLE]

Then the finite $N$ Lax matrix is the tridiagonal real symmetric matrix with

[TABLE]

$(L_{N})_{1,N}=(L_{N})_{N,1}=a_{N}$ because of periodic boundary conditions, and $(L_{N})_{i,j}=0$ otherwise. For infinite volume, the index is extended to $j\in\mathbb{Z}$ and the corresponding tridiagonal matrix is denote by $L$. We also introduce $L^{\mathrm{off}}$ as the off-diagonal part of $L$. The conserved fields of the Toda chain are then given by

[TABLE]

with the obvious density

[TABLE]

From the equations of motion one obtains the corresponding total currents [20],

[TABLE]

with density

[TABLE]

As in the non-integrable case, in addition there is the conserved stretch and its current

[TABLE]

The finite volume generalized Gibbs ensemble (GGE) is given by

[TABLE]

One should think of the potential $V$ as a finite power series with a strictly positive even leading coefficient. Physically this would correspond to a finite number of chemical potentials. But presumably a larger class of $V$ will work as well. Thermal equilibrium corresponds to $V(x)=2\beta(x-\tfrac{1}{2}u)^{2}$, compare eq. (2.7). However, for notational simplicity we set $P=\beta\tilde{P}$. The Toda free energy is defined through the infinite volume limit

[TABLE]

The Lax matrix has the eigenvalues $\lambda_{\ell}$ and eigenvectors $\psi_{\ell}(j)$, both depending on $N$,

[TABLE]

The observables of interest are then

[TABLE]

The first integrand above is the normalized density of states (DOS) of the Lax matrix with entries distributed according to GGE with parameters $P,V$. In the second expression there is an extra weight coming from the eigenvectors. In the limit $N\to\infty$ we expect both densities to converge to a limit,

[TABLE]

In fact, this limit should be self-averaging. By definition, $\rho_{Q}(x)\geq 0$ and $\int dx\rho_{Q}(x)=1$. But the weights in (3.12) have no definite sign and the same is true for $\rho_{J}(x)$. Denoting the finite volume GGE expectation by $\langle\cdot\rangle_{N,P,V}$, and $\langle\cdot\rangle_{P,V}$ its infinite volume limit, one arrives at

[TABLE]

for the conserved fields and correspondingly for the currents,

[TABLE]

(Existence of the thermodynamic limits on the left-hand side follows by extensivity.) As will be explained, variational type formulas for $F_{\mathrm{toda}}(P,V)$ and $\rho_{Q}(x)$ [20] are available. There is also a prediction for $\rho_{J}(x)$ based on the collision rate assumption [21, 22]. So our strategy will be to sample the quantities in (3.13) and to compare them with the theoretical predictions.

4 Free energy, random matrix model

We start from the free energy functional

[TABLE]

which has to be minimized under the constraint

[TABLE]

Since $\mathcal{F}_{P}^{\mathrm{MF}}$ is convex, the minimizer is unique and will be denoted by $\rho^{*}$. The $P$-dependence can be trivially shifted to the constraint by setting $\varrho(x)=P\rho(x)$ and defining

[TABLE]

to be minimized under the constraint

[TABLE]

Clearly, $P\mathcal{F}_{P}^{\mathrm{MF}}(P^{-1}\varrho)=\mathcal{F}(\varrho)-P\log P$ and $\varrho^{*}(x)=P\rho^{*}(x)$.

Introducing the Lagrange multiplier, $\mu$, the Euler-Lagrange equation for $\mathcal{F}$ reads

[TABLE]

Following common practice, one writes $\varrho_{\mu}=\mathrm{e}^{-\varepsilon}$ with pseudo-energies $\varepsilon$. Then (4.5) turns into the central identity

[TABLE]

Originally, in the case $V(x)=x^{2}$, this equation was obtained from the semiclassical limit of the TBA equation for the quantum Toda chain, hence the name classical TBA, although the Bethe ansatz appeared nowhere in our discussion.

For the GGE Toda free energy one finds that

[TABLE]

Equivalently

[TABLE]

The DOS, $\rho_{Q}$, is related to the functional derivative of the Toda free energy w.r.t. $V$. Using (4.3), (4.4) one obtains

[TABLE]

Let us define

[TABLE]

Then $\varrho^{*}(x,P)=\varrho_{\mu(P)}(x)$ and differentiating w.r.t. $P$,

[TABLE]

Finally one can differentiate TBA w.r.t. $\mu$ to obtain

[TABLE]

which replaces the $\mu$-differentiation by an integral equation.

The chemical potential $\mu(P)$ is convex down with maximum $\mu_{\mathrm{max}}$ at $P_{\mathrm{c}}$ and $\mu(P)\to-\infty$ for $P\to 1,\infty$. $P<P_{\mathrm{c}}$ is the low pressure and $P>P_{\mathrm{c}}$ is the high pressure branch. At $P=P_{\mathrm{c}}$ one has $\langle r_{0}\rangle_{P_{\mathrm{c}},V}=0$ and hence $q_{N}-q_{0}=\mathcal{O}(\sqrt{N})$. At low pressure the positions are approximately ordered, while at high pressure they are reverse ordered. Thus for $\mu\neq\mu_{\mathrm{max}}$ the integral equation (4.6) has two solutions [22]. An example is shown in Figure 1. As observed in [22], the numerical solution of (4.6) by iteration with fixed $\mu$ converges systematically to the low pressure branch, which should be considered an artifact of the iteration scheme. Instead, the two solutions can be reliably found by a Fokker-Planck or Dyson Brownian motion approach, which will be outlined below.

For $P\to\infty$, the entropy term in (4.1) can neglected and one arrives at the much studied variational problem for the ground state of the one-dimensional log gas. For $V(x)=x^{2}$ its solution is the Wigner semicircle law. Since the Lax matrix DOS is its derivative, one obtains $\rho_{Q}(x)=\pi^{-1}(1-x^{2})^{-\frac{1}{2}}$. In Fig. 1, the pressure is chosen to be $P=4.1$ and one observes that the DOS already starts to build up a singularity at the two edges $x=\pm 1$.

5 Dyson’s Brownian motion

We differentiate TBA w.r.t. to $x$ and obtain

[TABLE]

which can be viewed as the stationary solution of the nonlinear Fokker-Planck equation

[TABLE]

This suggests that the minimizer can be obtained from simulating many diffusing particles subject to an external potential $V$ and a logarithmic mean field interaction. Their stochastic differential equations read

[TABLE]

with $\{b_{j}(t),j=1,...,N\}$ a collection of independent standard Brownian motions. In the conventional Dyson’s Brownian motion the interaction is strong, i.e. $\mathcal{O}(1)$, and denoted by $\beta$. In this case, the time-stationary distribution is the much studied $\beta$-ensemble of random matrix theory, which refers to the energy dominated regime [25]. The variational principle (4.1) still holds, provided the entropy term is dropped. In our context the interaction is mean field and entropy and energy are both of order $N$. Abstractly, it is ensured that the dynamics (5.3) reaches a unique stationary state and for its empirical density one has the limiting behavior

[TABLE]

Thus, rather than to trying to numerically solve TBA, one can run the dynamics (5.3) until the stationary regime is reached. According to (5.4), sampling the positions of the particles in the steady state provides a good approximation to the solution of TBA.

Alternatively, one may integrate the Fokker-Planck PDE numerically. Starting from an arbitrary initial condition $\rho(x,0)$ such that $\int\rho(x,0)dx=1$ and $\rho(x,0)>0$, the solution is guaranteed to converge. In practice, we choose a standard Gaussian initial condition $\rho(x,0)=e^{-x^{2}/2}/\sqrt{2\pi}$; for all GGE ensembles we considered, convergence (defined in practice by $|\partial_{t}\rho|<10^{-5}$ for all $x$) is reached at $t\lesssim 5$.

6 GGE averaged currents

For the Toda lattice the collision rate assumption becomes

[TABLE]

In the case that

[TABLE]

the current DOS is conjectured to be equal to

[TABLE]

We test this conjecture numerically for a few different GGE states $(P,V)$. For each state, we first integrate the Fokker-Planck equation (5.2) until convergence to find the solution to the TBA equation (with a given pressure). The quantities on the RHS of (6.3) are obtained by numerically inverting the linear integral equations, see (4.12) and also Appendix 7. This provides the analytical prediction for $\rho_{J}(x)$, which is to be compared with direct simulations of the Lax matrix (of size $N=1024$), via Eq. (3.15). When $V(x)=ax^{2}/2$, the matrix elements $a_{j}$ and $b_{j}$ are independent and can be sampled directly, see Eq. (3.8). We also consider cases where $V(x)$ is not quadratic. In that case, we sample the GGE ensemble by Markov chain Monte Carlo. In either case, we measure $\rho_{J}(x)$ using (3.13) by exact diagonalization and averaging over many samples. The results for several GGE ensembles are shown in Figure 1. An excellent agreement is found in all tested cases. In particular, we exhibit the two solutions to the TBA equation with a fixed $\mu$, which describe distinct GGE ensembles, with opposite signs of $\langle r_{0}\rangle_{P,V}$. The conjecture for the current (6.3) appears to be valid for all GGE ensembles.

7 Conclusions

While we have tested only four different parameter values, the agreement between theory and numerical results is very convincing, and lends additional support to previous findings. Beyond the rigorously understood case of hard rods [26], the classical Toda chain appears to be the only integrable many-body system for which the collision rate assumption can be tested directly, for the full current density of states $\rho_{J}(x)$.

The dynamical assumption underlying generalized hydrodynamics is the propagation of local GGE states on physically reasonable space-time scales. Our discussion does not touch on this issue at all. Rather, if the dynamics of the Toda lattice propagates as a local GGE, then the local Lagrange multipliers have to be updated according to the kinetic equation based on (6.1).

Acknowledgments. We thank Benjamin Doyon and Joel E. Moore for helpful comments on the manuscript. We acknowledge support from the ERC synergy Grant UQUAM, the DOE grant DE-SC0019380 (XC) and the DARPA DRINQS program (VBB).

Appendix A Appendix: Some identities

To be self-contained, we list some well-known facts [2, 3]. Let us introduce the integral operator

[TABLE]

Then the TBA equations read

[TABLE]

Setting $\mathrm{e}^{-\varepsilon}=\varrho_{\mu}$, the dressing of a function $\psi$ is defined by

[TABLE]

On the right, $\varrho_{\mu}$ is regarded as a multiplication operator, $(\varrho_{\mu}\psi)(x)=\varrho_{\mu}(x)\psi(x)$. From (4.12) we conclude

[TABLE]

Here $[1]$ stands for the constant function $\psi(x)=1$ and similarly $[x]$ for the linear function $\psi(x)=x$.

The collision rate assumption can be expressed as

[TABLE]

whose solution can be written as

[TABLE]

To verify we start from (A.6) and use (A.4) as

[TABLE]

equivalently

[TABLE]

From (4.12) we note that the ratio equals $1+T\partial_{\mu}\rho_{\mu}$, thus yielding (A.5).