Generalized Gibbs Ensembles of the classical Toda chain

TL;DR

This paper derives the generalized free energy of the classical Toda chain using matrix models, linking it to log-gas systems, and discusses the evolution of local density of states via generalized hydrodynamics.

Contribution

It introduces a novel approach to compute the free energy of the Toda chain through the Dumitriu-Edelman matrix model, connecting integrable systems with log-gas models.

Findings

Derived the generalized free energy of the Toda chain.

Mapped the Toda chain to a one-dimensional log-gas system.

Identified the local density of states as an object evolving under generalized hydrodynamics.

Abstract

The Toda chain is the prime example of a classical integrable system with strictly local conservation laws. Relying on the Dumitriu-Edelman matrix model, we obtain the generalized free energy of the Toda chain and thereby establish a mapping to the one-dimensional log-gas with an interaction strength of order 1/N. The (deterministic) local density of states of the Lax matrix is identified as the object, which should evolve according to generalized hydrodynamics. In the current version missing factors of 2 are corrected.

Equations

$$H=\sum_{j=1}^{N}\big{(}\tfrac{1}{2}p_{j}^{2}+\mathrm{e}^{-r_{j}}\big{)},\quad r_{j}=q_{j+1}-q_{j},$$

$$\frac{d}{dt}{q}_j=p_j,\frac{d}{dt}{p}_j={\mathrm{e}}^{-r_{j-1}}-{\mathrm{e}}^{-r_j},$$

$$a_j={\mathrm{e}}^{-r_j/2},b_j=p_j.$$

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

$$(B_N{)}_{j,j}=0,(B_N{)}_{j,j+1}=-(B_N{)}_{j+1,j}=-a_j,$$

$$\frac{d}{dt}{L}_N=[\frac{1}{2}{B}_N,L_N].$$

$$Q^{[n],N}=\mathrm{tr}\big{[}(L_{N})^{n}\big{]}=\sum_{j=1}^{N}((L_{N})^{n})_{j,j}.$$

$$L_{j,j}=b_j,L_{j,j+1}=L_{j+1,j}=a_j,$$

$$Q_j^{[n]}=(L^n{)}_{j,j}$$

$$Q_{j}^{[1]}=p_{j},\quad\tfrac{1}{2}Q_{j}^{[2]}=\tfrac{1}{2}\big{(}p_{j}^{2}+\mathrm{e}^{-r_{j-1}}+\mathrm{e}^{-r_{j}}\big{)}.$$

$$\frac{d}{dt}{Q}_j^{[n]}=J_j^{[n]}-J_{j+1}^{[n]}.$$

$$\frac{d}{dt}L=[\frac{1}{2}B,L],$$

$$\frac{d}{dt}{L}^n=[\frac{1}{2}B,L^n]$$

$$\frac{d}{dt}{Q}_j^{[n]}=(B{L}^n-L^{n}B{)}_{j,j}=-a_j(L^n{)}_{j+1,j}+a_{j-1}(L^n{)}_{j,j-1}.$$

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

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

$$Z_{N}^{-1}\exp\Big{[}\sum_{n=0}^{\kappa}\mu_{n}Q^{[n]}\Big{]},$$

$$V(x)=-\sum _{n=0}^{\kappa }{\mu }_n{x}^n.$$

$$Z_{N}^{-1}\exp\big{(}-\mathrm{tr}[V(L_{N})]\big{)}.$$

$$Q^{[\mathrm{s}],N}=\sum _{j=1}^N{r}_j$$

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

$${\mu }_N^{(V,P)}({\mathrm{d}}^{N}r{\mathrm{d}}^{N}p)=(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.$$

$$\mathrm{tr}[V(L_N)]\ge c_{0}N+c_2\mathrm{tr}[(L_N{)}^2],$$

$$⟨r_0{⟩}_{\tilde{P}}=\log \beta -\frac{{\Gamma }^{\prime }(\beta \tilde{P})}{\Gamma (\beta \tilde{P})}.$$

$$Z_{\mathrm{toda},N}=Z_{L,N}(Z_{L,N}{)}^{-1}{\int }_{{\mathbb{R}}^{2N}}\prod _{j=1}^N\mathrm{d}{r}_j\mathrm{d}{p}_j{\mathrm{e}}^{-P{r}_j}{\mathrm{e}}^{-\frac{1}{2}\mathrm{tr}[(L_N{)}^2]}{\mathrm{e}}^{-\mathrm{tr}[\tilde{V}(L_N)]}.$$

$$Z_{L,N}=\int_{\mathbb{R}^{2N}}\prod_{j=1}^{N}\mathrm{d}r_{j}\mathrm{d}p_{j}\mathrm{e}^{-Pr_{j}}\mathrm{e}^{-\frac{1}{2}\mathrm{tr}[(L_{N})^{2}]}=\big{(}\sqrt{2\pi}\Gamma(P)\big{)}^{N}.$$

$$L_{1,1}={\xi }_{\mathrm{G}},L_{1,2}=\frac{1}{\sqrt{2}}{\chi }_{2P}.$$

$$F_{\mathrm{toda}}=-\underset{N\to \infty }{\lim }\frac{1}{N}\log Z_{\mathrm{toda},N}=-\log (\sqrt{2\pi }\Gamma (P))+F_{\mathrm{L}},$$

$$F_{\mathrm{L}}=-\lim_{N\to\infty}\tfrac{1}{N}\log\mathbb{E}_{L_{N}}\big{(}\mathrm{e}^{-\mathrm{tr}[\tilde{V}(L_{N})]}\big{)}.$$

$$F_{\mathrm{toda}}((1-\eta )V_0+\eta V_1,(1-\eta )P_0+\eta P_1)\le (1-\eta )F_{\mathrm{toda}}(V_0,P_0)+\eta F_{\mathrm{toda}}(V_1,P_1).$$

Full text

Generalized Gibbs Ensembles of the Classical

Toda Chain

Herbert Spohn

Zentrum Mathematik and Physik Department, TUM,

Boltzmannstraße 3, 85747 Garching, Germany, spohn@tum.de

Version 3 of this preprint is published online by Journal of Statistical Physics. Unfortunately, on the right hand side of the two equations in (2.27) a factor $2$ is missing. To properly correct, factors of $2$ dis- and reappear in various formulas throughout Sections 2 and 3. In this version the text, including references, remains untouched, only the factors of 2 are properly taken care of.

Abstract. The Toda chain is the prime example of a classical integrable system with strictly local conservation laws. Relying on the Dumitriu-Edelman matrix model, we obtain the generalized free energy of the Toda chain and thereby establish a mapping to the one-dimensional log-gas with an interaction strength of order $1/N$. The (deterministic) local density of states of the Lax matrix is identified as the object, which should evolve according to generalized hydrodynamics.

*It is my great pleasure to dedicate this article to Joel Lebowitz as teacher, as guide to yet unexplored scientific territories, and with gratitude for a lasting friendship. *

25.11.2019

1 Introduction

Fluid hydrodynamics is based on the notion of local equilibrium. If in the initial state the thermodynamic parameters are slowly varying in space, then they will vary also slowly in time. As a consequence the local equilibrium parameters satisfy the Euler equations which is a closed set of hyperbolic conservation laws. The structure of the Euler equations is universal, while the specific fluid properties enter through a model-dependent pressure. The central cornerstone of hydrodynamic theory is the validity of local microscopic conservation laws. The underlying microscopic dynamics could be either classical, or quantum, or stochastic. If the Navier-Stokes corrections are included, then the time range accessible to hydrodynamics encompasses even the approach to thermal equilibrium.

Over the recent years the hydrodynamic approach has been most successfully extended to integrable quantum chains. The main novel feature is an infinite number of conservation laws. Thus thermal equilibrium must be replaced by a generalized Gibbs ensemble (GGE), which depends on an infinite set of chemical potentials. However, to even write down the macroscopic conservation laws one has to compute the GGE average of the conserved fields and their currents. For the field averages a suitable extension of the Bethe ansatz usually suffices. The real difficulty consists then in finding out about the average currents, which for interacting systems cannot be written as derivatives of the GGE free energy. The generic structure of the currents was first derived in [1, 2], leading to a coupled set of equations called generalized hydrodynamics (GHD). Fully worked out examples have been accomplished, as the quantum sinh-Gordon model [1], the XXZ spin chain [2], the spin $\tfrac{1}{2}$ Hubbard model [3], and the Lieb-Liniger model [1] with a recent experimental confirmation [4].

With this background, one might wonder about classical integrable systems. The still instructive blue-print is a fluid of hard rods in one dimension [5, 6], for which domain wall initial states have been studied [7] and even the Navier-Stokes correction has been proved [8]. In the context of the Korteweg-de-Vries equation, Zakharov [9] proposed a hydrodynamic type theory. We refer to [10] for more details. For the relativistic sinh-Gordon model, also an integrable classical field theory, GHD is derived and studied in [11]. Besides their intrinsic interest, classical models are more accessible to molecular dynamics simulations. A system size of $10^{4}$ particles is standard, while quantum DMRG simulations would not go much beyond $10^{2}$ spins. Accurate simulations are already available for the Toda lattice [12], the Faddeev-Takhtajan spin chain [13], and a discrete integrable version of the sinh-Gordon model [11]. With classical models one has a better chance to test the range of validity of GHD.

In this note I will consider the classical Toda chain. My main novel observation is to connect the Dumitriu-Edelman matrix model with the Toda Lax matrix in thermal equilibrium. The generalized Gibbs free energy of the Toda chain turns out to be related to the $\beta$-ensembles of random matrix theory in the mean-field regime, for which the $\beta$-parameter scales as $1/N$. Thereby an exact variational formula is obtained for the density of states of the Lax matrix, when its matrix elements are distributed according to some GGE.

2 Toda lattice, local conservation laws, and GGE

Following the notation of M. Toda [14], Section 3, the hamiltonian of the Toda chain on sites $j=1,...,N$ is given by

[TABLE]

where periodic boundary conditions, $q_{N+j}=q_{j}+\ell$, are assumed. The positional increments, $r_{j}$, are the stretches. They tell us about the elongation/compression of the chain and can be of either sign. The equations of motion read

[TABLE]

and are viewed as a discrete nonlinear wave equation, just as the discrete nonlinear Schrödinger equation. The underlying lattice of labels is $\mathbb{Z}$, resp. some bounded interval of $\mathbb{Z}$. Alternatively, one can interpret (2.2) as labelled point particles moving on the real line, $\mathbb{R}$, see [15] for more details. We define

[TABLE]

Then the finite $N$ Lax matrix [16] 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. There is also the antisymmetric $B_{N}$ matrix, $B_{N}=-(B_{N})^{\mathrm{T}}$, with

[TABLE]

$(B_{N})_{1,N}=-(B_{N})_{N,1}=a_{N}$ because of periodic boundary conditions, and $(B_{N})_{i,j}=0$ otherwise. Then, under the Toda dynamics,

[TABLE]

Since $B_{N}$ is antisymmetric, the eigenvalues of $L_{N}$ are time-independent and the locally conserved fields $Q^{[n],N}$, $n=1,...,N$, of the Toda lattice can be written as

[TABLE]

Switching to the infinite lattice. As a generic result, the dynamics of the chain is well defined for initial conditions $\{q_{j},p_{j},j\in\mathbb{Z}\}$, which increase at infinity slower than exponential, see [17] for more precise statements. Let us introduce the two-sided infinite Lax matrix, $L$, by

[TABLE]

and $L_{i,j}=0$ otherwise, $i,j\in\mathbb{Z}$. From (2.7) we infer that the conserved fields have a density given by

[TABLE]

at site $j$. Note that $Q_{j}^{[n]}$ depends at most on the variables in $[j-n,...,j+n]$. Explicit expressions for $n=1,...,5$ are listed in [18]. In particular, for momentum and energy,

[TABLE]

Since there is an explicit expression for $Q_{j}^{[n]}$, one can determine the currents defined through

[TABLE]

Using that

[TABLE]

with $B$ the two-sided infinite version of $B_{N}$, one obtains

[TABLE]

and

[TABLE]

Since $L$ is symmetric, the second term is the left shift of the first term and hence one local version of the current is

[TABLE]

with $L^{\mathrm{off}}$ denoting the off-diagonal part of $L$. The total current at finite volume $N$ can be written as

[TABLE]

We turn to the generalized Gibbs ensembles (GGE) of the finite $N$ Toda chain, our main goal being their infinite volume free energy. Since the Toda lattice has many conserved fields, the canonical prescription is the introduce a chemical potential, $\mu_{n}$, for each field, i.e. to consider

[TABLE]

where $\kappa>0$ is assumed to be even and $\mu_{\kappa}<0$. If the conventional picture of thermalization for isolated systems still applies to integrable models, one would expect that in the long time limit the statistics of local observables would be provided by (2.17). Of course, the chemical potentials have to be chosen in such a way that the expectation values of the conserved fields in the initial state are matched. It will be convenient to have the short hand

[TABLE]

Then the GGE can be written as

[TABLE]

Physically our expression for the GGE misses the fact that the compression of the Toda lattice will also change on the hydrodynamic scale. In case of a generic anharmonic chain, an external pressure is enforced by, for example, fixing the left end and pushing/pulling the right end. From the equivalence of ensembles, see the discussion in [19, 20], it follows that the exponential of (2.17) has to be complemented by the stretch field

[TABLE]

with dual parameter $P$, which controls the pressure in the chain. Only for $P>0$ one has a well-defined GGE. For $P=0$, even worse for $P<0$, the chain would simply fall apart. From the equations of motion (2.2), obviously $r_{j}$ is locally conserved with corresponding current $-p_{j}$. Hence we complete the list of conservation laws by

[TABLE]

The finite volume GGE is then

[TABLE]

Since $\mu_{\kappa}<0$, there are constants $c_{0},c_{2}$, $c_{2}>0$, such that $V(x)\geq c_{0}+c_{2}x^{2}$ and hence

[TABLE]

which ensures a finite partition function.

Thermal equilibrium corresponds to the choice $V(x)=\tfrac{1}{2}\beta x^{2}$, $\beta>0$ the inverse temperature. [Below, $\beta$ will also be used as parameter of the log-gas. Given the context, no ambiguity should arise.] The physical pressure, $\tilde{P}$, is equated with the thermally averaged force between a nearest neighbor pair of particles, which then yields $P=\beta\tilde{P}$. For convenience $P$ is still called pressure. In thermal equilibrium the average stretch is computed as

[TABLE]

Hence there is a critical $\tilde{P}_{\mathrm{c}}$ such that $\langle r_{0}\rangle_{\tilde{P}_{\mathrm{c}}}=0$. For $\tilde{P}<\tilde{P}_{\mathrm{c}}$ the positions are predominantly ordered as $q_{j}<q_{j+1}$ and $\langle r_{0}\rangle_{\tilde{P}}\to\infty$ as $\tilde{P}\to 0$. On the other hand, for $\tilde{P}>\tilde{P}_{\mathrm{c}}$ the positions are predominantly reversely ordered and $\langle r_{0}\rangle_{\tilde{P}}\simeq-\log\tilde{P}$.

One might wonder, whether one could soften the somewhat stringent assumption (2.18). In fact, for later considerations only a confining potential $V$ with some smoothness will be required. Not to interrupt the flow of the argument we collect some technical remarks in the Appendix.

We split as $V(x)=\tfrac{1}{2}x^{2}+\tilde{V}(x)$ to obtain

[TABLE]

By an explicit computation

[TABLE]

Thus in (2.25) we average $\exp(-\mathrm{tr}[\tilde{V}(L_{N})])$ over a normalized reference measure, its expectation being denoted by $\mathbb{E}_{L_{N}}(\cdot)$. Under this measure $L_{N}$ becomes a symmetric random Jacobi matrix. Let $\xi_{\mathrm{G}}$ be a unit Gaussian random variable, probability density function (2\pi)^{-\frac{1}{2}}\exp\big{(}-x^{2}/2\big{)}, and let $\chi_{k}$ a chi-distributed random variable with parameter $k>0$, probability density function \big{(}2^{(k/2)-1}\Gamma(k/2)\big{)}^{-1}x^{k-1}\exp\big{(}-x^{2}/2\big{)} for $x>0$ and zero otherwise. Then $\{(L_{N})_{j,j},j=1,...,N\}$ and $\{(L_{N})_{j,j+1},j=1,...,N\}$ are two independent families of i.i.d random variables, where in distribution

[TABLE]

The free energy of the Toda chain is now expressed as

[TABLE]

where

[TABLE]

The Toda free energy depends on the potential $V$ and the pressure $P$, and likewise $F_{\mathrm{L}}$. Using the standard Hölder inequality for the partition function, one obtains for $0\leq\eta\leq 1$

[TABLE]

Thus as a function of the intensive parameters $F_{\mathrm{toda}}$ is convex down, while being convex up as function of the extensive fields. In particular, for fixed $V$, $F_{\mathrm{toda}}(P)$ is convex down and $\partial_{P}F_{\mathrm{toda}}(P)=\langle r_{0}\rangle_{P}$ is decreasing, consistent with the explicit case (2.24). To keep the notation slim, the dependence on $V,P$ is not displayed, except where needed from the context.

The expression (2.29) looks like an unsurmountable obstacle. To my own surprise random matrix theory turns out to be the supportive tool.

3 Random matrix model, free energy

Dumitriu and Edelman [21] introduced a Jacobi matrix model which has the same eigenvalue distribution as the $\beta$-ensemble from random matrix theory [22]. Their $N\times N$ matrix model is

[TABLE]

with $(T_{N,\beta})_{i,j}=0$ otherwise, $i,j=1,...,N$. The random variables $\{\xi_{\mathrm{G},j},\chi_{(N-j)\beta},j=1,...,N\}$ are independent. The average with respect to $T_{N,\beta}$ is denoted by $\mathbb{E}_{T_{N,\beta}}(\cdot)$.

Proposition 3.1. [21] Let $\{\lambda_{1},...,\lambda_{N}\}$ be the eigenvalues of $T_{N,\beta}$. Then their joint probability density function is given by the $\beta$-ensemble of random matrix theory, namely

[TABLE]

with partition function

[TABLE]

To establish the connection with the Lax matrix we set $\beta=2P/N$. Then locally close to $\lfloor(1-u)N\rfloor$, $0\leq u\leq 1$, the matrix $T_{N,2P/N}$ will look like $L_{uP}$. Since the parameter of the $\chi$ random variables changes on scale $1/N$, the corresponding free energies merely add up, see Appendix. Thus, setting

[TABLE]

one concludes that

[TABLE]

Relying on Proposition 3.1,

[TABLE]

where the quadratic term in (3.2) recombines with $\tilde{V}$ to arrive at $V$. Thus the large $N$ limit of the properly scaled Dumitriu-Edelman matrix is reduced to the mean-field limit for a conventional statistical mechanics model. Its central object is the empirical density

[TABLE]

As $N\to\infty$ the sequence $\rho_{N}(\mathrm{d}x)$ converges to the almost sure limit, $\rho^{*}(x)\mathrm{d}x$, which is determined by a variational formula. For a detailed and instructive discussion we refer to [23], Section 2.3. Note that in the following discussion $V$ remains fixed.

More precisely one defines the free energy functional

[TABLE]

The first two terms are the obvious continuum limit of the corresponding expressions in (3). The entropy term results from the product reference measure. $\rho^{*}$ is then the minimizer of $\mathcal{F}_{P}^{\mathrm{MF}}$ under the constraint

[TABLE]

and the limiting free energy associated to (3) is given by $\mathcal{F}_{P}^{\mathrm{MF}}(\rho^{*})$. Since $\mathcal{F}_{P}^{\mathrm{MF}}$ is convex, the minimizer is unique. Then

[TABLE]

Hence

[TABLE]

Changing variables and differentiating with respect to $P$ yields

[TABLE]

$P$ appears only in the strength of the quadratic interaction term. Thus it turns out to be more convenient to introduce $\varrho=P\rho$ and the modified free energy functional

[TABLE]

Then $P\mathcal{F}_{P}^{\mathrm{MF}}(P^{-1}\varrho)=\mathcal{F}(\varrho)-P\log P$. The $P$-dependence has now been moved to the constraint as

[TABLE]

We denote the minimizer of $\mathcal{F}$ by $\varrho^{*}(P)$. Then

[TABLE]

To remove the constraint one introduces the Lagrange multiplier $\lambda$ and sets

[TABLE]

A minimizer of $\mathcal{F}_{\lambda}(\varrho)$ is denoted by $\varrho_{\lambda}^{*}$ as determined by the solution of the Euler-Lagrange equation

[TABLE]

$\lambda$ has to be adjusted such that

[TABLE]

To obtain the Toda free energy, we differentiate as

[TABLE]

Thus integrating (3.17) against $\partial_{P}\varrho^{*}(P)$ one obtains

[TABLE]

From (3.15) it then follows that

[TABLE]

As mentioned already, $F_{\mathrm{toda}}(P)$, hence $\lambda(P)$, is convex down and $\lambda^{\prime}(P)=\langle r_{0}\rangle_{P}$ is monotone deceasing. For $P\to 0$ the stretch becomes large and $\lambda^{\prime}(P)\to\infty$. On the other hand for, $P\to\infty$ the average stretch turns negative and $\lambda^{\prime}(P)\to-\infty$. Hence there is a unique pressure, $P_{\mathrm{c}}$, such that $\lambda^{\prime}(P_{\mathrm{c}})=0$. The solution to (3.17) has exactly two branches, labelled by $\lambda,-$ for $0<P<P_{\mathrm{c}}$ and $\lambda,+$ for $P_{\mathrm{c}}<P$. More precisely (3.18) reads

[TABLE]

which defines the two branches $P_{\pm}(\lambda)$ with $\lambda<\lambda(P_{\mathrm{c}})$. In particular

[TABLE]

Our interest are the Toda GGE averages

[TABLE]

where $\langle\cdot\rangle$ stands for the infinite volume GGE average. These expectations can be computed through first order derivatives of the Toda free energy, starting from the expression (3.15). For the average of the $n$-th conserved field one obtains

[TABLE]

Therefore we first introduce the linearization of $\varrho^{*}$ as

[TABLE]

Then

[TABLE]

where we used that $\int\mathrm{d}x\varrho^{*\{n\}}(x)=0$. Integrating the Euler-Lagrange equation (3.17) at $\lambda=\lambda(P)$ against $\varrho^{*\{n\}}$, the terms on the right of (3) vanish and

[TABLE]

For the average stretch one differentiates (3.21) to obtain

[TABLE]

By the same procedure one can work out the susceptibility matrix

[TABLE]

We note that $\varrho^{*\{n\}}$ satisfies the linear equation

[TABLE]

Then

[TABLE]

We return to the Lax matrix $L_{N}$ with eigenvalues $\ell_{1},...,\ell_{N}$ and hence with empirical density

[TABLE]

From (3.28) one deduces that on average

[TABLE]

As to be discussed in Sect. 4, when integrated against some smooth test function, the convergence is actually almost surely. Thus $\partial_{P}\varrho^{*}(x,P)$ is the deterministic density of states (DOS) of the Lax matrix, whose matrix elements are distributed according to the GGE at parameters $V,P$. In particular, $\partial_{P}\varrho^{*}(x,P)\geq 0$, which cannot be deduced directly from its definition. On the hydrodynamic scale, by assumption, the (deterministic) local DOS is initially changing slowly in space and, because of the conservation law, this feature is expected to persist as slow variation in space-time. However $\partial_{P}\varrho^{*}(x,P)$ does not carry any information about the average stretch $\lambda^{\prime}(P)$. Another option, comprehending all data, would be the two branches of $\partial_{\lambda}\varrho_{\lambda,\pm}^{*}$. The corresponding pressure is defined through (3.22). Since

[TABLE]

“$-$” referring to $0<P<P_{\mathrm{c}}$ and “$+$” referring to $P_{\mathrm{c}}<P$, one concludes that

[TABLE]

with $\langle r_{0}\rangle_{P}\partial_{\lambda}\varrho_{\lambda,\pm}^{*}(x)\geq 0$ and normalized to 1.

Equation (3.21), and its consequences (3.28), (3.29), (3.34) are the main result of our contribution. For the Toda chain, they establish exact expressions for the GGE free energy and the GGE average of the conserved fields.

4 Mean-field limit for Dyson’s Brownian motion

While in the context of the Toda lattice only the static mean-field limit is of relevance, the dynamic approach to this limit turns out to be instructive. In particular, additional identities will be generated. We thus consider the stochastic particle system on $\mathbb{R}$ governed by

[TABLE]

with $\{b_{j}(t),j=1,...,N\}$ a collection of independent standard Brownian motions. This is Dyson’s Brownian motion in an external potential $V$. For easier comparison we have chosen the more neutral $\alpha$, instead of $P$. The interaction has strength $1/N$, which corresponds to a standard mean-field limit. We assume an initial measure such that for all smooth test functions, $f$, the initial empirical density has the almost sure limit $\rho_{0}$,

[TABLE]

Then also for $t>0$

[TABLE]

almost surely, where $\rho(x,t)$ is the solution to the mean-field equation

[TABLE]

with initial condition $\rho(x,0)=\rho_{0}(x)$.

The unique stationary probability measure of (4.1) is given by $\mu_{\alpha,N}^{\mathrm{MF}}$ with density

[TABLE]

The corresponding empirical density converges to the stationary solution of (4.4). More precisely, under $\mu_{\alpha,N}^{\mathrm{MF}}$, for all smooth test functions $f$,

[TABLE]

almost surely, where $\rho_{\mathrm{s}}$ satisfies

[TABLE]

Point-wise this equation reads

[TABLE]

Since Dyson’s Brownian motion is time-reversible, stationarity is equally defined by the large round bracket in (4.8) to vanish. The respective free energy (large deviation) functional is (3.8), setting $P=\alpha$. Its minimizer $\rho^{*}$ satisfies a corresponding Euler-Lagrange equation. Differentiating this equation with respect to $x$ and then multiplying by $\rho^{*}$ yields

[TABLE]

Thus we conclude $\rho_{\mathrm{s}}=\rho^{*}$.

Under general conditions on $V$, the law of large numbers (4.3), (4.4) has been proved by Cépa and Lépingle [24]. Much studied is also the strong coupling regime [25], for which $2\alpha/N$ in (4.1) is substituted by $\beta>0$. The two regimes have a very different $N$-dependence. In mean-field the free energy is of order $N$, the empirical density is strictly positive with size of order 1, and the fluctuations are of order $\sqrt{N}$. On the other side, for fixed $\beta$ there is a strong repulsion between particles. Thus the free energy is of order $N^{2}$ and the empirical density is extended over a region of size $\sqrt{N}$, while the shape function is generically of compact support. Hence the diffusion term is of order $1/N$ and in the variational formula the entropy term may be dropped. Particles can hardly move any more and fluctuations are of order 1. However the large $\alpha$ limit of $\rho_{\mathrm{s}}$ seems to match with the strong coupling steady state in the limit of small $\beta$.

5 Toda chain in thermal equilibrium

Thermal equilibrium corresponds to the special choice $V(x)=\tfrac{1}{2}x^{2}$, in which case the stationary density $\rho_{\mathrm{s}}$ is known explicitly and given by

[TABLE]

where

[TABLE]

Physically, one should introduce the proper inverse temperature and pressure. This will be done at the end of this section.

Opper [26] derived (3.17) and (3.18) for $V(x)=\tfrac{1}{2}x^{2}$ through the semi-classical limit of the quantum Toda chain and obtained the exact solution (5.1), (5.2), see also the earlier work [27] and the discussion in [14], Section 6.6. Of course, the connection to random matrices was not known yet. The more recent investigation [28] starts from a modified Dyson’s Brownian motion, such that in a particular limit one arrives at the mean-field interaction. As a common feature of these studies the Stieltjes transform of (4.8) is shown to satisfy a local quadratic equation. It can be linearized to arrive at a linear second order differential equation, the Weber equation, which then leads to the explicit solution (5.1), (5.2). For small $\alpha$, $\rho_{\mathrm{s}}$ is close to the unit Gaussian, while $\rho_{\mathrm{s}}$ converges to the Wigner semi-circle law in the limit $\alpha\to\infty$. A related study has been carried out for $\beta$-Wishart ensembles [29].

Note that \frac{d}{d\alpha}\big{(}\alpha\rho_{\mathrm{s}}\big{)} is the DOS of the two-sided Lax matrix, $L$, defined through (2.27) in the two-sided limit $N\to\infty$. $L$ can be viewed as a random Jacobi matrix and from this perspective Duy and Shirai [30] and Duy [31] independently employ the connection to the Dumitriu and Edelman theorem to establish the exact DOS.

Choosing $f$ in (4.7) to be a monomial one obtains a recursion relation for the moments of $\rho_{\mathrm{s}}$. The odd moments vanish and for the even moments, $c_{n}=\int_{\mathbb{R}}\mathrm{d}x\rho_{\mathrm{s}}(x)x^{2n}$,

[TABLE]

This recursion relation is stated already in [30].

The large $\alpha$ behavior can read off from (2.7) by expanding \mathrm{tr}\big{[}(L_{N})^{n}\big{]}. Each term corresponds to a random walk, step size $0,\pm 1$, starting and ending at site $j$, and weight resulting from the thermal average. For large $\alpha$ to leading order only paths with no step of size [math] count and hence

[TABLE]

On the right side one observes the Catalan numbers. Hence $\rho_{\mathrm{s}}$ is the normalized Wigner semi-circle probability distribution function

[TABLE]

and the Lax matrix DOS equals \big{(}\pi\sqrt{4\alpha-x^{2}}\big{)}^{-1} with $x^{2}\leq 4\alpha$.

As to be explained in Section 6 below, also the Gaussian fluctuations are of interest in connection with hydrodynamics linearized around thermal equilibrium. The CLT for the Lax matrix is studied in [32] with a somewhat different perspective. Our goal is to have reasonably explicit expressions for the covariance matrix.

We start from stationary Dyson’s Brownian motion (4.1) with $V(x)=\tfrac{1}{2}x^{2}$ and define the fluctuation field

[TABLE]

for some smooth test function $f$. Then, in distribution,

[TABLE]

with $\phi(f,t)$ an infinite-dimensional stationary Ornstein-Uhlenbeck process governed by

[TABLE]

where $\xi$ denotes standard space-time white noise. The stationary covariance $C$, with kernel $C(x,y)$, is then defined through the equation

[TABLE]

where

[TABLE]

and $\langle\cdot;\cdot\rangle$ denotes the usual $L^{2}$ inner product. A corresponding result for strong coupling regime has been proved in [25]. The martingale identities used there seem to extend to our case.

For $\alpha=0$ one confirms that the solution to (5.10) is given by

[TABLE]

with $\rho_{\mathrm{G}}$ the unit Gaussian. For $\alpha>0$ the most instructive choice for the basis of functions seems to be simply the powers $x^{n},n\geq 0$. We set

[TABLE]

and note the boundary condition $C_{0,n}=0=C_{n,0}$ to be used below. By reflection symmetry of (3), $C_{m,n}=0$ for odd $m+n$. Denoting by $\langle\cdot\rangle_{\rho_{\mathrm{s}}}$ the average with respect to $\rho_{\mathrm{s}}$ and by $\rho_{\mathrm{s}}$ being even, for $n>1$ the action of $D$ reads

[TABLE]

which defines the matrix $(\Lambda^{\mathrm{T}})_{m,n}$. In its explicit form

[TABLE]

while all other matrix elements vanish, $m,n\geq 1$. Then the stationary condition (5.9) translates into

[TABLE]

One computes $C_{1,1}=\langle x^{0}\rangle_{\rho_{\mathrm{s}}}=1$. Assume that $C_{m,n}$ is known for $1\leq m,n\leq\kappa$. Then $C_{m,\kappa+1}$, $m=1,...,\kappa+1$, can be obtained iteratively through (5.15) upon starting from the boundary. The remaining matrix elements from the set $1\leq m,n\leq\kappa+1$ follow by symmetry. Hence, not only the mean, but also the covariance matrix allows for a recursive computational scheme.

Finally, we extend our results for the Toda lattice at thermal equilibrium with inverse temperature $\beta$ and physical pressure $P$. For this purpose we have to adjust our notation to a more explicit form. We set $\alpha=\beta P$ and add the $\alpha$- and $\beta,P$-dependence in the obvious way as $\rho_{\mathrm{s},\alpha}$, $L(P,\beta)$. Let us recall the random Lax matrix

[TABLE]

$(L_{\alpha})_{i,j}=0$ otherwise, $i,j\in\mathbb{Z}$, and all random matrix elements being independent. Now, with $\tilde{H}=H+P\sum_{j}r_{j}$,

[TABLE]

which defines the thermal Lax matrix $L(P,\beta)$. In other words,

[TABLE]

in distribution. Therefore the DOS scales as

[TABLE]

which for the moments implies

[TABLE]

The stretch average is treated separately by setting

[TABLE]

and

[TABLE]

6 Linearized hydrodynamics

As discussed in [33], for writing down the linearized Euler equation two matrices are required. In principle the scheme works for any GGE, but in our context we linearize around thermal equilibrium, hence the thermal average $\langle\cdot\rangle_{P,\beta}$ appears with $\beta$ the inverse temperature and $P$ the physical pressure, while $\beta P$ is the Lagrange multiplyer for $r_{j}$. The first matrix is the static susceptibility defined by

[TABLE]

where $\langle\cdot;\cdot\rangle_{P,\beta}$ stands for the second cumulant. Note that by strict locality the sum consists of finitely many terms only. The second matrix is the cross correlation between fields and currents,

[TABLE]

In the limit of zero lattice spacing, the lattice of labels, the linearized Euler equations read

[TABLE]

because of the transformation from intensive to extensive variables. In passing we recall that the matrix $\mathsf{B}\mathsf{C}^{-1}\mathsf{B}$ is the Drude weight of the Toda chain in thermal equilibrium, in the sense that \big{(}\mathsf{B}\mathsf{C}^{-1}\mathsf{B}\big{)}_{m,n} equals the long time asymptotics of the correlation between the $m$-th and $n$-th total current [34, 33, 35].

Let us deviate for a moment to explain why, despite its asymmetric looking definition, one finds that $\mathsf{B}$ is actually symmetric, see also [20]. For this purpose first note that for an initial GGE the process $Q^{[n]}_{j}(t)$ is space-time stationary. Hence

[TABLE]

which upon summing over $j$ and differentiating at $t=0$ turns into

[TABLE]

Thus $\mathsf{B}=\mathsf{B}^{\mathrm{T}}$ by partial summation.

In connection with the molecular dynamics [12], of particular interest is the field-field time correlation

[TABLE]

For example, up to factors of $2$, the matrix elements $m,n=\mathrm{s},1,2$ would be stretch, momentum, and energy correlations, including their cross-correlations. In leading order one would expect that such correlations are equal to the solution to (6.3) with random initial conditions, whose covariance equals

[TABLE]

The solution is given by

[TABLE]

Since $\mathsf{A}\mathsf{C}=\mathsf{B}$, the more symmetric form reads

[TABLE]

which admits the spectral resolution $\rho_{m,n}(\mathrm{d}\lambda)$ such that

[TABLE]

In particular $\mathsf{S}_{m,n}(x,t)$ scales exactly as $t^{-1}f(x/t)$. On the other hand the true microscopic correlations attain their scaling behavior only for large $j,t$ both of the same order and

[TABLE]

with $x$ the continuum version of $j$. Molecular dynamics simulations of the lowest order equilibrium time correlations are available for parameters $N=1024$, $\beta=1$, $P=1$ at times 200, 300. The average is over $10^{6}-10^{7}$ initial thermal samples and ballistic scaling is well satisfied [12].

Some general properties of $\mathsf{C},\mathsf{B}$ can be obtained from time-reversal $p_{j}\leadsto-p_{j}$ for all $j$, denoted by $\mathcal{R}$. The equilibrium measure is even under $\mathcal{R}$. By inspection $\mathcal{R}Q^{[n]}_{0}=(-1)^{n}Q^{[n]}_{0}$ and similarly for the currents $\mathcal{R}J^{[n]}_{0}=(-1)^{n+1}J^{[n]}_{0}$. Hence, for all $m,n=\mathrm{s},1,2,...\,$, $\mathrm{s}$ considered as even,

[TABLE]

and

[TABLE]

We call matrices with the left property in (6.13) of type $\mathcal{A}_{0}$ and with the right property of type $\mathcal{A}_{1}$. One easily verifies that if $C_{1},C_{2}\in\mathcal{A}_{0}$, then $C_{1}C_{2}\in\mathcal{A}_{0}$. Also if $C_{1}\in\mathcal{A}_{0}$ and $B_{1}\in\mathcal{A}_{1}$, then $A_{1}B_{1},B_{1}A_{1}\in\mathcal{A}_{0}$. Hence $\mathsf{C}^{-\frac{1}{2}}\mathsf{B}\mathsf{C}^{-\frac{1}{2}}\in\mathcal{A}_{0}$. This property ensures that $\rho_{m,n}(\mathrm{d}\lambda)=\rho_{m,n}(-\mathrm{d}\lambda)$ and thus $\mathsf{S}(x,t)=\mathsf{S}(-x,t)$.

7 Conclusions

Beyond the Toda chain there are several classical integrable models which admit a Lax pair. Thus it would be of interest whether similar ideas as in this note could be used to construct GGEs. For one-dimensional continuum models, as the sinh-Gordon model and the nonlinear Schrödinger equation, because of ultraviolet divergencies, the GGE tends to be a singular object. But rather commonly there are lattice versions which might be more approachable. Possibly there is some generic structure of which the Toda chain is just one example.

If in the Euler-Lagrange equation (3.17) one substitutes $\varrho_{\lambda}^{*}(x)=\exp[-\varepsilon(x)]$, then

[TABLE]

which can be viewed as the semi-classical limit of the thermodynamic Bethe ansatz (TBA) for the quantum Toda chain. In fact, for $V(x)=\tfrac{1}{2}x^{2}$, this is how (7.1) has been derived first [27, 26], see also [36].

While writing up my results, Benjamin Doyon accomplished a derivation of the Toda GHD [15]. Doyon starts from the fluid picture of Toda particles moving on the real line, obtains (7.1) using inverse scattering theory, and thereby confirms that $\partial_{\lambda}\varrho_{\lambda}^{*}$ can be interpreted as density of quasi-particles. GHD is written for both, lattice and fluid, coordinate frames. Independently, GHD is derived through a semi-classical limit starting from the quantum Toda chain, also including a discussion of quantum corrections [37]. Both contributions state the same definite prediction for the DOS linked to the currents. Comparing with (2.16), it differs from the Lax matrix DOS only by the extra factor $L_{N}^{\mathrm{off}}$, which means that now some information on the eigenvectors of the Lax matrix is required. In the more recent study [38], the DOS of the Lax matrix is sampled by Monte-Carlo techniques for quadratic and double well potentials. One finds very good agreement with numerically iterating (4.4) and letting the solution relax to the steady state. Also, the prediction for the GGE average currents is tested. While an analytical argument is still missing, the simulations are strongly favoring the conjecture.

Acknowledgements. I am grateful to Duy Trinh and Simone Warzel for most helpful comments at an early stage of this project. I benefitted from most stimulating exchanges with Benjamin Doyon and thank Christian Mendl for a careful reading of a first draft. When completing my work, Xiangyu Cao and Vir Bulchandani informed me about their study of the Toda lattice with domain wall initial conditions. In particular, they pointed out to me the two-valuedness of $P(\lambda)$.

Appendix A Appendix: infinite volume limit, linear ramp

We use Flaschka variables and restate the issue, omitting irrelevant constants. The variables are $\{a_{j}\in\mathbb{R}_{+},b_{j}\in\mathbb{R},j=1,...,N\}$. Then the Toda GGE reads

[TABLE]

and the linearly ramped Toda measure is given by

[TABLE]

Up to inverting the ramp, normalization, and a boundary term, the ramped Toda measure agrees with (3.1), (3).

We will establish three properties.

(i) The following limit exists and, up to constants, agrees with $F_{\mathrm{toda}}(P)$,

[TABLE]

(ii) The infinite volume limit of $\mu_{\mathrm{to},N}$, shifted by $-\tfrac{1}{2}N$, exists. The limit measure, $\mu_{\mathrm{to}}$, defined on $(\mathbb{R}_{+}\times\mathbb{R})^{\mathbb{Z}}$ is stationary under spatial shifts and exponentially mixing.

(iii) For the linearly ramped free energy one obtains

[TABLE]

Our proof relies on transfer matrix techniques. We set $N=\kappa M$, $M$ integer, and denote $\mathsf{a}_{m}=(a_{\kappa(m-1)+1},...,a_{\kappa(m-1)+\kappa})$, $\mathsf{b}_{m}=(b_{\kappa(m-1)+1},...,b_{\kappa(m-1)+\kappa})$, $m=1,...,M$. We break the system into blocks of length $\kappa$. Then the interaction is only between neighboring blocks and, denoting by $\mathsf{a},\mathsf{b}$ the variables in block $m$ and by $\mathsf{a}^{\prime},\mathsf{b}^{\prime}$ the variables in block $m+1$, one can choose a transfer matrix $\mathcal{T}_{P}(\mathsf{a},\mathsf{b}|\mathsf{a}^{\prime},\mathsf{b}^{\prime})$ such that

[TABLE]

where $\mathrm{tr}_{\mathcal{H}}$ refers to the trace in the Hilbert space \mathcal{H}=L^{2}\big{(}(\mathbb{R}_{+}\times\mathbb{R})^{\tau},\mathrm{d}^{\tau}\mathsf{a}\mathrm{d}^{\tau}\mathsf{b}\big{)}. $\mathcal{T}_{P}>0$ pointwise, $\mathcal{T}_{P}$ can be choosen to be symmetric, and $\mathcal{T}_{P}$ is trace class. By the Perron-Frobenius theorem, $\mathcal{T}_{P}$ has thus a non-degenerate, strictly positive eigenvalue, $\lambda_{\mathrm{max}}(P)$, which is separated by a gap from the rest of the spectrum, and

[TABLE]

By an extension of the argument, expectation values of local observables have an infinite volume limit. Exponential mixing is a consequence of the spectral gap.

For the ramped partition function the transfer matrix becomes now $m$-dependent. By perturbation theory, for the transfer matrix from block $m$ to block $m+1$ the maximal eigenvalue equals $\lambda_{\mathrm{max}}(Pm/M)$ up to an error of order $1/N$. Also the spectral gap is at most shifted by order $1/N$. Hence, for large $N$, up to an error of order $\mathrm{e}^{-\gamma N}$ with the constant $\gamma>0$ related to the spectral gap, one arrives at

[TABLE]

which in the limit $N\to\infty$ converges to the right hand side of (A.4).

Our proof relies heavily on the assumption (2.18) assuring a finite power series. Once having moved to the mean-field limit, any confining potential with some minimal smoothness would do.