Hydrodynamic Scales of Integrable Many-Particle Systems

Herbert Spohn

Departments of Mathematics and Physics, Technical University Munich,

Boltzmannstr. 3, 85747 Garching, Germany

December 6, 2023

To my grandson Lio Spohn for his continuing support

###### Contents

  1. 1 Overview
  2. 2 Dynamics of the classical Toda lattice
    1. 2.1 Locally conserved fields and their currents
    2. 2.2 Action-angle variables, notions of integrability
      1. 2.2.1 Conventional integrability of the Toda lattice
      2. 2.2.2 Hydrodynamic perspective on integrability
    3. 2.3 Scattering theory
      1. 2.3.1 Hard rod fluid
      2. 2.3.2 Two-particle Toda lattice
      3. 2.3.3 $N$-particle Toda lattice
  3. 3 Static properties
    1. 3.1 Generalized Gibbs ensembles
    2. 3.2 Lax matrix filter and local GGEs
    3. 3.3 Generalized free energy
    4. 3.4 Lax density of states, TBA equation
      1. 3.4.1 Thermal equilibrium
    5. 3.5 Mean-field techniques
  4. 4 Dyson Brownian motion
    1. 4.1 Macroscopic equation, law of large numbers
    2. 4.2 Fluctuation theory
  5. 5 Hydrodynamics for hard rods
    1. 5.1 Hard rod fluid
    2. 5.2 Hard rod lattice
    3. 5.3 TBA, collision rate ansatz
  6. 6 Equations of generalized hydrodynamics
    1. 6.1 Average currents
    2. 6.2 Hydrodynamic equations
  7. 7 Linearized hydrodynamics and GGE dynamical correlations
    1. 7.1 Equilibrium spacetime correlations for nonintegrable chains
    2. 7.2 GGE spacetime correlations for the Toda lattice
  8. 8 Domain wall initial states
  9. 9 Toda fluid
    1. 9.1 Euler equations
    2. 9.2 Generalized free energy - again
  10. 10 Hydrodynamics of soliton gases
    1. 10.1 Soliton gas of the Korteweg-de Vries equation
    2. 10.2 Soliton gas of the Toda lattice
    3. 10.3 Comparing soliton- and particle-based hydrodynamics
  11. 11 Calogero models
    1. 11.1 Hyperbolic Calogero model, charges and currents
    2. 11.2 Scattering coordinates
    3. 11.3 Generalized free energy
    4. 11.4 Hydrodynamic equations
    5. 11.5 Classical Bethe equations
    6. 11.6 Trigonometric Calogero-Moser model
  12. 12 Discretized nonlinear Schrödinger equation
    1. 12.1 Continuum wave equations
    2. 12.2 Ablowitz-Ladik discretization
    3. 12.3 Circular random matrices with pressure ramp
    4. 12.4 Average currents
    5. 12.5 Hydrodynamic equations
    6. 12.6 Modified Korteweg-de Vries equation
  13. 13 Hydrodynamics for the Lieb-Liniger $\delta$-Bose gas
    1. 13.1 Bethe ansatz
    2. 13.2 Bethe root densities, free energy, TBA equations
    3. 13.3 Charge currents, hydrodynamic equations
    4. 13.4 Comparing Toda fluid and $\delta$-Bose gas, generic TBA structure
    5. 13.5 Gaudin matrix
  14. 14 Quantum Toda lattice
    1. 14.1 Integrability, monodromy matrix
    2. 14.2 Spectral properties
    3. 14.3 GGE and hydrodynamics
  15. 15 Beyond the Euler time scale
    1. 15.1 General framework
    2. 15.2 Nonintegrable chains
    3. 15.3 Navier-Stokes equations

Preface

In August 1974 for the first time I attended a summer school, which happened
to be number three in the series on “Fundamental Problems in Statistical
Mechanics”. The school took place at the Agricultural University of
Wageningen. We were approximately 70 participants from 25 countries. The
school lasted nearly three weeks with 11 lecture courses, each four hours
long, and various more specialized seminars. As to be expected, present were
the big shots, as Eddy Cohen who with moderate success tried to slow down the
speaker through posing questions. Besides lectures, the truly exciting part of
the school was to meet fellow youngsters who had similar interests and
struggled more or less with the same difficulties.

At the time, critical phenomena and the just invented RG methods were the
overwhelming topic. Fortunately the Dutch physics community has a long
tradition in Statistical Mechanics and therefore a wide range of topics were
covered. I vividly recall the lectures by Nico van Kampen on “Stochastic
differential equations”. Joe Ford lectured on “The statistical mechanics of
classical analytic dynamics”, the early days of deterministic chaos. In fact,
at the end of his lectures Joe mentioned the Toda lattice describing the just
confirmed integrability. Perhaps I should have had listened with more care.
For me the strongest impact had the lectures presented by Piet Kasteleyn on
“Exactly solvable lattice models”, explaining the fascinating link between
equilibrium statistical mechanics and models from quantum many-body physics.

My second encounter with integrable systems is related to the study of Dyson
Brownian motion, which is an integrable stochastic particle system. At the
time, to me the model was an intriguing example for the hydrodynamics of a
many-particle system with long range forces. The third encounter was triggered
by the KPZ revolution, which brought me in contact with other corners of
integrable systems. Around 2016 I first learned about the activities
investigating the hydrodynamic scales for integrable quantum many-body
systems. I could not resist. Of course, major insights had been accomplished
already. But, apparently, classical integrable many-particle systems were in a
state of dormancy. This is how my enterprise got started.

During the ongoing project, I had many insightful comments and good advice.
Gratefully acknowledged are Mark Adler, Amol Aggarwal, Vir Bulchandani,
Xiangyu Cao, Kedar Damle, Avijit Das, Percy Deift, Jacopo De Nardis, Atharv
Deokule, Abhishek Dhar, Maurizio Fagotti, Pablo Ferrari, Patrik Ferrari,
Chiara Franceschini, Tamara Grava, Alice Guionnet, David Huse, Thomas
Kappeler, Karol Kozlowski, Thomas Kriecherbauer, Manas Kulkarni, Anupam Kundu,
Aritra Kundu, Gaultier Lambert, Joel Lebowitz, Guido Mazzuca, Ken McLaughlin,
Christian Mendl, Pierre van Moerbecke, Joel Moore, Fumihiko Nakano, Neil
O’Connell, Stefano Olla, Lorenzo Piroli, Balázs Pozsgay, Michael Prähofer,
Sylvain Prolhac, Tomaz Prosen, Keiji Saito, Makiko Sasada, Tomohiro Sasamoto,
Naoto Shiraishi, Jörg Teschner, Khan Duy Trinh, Simone Warzel, and Takato
Yoshimura.

Special thanks are due to Benjamin Doyon. Our encounter at Pont-à-Mousson is
well remembered.

When working on the manuscript I had the opportunity to stay for extended time
spans at the Mathematical Research Institute at Berkeley, the Galileo Galilei
Institute at Firenze, the Newton Institute at Cambridge, and the International
Center for Theoretical Sciences at Bengaluru, in chronological order. I am
most thankful for such generous invitations.

##  1 Overview

Hydrodynamics is based on the observation that the motion of a large assembly
of strongly interacting particles is constrained by local conservation laws.
As a result, local equilibrium is established over an initial time span to be
followed by a much longer time window when local equilibrium parameters are
governed by the hydrodynamic evolution equations. The initial time span could
shrink to microscopic times when the system starts out already in local
equilibrium. It is a matter of fact that a vast amount of interesting physics
is covered by the hydrodynamic approach. Historically the best known example
are simple fluids for which hydrodynamics is synonymous with fluid dynamics.

Already for simple fluids the hydrodynamic approach carries the seed for
further extensions, since the equilibrium phase diagram is richly structured.
Most common is the occurrence of a liquid-gas phase transition. This discrete
order parameter has now to be added as a further parameter characterizing
local equilibrium. For example, gas and fluid phase may spatially coexist and
the respective interface is then an additional slow degree of freedom, to be
included in the macroscopic dynamics.

Close to critical points, conventional hydrodynamics has to be augmented by
more refined theories. At lower temperatures, generically a solid phase
stabilizes. Due to slow relaxation of solids, for the dynamics of the solid-
gas interface mostly non-hydrodynamic modelling is used. Bosonic particles at
low temperatures will form a condensate. One then employs a hydrodynamic two-
fluid model, which governs the superfluid interacting with the normal fluid.
Going beyond short-range interactions, magneto-hydrodynamics describes the
motion of fluids made up of charged particles, also including the Maxwell
field as additional dynamical degrees of freedom. Relativistic hydrodynamics
becomes relevant for extreme events such as the formation of superdense
neutron stars, relativistic jets, and Gamma ray bursts. Each topic mentioned
is part of a vast enterprise with ongoing research.

Approximately seven years ago a novel item was added to our list under the
name of generalized hydrodynamics (GHD). Physically perhaps not as far
reaching as other areas mentioned, generalized hydrodynamics relies on an
amazing twist. The novel topic is concerned with integrable many-particle
models for which the number of conserved fields is proportional to system
size, in sharp contrast to the models listed before which have only a few
conserved fields, for example number, momentum, energy, plus broken symmetries
in case of simple fluids. At first glance the mere idea of a hydrodynamic
description of the time evolution of such an integrable system sounds like an
intrinsic contradiction. After all, establishing local equilibrium relies on
chaotic dynamics which is just the opposite of integrability. But the huge
number of degrees of freedom helps. Since integrable many-particle systems
have an extensive number of local conservation laws, local equilibrium must
now be characterized by a correspondingly large number of chemical potentials.
It is this feature which is called “generalized”. In the limit of infinite
system size, the hydrodynamic fields are labelled by a parameter taking
integer values, $n=0,1,\ldots$, or possibly by more complicated labelling
schemes. As a consequence, writing down the coupled set of hyperbolic
conservation laws is already a major obstacle.

The notion of generalized Gibbs ensemble (GGE) was introduced earlier and
studied systematically in a related context, known as quantum quench. But the
issue is generic. One starts from a spatially homogeneous random state and
wants to identify the random state reached after a long time. More physically,
one prepares a homogeneous state of a particular hamiltonian dynamics and then
abruptly changes the dynamics (the quench). If the quench dynamics is not
integrable, generically one expects the system to thermalize with parameters
determined by the conserved fields when averaged over the initial state. But
an integrable system has many conserved fields and the final state will depend
on an extensive set of parameters. Such asymptotic states are called GGE.

For many-particle systems to be integrable requires fine-tuned interactions.
Nevertheless the list is not so short. The first and still much studied model
is the Lieb-Liniger $\delta$-Bose gas from 1963. The Toda lattice was
discovered in 1967 and its integrability being firmly established seven years
later through the construction of a Lax matrix. Further examples are the XXZ
spin chain, the one-dimensional spin-$\tfrac{1}{2}$ Fermi-Hubbard model, the
classical particle models of Calogero, and continuum wave equations as
Korteweg-de Vries, nonlinear Schrödinger, and sinh-Gordon. In fact, the
central goal of our notes is to argue that

On a hydrodynamic scale all integrable many-particle systems are structurally
alike.

Given the diversity of microscopic models such a claim is surprisingly bold.
On the other hand, as to be discussed, the route to tackle the hydrodynamic
scale will depend on the specific model. The precise meaning of our claim will
unfold. But to provide at least a very preliminary glimpse, in all models
under study the two-body scattering shift will be a crucial piece of the
hydrodynamic description.

As familiar from fluids, Euler equations refer to the ballistic scale, which
is characterized by space and time to be of same order of magnitude. Formally,
entropy is locally conserved. Transport properties arise at longer diffusive
time scales and are included in the hydrodynamic equations through the Navier-
Stokes correction. For integrable systems the same distinctions apply, at
least in principle. This topic will be briefly touched upon in the very last
chapter of our notes. Otherwise, hydrodynamics is understood as ballistic
Euler type scaling.

From a broader perspective, in deriving the equations governing the motion on
the hydrodynamic scale one faces several difficulties.
(i) For a given system, the local conservation laws have to be listed in terms
of which GGEs can be constructed. In a somewhat vague sense, this list has to
be complete, since hydrodynamically all fields are expected to be coupled to
each other.
(ii) The structure of the generalized free energy has to be understood,
including its first order derivatives which are linked to the GGE averaged
conserved fields.
(iii) To complete the hydrodynamic equations, one has to know the GGE averaged
currents as a functional of the GGE averaged conserved fields.

My exposition is not a review, even though much of the relevant literature
will be cited. To establish a guiding backbone, the classical Toda lattice is
discussed in considerable detail. Particularly introduced are two distinct
strategies (1) a closed system with a linearly varying pressure and (2) the
canonical transformation to scattering coordinates. The first method will also
be applied to the Ablowitz-Ladik discretized nonlinear Schrödinger equation
and the second one serves well for the Calogero fluid. The key quantum models
accounted for will be the Lieb-Liniger $\delta$-Bose gas and the quantized
Toda lattice, for both models relying on the Bethe ansatz as strategy.

Let me refrain from further comments on the content of my notes and rather
turn to some remarks on the history of the subject. The one-dimensional system
of classical hard rods was studied around 1970. The system is integrable,
since in a collision momenta are merely exchanged. Conserved is any one-
particle sum function depending only on the momenta. Due to the hard core the
hydrodynamic fields are nonlinearly coupled, in sharp contrast to an ideal
gas. Jerry Percus first derived the hydrodynamic equations. In the 1980ies
Roland Dobrushin and collaborators analyzed in much greater detail the time
evolution of hard rods. He well understood the hydrodynamic perspective,
including the issue of Navier-Stokes corrections. But at the time no tools
were available for handling more intricate models. In retrospect, the true
simplification of hard rods is a two-particle scattering shift which is
independent of the incoming quasiparticle velocities.

Another early line of research concerns the Korteweg-de Vries equation which
is an integrable nonlinear wave equation accessible through the inverse
scattering transform. In the mid 1990ies Vladimir E. Zhakarov studied a low
density gas of solitons and derived the respective kinetic equation for the
spacetime dependence of the soliton counting function. The extension to a
dense soliton gas, to say the respective hydrodynamic equations, has been
obtained by Gennady El in 2003.

Generalized hydrodynamics as a systematic research activity relies on a
breakthrough advance in 2016 independently by the two groups: B. Bertini, M.
Collura, J. De Nardis, M. Fagotti and O.A. Castro-Alvaredo, B. Doyon, T.
Yoshimura. They discovered a general scheme of how to write down the average
currents, thereby covering classical field theories and quantum many-body
systems. Only with such an input the equations of generalized hydrodynamics
could be written with confidence, herewith opening the door to applications of
physical interest. Such detailed studies strongly support our claim that on a
hydrodynamic scale all integrable many-particle models look alike.

In condensed matter physics the notion quantum many-body is widely used. Many-
particle system is more natural for models from classical mechanics. Both
notions are employed interchangeably and comprise also classical and quantum
field theories.

How the text is structured. Our material is arranged in 15 Chapters with
sections varying in number. Longer subsections might be separated by boldface
headers. No quotations are provided in the main text. Instead, at the end of
each chapter, one finds ”Notes and references”, which roughly means a
bibliography with extended comments. In addition, there are Inserts separated
from the main text by $\blackdiamond\blackdiamond$ Header. $\cdots$
$\blackdiamond\blackdiamond$. Typically an Insert deals with a closely related
topic, which however can be touched upon only superficially. Also some more
technical derivations have been shifted to Inserts. The idea is that at first
reading an Insert can be skipped, except for conventions on notation.

Notes and references

Overview

The Proceedings of the Wageningen summer school have been edited by E.D.G.
Cohen [58]. My work on Dyson Brownian motion is published in Spohn [265]. The
KPZ revolution is covered in many articles from which only the overviews
Corwin [63], Quastel and Spohn [241], Spohn [268] and Takeuchi [279] are
quoted.

Section 1

A useful account on the developments prior to generalized hydrodynamics can be
found in a special volume on “Quantum Integrability in Out-of-Equilibrium
Systems” edited by Calabrese et al. [43]. Its central theme are quantum
quenches starting from a spatially homogeneous initial state.

The notion “scattering shift” is convenient but less widely used. It refers to
the fact that, when two particles undergo a scattering motion, each trajectory
is asymptotically of the form $v_{j}t+\phi_{j}^{\pm}$, $j=1,2$, as
$t\to\pm\infty$, i.e. free motion linear in time and on top a constant
displacement as first order correction, which is the scattering shift. In
quantum mechanical two-body scattering, the wave function for the relative
motion is asymptotically of the form
$\exp\big{[}\mathrm{i}(-kx+\tfrac{1}{2}k^{2}t+\theta(k))\big{]}$. Then a
narrow wave packet, centered at $k_{0}$ in momentum space, travels in physical
space with velocity $k_{0}$ and is displaced by $\theta^{\prime}(k_{0})$.
$\theta$ is the phase shift, while its derivative, $\theta^{\prime}$, is the
scattering shift.

Percus [230] wrote down what he called a kinetic equation. In his prior work,
in collaboration with Lebowitz and Sykes [190], the exact spacetime two-point
function of hard rods in thermal equilibrium was obtained. Boldrighini,
Dobrushin, and Suhov [25] prove, under fairly general assumptions on the
initial probability measure, the validity of the Euler equation for a system
of hard rods under ballistic scaling. They also established that the
hydrodynamic equation has smooth solutions. Dobrushin [81] is his vision on
hydrodynamic limits. In the context of the Korteweg-de Vries equation,
Zakharov [299] studied a low density gas of solitons with Poisson distributed
centers and statistically independent soliton velocities. He argued for a
Boltzmann type kinetic equation. The hydrodynamic extension to a dense soliton
fluid has been accomplished in El [101], El and Kamchatnov [103], see also El
et al. [104], Carbone et al. [51]. Soliton-based hydrodynamics will be
explained in Chapter 10 with more references to be added. The upswing of
generalized hydrodynamics can be traced back to the two independent seminal
contributions Castro-Alvaredo et al. [53] and Bertini et al.[22], in which a
general scheme for the computation of GGE averaged currents is presented, GGE
being the standard acronym for “generalized Gibbs ensemble”. Thus Euler type
equations could be written down based on convincing theoretical reasoning.

##  2 Dynamics of the classical Toda lattice

As well known at the time when Morikazu Toda started his studies on the
lattice with exponential interactions, shallow water waves in long channels
have peculiar dynamical properties. One observes solitary waves and soliton
collisions. In the latter, two incoming solitons emerge with their original
shape after an intricate dynamical process. The Korteweg-de Vries (KdV)
equation, a one-dimensional nonlinear wave equation, provides an accurate
theoretical description of these phenomena. In 1967 Toda investigated whether
discrete wave equations also might have solitary type dynamics. A point in
case is the standard wave equation and the harmonic lattice as its
discretization, both of them linear equations. With ingenious insight Toda
discovered that a lattice with exponential interactions, later known as Toda
lattice, exhibits the same dynamical features as the KdV equation. It took
another seven years until it had been firmly established that the $N$-particle
Toda lattice is indeed integrable with $N+1$ conservation laws.

In dimensionless form, the hamiltonian of the Toda lattice reads

$H_{\mathrm{to}}=\sum_{j\in\mathbb{Z}}\big{(}\tfrac{1}{2}p_{j}^{2}+\mathrm{e}^{-(q_{j+1}-q_{j})}\big{)},$
(2.1)

where $j\in\mathbb{Z}$ is the particle label and $(q_{j},p_{j})$ are position
and momentum of the $j$-th particle. One could introduce a particle mass, $m$,
coupling strength, $g$, and decay parameter, $\gamma$, as
$(p_{j}^{2}/2m)+g\exp(-\gamma(q_{j+1}-q_{j}))$ with $m,g,\gamma>0$. But
through rescaling spacetime the standard form (2.1) is recovered. The Toda
chain has no free parameters. Following Newton the equations of motion are

$\frac{\mathrm{d}}{\mathrm{d}t}q_{j}=p_{j},\quad\frac{\mathrm{d}}{\mathrm{d}t}p_{j}=\mathrm{e}^{-(q_{j}-q_{j-1})}-\mathrm{e}^{-(q_{j+1}-q_{j})},\quad
j\in\mathbb{Z}.$ (2.2)

Physically this equation can be viewed in two different ways. (i) The
displacements $q_{j}(t),j\in\mathbb{Z}$, are regarded as the lattice
discretization of a continuum wave field $q(x,t)$, $x\in\mathbb{R}$. We call
this the lattice, or field theory, picture. (ii) The fluid picture is to
literally view $q_{j}(t),j=1,\ldots,N$, as positions of particles moving on
the real line. However they do not interact pairwise as would be the case for
a real fluid. The hamiltonian is not invariant under relabelling of particles.
Toda was mostly thinking of a lattice discretization. But the fluid picture is
easier to visualize. Of course, there is only a single set of equations of
motion and one can switch back and forth between the two options.

At this stage, a review of the vast research on the Toda lattice can neither
be supplied nor is it intended. We would have to refer to original research
articles, monographs, reviews, and textbooks. However, it can be safely
summarized that almost exclusively problems have been studied for which
physically the chain is at zero temperature. Examples are multi-soliton
solutions and the spatial spreading of local perturbations of an initially
periodic particle configuration. In contrast, our focus are random initial
data with an energy proportional to system size and away from a ground state
energy. A paradigmatic set-up would be the thermal state at some non-zero
temperature. For such an enterprise novel techniques are required. In
particular, the issue of large system size has to be properly understood.

Figure 1: Displayed are the trajectories, $q_{j}(t)$, of 16 Toda particles
governed by $H_{\mathrm{cell},16}$ in a box of size $\ell=20$. Very roughly
the stretch $\nu=1.25$ and inverse temperature $\beta=1$, which corresponds to
medium pressure $P=0.7$. Clearly visible are two-particle collisions, but also
more complicated intertwined structures develop. The trajectory of the 9-th
particle is shown in magenta. Quasiparticles are introduced by the condition
to approximately maintain their velocity, except for collisions. One choice
for a quasiparticle trajectory with the same starting point as the 9-th
particle is shown in red. From [271].

As illustration in Figure 1 we show a numerical simulation of the particle
dynamics and note that the widely used visualization as a sequence of two-body
collision seems to be of limited value.

###  2.1 Locally conserved fields and their currents

Our first task is to elucidate the integrable structure of the Toda lattice.
For this purpose we introduce the stretch

$r_{j}=q_{j+1}-q_{j},$ (2.3)

also the free distance, or free volume, between particles $j$ and $j+1$, and
the Flaschka variables

$a_{j}=\mathrm{e}^{-r_{j}/2},\quad b_{j}=p_{j}.$ (2.4)

The stretch can have either sign, while $a_{j}>0$. The $a$’s and $b$’s are
conventional notation, but we will avoid the duplication of symbols by using
only the momentum $p_{j}$. In terms of these variables the equations of motion
read

$\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}.$
(2.5)

Hence $a_{j}$ couples to the right neighbor and $p_{j}$ to the left one. In
principle, we could have set $a_{j}=\tau\mathrm{e}^{-r_{j}/2}$ and $b_{j}=\tau
p_{j}$, which amounts to a mere time change. Flaschka picked
$\tau=\tfrac{1}{2}$. Below we will explain why $\tau=1$ is singled out in our
context.

For the purpose of thermodynamics, one first considers the lattice
$[1,\ldots,N]$ with periodic boundary conditions. In Flaschka variables they
amount to the obvious condition

$a_{0}=a_{N},\quad p_{N+1}=p_{1}.$ (2.6)

This choice is also called the closed or periodic chain, for which the phase
space is $\Gamma_{\\!N}^{\circ}=(\mathbb{R}_{+}\times\mathbb{R})^{N}$. Note
that

$\frac{d}{dt}\sum_{j=1}^{N}\log a_{j}^{2}=0,$ (2.7)

since carrying out the time derivative yields a telescoping sum which vanishes
by (2.6). Hence

$\sum_{j=1}^{N}r_{j}(t)=\ell,$ (2.8)

with some constant $\ell$, which can have either sign. Under this constraint,
an equivalent description is to consider the infinite Toda chain and to impose
the initial conditions

$q_{j+N}=q_{j}+\ell,\quad p_{j+N}=p_{j},$ (2.9)

for some $\ell\in\mathbb{R}$ and all $j\in\mathbb{Z}$, which then holds at any
time. Cutting the real line into cells, each of size $\ell$, in every cell
there are $N$ particles and their dynamics is governed by the cell hamiltonian

$H_{\mathrm{cell},N}=\sum_{j=1}^{N}\tfrac{1}{2}p_{j}^{2}+\sum_{j=1}^{N-1}\mathrm{e}^{-(q_{j+1}-q_{j})}+\mathrm{e}^{-(\ell-(q_{N}-q_{1}))}.$
(2.10)

The positions are unconstrained. Total momentum is conserved and, under the
dynamics generated by $H_{\mathrm{cell},N}$,

$\sum_{j=1}^{N}q_{j}(t)=\sum_{j=1}^{N}q_{j}(0)+t\sum_{j=1}^{N}p_{j}(0).$
(2.11)

The center of mass moves with constant velocity. The internal degrees of
freedom are the stretches $r_{1}(t),\ldots,r_{N}(t)$ subject to the time-
independent constraint (2.8). The stretches move in a potential which
increases exponentially in all directions.

In the literature, periodic boundary conditions are often stated as
$q_{N+1}=q_{1}$, which corresponds to the special case $\ell=0$. Physically
the parameter $\ell$ is of crucial importance because through it the stretch
per particle, $\ell/N$, is controlled. In the fluid picture the unit cell has
length $\ell$ and contains $N$ particles. The physical particle density is
$N/|\ell|$. But often it is more natural to work with the signed particle
density $N/\ell$, which can be negative.

Out of the Flaschka variables one forms the tridiagonal Lax matrix, $L_{N}$,

$L_{N}=\begin{pmatrix}p_{1}&a_{1}&0&\cdots&a_{N}\\\
a_{1}&p_{2}&a_{2}&\ddots&0\\\ 0&a_{2}&p_{3}&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&a_{N-1}\\\ a_{N}&0&\cdots&a_{N-1}&p_{N}\\\
\end{pmatrix},$ (2.12)

for $N\geq 3$, and its partner matrix

$\qquad B_{N}=\frac{1}{2}\begin{pmatrix}0&-a_{1}&0&\cdots&a_{N}\\\
a_{1}&0&-a_{2}&\ddots&0\\\ 0&a_{2}&0&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&-a_{N-1}\\\ -a_{N}&0&\cdots&a_{N-1}&0\\\
\end{pmatrix}.$ (2.13)

$L_{N},B_{N}$ is called a Lax pair. While $L$ is the common notation for a Lax
matrix, the partner matrix is also denoted by $A$ and $M$. Here we follow the
convention of Toda Toda [290]. Clearly, $L_{N}$ is symmetric and $B_{N}$ skew
symmetric, $(L_{N})^{\mathrm{T}}=L_{N}$, $(B_{N})^{\mathrm{T}}=-B_{N}$, with T
denoting the transpose of a matrix. Later on for the adjoint of an operator
also the more common ∗ will be used. From the equations of motion (2.5), one
verifies that

$\frac{d}{dt}L_{N}=[B_{N},L_{N}]$ (2.14)

with $[\cdot,\cdot]$ denoting the commutator,
$[B_{N},L_{N}]=B_{N}L_{N}-L_{N}B_{N}$. Since $B_{N}$ is skew symmetric,
$L_{N}(t)$ is isospectral to $L_{N}=L_{N}(0)$. Thus the eigenvalues of $L_{N}$
are conserved. Actually any matrix $B_{N}$ would do, because

$\frac{d}{dt}\mathrm{tr}\big{[}(L_{N})^{n}\big{]}=\sum_{j=0}^{n-1}\mathrm{tr}\big{[}(L_{N})^{j}[B_{N},L_{N}](L_{N})^{n-1-j}\big{]}=\mathrm{tr}\big{[}B_{N},(L_{N})^{n}\big{]}=0$
(2.15)

for all $n$.

$\blackdiamond\blackdiamond$ Vector notation. In our text various $N$-vectors
will appear. The standard notation $x\in\mathbb{R}^{N}$,
$x=(x_{1},\ldots,x_{N})$, is adopted. The $N$-dimensional volume element is
denoted by $\mathrm{d}^{N}x$. Also $x\in\mathbb{R}$ will be used. The
distinction should be obvious from the context. $\blackdiamond\blackdiamond$

$\blackdiamond\blackdiamond$ Notation for time-dependence. For time-dependent
quantities, as $X_{j}(t)$, we use the convention $X_{j}(0)=X_{j}$ and refer to
$X_{j}$ as time-zero field. While this notation is convenient, it might be
ambiguous. For example in (2.2) we should have used $q_{j}(t)$ with initial
condition $q_{j}$. We anticipate that the exact meaning will be clear from the
context. $\blackdiamond\blackdiamond$

$\blackdiamond\blackdiamond$ Phase spaces. It is recommendable to keep track
of phase spaces. We will use $\Gamma$ as a generic symbol and $\Gamma_{n}$ to
indicate a phase space of dimension $2n$. More specifically the notation is
$\Gamma_{\\!N}=\mathbb{R}^{N}\times\mathbb{R}^{N}$,
$\Gamma_{\\!N}^{\circ}=\mathbb{R}_{+}^{N}\times\mathbb{R}^{N}$,
$\Gamma_{\\!N}^{\diamond}=\mathbb{R}_{+}^{(N-1)}\times\mathbb{R}^{N}$, and
$\Gamma_{\\!N}^{\triangleright}=\mathbb{W}_{N}\times\mathbb{R}^{N}$ with Weyl
chamber $\mathbb{W}_{N}=\\{q_{1}<\ldots<q_{N}\\}$.
$\blackdiamond\blackdiamond$

Let us write the eigenvalue problem for $L_{N}$ as

$L_{N}\psi_{\alpha}=\lambda_{\alpha}\psi_{\alpha}$ (2.16)

with $\alpha=1,\dots,N$. Then $\lambda_{\alpha}(a,p)$ is some function on
phase space which does not change under the Toda time evolution. However
$\lambda_{\alpha}$ is a highly nonlocal function, in general. For example,
considering the dependence of $\lambda_{\alpha}$ only on $p_{1}$ and
$p_{N/2}$, this will not split into a sum as $g(p_{1})+\tilde{g}(p_{N/2})$
even approximately. Physically more relevant are local conservation laws. For
the Toda lattice they are easily obtained through forming the trace as

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

The second identity confirms that $Q^{[n],N}$ is conserved and the fourth
identity that the density $Q^{[n],N}_{j}$ is local. Indeed, taking $n<N/2$ and
expanding as

$Q^{[n],N}_{j}=\sum_{j_{1}=1}^{N}\dots\sum_{j_{n-1}=1}^{N}(L_{N})_{j,j_{1}}(L_{N})_{j_{1},j_{2}}\dots(L_{N})_{j_{n-1},j},$
(2.18)

the density $Q^{[n],N}_{j}$ depends only on the variables
$\\{a_{j-i},p_{j-i},\ldots,a_{j+i},p_{j+i},i=0,\ldots,n-1\\}$, modulo $N$.

The two-sided infinite volume limit of $L_{N}$ is the tridiagonal Lax matrix
$L$, and correspondingly the partner matrix $B$, which are now operators
acting on the Hilbert space $\ell_{2}(\mathbb{Z})$ of square-integrable two-
sided sequences over the lattice $\mathbb{Z}$. The Lax pair still satisfies

$\frac{d}{dt}L=[B,L].$ (2.19)

Of course $\mathrm{tr}[L^{n}]$ makes no sense, literally. But the infinite
volume density

$Q^{[n]}_{j}=L_{j,j}$ (2.20)

is well defined. $Q^{[n]}_{j}$ is a finite polynomial in the variables
$\\{a_{i},p_{i},|i-j|\leq n-1\\}$, which can be grasped more easily by using
the random walk expansion derived from the $N=\infty$ version of (2.18). The
walk is on $\mathbb{Z}$ with $n$ steps of step size $0,\pm 1$, starting and
ending at $j$. A step from $i$ to $i$ carries the variable $p_{i}$ and the
step either from $i$ to $i+1$ or from $i+1$ to $i$ carries the variable
$a_{i}$. For a given walk one forms the product along the path, which is a
monomial of degree $n$, compare with Figure 2. $Q^{[n]}_{j}$ is then obtained
by summing over all admissible walks. By translation invariance of the model,
$Q^{[n]}_{j}$ and $Q^{[n]}_{j+i}$ are identical polynomials, except that the
particle label is shifted by $i$. Just as for the hamiltonian (2.1), formally
we still write

$Q^{[n]}=\sum_{j\in\mathbb{Z}}Q^{[n]}_{j},$ (2.21)

where the superscript $[n]$ ranges over positive integers.

Figure 2: A random walk path contributing to $Q^{[6]}_{j}$ and, according to
the rules, carrying the weight
$a_{j}p_{j+1}a_{j}a_{j-1}a_{j-1}p_{j}=a_{j-1}^{2}a_{j}^{2}p_{j}p_{j+1}$.

As remarked already, in addition the stretch is locally conserved, which
carries the label $0$,

$Q^{[0]}_{j}=r_{j},\qquad Q^{[0],N}=\sum_{j=1}^{N}r_{j}.$ (2.22)

The three lowest order fields have an immediate physical interpretation as
stretch, momentum, and energy density,

$Q^{[0]}_{j}=r_{j},\quad
Q^{[1]}_{j}=p_{j},\quad\tfrac{1}{2}Q^{[2]}_{j}=\tfrac{1}{2}\big{(}p_{j}^{2}+a_{j}^{2}+a_{j-1}^{2}\big{)}.$
(2.23)

Obviously, there cannot be physical names for higher orders. In the community
of quantum integrable systems, $Q^{[n]}$ is called the $n$-th conserved charge
or merely the $n$-th charge, which serves as a concise notion but carries no
specific physical meaning. We will use local conserved field and local charge
interchangeably.

Since the Toda hamiltonian has a local energy density, any locally conserved
field must satisfy a continuity equation, in other words the lattice version
of a local conservation law. We consider the infinite lattice and compute

$\frac{d}{dt}Q_{j}^{[n]}=(BL^{n}-L^{n}B)_{j,j}=a_{j-1}(L^{n})_{j,j-1}-a_{j}(L^{n})_{j+1,j}=J_{j}^{[n]}-J_{j+1}^{[n]}.$
(2.24)

Hence $J_{j+1}^{[n]}$ is the current of the $n$-th conserved field from $j$ to
$j+1$ and $J_{j}^{[n]}$ the current from $j-1$ to $j$. Defining the lower
triangular matrix $L^{\scriptscriptstyle\downarrow}$ by
$(L^{\scriptscriptstyle\downarrow})_{j+1,j}=a_{j}$ for all $j$ and
$(L^{\scriptscriptstyle\downarrow})_{i,j}=0$ otherwise, a more concise
expression is

$J_{j}^{[n]}=(L^{n}L^{\scriptscriptstyle\downarrow})_{j,j},$ (2.25)

$n=1,2,\ldots\,$. For the stretch

$\frac{d}{dt}Q_{j}^{[0]}=-p_{j}+p_{j+1},\quad J_{j}^{[0]}=-p_{j},$ (2.26)

and hence

$J_{j}^{[0]}=-Q_{j}^{[1]}.$ (2.27)

This innocent looking equation will have surprising consequences.

$\blackdiamond\blackdiamond$ Ambiguity of densities. In the way just
presented, the densities, both for field and current, seem to be unique. This
is not the case, however. We have made a particular choice which will be
useful when investigating the hydrodynamic scale of the Toda lattice. A
further common scheme is to require the density to depend on a minimal number
of lattice sites. In terms of the random walk representation described above,
this amounts to a density $\tilde{Q}_{j}^{[n],N}$ defined by summing over all
closed admissible paths with minimum equal to $j$. As an example, for $n=4$
the minimal version of the density is given by

$\tilde{Q}_{j}^{[4]}=p_{j}^{4}+4(p_{j}^{2}+p_{j}p_{j+1}+p_{j+1}^{2})a_{j}^{2}+2a_{j}^{4}+4a_{j}^{2}a_{j+1}^{2}.$
(2.28)

Rather than trying to dwell on generalities, we illustrate the issue by
considering the energy density, $n=2$. From (2.18) we have

$Q_{j}^{[2]}=p_{j}^{2}+a_{j-1}^{2}+a_{j}^{2}$ (2.29)

with the current density

$J_{j}^{[2]}=a_{j-1}^{2}(p_{j-1}+p_{j}).$ (2.30)

The minimal version of the energy density would be

$\tilde{Q}_{j}^{[2]}=p_{j}^{2}+2a_{j}^{2},$ (2.31)

having the current density

$\tilde{J}_{j}^{[2]}=2a_{j-1}^{2}p_{j}.$ (2.32)

At infinite volume, we consider the spatial sums
$\sum_{j=1}^{N}Q_{j}^{[2]}=Q^{[2],N\infty}$ and
$\sum_{j=1}^{N}\tilde{Q}_{j}^{[2]}=\tilde{Q}^{[2],N\infty}$. They differ only
by a boundary term and hence the spatial average
$N^{-1}\big{(}Q^{[2],N\infty}-\tilde{Q}^{[2],N\infty}\big{)}=\mathcal{O}(1/N)$.
On the other hand, the corresponding total currents, $J^{[2],N\infty}$ and
$\tilde{J}^{[2],N\infty}$, differ by $\mathcal{O}(N)$, but

$J_{j}^{[2]}-\tilde{J}_{j}^{[2]}=\frac{d}{dt}(a_{j-1})^{2}.$ (2.33)

The difference is a total time derivative and thus vanishes when averaged over
a time-stationary probability measure, e.g. thermal average. As conclusion,
while there is ambiguity on the microscopic scale, upon averaging over large
spacetime cells this amounts to only small surface type correction terms. In
particular, the hydrodynamic equations for the Toda lattice do not depend on
the particular choice of microscopic densities. $\blackdiamond\blackdiamond$

###  2.2 Action-angle variables, notions of integrability

####  2.2.1 Conventional integrability of the Toda lattice

A cornerstone of hamiltonian dynamics is the abstract characterization of
integrable systems. Just to recall, given is some hamiltonian, $H$, on a phase
space $\Gamma_{n}$ of dimension $2n$. The dynamics generated by $H$ is called
integrable, if there are $n$ differentiable functions, , the action variables,
on phase space which have the following properties: (i) They are conserved,
which means that the Poisson brackets $\\{I_{j},H\\}=0$. (ii) They are in
involution, i.e. $\\{I_{i},I_{j}\\}=0$ for all $i,j=1,\ldots,n$. (iii)
$I_{1},\ldots,I_{n}$ span a $n$-dimensional hyper-surface in $\Gamma_{n}$,
which is compact and connected. In particularly, no scattering orbits are
permitted. The hypersurface is assumed to be invariant under the hamiltonian
flow generated by $I_{j}$, for all $j$. So to speak, as regards to the
dynamics generated by $I_{j}$, the hypersurface has is no boundary. Then the
Arnold-Liouville theorem states that there exists a canonical transformation
to action variables $(I_{1},\ldots,I_{n})$ and the canonically conjugate angle
variables $(\vartheta_{1},\ldots,\vartheta_{n})\in\mathbb{T}^{n}$,
$\mathbb{T}^{n}$ the $n$-dimensional torus, such that the transformed
hamiltonian $\tilde{H}$ depends only on $I$. In these variables the dynamics
trivializes as

$\frac{d}{dt}\vartheta_{j}=\omega_{j},\qquad\omega_{j}=\partial_{I_{j}}\tilde{H}(I)$
(2.34)

for $j=1,\ldots,n$. The angle $\vartheta_{j}$ moves on the unit torus with
frequency $\omega_{j}$.

This characterization applies also to the closed Toda chain. We consider a
lattice of $N$ sites, the phase space
$\Gamma_{\\!N}^{\circ}=(\mathbb{R}_{+}\times\mathbb{R})^{N}$, and the
evolution (2.5) in terms of the Flaschka variables $a=(a_{1},\ldots,a_{N})$
and momenta $p=(p_{1},\ldots,p_{N})$. These are not canonical variables.
However, instead of the usual Poisson bracket, one can introduce a nonstandard
Poisson bracket by first defining the $N\times N$ matrix

$\qquad A_{N}=\frac{1}{2}\begin{pmatrix}-a_{1}&0&\cdots&0&a_{N}\\\
a_{1}&-a_{2}&0&\ddots&0\\\ 0&a_{2}&-a_{3}&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&0\\\ 0&\cdots&0&a_{N-1}&-a_{N}\end{pmatrix}$
(2.35)

and adjusting the Poisson bracket to

$\\{f,g\\}=\langle\nabla_{p}f,A_{N}\nabla_{a}g\rangle-\langle\nabla_{p}g,A_{N}\nabla_{a}f\rangle$
(2.36)

with $\langle\cdot,\cdot\rangle$ denoting the inner product in
$\mathbb{R}^{N}$. For the usual Poisson bracket, $A_{N}$ would be the identity
matrix. In the Flaschka variables

$H_{\mathrm{to},N}=\sum_{j=1}^{N}\big{(}\tfrac{1}{2}p_{j}^{2}+a_{j}^{2}\big{)}$
(2.37)

and the equations of motion (2.5) can be written in hamiltonian form as

$\frac{d}{dt}p_{j}=\\{p_{j},H_{\mathrm{to},N}\\},\qquad\frac{d}{dt}a_{j}=\\{a_{j},H_{\mathrm{to},N}\\}.$
(2.38)

The matrix $A_{N}$ is of rank $N-1$, the oblique projection for the eigenvalue
$0$ being

$|a_{1}^{-1},a_{2}^{-1},\ldots,a_{N}^{-1}\rangle\langle 1,1,\ldots,1|,$ (2.39)

i.e. the corresponding left eigenvector of $A_{N}$ equals
$-2\nabla_{a}Q^{[0],N}$ and the right one $\nabla_{p}Q^{[1],N}$. The Poisson
structure (2.36) is degenerate. But it can be turned non-degenerate simply by
fixing the two conservation laws as $Q^{[0],N}=c_{0}$, $Q^{[1],N}=c_{1}$ with
an arbitrary choice of the real parameters $c_{0},c_{1}$. The new phase space
becomes $\Gamma_{\\!{N-1}}^{\circ}$ and the dynamical evolution equations
involve only the variables $(a_{1},\ldots,a_{N-1},p_{1},\ldots,p_{N-1})$,
which are a hamiltonian system with a non-degenerate Poisson bracket
structure. The phase space for the action-angle variables can be chosen as
$\Gamma_{N-1}$ with coordinates $(x_{1},\ldots,x_{N-1},y_{1},\ldots,y_{N-1})$
and corresponding action-angle variables $(\vartheta,I)$, with
$I_{j}\in\mathbb{R}_{+},\vartheta_{j}\in[0,2\pi]$, $j=1,\ldots,N-1$, as

$x_{j}=\sqrt{I_{j}}\cos\vartheta_{j},\quad
y_{j}=\sqrt{I_{j}}\sin\vartheta_{j},$ (2.40)

which are known as global Birkhoff coordinates. As proved in 2008 by A.
Henrici and T. Kappeler, there is a canonical transformation
$\Phi:\Gamma_{\\!{N-1}}^{\circ}\to\Gamma_{N-1}$ such that the transformed
hamiltonian, $H_{\mathrm{aa}}$, depends only on the action variables,
$H_{\mathrm{aa}}=H_{\mathrm{aa}}(I_{1},\ldots,I_{N-1})$. In fact, as for us
crucial property, $H_{\mathrm{aa}}$ is a strictly convex, real-analytic
function. This means that the phases
$\omega_{j}=\partial_{I_{j}}H_{\mathrm{aa}}$ are incommensurate Lebesgue
almost surely. In other words, $H_{\mathrm{aa}}$ has no linear pieces, as
would be the case for a system of harmonic oscillators. The Toda lattice is
phase-mixing: starting with some probability measure on $\Gamma_{N-1}^{\circ}$
with a continuous density function, in the long time limit the density will
become uniform on almost every torus of dimension $N-1$ with an amplitude
computed from the initial density. The “almost every” is required because of
tori with commensurate frequencies, which however form a set of Lebesgue
measure zero.

$\blackdiamond\blackdiamond$ Pitfalls of classical integrability. We consider
particles on the real line governed by the standard hamiltonian

$H_{\mathrm{mec},N}=\sum_{j=1}^{N}\tfrac{1}{2}p_{j}^{2}+\sum_{i<j=1}^{N}V_{\mathrm{mec}}(q_{i}-q_{j}).$
(2.41)

The mechanical interaction potential is assumed to be even,
$V_{\mathrm{mec}}(x)=V_{\mathrm{mec}}(-x)$, and repulsive,
$V_{\mathrm{mec}}^{\prime}(x)<0$ for $x>0$. Furthermore the potential decays
at infinity such that $|V_{\mathrm{mec}}(x)|<x^{-\gamma}$ for large $x$ with
some $\gamma>1$ and diverges at the origin, thereby ensuring that particles do
not cross. Hence the phase space equals $\Gamma_{\\!N}^{\triangleright}$.
Under such assumptions it is proved that asymptotic momenta exist,

$\lim_{t\to\pm\infty}p_{j}(t)=p_{j}^{\pm},$ (2.42)

and also the respective scattering shifts

$\lim_{t\to\pm\infty}q_{j}(t)-p_{j}^{\pm}t=\phi_{j}^{\pm},\qquad\phi_{j}^{\pm}\in\mathbb{R},$
(2.43)

as functions of the initial conditions. This defines the scattering map
$\Phi^{-1}:(q,p)\mapsto(p^{\pm},\phi^{\pm})$. The asymptotic momenta are
ordered as $p_{1}^{+}<\ldots<p_{N}^{+}$ and the scattering map is one-to-one
on $\Gamma_{\\!N}^{\triangleright}$. The scattering map is canonical, which
implies that Poisson brackets are conserved. In particular,

$\\{p_{i}^{+},p_{j}^{+}\\}=0$ (2.44)

and

$H_{\mathrm{mec},N}\circ\Phi=\sum_{j=1}^{N}\tfrac{1}{2}(p_{j}^{+})^{2}.$
(2.45)

For the past asymptotic momenta are anti-ordered and the corresponding
properties hold as well.

Clearly, the family $\\{p_{j}^{+},j=1,\ldots,N\\}$ is in involution and
conserved. As premature reaction the mechanical system could be classified as
integrable. However, the Arnold-Liouville theorem does not apply. Instead of
quasi-periodic motion on tori, according to (2.45), in action-angle variables
the dynamics trivializes as $\phi^{+}_{j}+p_{j}^{+}t$. For periodic boundary
conditions integrability can be tested through molecular dynamics simulations.
For a generic choice of $V_{\mathrm{mec}}$ the particle system will thermalize
in the long time limit. Such dynamics is very different from the one of the
Toda fluid on a ring. $N$-body integrability requires fine-tuning of the
interaction potential. One option would be to ask for a Lax matrix. Then,
under our conditions, it is known that the only choices are
$V_{\mathrm{mec}}(x)=\sinh^{-2}(x)$ and $V_{\mathrm{mec}}(x)=x^{-2}$, see
Chapter 11 for more details. But this option is ad hoc without link to
standard definitions. In the following subsection we will argue that for
hydrodynamic purposes the natural defining property is a quasilocal density
for the conserved fields. $\blackdiamond\blackdiamond$

####  2.2.2 Hydrodynamic perspective on integrability

The Toda lattice is a very peculiar dynamical system in the sense that it is
integrable for every system size $N$, which we call integrable many-particle
or integrable many-body. Now the large $N$ limit is in focus and from a
physics perspective the conventional definition of integrability might have to
be reconsidered. This is even more urgent, since the naive extension of
classical integrability to quantum systems fails. A further desideratum would
be a notion referring directly to the infinite lattice. For hydrodynamics the
central building blocks are local conservation laws, to be more precise the
hamiltonian and the conservation laws are constructed from a strictly local
density. We thus propose to call an infinitely extended system nonintegrable,
if it admits only a few strictly local conservation laws. The system is called
integrable, if it possesses an infinite number of linearly independent local
conservation laws.

Starting from a strictly local density supported on an interval of $m$ sites,
the property to be the density of a conservation law refers to a phase space
of dimension $2(n+\ell)$, in case the hamiltonian density is supported on
$\ell$ sites. The condition of being in involution is no longer mentioned, but
it seems to hold in concrete examples, possibly after first adjusting either
the classical phase space or the Poisson bracket. Note that local conservation
laws have a linear structure, in the sense that the sum of two local
conservation laws is again a local conservation law.

For the Toda lattice the Lax matrix is tridiagonal implying strictly local
conservation laws. Another integrable model for a fluid is the Calogero system
with a $1/\sinh^{2}$ pair potential. This potential decays exponentially and
the Lax matrix is fully occupied. Hence the conserved fields of the Calogero
fluid cannot be strictly local. Physically, the hydrodynamic scale has to
include quasilocal fields with densities having exponential tails. However,
depending on the model, the precise borderline could be a subtle issue.

$\blackdiamond\blackdiamond$ Toda local conservation laws. The Toda chain is
integrable in the hydrodynamic sense with densities of the locally conserved
fields stated in (2.20) and (2.22). However, as a stronger property, one would
like to establish that there are no further local conservation laws. Currently
this is a conjecture and more studies would be needed. Still, a precise
formulation is worthwhile. As before, periodic boundary conditions are
understood. We assume some general density function $f$ of support $\kappa$,
in other words $f(a_{1},p_{1},\ldots,a_{\kappa},p_{\kappa})$. Then the shifted
densities are
$f_{j}(a,p)=f(a_{j+1},p_{j+1},\ldots,a_{j+\kappa},p_{j+\kappa})$, $f_{0}=f$,
and the conditions for being a local conservation law read

$Q^{[f],N}=\sum_{j=1}^{N}f_{j},\qquad\\{Q^{[f],N},H_{\mathrm{to},N}\\}=\sum_{j=1}^{N}\\{f_{j},H_{\mathrm{to},N}\\}=0$
(2.46)

for $N>2\kappa$. If so, each of the Poisson brackets is a local function. By
translation invariance the sum has to be telescoping and thus necessarily
there exists a current function, $J_{j}$, of support of size $\kappa+2$, such
that $\\{f_{j},H_{\mathrm{to},N}\\}=J_{j}-J_{j+1}$. According to the already
proven conventional integrability there must be some function, $G$, such that

$Q^{[f],N}=G\big{(}Q^{[0],N},\ldots,Q^{[N],N}\big{)}.$ (2.47)

Our conjecture claims that $G$ is necessarily linear.
Conjecture: For fixed $\kappa$ and $N$ sufficiently large, there exists
coefficients $c_{0},c_{1},\ldots,c_{\kappa}$ such that

$Q^{[f],N}=\sum_{m=0}^{\kappa}c_{m}Q^{[m],N}.$ (2.48)

One argument in favor of the conjecture comes from a simple observation.
Consider some locally conserved field, $Q^{[n],N}$. Then $(Q^{[n],N})^{2}$ is
also conserved, but no longer local. The condition of locality should be
strong enough to force a linear function in (2.47).

$\blackdiamond\blackdiamond$

The reader might find the evidence for linking integrability and conservation
laws not particularly convincing. Agreed, but it should be noted that our
definition translates one-to-one to quantum spin chains, and also continuum
quantum models. In the former case there are two concrete results strongly
supporting the hydrodynamic notion of integrability. Considered is the XYZ
spin chain with couplings $J_{x},J_{y},J_{z}$ and external magnetic field,
$h$, pointing in the $z$-direction. In terms of the Pauli spin-$\tfrac{1}{2}$
matrices, $\sigma^{x},\sigma^{y},\sigma^{z}$, the hamiltonian reads

$H_{\mathrm{XYZ}}=\sum_{j\in\mathbb{Z}}\big{(}J_{x}\sigma_{j}^{x}\sigma_{j+1}^{x}+J_{y}\sigma^{y}_{j}\sigma_{j+1}^{y}+J_{z}\sigma^{z}_{j}\sigma^{z}_{j+1}-h\sigma_{j}^{z}\big{)}.$
(2.49)

This model is integrable for $h=0$ and for $h\neq 0$ in case $J_{x}=J_{y}$,
the XXZ model with external field. The model is expected to be non-integrable
for any other choice of parameters. Now, for the case of non-zero coupling
constants and $h\neq 0$, $J_{x}\neq J_{y}$, it is proved that there is only a
single local conservation law, namely the hamiltonian $H_{\mathrm{XYZ}}$
itself. This is, so to speak, the fully chaotic case. Only energy is
transported and the chain thermalizes, in the sense that expectations of local
observables converge to the thermal average in the long time limit.

On the other hand for $h=0$ and nonvanishing couplings the model is
integrable. Most effectively the strictly local conserved charges are computed
through the boost operator. The second result states that, if the coupling
constants are non-vanishing, then any local conservation law is a finite
linear combination of the already known conservation laws. This is the precise
analogue of our conjecture. However, around 2014, for the XXZ chain with
magnetic field it was discovered that our notion of locality is indeed too
restrictive. The XXZ chain possesses in addition quasilocal charges which have
to be included in a hydrodynamic description, see Notes for details.

In passing we note that the XXX chain at $h=0$ is integrable and all three
spin components, $\sigma_{j}^{x},\sigma_{j}^{y},\sigma_{j}^{z}$, are
conserved. However they do not commute with each other. There is still a tower
of local charges, $Q_{2},Q_{3},\ldots$ in the usual notation, which commute
with each other and with each spin-component. Strictly speaking the involution
property is violated.

Our discussion raises some difficult issues. From the available evidence,
there is a dichotomy, either a few conservation laws or infinitely many. One
does not know whether there is a deep reason behind or merely reflects the
limited class of models studied. Personally I believe in the first option. In
this context, particularly intriguing are nearly integrable systems. For
finite $N$, the KAM theorem provides information on the stability of the
invariant tori, confirming the coexistence of chaotic and integrable regions
in phase space. But in our definition the limit $N\to\infty$ is taken first
and integrable regions might be rare.

###  2.3 Scattering theory

In hydrodynamics the system is confined and thus interactions persist without
interruption, which in the long time limit then leads to some sort of
statistical equilibrium. A dynamically distinct set-up is scattering: in the
distant past particles are in the incoming configuration, for which they are
far apart and do not interact. A time span of multiple collision processes
follows. In the far future the outgoing particles move freely again. If the
interaction potential is repulsive, no bound states can be formed. To
illustrate the special features of scattering for a integrable many-body
systems, we first discuss a fluid consisting of hard rods, which will serve as
an instructive example also later on.

####  2.3.1 Hard rod fluid

We consider $N$ hard rods, rod length $\mathsfit{a}>0$, moving on the real
line. The hamiltonian reads

$H_{\mathrm{hr},N}^{\diamond}=\sum_{j=1}^{N}\tfrac{1}{2}p_{j}^{2}+\sum_{j=1}^{N-1}V_{\mathrm{hr}}(q_{j+1}-q_{j}),$
(2.50)

with the hard rod potential $V_{\mathrm{hr}}(x)=\infty$ for
$|x|<\mathsfit{a}/2$ and $V_{\mathrm{hr}}(x)=0$ for $|x|\geq\mathsfit{a}/2$.
Hard rods collide elastically with their two neighbors, except for the border
particles with only one neighbor. Obviously, the system is integrable with
one-particle sum functions, $\sum_{j=1}^{N}f(p_{j})$, being conserved.

Since particles cannot cross, we order $q_{1}(t)<q_{2}(t)<\ldots<q_{N}(t)$.
For sufficiently long times towards future and past,

$q_{j}(t)=p_{j}^{+}t+\phi_{j}^{+},\quad t\to\infty,\qquad
q_{j}(t)=p_{j}^{-}t+\phi_{j}^{-},\quad t\to-\infty,$ (2.51)

where $\phi_{j}^{+}$ is the forward in time and $\phi_{j}^{-}$ the backward
scattering shift. When comparing with the point dynamics, $\mathsfit{a}\to 0$,
one concludes

$p_{j}^{+}=p_{N-j+1}^{-},\qquad p_{N}^{-}<\ldots<p_{1}^{-},\quad
p_{1}^{+}<\ldots<p_{N}^{+}.$ (2.52)

The relative scattering shift of particle $j$, $\kappa_{j}$, is defined by the
deviation relative to the point particle dynamics. Hence

$\kappa_{j}=\phi_{(}N-j+1)^{+}-\phi_{j}^{-}.$ (2.53)

Just considering the intersection points arising from $N$ straight spacetime
lines, one concludes

$\kappa_{j}=\sum_{1\leq i\leq N,i\neq
j}\mathrm{sgn}(p_{j}^{-}-p_{i}^{-})\phi_{\mathrm{hr}}(p_{j}^{-}-p_{i}^{-}),\qquad\phi_{\mathrm{hr}}(w)=-\mathsfit{a},$
(2.54)

when written for incoming momenta. A similar expression holds for outgoing
momenta. As a hallmark of integrable many-body systems, the relative
scattering shift is the properly signed sum of two-particle scattering shifts.

More intuitive is the notion of a quasiparticle, which maintains its velocity
through a collision. Since for hard rods the collision time vanishes, a
quasiparticle moves along a straight line interrupted by jumps of size
$\mathsfit{a}$ either to the right or left. Quasiparticles are ordered
increasingly according to the ingoing particle configuration. Then the
relative scattering shift $\kappa_{j}$ is the accumulated spatial shift of the
$j$-th quasiparticle.

Our definitions involve sign conventions, which vary from author to author. As
illustrated in Figure 3, by our rules, a negative scattering shift means that
through a collision the two incoming particles get pushed apart relative to
the free particle motion, while for a positive scattering shift they get
pulled closer. Hence the two-particle scattering shift of the hard rod fluid
equals $-\mathsfit{a}$, which corresponds to the relative scattering shift of
particle 1 in case incoming momenta are ordered as $p_{2}<p_{1}$.

Figure 3: Two-particle scattering. The numbers refer to particle labels. The
red trajectory is quasiparticle 1 and blue quasiparticle 2. Quasiparticle 1 is
shifted by $\phi_{12}$, negative in the example, and quasiparticle 2 by
$-\phi_{12}$.

For physical rods the length is positive. But our rules make perfectly sense
also for negative $\mathsfit{a}$. The two rods pass through each other until
they reach the now negative distance $\mathsfit{a}$. At that moment the
velocities are exchanged. Since all rods have the same length, no ambiguities
arise. Without much ado, the assumption $\mathsfit{a}>0$ will be dropped.
However, in our verbal explanations we have to be more cautious.

####  2.3.2 Two-particle Toda lattice

Scattering is physically more intuitive in the fluid picture, i.e. particles
move on the real line, also called the open chain. For two particles, the
equations of motion are

$\ddot{q}_{1}(t)=-\mathrm{e}^{q_{1}-q_{2}},\quad\ddot{q}_{2}(t)=\mathrm{e}^{q_{1}-q_{2}}.$
(2.55)

The relative motion, $q_{2}(t)-q_{1}(t)$, corresponds to a single particle
subject to the potential $2\mathrm{e}^{-x}$. We impose the asymptotic
conditions $p_{1}(-\infty)=p_{1}^{-}$ and $p_{2}(-\infty)=p^{-}_{2}$ with
$p_{2}^{-}<p_{1}^{-}$. Adjusting the initial time such that
${q}_{2}(0)-{q}_{1}(0)$ is at the turning point, i.e.
$\dot{q}_{2}(0)-\dot{q}_{1}(0)=0$, the solution to (2.55) becomes

$\displaystyle
q_{1}(t)=\tfrac{1}{2}(p^{-}_{1}+p^{-}_{2})t-\log\big{(}\gamma^{-1}\cosh(\gamma
t)\big{)},$ $\displaystyle
q_{2}(t)=\tfrac{1}{2}(p^{-}_{1}+p^{-}_{2})t+\log\big{(}\gamma^{-1}\cosh(\gamma
t)\big{)},$ (2.56)

with $\gamma=\tfrac{1}{2}(p^{-}_{1}-p^{-}_{2})>0$. The large time asymptotics
is given by

$q_{1}(t)=\begin{cases}p^{-}_{1}t+\log|p^{-}_{1}-p^{-}_{2}|,\\\\[4.30554pt]
p^{-}_{2}t+\log|p^{-}_{1}-p^{-}_{2}|,\end{cases}q_{2}(t)=\begin{cases}p_{2}^{-}t-\log|p^{-}_{1}-p^{-}_{2}|,&\quad
t\to-\infty,\\\\[4.30554pt] p^{-}_{1}t-\ \log|p^{-}_{1}-p^{-}_{2}|,&\quad
t\to\infty.\end{cases}$ (2.57)

Since quasiparticle one is shifted by $2\log|p^{-}_{1}-p^{-}_{2}|$ and
quasiparticle two by $-2\log|p^{-}_{1}-p^{-}_{2}|$, we conclude that the Toda
two-particle relative scattering shift is given by

$2\log|p^{-}_{1}-p^{-}_{2}|=\phi_{\mathrm{to}}(p^{-}_{1}-p^{-}_{2}).$ (2.58)

The scattering shift has no definite sign. For $|p^{-}_{1}-p^{-}_{2}|=1$ the
scattering shift vanishes. For $|p^{-}_{1}-p^{-}_{2}|<1$ the scattering shift
is negative, just as for hard rods. The trajectories of the two Toda particles
look similar to the ones of hard rods, but the hard rod zero collision time is
smeared to an exponential with rate $\gamma$. For $|p^{-}_{1}-p^{-}_{2}|>1$
the scattering shift is positive. The trajectories spatially cross each other,
still approaching their asymptotic motion exponentially fast.

####  2.3.3 $N$-particle Toda lattice

. For $N$ particles the hamiltonian of the open chain reads

$H^{\diamond}_{\mathrm{to},N}=\sum_{j=1}^{N}\tfrac{1}{2}p_{j}^{2}+\sum_{j=1}^{N-1}\mathrm{e}^{-(q_{j+1}-q_{j})},$
(2.59)

where the superscript ⋄ is used to indicate the open chain. In the Flaschka
variables the equations of motion become

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

with the boundary conditions $a_{0}=0,a_{N}=0$. The Lax matrix
$L^{\diamond}_{N}$ equals $L_{N}$ except for $a_{N}=0$, correspondingly for
the partner matrix $B_{N}^{\diamond}$. The time evolution is encoded as

$\frac{d}{dt}L_{N}^{\diamond}=[B_{N}^{\diamond},L_{N}^{\diamond}].$ (2.61)

In general, modifying boundary conditions is likely to break integrability.
But the open Toda chain is still integrable.

Figure 4: Scattering shifts $\kappa_{j}$. The numbers refer to quasiparticle
labels. In the left part the sequence of collisions is
$\phi_{23},\phi_{13},\phi_{12}$ and in the right part
$\phi_{12},\phi_{13},\phi_{23}$. Hence $\kappa_{1}=\phi_{12}+\phi_{13}$,
$\kappa_{2}=-\phi_{12}+\phi_{23}$, and $\kappa_{3}=-\phi_{23}-\phi_{13}$,
independently of the particular sequence.

The time zero phase point is denoted by $(q,p)\in\Gamma_{\\!N}$. Since
particles repel, for sufficiently large $t$ one has $q_{1}(t)<\ldots<q_{N}(t)$
and $p_{1}(t)<\ldots<p_{N}(t)$. We order the eigenvalues of $L_{N}^{\diamond}$
as $\lambda_{1}<\ldots<\lambda_{N}$, which defines the Weyl chamber
$\mathbb{W}_{N}$. Then, in the limit $t\to\infty$,

$q_{j}(t)\simeq\lambda_{N-j+1}t+\phi_{j}^{+}$ (2.62)

with the forward scattering shift $\phi^{+}_{j}$. Correspondingly in the past,

$q_{j}(t)\simeq\lambda_{j}t+\phi_{j}^{-}$ (2.63)

for $t\to-\infty$. The relative scattering shift is given by

$\kappa_{j}=\sum_{1\leq i\leq N,i\neq
j}\mathrm{sgn}(\lambda_{j}-\lambda_{i})\phi_{\mathrm{to}}(\lambda_{j}-\lambda_{i}),$
(2.64)

compare with (2.54) for which the same conventions are used. In view of Figure
1 this result is very surprising. Despite the intricate pattern of multiple
collisions, at very long times particles manage to have a scattering shift
which is the weighted sum of two-particle scattering shifts. As one
consequence of (2.64), the scattering shift does not depend on the order of
collisions. For quantum mechanical many-body systems this property is known as
Yang-Baxter relation, see Chapter 14, and most commonly illustrated for three
particles as in Figure 4.

We return to the scattering map $\Phi^{-1}$ defined through (2.62) and for
notational simplicity denote $\phi^{+}$ by $\phi$. As limit of canonical
transformations the scattering map $\Phi^{-1}:(q,p)\mapsto(\lambda,\phi)$ is
symplectic, which means that the asymptotic momenta and the scattering shifts
are canonical coordinates. In other words, as functions on $\Gamma_{\\!N}$ the
Poisson brackets read

$\\{\lambda_{i},\lambda_{j}\\}=0,\quad\\{\phi_{i},\phi_{j}\\}=0,\quad\\{\phi_{i},\lambda_{j}\\}=\delta_{ij}.$
(2.65)

The variables $(\lambda,\phi)$ are the scattering analogue of action-angle
variables. Only the angles vary over $\mathbb{R}$, rather than taking values
on a torus. To make this distinction also verbally one should speak of
variables for asymptotic momenta and scattering shifts. But such practice
becomes unwieldy and the common usage is scattering coordinates and more
specifically action-angle variables, keeping in mind that a scattering
situation is discussed.

Somewhat unexpectedly, for the Toda lattice the scattering map can be made
explicit. The formulas simplify by considering the inverse transformation. We
define the map $\Phi:\Gamma_{\\!N}^{\triangleright}\to\Gamma_{\\!N}$, i.e.
$\Phi:(\lambda,\phi)\mapsto(q,p)$ and
$\Gamma_{\\!N}^{\triangleright}=\mathbb{W}_{N}\times\mathbb{R}^{N}$. The map
$\Phi$ is one-to-one, holomorphic, and given by

$q_{j}=\log\big{(}\sigma_{N+1-j}/\sigma_{N-j}\big{)},\qquad
p_{j}=\big{(}\dot{\sigma}_{N+1-j}/\sigma_{N+1-j}\big{)}-\big{(}\dot{\sigma}_{N-j}/\sigma_{N-j}\big{)}.$
(2.66)

Here $\sigma_{0}=1$ and

$\sigma_{k}=\sum_{|I|=k}\exp\Big{(}\sum_{i\in I}\phi_{i}\Big{)}\prod_{i\in
I,j\notin I}|\lambda_{i}-\lambda_{j}|^{-1}$ (2.67)

for $k=1,\ldots,N$, where the first sum is over all subsets
$I\subset\\{1,\ldots,N\\}$ of cardinality $|I|=k$. The symbol
$\dot{\sigma}_{k}$ refers to the Poisson bracket with the transformed
hamiltonian,

$\dot{\sigma}_{k}=\big{\\{}\sigma_{k},\sum_{j=1}^{N}\tfrac{1}{2}\lambda_{j}^{2}\big{\\}}.$
(2.68)

Later on, we will use these expressions to compute the generalized free
energy, see Section 9.2.

Notes and references

Section 2.0

The ground breaking discoveries are Toda [287, 288]. The second edition of the
book Toda’s “Theory of Nonlinear Lattices” [290] is still the most complete
account up to 1989. Faddeev and Takhtajan [108] is a widely used standard
monograph on classical integrable systems. A somewhat more elementary
introduction is Arutyunov [1]. Specifically the Toda lattice is reviewed by
Krüger and Teschl [180]. Closer to hydrodynamics is the study of a particular
shock problem by Venakides et al. [293]. The fifty years anniversary volume
edited by Bazhanov et al. [19] provides a glimpse on research in vastly
diverse directions.

Section 2.1

Based on explicit soliton solutions, Toda conjectured integrability. For the
case of three particles, Ford et al. [117] obtained very supporting numerical
Poincaré plots. The integrals of motion in full generality were obtained by
Hénon [149] by a tricky enumerative argument. Hénon was worried about
locality. Flaschka [113] had the advantage of working at the Courant
Institute, at which Peter Lax [189] introduced his matrix in the context of
the Korteweg-de Vries equation, see Section 10.1. Once the Lax matrix had been
discovered, locally conserved fields are easily constructed. Independently,
the Lax matrix for the Toda lattice has been reported by Manakov [197].
Apparently, at the time and later on, currents were hardly in focus, one
exception being Shastry and Young [256], who study the energy Drude weight in
thermal equilibrium.

Section 2.2

To find out the canonically conjugate angles is a much more technical
enterprise, which was accomplished in a series of papers by Henrici and
Kappeler [150, 151, 152], starting from an early proposal by Flaschka and
McLaughlin [115], see Ferguson et al. [110] for complimentary aspects. The
pitfalls are also discussed in Ruijsenaars [250], Section 5. Lucid proofs of
the stated properties are presented by Hubacher [156]. As a further result,
the additivity (2.64) is deduced when assuming the property
$p_{j}^{+}=p_{N+1-j}^{-}$. In spirit this result ensures the existence of
quasiparticles. However, the abstract proof does not tell us for which
interaction potentials the assumption holds.

For quantum systems, the notion of integrability is controversial, since the
naive transcription of the classical notion would mean that every
eigenprojection of the hamiltonian is conserved, in itself not such a helpful
observation. We refer to Caux and Mossel [55] for an exhaustive discussion.
The link between integrability and local conservation laws has been mostly
pushed by the quantum community, see Grabowski and Mathieu [136, 137] for
early work. For classical systems, in general, this avenue still needs to be
further developed. The mentioned results for the XYZ chain are prototypical
for what one would like to achieve. The nonintegrable case is a result of
Shiraishi [258]. The integrable case, $h=0$, has been studied already by
Grabowski and Mathieu [136, 137] with recent progress Nozawa and Fukai [221].
Considering the XXZ chain with parameters $J_{x}=J_{y}=1$, $J_{z}=\Delta$ and
using only strictly local charges of the spin chain, one computes the spin
Drude weight at zero magnetization, i.e. the persistent spin current, by using
the Mazur formula. By spin inversion symmetry this weight turns out to be
identically $0$. On the other hand for $0\leq\Delta<1$, numerical evidence and
exact steady state results for the boundary driven chain indicate that the
Drude weight does not vanish. The puzzle is resolved by the construction of
quasilocal conserved charges in Mierzejewski et al. [208], see also related
work Ilievski et al. [162, 163], for a broader perspective Doyon [86], and the
recent review Ilievski [159]. The Drude weight is nowhere continuous in its
dependence on $\Delta$.

The relation between integrability and conservation laws has been investigated
also in the context of $1+1$-dimensional quantum field theories. Claimed is
indeed a dichotomy, in the following sense: If beyond the conservation of
energy and momentum there is a single higher order locally conserved charge,
then necessarily the theory is integrable, in the sense of possessing
infinitely many conservation laws, see Coleman and Mandula [59], Iagolnitzer
[157, 158], Parke [229], and Doyon [85].

Section 2.3

In a beautiful piece of analysis Moser [214] proves the scattering shift for
the $N$-particle Toda lattice. An account of his work can be found in Toda
[290]. From a different perspective, a more recent discussion are the notes of
Deift et al. [71] from his course at the Courant Institute in Spring 2019. The
action-angle map $\Phi$ was established by Ruijsenaars [248], where also the
stated properties are proved. Asymptotic momenta equal eigenvalues of the Lax
matrix. In spirit, there is also a corresponding algebraic identity for the
scattering shifts, which has been constructed in analogy to more accessible
models. But to check directly the validity of the Poisson bracket relations in
(2.65) seems to be completely out of reach. Thus as a major difficulty, one
first has to establish agreement between algebraic approach and scattering
map. A more recent point of view is developed by Fehér [109].

The distinction between integrable and nonintegrable scattering has been
studied in detail for one-particle systems. An example is the four hill
potential
$V_{\mathrm{fh}}(q_{1},q_{2})=q_{1}^{2}q_{2}^{2}\exp\big{(}-q_{1}^{2}-q_{2}^{2}\big{)}$.
Chaotic scattering is reviewed by Seoane and Sanjuán [255].

##  3 Static properties

Considering the Toda lattice with periodic boundary conditions and a highly
excited initial condition $q,p$, one would expect that in the long time limit
a statistically stationary state is reached. For a simple fluid, this state
would be thermal equilibrium. But the motion of Toda particles is highly
constrained through the conservation laws. Still, there are lots and lots of
random like collisions. Following Boltzmann, a natural guess for the
statistically stationary state is a generalized microcanonical ensemble,
namely the uniform measure on the $(N-1)$-dimensional torus $\mathbb{T}^{N-1}$
at fixed values of Lax eigenvalues $\lambda_{1},\ldots,\lambda_{N}$ and of
$Q^{[0],N}$, see the discussion in the beginning of Section 2.2. The
corresponding thermodynamics thus depends on $N+1$ extensive parameters. As
for simple fluids, the first step towards hydrodynamic equations is a study of
such generalized thermodynamics.

###  3.1 Generalized Gibbs ensembles

For fixed number of lattice sites, the phase space is $(r,p)\in\Gamma_{\\!N}$
with the a priori weight

$\prod_{j=1}^{N}\mathrm{d}r_{j}\mathrm{d}p_{j}\delta\Big{(}\sum_{j=1}^{N}r_{j}-\ell\Big{)}$
(3.1)

for some $\ell\in\mathbb{R}$, where we included already the microcanonical
constraint resulting from the boundary conditions (2.9). This measure is
invariant under the flow generated by Eq. (2.5). Since particles are
distinguishable, there is no factor of $1/N!$ in front. The remaining
conserved fields are taken into account through the grand canonical type
Boltzmann weight

$\exp\\!\Big{(}-\sum_{n=1}^{N}\mu_{n}Q^{[n],N}\Big{)}.$ (3.2)

Here $\mu_{1},\ldots,\mu_{N}$ are the intensive parameters. Only the low order
ones have a physical interpretation, specifically $\mu_{2}=\tfrac{1}{2}\beta$
with $\beta$ the inverse temperature and $\mu_{1}$ as control parameter for
the average total momentum. As common in Statistical Mechanics we invoke the
equivalence of ensembles to lift the delta constraint by the substitution

$\delta\big{(}Q^{[0],N}-\ell\big{)}\quad\Rightarrow\quad\exp\\!\big{(}-PQ^{[0],N}\big{)}.$
(3.3)

Such kind of equivalence has been extensively studied in rigorous statistical
mechanics and presumably some of the techniques can be used also for the Toda
lattice. Along with other items, we have to leave this problem for future
studies. For a general anharmonic chain the physical pressure, $\mathfrak{p}$,
is defined as the average force between neighboring particles in thermal
equilibrium. Using a simple integration by parts, one obtains the relation
$\mathfrak{p}=\beta^{-1}P$. We still refer to $P$ as pressure, since it is the
thermodynamic dual of the stretch. To have an integrable Boltzmann weight,
$P>0$ is required. In combination the generalized Gibbs ensemble (GGE) is
defined through

$\prod_{j=1}^{N}\mathrm{d}r_{j}\mathrm{d}p_{j}\exp\\!\Big{(}-PQ^{[0],N}-\sum_{n=1}^{N}\mu_{n}Q^{[n],N}\Big{)},$
(3.4)

which still has to be normalized. The GGE is invariant under the Toda
dynamical flow.

To study properties of GGE the first natural step is to transform (3.4) to
Flaschka variables. For conciseness we introduce

$V(w)=\sum_{n=1}^{\infty}\mu_{n}w^{n}.$ (3.5)

The chemical potentials, $\mu_{n}$, are assumed to be independent of $N$. Then
the transformed density reads

$\exp\\!\big{(}-\mathrm{tr}[V(L_{N})]\big{)}\prod_{j=1}^{N}\mathrm{d}p_{j}\prod_{j=1}^{N}\mathrm{d}a_{j}\frac{2}{a_{j}}(a_{j})^{2P},$
(3.6)

which is defined on the phase space of the Flaschka variables, i.e. on
$\Gamma_{\\!N}^{\circ}$.

$\blackdiamond\blackdiamond$ Confining potential. A many-particle system is
defined by a particular interaction potential, which throughout will be
denoted by $V_{\bullet}$. For example $V_{\mathrm{to}}$ is the interaction
potential of the Toda lattice and $V_{\mathrm{li}}$ the $\delta$-potential of
the Lieb-Liniger model. The potential (3.5) could be called a generalized
chemical potential, a terminology which however would entirely miss the
central role of $V$. We wrote $V$ in terms of a power series. But there is no
compelling reason to do so. $V$ is simply a rather generic function on
$\mathbb{R}$. In the context of GGE, its purpose is to properly confine the
eigenvalues of $L_{N}$. We thus assume that $V$ is continuous and bounded
linearly from below as $V(w)\geq c_{0}+c_{1}|w|$ with $c_{1}>0$. Based on such
reasoning $V$ is called confining potential. It carries no relation to the
interaction potential. Rather $V$ should be viewed as a thermodynamic
variable, which so to speak labels the GGEs. The quadratic confining
potential, $V(w)=\tfrac{1}{2}\beta w^{2}$, corresponds to thermal equilibrium.
$\blackdiamond\blackdiamond$

$\blackdiamond\blackdiamond$ Infinite volume limit, exponential mixing,
equivalence of ensembles. As for other Gibbs measures, one might want to know
about the existence of the infinite volume limit for the normalized sequence
of measures in (3.6), the limit measure being independent of boundary
conditions, and a bound on the decay of correlations. If the confining
potential is a finite polynomial with even strictly positive leading term, say
$\kappa$, then the confining potential is bounded from below and such
properties can be answered by using transfer matrix techniques. In the
language of statistical mechanics, $Q^{[n]}$ has a range of size $n$ and hence
$\mathrm{tr}[V(L_{N})]$ has range $\kappa$. One cuts $[1,\ldots,N]$ in blocks
of size $\kappa$. The density in (3.6) can then be written as an
$(N/\kappa)$-fold power of the transfer matrix. This is just like the familiar
case of the one-dimensional Ising model, in which case the transfer matrix is
a $2\times 2$ matrix. For the Toda lattice the transfer matrix is given by an
integral kernel with arguments in $\Gamma_{\kappa}^{\circ}$. By the Perron-
Frobenius theorem, the transfer matrix has a unique maximal eigenvalue, which
is separated by a gap from the rest of the spectrum. With this input, one
concludes that there is a unique limit measure. In one dimension, phase
transitions would occur only if the interaction potential has a decay slower
than (range)-2, much slower than the case under consideration here. The
spectral gap also ensures exponential decay of correlations. If
$\kappa=\infty$ such methods fail completely. Other techniques will have to be
developed, see Notes.

A presumably more delicate issue is the equivalence of ensembles. This refers
to replacing the sharp constraint $\delta\big{(}H_{N}-\mathsfit{e}N\big{)}$ by
$\exp\big{(}-\beta H_{N}\big{)}$. On the thermodynamic level, in the limit
$N\to\infty$, this amounts to a Legendre transform between $\mathsfit{e}$ and
$\beta$. As a stronger property the statistics of local observables is the
same provided $\mathsfit{e}$ and $\beta$ are related according to
thermodynamics. For the Toda lattice we switched from $\nu N$ to $P$, a step
which is expected to be accessible to current techniques. More delicate is the
closed Toda chain as discussed in Section 2.2. In principle one should fix all
action variables. The integral over the tori is trivial and one is left with
$\sum_{j=1}^{N}V(\lambda_{j})$. Abstractly this sum depends only on $I$, but a
more concrete characterization does not seem to be available. Thus a fully
microcanonical approach will not be pursued any further.
$\blackdiamond\blackdiamond$

###  3.2 Lax matrix filter and local GGEs

This section is somewhat premature, since so far the only model discussed is
the Toda lattice. Still, before entering in computational details, we should
explain the physics underlying the notion of local GGE’s as the central theme
of hydrodynamic scales. Some of the material will be explained in greater
detail in the following.

Let us start from the example of a classical ideal gas in one dimension, which
consists of many particles moving along straight lines as
$q_{j}(t)=q_{j}+p_{j}t$ for initial conditions $(q,p)$. The momenta are
conserved and hence independent of time. For general random initial data, one
introduces the one-particle distribution function

$\rho_{1}(x,t;w)=\big{\langle}\sum_{j}\delta(x-q_{j}(t))\delta(w-p_{j})\big{\rangle}_{\mathrm{ini}}$
(3.7)

average over the initial probability density function. For an ideal gas,
$\rho_{1}(x,t;w)=\rho_{1}(x-wt,0;w)$.

For the hydrodynamic scale of the ideal gas we will use the notation
$\rho_{\mathrm{Q}}(x,t;w)$ which superficially looks rather similar to
$\rho_{1}$. Here $(x,t)$ refers to a macroscopic spacetime point, which is the
center of a microscopic cell of size $\ell$ containing $N(\ell)$ particles.
$N(\ell)/\ell$ is of order 1 and each cell contains a large number of
particles. The index $\mathrm{Q}$ should be a reminder that
$\rho_{\mathrm{Q}}(x,t;w)$ carries the information on the average conserved
fields in the considered cell. $\rho_{\mathrm{Q}}(x,t;w)$ has a double
meaning. We define the empirical velocity distribution function through

$\frac{1}{\ell}\sum_{j=1}^{N(\ell)}\delta(w-p_{j})\simeq\rho_{\mathrm{Q}}(x,t;w),$
(3.8)

where the sum is over all particles in the cell centered at $(x,t)$. The left
hand side is random, but self-averaging in the sense that fluctuations vanish
for large $N(\ell)$. In the limit, $\rho_{\mathrm{Q}}(x,t;w)$ is nonrandom and
its integral with respect to $w$ is the macroscopic density $\bar{\rho}(x,t)$.
In addition $\rho_{\mathrm{Q}}(x,t;w)$ uniquely characterizes the local GGE
governing the statistics of particles close to $(x,t)$. For an ideal gas
velocities are independent with common probability density function
$\rho_{\mathrm{Q}}(x,t;w)/\bar{\rho}(x,t)$. The positional distribution is
independent of velocities, more precisely, a Poisson point process with
density $\bar{\rho}(x,t)$, which means that inter-particle distances are
independent and exponentially distributed. Note that in the limit of infinite
scale separation, the local GGE lives on the entire line and is translation
invariant. The construction (3.8) holds also for a hard rod fluid. Only now
the local GGE has a positional distribution of particles satisfying the hard
core constraint.

The at first sight truly surprising claim is that for classical integrable
many-body systems the definition in (3.8) can still be used provided the
momenta $p_{j}$ are replaced by the eigenvalues $\lambda_{j}$ of the local Lax
matrix $L(x,t)$. This is what we mean by a Lax matrix filter. Local positions
and momenta are noisy. But by inserting these data in the local Lax matrix and
determining the local density of states (DOS), magically, one filters the
slowly varying degrees of freedom. Compared with the ideal gas, rather than
sampling local velocities one has to sample the eigenvalues of the local Lax
matrix. More precisely, the Lax matrix is random under the local GGE and has
an empirical DOS according to

$\frac{1}{\ell}\sum_{j=1}^{N(\ell)}\delta(w-\lambda_{j}).$ (3.9)

The sum is over eigenvalues of the Lax matrix constructed from the particular
fluid cell under consideration. As valid in great generality, the DOS is self-
averaging and thus has a deterministic limit, as before denoted by
$\rho_{\mathrm{Q}}(x,t;w)$, but now referring to the Lax DOS. The function
$\rho_{\mathrm{Q}}(x,t;w)$ uniquely determines the underlying local GGE in the
cell $(x,t)$. The sought for hydrodynamic equations are an evolution equation
for $\rho_{\mathrm{Q}}(x,t;w)$, which is of conservation type

$\partial_{t}\rho_{\mathrm{Q}}(x,t;w)+\partial_{x}\rho_{\mathrm{J}}(x,t;w)=0.$
(3.10)

A major task will be to figure out the $\rho_{\mathrm{J}}(x,t;w)$ as a
functional of $\rho_{\mathrm{Q}}(x,t;w)$ at fixed fluid cell $(x,t)$. Concrete
examples will come. For some models the Lax matrix is unitary and eigenvalues
lie on the unit circle. Only through a more detailed analysis such properties
can figured out.

###  3.3 Generalized free energy

The Toda partition function is defined by

$Z_{\mathrm{to},N}(P,V)=\int_{\Gamma_{\\!N}^{\circ}}\prod_{j=1}^{N}\mathrm{d}p_{j}\prod_{j=1}^{N}\mathrm{d}a_{j}\frac{2}{a_{j}}(a_{j})^{2P}\exp\\!\big{(}-\mathrm{tr}[V(L_{N})]\big{)}.$
(3.11)

Accordingly, normalizing the expression in (3.6), one arrives at the
probability measure

$\mu^{\mathrm{GGE}}_{P,V,N}(\mathrm{d}^{N}a\,\mathrm{d}^{N}p)=\frac{1}{Z_{\mathrm{to},N}(P,V)}\prod_{j=1}^{N}\mathrm{d}p_{j}\prod_{j=1}^{N}\mathrm{d}a_{j}\frac{2}{a_{j}}(a_{j})^{2P}\exp\\!\big{(}-\mathrm{tr}[V(L_{N})]\big{)}.$
(3.12)

This measure is time-stationary under the dynamics (2.5), since the a priori
measure and the eigenvalues of $L_{N}$ do not change in time. The expectations
of $\mu^{\mathrm{GGE}}_{P,V,N}$ will be denoted by
$\langle\cdot\rangle_{P,V,N}$. As central thermodynamic object, the free
energy per lattice site is defined by

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

For hydrodynamics, a crucial input is the GGE average of the conserved fields.
They can be computed as derivative with respect to $P$ and as variational
derivative with respect to $V$ of the free energy. In terms of the eigenvalues
of the Lax matrix, the averages can be written as

$\frac{1}{N}\langle
Q^{[n],N}\rangle_{P,V,N}=\frac{1}{N}\big{\langle}\sum_{j=1}^{N}(\lambda_{j})^{n}\big{\rangle}_{P,V,N}=\int_{\mathbb{R}}\mathrm{d}ww^{n}\langle\rho_{\mathrm{Q},N}(w)\rangle_{P,V,N},$
(3.14)

where

$\rho_{\mathrm{Q},N}(w)=\frac{1}{N}\sum_{j=1}^{N}\delta(w-\lambda_{j})$ (3.15)

is the empirical density of states of the Lax matrix $L_{N}$. Here empirical
refers to the fact that the DOS is defined for every collection of eigenvalues
$\\{\lambda_{1},\ldots,\lambda_{N}\\}$. Thus $\rho_{\mathrm{Q},N}(w)$ is a
random function under $\mu^{\mathrm{GGE}}_{P,V,N}$. More properly,
$\rho_{\mathrm{Q},N}$ is a random probability measure supported on $N$ points
each with weight $1/N$.

This observation suggests a novel perspective. The Lax matrix becomes a random
matrix under $\mu^{\mathrm{GGE}}_{P,V,N}$. Thermal equilibrium is particularly
simple. Since $\mathrm{tr}[(L_{N})^{2}]=\sum_{j=1}^{N}(p_{j}^{2}+2a_{j}^{2})$,
the diagonal and off-diagonal matrix elements of $L_{N}$ are families of
independent identically distributed (i.i.d.) random variables. For all other
GGEs the matrix elements are correlated. The DOS encodes the complete
statistical information on the conserved fields. In fact
$\rho_{\mathrm{Q},N}(w)$ is self-averaging with fluctuations of order
$1/\sqrt{N}$ and the limit

$\lim_{N\to\infty}\rho_{\mathrm{Q},N}(w)=\rho_{\mathrm{Q}}(w)$ (3.16)

exists with some nonrandom limiting density $\rho_{\mathrm{Q}}$, which
generically is a smooth function.

$\blackdiamond\blackdiamond$ The Dumitriu-Edelman identity. In 2002 Dumitriu
and Edelman studied the $\beta$-ensembles of random matrix theory. To ease a
comparison, I describe their result in the original notation, in which some
symbols will reappear with a different meaning. The proper translation will be
obvious, however. Their starting point is a $n\times n$, symmetric,
tridiagonal matrix $T$ with real matrix elements $T_{j,j}=a_{j}$,
$T_{j,j+1}=T_{j+1,j}=b_{j}>0$, and zero otherwise, compare with $L^{\diamond}$
from (2.61), in particularly,
$(a,b)\in\Gamma^{\diamond}_{\\!n}=\mathbb{R}^{n}\times\mathbb{R}_{+}^{n-1}$.
The eigenvalues of $T$ are ordered as $\lambda_{1}<\ldots.<\lambda_{n}$ and
the first component of an eigenvector is denoted by $\psi_{j}(1)=q_{j}>0$ with
$|q|=1$ imposed. We set $\mathrm{d}^{n}a=\prod_{j=1}^{n}\mathrm{d}a_{j}$,
$\mathrm{d}^{n}\lambda=\prod_{j=1}^{n}\mathrm{d}\lambda_{j}$ and
$\mathrm{d}^{n}q$ the surface element of the unit $n$-sphere. As a general
fact of Jacobi matrices with positive off-diagonal matrix elements, there is a
one-to-one and onto map $\Psi:(a,b)\mapsto(\lambda,q)$. Dumitriu and Edelman
managed to obtain some information on the Jacobian of $\Psi$. Their identity
has a free parameter $\beta>0$ and reads

$\mathrm{e}^{-\mathrm{tr}[V(T)]}\,\mathrm{d}^{n}a\prod_{j=1}^{n-1}\frac{2}{b_{j}}(b_{j})^{\beta
j}\mathrm{d}b_{j}=\exp\Big{(}-\sum_{j=1}^{n}V(\lambda_{j})\Big{)}\big{(}n!\tilde{\zeta}_{n}(\beta)\Delta(\lambda)^{\beta}\mathrm{d}^{n}\lambda\big{)}\Big{(}c_{q}^{\beta}\prod_{j=1}^{n}(q_{j})^{\beta-1}\mathrm{d}^{n}q\Big{)}.$
(3.17)

Here

$c_{q}^{\beta}=2^{n-1}\Gamma(\tfrac{1}{2}\beta
n)\Gamma(\tfrac{1}{2}\beta)^{-n}$ (3.18)

normalizes the third factor on the left to 1. In the second factor,
$\Delta(\lambda)$ is the Vandermonde determinant

$\Delta(\lambda)=\prod_{1\leq i<j\leq n}(\lambda_{j}-\lambda_{i}),$ (3.19)

while

$\tilde{\zeta}_{n}(\beta)=\Gamma(\tfrac{1}{2}\beta
n)^{-1}\Gamma(1+\tfrac{1}{2}\beta)^{n}\prod_{j=1}^{n}\frac{\Gamma(\tfrac{1}{2}\beta
j)}{\Gamma(1+\tfrac{1}{2}\beta j)}$ (3.20)

is the proportionality constant. $\blackdiamond\blackdiamond$

In the Dumitriu-Edelman identity (3.17) we substitute $n\leadsto N$,
$a_{j}\leadsto p_{j}$, and $b_{j}\leadsto a_{j}$. Since the scale parameter
mentioned below Eq. (2.4) has been set to $\tau=1$, one notes that
$T=L_{N}^{\diamond}$. Therefore integrating both side of (3.17) over the
entire phase space, first making the choice

$\beta=\frac{2P}{N},$ (3.21)

yields

$\displaystyle
Z_{\mathrm{de},N}(P,V)=\int_{\Gamma^{\diamond}_{\\!N}}\exp\\!\big{(}-\mathrm{tr}[V(L^{\diamond}_{N})]\big{)}\prod_{j=1}^{N}\mathrm{d}p_{j}\prod_{j=1}^{N-1}\mathrm{d}a_{j}\frac{2}{a_{j}}(a_{j})^{2(j/N)P}$
$\displaystyle\hskip
18.0pt=\zeta_{N}(P)\int_{\mathbb{R}^{N}}\mathrm{d}^{N}\lambda\exp\\!\Big{(}-\sum_{j=1}^{N}V(\lambda_{j})+\frac{P}{N}\sum_{i,j=1,i\neq
j}^{N}\log|\lambda_{i}-\lambda_{j}|\Big{)}$ (3.22)

with prefactor

$\tilde{\zeta}_{N}(2P/N)=\zeta_{N}(P)=\Gamma(P)^{-1}\Gamma(1+\tfrac{P}{N})^{N}\prod_{j=1}^{N}\frac{\Gamma(\tfrac{j}{N})}{\Gamma(1+\tfrac{j}{N})}.$
(3.23)

The integration over eigenvalues is not ordered, which takes care of the
factor $N!$. Also, the integration over the $N$-sphere is normalized to $1$.

The term on the right side of (3.3) requires more explanations. For the
normalization one obtains

$\lim_{N\to\infty}-\frac{1}{N}\log\zeta_{N}(P)=\log P-1.$ (3.24)

Otherwise the partition function is the one of the repulsive one-dimensional
log gas. $V$ turns out to be the confining potential of the log gas, which is
the real reason for our original choice of name. In the standard log gas the
interaction strength is of order $1$, which implies that the free energy is
dominated by the energy term. But in our case the interaction strength is
$1/N$, which is the standard mean-field scaling. Such a problem can be handled
through the study of a free energy functional. To distinguish from (3.15), one
introduces

$\varrho_{N}(w)=\frac{1}{N}\sum_{j=1}^{N}\delta(w-\lambda_{j})$ (3.25)

with $\lambda$ now referring to (3.3). Except for the diagonal contribution,
the integrand of (3.3) can be written as

$\exp\Big{[}-N\Big{(}\int_{\mathbb{R}}\mathrm{d}w\varrho_{N}(w)V(w)-P\int_{\mathbb{R}}\mathrm{d}w\int_{\mathbb{R}}\mathrm{d}w^{\prime}\varrho_{N}(w)\varrho_{N}(w^{\prime})\log|w-w^{\prime}|\Big{)}\Big{]}.$
(3.26)

For given $\varrho_{N}$ the corresponding volume element of
$\mathrm{d}^{N}\lambda$ is approximately

$\exp\Big{[}-N\int_{\mathbb{R}}\mathrm{d}w\varrho_{N}(w)\log\varrho_{N}(w)\Big{]}.$
(3.27)

The large $N$ limit of the partition function is then determined by the mean-
field free energy functional

$\mathcal{F}_{\mathrm{de}}^{\circ}(\varrho)=\int_{\mathbb{R}}\mathrm{d}w\varrho(w)\Big{(}V(w)+\log\varrho(w)+\log
P-P\int_{\mathbb{R}}\mathrm{d}w^{\prime}\varrho(w^{\prime})\log|w-w^{\prime}|\Big{)}.$
(3.28)

The actual free energy is obtained by minimizing over all $\varrho$ with
$\varrho\geq 0$ and $\int_{\mathbb{R}}\mathrm{d}w\varrho(w)=1$. As will be
discussed below, there is a unique minimizer, $\varrho^{\star}$, and thus

$\lim_{N\to\infty}-\frac{1}{N}\log
Z_{\mathrm{de},N}(P,V)=F_{\mathrm{de}}(P,V)=\mathcal{F}_{\mathrm{de}}^{\circ}(\varrho^{\star})-1.$
(3.29)

To obtain the Toda free energy, we note that the Dumitriu-Edelman partition
function has a pressure changing linearly with slope $2P/N$. This is not
exactly what is required, since the pressure is constant for the Toda chain.
But in a large segment $\ell$, still with size $\ell\ll N$, the pressure is
constant to a very good approximation. Since GGEs have good spatial mixing
properties, local free energies merely add up and for the Dumitriu-Edelman
free energy one concludes

$F_{\mathrm{de}}(P,V)=\int_{0}^{1}\mathrm{d}uF_{\mathrm{to}}(uP,V)$ (3.30)

and hence

$F_{\mathrm{to}}(P,V)=\partial_{P}(PF_{\mathrm{de}}(P,V)).$ (3.31)

While the just presented derivation of $F_{\mathrm{de}}(P,V)$ is fairly
standard, one may wonder about the missing steps. A poor man’s version will be
explained in Section 3.5. Besides, the topic has been studied extensively with
methods covering the case of interest. Less standard is the linear pressure
ramp leading to the identity (3.30).

In the free energy functional (3.28) the quadratic term has the kernel
$\log|w-w^{\prime}|$ that resulted from the Dumitriu-Edelman change of
coordinates. The same expression came up already in Section 2.3 in the context
of scattering theory, which might be considered as purely accidental. In fact,
we touched upon a generic feature of generalized free energies for integrable
models. For the Toda lattice the connection to scattering theory will be
elucidated in Section 9.2.

It turns out to be more convenient to absorb $P$ into $\varrho$ by setting
$\rho=P\varrho$. Then
$P\mathcal{F}_{\mathrm{de}}^{\circ}(P^{-1}\rho)=\mathcal{F}_{\mathrm{de}}(\rho)$
with the transformed free energy functional

$\mathcal{F}_{\mathrm{de}}(\rho)=\int_{\mathbb{R}}\mathrm{d}w\rho(w)\Big{(}V(w)+\log\rho(w)-\int_{\mathbb{R}}\mathrm{d}w^{\prime}\rho(w^{\prime})\log|w-w^{\prime}|\Big{)}.$
(3.32)

$\mathcal{F}_{\mathrm{de}}$ has to be minimized under the constraint

$\rho(w)\geq 0,\quad\int_{\mathbb{R}}\mathrm{d}w\rho(w)=P$ (3.33)

with minimizer denoted by $\rho^{\star}$. Then

$F_{\mathrm{to}}(P,V)=\partial_{P}\mathcal{F}_{\mathrm{de}}(\rho^{\star})-1.$
(3.34)

The constraint (3.33) is removed by introducing the Lagrange multiplier $\mu$
as

$\mathcal{F}_{\mathrm{de}}^{\bullet}(\rho)=\mathcal{F}_{\mathrm{de}}(\rho)-\mu\int_{\mathbb{R}}\mathrm{d}w\rho(w).$
(3.35)

A minimizer of $\mathcal{F}_{\mathrm{de}}^{\bullet}(\rho)$ is denoted by
$\rho_{\mathsf{n}}$ and by $\rho_{\mathsf{n},\mu}$ when keeping track of the
$\mu$-dependence. It is determined as solution of the Euler-Lagrange equation

$V(w)+\log\rho_{\mathsf{n},\mu}(w)-\mu-2\int_{\mathbb{R}}\mathrm{d}w^{\prime}\rho_{\mathsf{n},\mu}(w^{\prime})\log|w-w^{\prime}|=0.$
(3.36)

The Lagrange parameter $\mu$ has to be adjusted such that

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

In fact, it will be more convenient to work directly with $\rho_{\mathsf{n}}$.
Note that the Lagrange parameter $\mu$ amounts to shifting the confining
potential as $V-\mu$. Thus $\mu$ could be incorporated in the definition (3.5)
of $V$ as $\mu_{0}=-\mu$. The minus sign is a standard convention for the
chemical potential dual to the particle number.

To obtain the Toda free energy, we differentiate as

$\partial_{P}\mathcal{F}_{\mathrm{de}}(\rho^{\star})=\int_{\mathbb{R}}\mathrm{d}w\partial_{P}\rho^{\star}(w)\Big{(}V(w)+\log\rho^{\star}(w)-2\int_{\mathbb{R}}\mathrm{d}w^{\prime}\rho^{\star}(w^{\prime})\log|w-w^{\prime}|\Big{)}+1.$
(3.38)

Integrating (3.36) against $\partial_{P}\rho^{\star}$ one arrives at

$\partial_{P}\mathcal{F}_{\mathrm{de}}(\rho_{\mathsf{n},\mu})=\mu+1$ (3.39)

and thus

$F_{\mathrm{to}}(P,V)=\mu(P,V).$ (3.40)

Sharing with other integrable models, the Toda lattice has the property that
its free energy is determined by a variational problem for densities over
$\mathbb{R}$, in our case normalized to $P$.

$\blackdiamond\blackdiamond$ Densities. Our generic symbol for a density is
$\rho$, resp. $\varrho$. We use $\varrho$ for normalized densities,
$\int\mathrm{d}w\varrho(w)=1$, while $\rho$ has a context dependent
normalization. For Toda and other integrable systems the densities are defined
on $\mathbb{R}$. Later on we will also encounter models where the densities
live on the unit circle. Density does not refer to position space. The
appropriate picture is a density in distorted momentum space. Since there are
several densities, one has to distinguish them by a label which in our
notation appears through a lower index, as in $\rho_{\mathrm{Q}}$ and
$\rho_{\mathrm{J}}$. We introduced already the density $\rho_{\mathsf{n}}$,
also called number density. But it will turn out to be convenient to introduce
the further densities $\rho_{\mathsf{p}}$, $\rho_{\mathsf{h}}$,
$\rho_{\mathsf{s}}$, called particle, hole, space density. The label is now in
serif to avoid confusion with arguments. All densities are functionals of $V$,
a dependence which is mostly suppressed in our notation. For example
$\rho_{\mathsf{n}}$ functionally depends on $V-\mu$, hence it is a function of
$\mu$ and written as $\rho_{\mathsf{n},\mu}$. $\rho_{\mathsf{n}}$ can be
regarded also as a function of $P$, since $\mu=\mu(P)$ for fixed $V$.

$\blackdiamond\blackdiamond$

###  3.4 Lax density of states, TBA equation

To obtain the hydrodynamic equations, GGE averages of the conserved fields are
required. As common practice in statistical mechanics, they are defined in the
infinite volume limit. For this purpose we adopt the volume $[-N,\ldots,N]$
with periodic boundary conditions and adjust our notation by using the label
$2N+1$ instead of $N$. Then, using translation invariance, in the infinite
volume limit

$\lim_{N\to\infty}\frac{1}{2N+1}\langle
Q^{[n],2N+1}\rangle_{P,V,2N+1}=\lim_{N\to\infty}\langle
Q^{[n],2N+1}_{0}\rangle_{P,V,2N+1}=\langle Q^{[n]}_{0}\rangle_{P,V}.$ (3.41)

On the right hand side $Q^{[n]}_{0}$ is a local function and
$\langle\cdot\rangle_{P,V}$ refers to the average with respect to the infinite
volume GGE at parameters $P,V$. Since in our context boundary terms should be
negligible, we assume this measure to be well-defined and independent of
boundary conditions. More pragmatically, only for particular observables, as
$Q^{[n]}_{0}$ and $Q^{[n]}_{0}Q^{[m]}_{j}$, the infinite volume average has to
exist.

To determine $\langle Q^{[n]}_{0}\rangle_{P,V}$, one can start from the
microscopic definition above and use that $Q^{[n],N}$ depends only on the
eigenvalues of the Lax matrix. The other method, employed here, is to simply
differentiate the infinite volume free energy. We start with $n=0$ and note
that the average stretch

$\nu=\langle
Q^{[0]}_{0}\rangle_{P,V}=\partial_{P}F_{\mathrm{to}}(P,V)=\partial_{P}\mu(P,V)=\Big{(}\int_{\mathbb{R}}\mathrm{d}w\partial_{\mu}\rho_{\mathsf{n}}(w)\Big{)}^{-1},$
(3.42)

where the last equality results from differentiating Eq. (3.37) as

$1=(\int\mathrm{d}w\partial_{\mu}\rho_{\mathsf{n},\mu}(w))\partial_{P}\mu(P).$
(3.43)

For $n\geq 1$ we perturb $V$ as $V_{\kappa}(w)=V(w)+\kappa w^{n}$ and
differentiate the free energy at $\kappa=0$. Then

$\langle
Q^{[n]}_{0}\rangle_{P,V}=\partial_{\kappa}F_{\mathrm{to}}(P,V_{\kappa})\big{|}_{\kappa=0}=\partial_{P}\partial_{\kappa}\mathcal{F}_{\mathrm{de}}(\rho^{\star}(P,V_{\kappa}))\big{|}_{\kappa=0}$
(3.44)

and, first introducing the linearization of $\rho^{\star}$ as

$\rho{{}^{\star}}^{\prime}=\partial_{\kappa}\rho^{\star}(P,V_{\kappa})\big{|}_{\kappa=0},$
(3.45)

one obtains

$\displaystyle\partial_{\kappa}\mathcal{F}_{\mathrm{de}}(\rho^{\star}(P,V_{\kappa}))\big{|}_{\kappa=0}=\int_{\mathbb{R}}\mathrm{d}w\rho^{\star}(w,P,V)w^{n}+\int_{\mathbb{R}}\mathrm{d}w\rho{{}^{\star}}^{\prime}(w)\Big{(}V(w)+\log\rho^{\star}(w,P,V)$
$\displaystyle\hskip
80.0pt-2\int_{\mathbb{R}}\mathrm{d}w^{\prime}\rho^{*}(w^{\prime},P,V)\log|w-w^{\prime}|\Big{)}$
(3.46)

using that $\int\mathrm{d}w\rho{{}^{\star}}^{\prime}(w)=0$. Integrating the
Euler-Lagrange equation (3.36) at $\mu=\mu(P)$ against
$\rho{{}^{\star}}^{\prime}$, the terms on the right side of (3.45) vanish and

$\langle
Q^{[n]}_{0}\rangle_{P,V}=\int_{\mathbb{R}}\mathrm{d}w\partial_{P}\rho^{\star}(w,P,V)w^{n}.$
(3.47)

Thus the Lax DOS is given by

$\rho_{\mathrm{Q}}(w)=\partial_{P}\rho^{\star}(w).$ (3.48)

Naively one might have guessed that the Lax DOS equals $\varrho^{\star}$. But
the slow linear variation of the pressure in the Dumitriu-Edelman identity
amounts to a slightly deviating result.

In the literature the Euler-Lagrange equation (3.36) is written differently by
formally defining a Boltzmann weight as

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

with quasi-energy $\varepsilon(w)$. Then

$\varepsilon(w)=V(w)-\mu-2\int_{\mathbb{R}}\mathrm{d}w^{\prime}\log|w-w^{\prime}|\mathrm{e}^{-\varepsilon(w^{\prime})},$
(3.50)

In addition one also introduces the particle density through

$\rho_{\mathsf{p},\mu}(w)=\partial_{\mu}\rho_{\mathsf{n},\mu}(w),$ (3.51)

at the moment just a convenient terminology.

The structure uncovered is familiar from the Yang-Yang thermodynamics of the
Lieb-Liniger $\delta$-Bose gas, which is an integrable quantum many-body
system and solved by Bethe ansatz. For quantum integrable systems the analogue
of (3.50) is called TBA (thermodynamic Bethe ansatz) equation. We will call
(3.50) classical TBA equation or simply TBA, despite the fact that no Bethe
ansatz had been used in its derivation. Some patience is required to fully
appreciate this analogy. Further evidence will be accumulated from an
alternative route based on scattering coordinates, see Section 9.2 for the
Toda lattice and Section 11.3 for the Calogero fluid. The quantum side of the
analogy will be covered in Chapter 13 for the $\delta$-Bose gas and in Chapter
14 for the quantum Toda lattice.

Later on we will use some identities based on TBA. We collect them here,
together with introducing standard notations. The Hilbert space of square
integrable functions on the real line is denoted by
$L^{2}(\mathbb{R},\mathrm{d}w)$ with scalar product

$\langle f,g\rangle=\int_{\mathbb{R}}\mathrm{d}wf(w)^{*}g(w).$ (3.52)

We will work mostly with real functions and then the complex conjugation in
(3.52) can be omitted. There will be many integrals over $\mathbb{R}$ and a
convenient shorthand is simply

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

To distinguish, an average over some probability measure is denoted by
$\langle\cdot\rangle_{P,V,N}$, carrying suitable subscripts. Starting from the
$Q^{[n]}$’s, so far a discrete basis has been used. Obviously any linear
combination of conserved fields is still conserved and, as in other linear
problems, the choice of basis is an important consideration. From the
viewpoint of Lax DOS, the label $n$ corresponds to the monomial $w^{n}$, which
will continued to be used and is denoted by

$\varsigma_{n}(w)=w^{n},$ (3.54)

including $\varsigma_{0}(w)=1$. More generally, the set of basis functions
will depend on the particular integrable model under consideration.

Let us define the integral operator

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

with $w\in\mathbb{R}$. Then the TBA equation can be rewritten as

$\varepsilon(w)=V(w)-\mu-(T\mathrm{e}^{-\varepsilon})(w).$ (3.56)

In addition one introduces the dressing of a real-valued function $f$ through

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

where $\rho_{\mathsf{n}}$ is regarded as multiplication operator, i.e.
$(\rho_{\mathsf{n}}f)(w)=\rho_{\mathsf{n}}(w)f(w)$. With our improved
notation, the Lax DOS (3.48) can be written as

$\rho_{\mathrm{Q}}=\partial_{P}\rho_{\mathsf{n}}=(\partial_{P}\mu)\partial_{\mu}\rho_{\mathsf{n}}=\nu\rho_{\mathsf{p}},$
(3.58)

where (3.51) has been used. Since $\rho_{\mathrm{Q}}$ is normalized,

$\nu\langle\rho_{\mathsf{p}}\rangle=1.$ (3.59)

Physically the central objects of the theory are $\nu$ and
$\nu\rho_{\mathsf{p}}$, since they encode the GGE average of the conserved
fields. Differentiating TBA with respect to $\mu$ we conclude

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

which expresses $\rho_{\mathsf{p}}$ as a functional of $\rho_{\mathsf{n}}$.
This relation can be inverted to yield

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

For later purposes, we also state the definition

$q_{n}=\langle
Q^{[n]}_{0}\rangle_{P,V}=\nu\langle\rho_{\mathsf{p}}\varsigma_{n}\rangle,$
(3.62)

$n\geq 1$, not to be confused with a particle position.

Identities as (3.56) to (3.61) will reappear in other, either classical or
quantum, integrable systems. Due to their wide use, these identities are
referred to as TBA formalism which in a specific way reflects the underlying
free energy functional.

$\blackdiamond\blackdiamond$ Uniqueness of solutions of the TBA equation. For
the Toda lattice at given $\mu$ the TBA equation has two solutions. At first
glance this looks surprising. In fact, in a standard numerical solution scheme
one follows a particular branch, say starting from small $P$, and encounters
an end-point at which instabilities arise. So, some explanations are in
demand.

Firstly $\mu(P)$ is concave and its derivative, $\nu(P)$, is strictly
decreasing. For example, in the case of thermal equilibrium
$\mu(P)=\log\sqrt{\beta/2\pi}+P\log\beta-\log\Gamma(P)$, which has a single
maximum at $P=P_{\mathrm{c}}$. The physics is rather obvious. At very small
$P$ the average stretch is huge and diverges as $P\to 0$. By increasing
pressure the stretch is decreased. Since there is no hard core, increasing $P$
even further the stretch becomes negative. In physical space, for small $P$,
up to small random errors, the labelling of particles is increasing. But at
large $P$ the labelling is reversed. At $P_{\mathrm{c}}$ the stretch vanishes
and the typical distance between particles with adjacent index is of order
$1/\sqrt{N}$. The function inverse to $\mu(P)$ has two branches, meaning that
for given $\mu$ there are two values of $P$.

Now considering the densities, by construction $\rho_{\mathsf{n}}\geq 0$,
$\nu\rho_{\mathsf{p}}\geq 0$, $\nu\langle\rho_{\mathsf{p}}\rangle=1$.
$\rho_{\mathsf{n},\mu(P)}$ is pointwise increasing in $P$ and varies smoothly
through $P_{\mathrm{c}}$, so does $\nu\rho_{\mathsf{p},\mu(P)}$. On the other
hand, $\rho_{\mathsf{p},\mu(P)}(w)$ diverges to $+\infty$ as $P$ approaches
$P_{\mathrm{c}}$ from the left, globally flips to $-\infty$ at
$P_{\mathrm{c}}$, and then flattens out as $P\to\infty$.
$\blackdiamond\blackdiamond$

####  3.4.1 Thermal equilibrium

. Thermal equilibrium corresponds to the quadratic confining potential
$V(w)=\tfrac{1}{2}\beta w^{2}$ with $\beta$ the inverse temperature. Only for
this particular case the diagonal entries of the Lax matrix,
$\\{p_{j},j\in\mathbb{Z}\\}$, are independent with $p_{j}$ a Gaussian random
variable of mean zero and variance $\beta^{-1}$. Hence
$\langle(p_{j})^{n}\rangle_{P,\beta}=0$ for odd $n$ and
$\langle(p_{j})^{n}\rangle_{P,\beta}=(n-1)!!(\beta)^{-n/2}$ for even $n$. The
off-diagonal entries, $\\{a_{j},j\in\mathbb{Z}\\}$, are also independent with
$a_{j}$ a $\chi$ distributed random variable with parameter $2P$. In
particular for the even moments
$\langle(a_{j})^{2n}\rangle_{P,\beta}=P(P+1)\ldots(P+n-1)$, $n=1,2,\ldots\,$.
Due to independence, the free energy of the chain is easily computed with the
result

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

However to figure out the entire DOS requires the TBA machinery.

We start from the Euler-Lagrange equation for the free energy (3.28), set
$V(w)=\tfrac{1}{2}\beta w^{2}$, differentiate with respect to $w$, and
multiply the resulting expression by $\varrho^{\star}$. Then

$(\beta
w+\partial_{w})\varrho^{\star}(w)-2P\int_{\mathbb{R}}\mathrm{d}w^{\prime}\frac{1}{w-w^{\prime}}\varrho^{\star}(w)\varrho^{\star}(w^{\prime})=0.$
(3.64)

Note that $\beta$ scales by setting

$\varrho^{\star}_{P,\beta}(w)=\sqrt{\beta}\varrho^{\star}_{P,1}(\sqrt{\beta}w).$
(3.65)

Hence, for simplicity, we set $\beta=1$, omit the explicit dependence on $P$,
and denote by $\varrho^{\star}$ the solution to (3.64).

Figure 5: For the Toda lattice shown is the current DOS $\rho_{\mathrm{J}}$
and as inset the corresponding DOS $\rho_{\mathrm{Q}}$. For panels (a), (b),
(c) the Toda lattice is in thermal equilibrium with $\beta=1$ and
$(\nu,P)=(7.64,0.138),\,(-1.28,4.10),\,(0,1.461)$, respectively. Panel (c)
corresponds to the critical pressure $P_{\mathrm{c}}$, while (a) and (b) have
the same value of $\mu$, thereby confirming that the TBA equation has two
distinct solutions. Panel (d) displays the same simulation for a quartic
confining potential at $P=2$, $\beta=1$. For the DOS, sampled are the
eigenvalues of the Lax matrix at $N=1024$ (red dots) and compared with the
numerical solution of the TBA equation as based on the nonlinear Fokker-Planck
equation (4.11). For the current DOS, the expression (2.25) is sampled from
the random Lax matrix and compared with $\rho_{\mathrm{J}}$ computed through
the TBA formalism, see Eq. (6.18). From Cao et al.[50].

Taking the Stieltjes transform,

$g(z)=\int_{\mathbb{R}}\mathrm{d}w\varrho{{}^{\star}}(w)\frac{1}{w-z},$ (3.66)

yields the equation

$zg(z)+\frac{d}{dz}g(z)+Pg(z)^{2}=-1.$ (3.67)

Setting $g(z)=u^{\prime}(z)/u(z)$, Eq. (3.67) transforms to the linear second
order differential equation

$u^{\prime\prime}(z)+zu^{\prime}(z)+Pu(z)=0.$ (3.68)

Since $\varrho{{}^{*}}$ is a probability density, the $|z|\to\infty$
asymptotics,

$u(z)\simeq\frac{c_{1}}{z^{P}},$ (3.69)

follows. Finally changing to the function $u(z)=\mathrm{e}^{z^{2}/4}y(z)$ one
arrives at the Schrödinger type equation

$y^{\prime\prime}(z)+\big{(}P-\tfrac{1}{2}-\tfrac{1}{4}z^{2}\big{)}y(z)=0,$
(3.70)

which can be solved in terms of parabolic cylinder functions. The appropriate
linear combination is determined by the asymptotic condition (3.69). Somewhat
unusually, the Stieltjes transform can be still inverted and yields the fairly
explicit expression

$\varrho{{}^{\star}}(w)=\frac{\Gamma(P)\mathrm{e}^{-w^{2}/2}}{\sqrt{2\pi}P|\hat{D}_{P}(w)|^{2}},\qquad\hat{D}_{P}(w)=\int_{0}^{\infty}\mathrm{d}t\,t^{P-1}\mathrm{e}^{\mathrm{i}wt}\mathrm{e}^{-\frac{1}{2}t^{2}}.$
(3.71)

Particular examples are shown in Figure 5, which however are obtained from
numerically solving a nonlinear Fokker-Planck equation rather than using
(3.71), see Chapter 4.

We reintroduce the dependence on $P,\beta$. For small $P$, the Lax off-
diagonal matrix elements $a_{j}\to 0$ and hence in leading order
$\varrho_{P,\beta}^{\star}(w)=(2\pi)^{-1/2}\exp[-\tfrac{1}{2}\beta w^{2}]$. On
the other hand, one can integrate (3.64) against $\varsigma_{n}(w)$ to obtain
a recursion relation for the even moments of $\varrho_{P,\beta}^{\star}$,
$c_{n}=\langle\varrho_{P,\beta}^{\star}(w)w^{2n}\rangle$,

$\beta c_{n}=(2n-1)c_{n-1}+\mathfrak{p}\beta\sum_{j=0}^{n-1}c_{n-1-j}c_{j},$
(3.72)

with $n=1,2,\ldots\,$ and $c_{0}=1$, where we switched to the physical
pressure $\mathfrak{p}=P/\beta$. For low temperatures, $\beta\to\infty$, the
term $(2n-1)c_{n-1}(\beta)$ can be neglected implying the asymptotic result
$c_{n}(\beta)\to c_{n}(\infty)$ with

$c_{n}(\infty)=\mathfrak{p}^{n}\frac{1}{n+1}\begin{pmatrix}2n\\\
n\end{pmatrix}.$ (3.73)

On the right side one notes the Catalan numbers. Hence
$\varrho_{\mathfrak{p}\infty,\infty}^{\star}$ is the normalized Wigner semi-
circle probability density function

$\varrho_{\mathfrak{p}\infty,\infty}^{\star}(w)=\frac{1}{2\pi\mathfrak{p}}\sqrt{4\mathfrak{p}-w^{2}}$
(3.74)

for $|w|\leq 2\sqrt{\mathfrak{p}}$. To obtain the Lax DOS one still has to act
with the operator $\partial_{P}P$. The low pressure Gaussian does not change.
For low temperatures. upon applying the operator
$\partial_{\mathfrak{p}}\mathfrak{p}$ one concludes the convergence

$\lim_{\beta\to\infty,P=\mathfrak{p}\beta}\rho_{\mathrm{Q}}(w)=\frac{1}{\pi\sqrt{4\mathfrak{p}-w^{2}}}$
(3.75)

for $|w|\leq 2\sqrt{\mathfrak{p}}$ in the sense of convergence of moments. The
Lax eigenvalues concentrate close to the borders $\pm 2\sqrt{\mathfrak{p}}$.
In probability theory the density (3.75) is known as centered arcsine law.

This convergence can also be obtained from the random walk expansion of
$\langle Q^{[n]}_{0}\rangle_{\mathfrak{p}\beta,\beta}$. Then every path with
at least one horizontal step carries a weight with at least one factor of
$\beta^{-1}$. Thus in the limit one has to count only the number of simple
random walks with $n$ steps starting at $j=0$ and ending at $j=0$, each step
carrying the weight $\mathfrak{p}^{1/2}$. This number is nonzero only for even
$n$, in which case $\langle
Q^{[2n]}_{0}\rangle_{\mathfrak{p}\infty,\infty}=(n+1)c_{n(\infty)}$.

###  3.5 Mean-field techniques

Not to interrupt the line of arguments, we only claimed that the free energy
for the Dumitriu-Edelman partition function can be obtained by minimizing the
free energy functional $\mathcal{F}_{\mathrm{de}}^{\circ}$. In fact, in all
integrable models studied so far, both classical and quantum, there is such a
variational characterization for the generalized free energy. Thus it is
worthwhile to elucidate the main characteristics independently from the
specific applications we are interested in. From a probabilistic perspective
the problem belongs to the theory of large deviations, as pioneered by Donsker
and Varadhan. For our purposes a more simple-minded statistical mechanics
approach will do. But beware: we only explain the main steps. More space would
be needed for a complete proof in particular cases.

We will outline what mathematicians call a proof by compactness. A simple
example would be a sequence of real numbers $y_{1},y_{2},\ldots$ with
$|y_{j}|\leq a$, the condition of compactness. One wants to determine the
large $j$ limit $y_{j}\to y_{\infty}$. Due to the uniform bound, one can
choose a subsequence $\tilde{y}_{j}$ such that
$\tilde{y}_{j}\to\tilde{y}_{\infty}$ for $j\to\infty$. The crucial input is
some information on $\tilde{y}_{\infty}$. In our context a natural example
would be to have some function $g:[-a,a]\to\mathbb{R}$ with the property that
$g(\tilde{y}_{\infty})=0$. Assuming in addition that $g$ has a unique zero,
$\bar{y}$, one concludes that $\lim_{j\to\infty}\tilde{y}_{j}=\bar{y}$. Since
this argument applies to any convergent subsequence, one arrives at
$\lim_{j\to\infty}y_{j}=\bar{y}$. To obtain some information on the rate of
convergence would require additional considerations.

We turn to the mean-field type problem by considering $N$ “spins”
$x_{1},\ldots,x_{N}$ with $x_{j}\in\mathbb{R}$. The “energy” of a spin
configuration is taken to be of the particular form

$H_{N}(x)=\sum_{j=1}^{N}V(x_{j})+\frac{1}{N}\sum_{i,j=1}^{N}W(x_{i},x_{j}).$
(3.76)

Thereby one defines partition function and Gibbs probability density relative
to $\mathrm{d}^{N}x$,

$Z_{N}=\int_{\mathbb{R}^{N}}\mathrm{d}^{N}x\exp\big{(}-H_{N}(x)\big{)},\qquad\mu^{N}(x)=Z_{N}^{-1}\exp\big{(}-H_{N}(x)\big{)}.$
(3.77)

The superscript $N$ refers to a probability density on $\mathbb{R}^{N}$ with
standard definitions for free energy, energy, and entropy,

$F(\mu^{N})=-\log Z_{N},\qquad E(\mu^{N})=\langle
H_{N}\rangle_{\mu^{N}},\qquad
S(\mu^{N})=-\int_{\mathbb{R}^{N}}\mathrm{d}^{N}x\mu^{N}\log\mu^{N}.$ (3.78)

As before, the main goal is the large $N$ limit of the empirical density

$\varrho_{N}(w)=\frac{1}{N}\sum_{j=1}^{N}\delta(w-x_{j}).$ (3.79)

The energy of our model is very special in the sense that it is merely a
functional of the empirical density,

$H_{N}=N\Big{(}\int_{\mathbb{R}}\mathrm{d}w\varrho_{N}(w)\big{(}V(w)+\int_{\mathbb{R}}\mathrm{d}w^{\prime}\varrho_{N}(w^{\prime})W(w,w^{\prime})\big{)}\Big{)}.$
(3.80)

The quadratic functional serves here only as illustration. More complicated
functionals will appear in Chapters 9 and 11. Of course the potential $V$ has
to be sufficiently confining. In a first run the kernel
$W(w,w^{\prime})=W(w^{\prime},w)$ could taken to be bounded. As central
observation, $\mu^{N}$ is invariant under permutation of indices. As a
striking consequence, for example, the two-point function $\langle
f_{1}(x_{i})f_{2}(x_{j})\rangle_{\mu^{N}}$ does not depend on the choice of
$(i,j)$. Clearly, such a large symmetry group will force special properties.

The mentioned compactness can be ensured by bounds as
$Z_{N}\leq\mathrm{e}^{c_{0}N}$. Thus there exists a subsequence, denoted again
as $\mu^{N}\to\mu_{\infty}$, where the limit is a probability measure on
$\mathbb{R}^{\mathbb{N}}$. Now by a theorem of Hewitt and Savage, in a
simplified context earlier studied by De Finetti, the limit measure must be a
convex combination of product probability measures. A pure product measure is
written as

$\mu_{\varrho}=\varrho\times\varrho\times\ldots$ (3.81)

with $\varrho\geq 0$ and $\langle\varrho\rangle=1$. Each factor has finite
entropy and hence is of the form $\varrho(w)\mathrm{d}w$ with finite

$\mathsfit{S}(\varrho)=-\int_{\mathbb{R}}\mathrm{d}w\varrho(w)\log\varrho(w).$
(3.82)

The entropy and energy per site equal

$s(\mu_{\varrho})=\mathsfit{S}(\varrho),\qquad e(\mu_{\varrho})=\langle
V,\varrho\rangle+\langle\varrho,W\varrho\rangle=\mathsfit{E}(\varrho).$ (3.83)

A convex combination of product measures is written as

$\mu=\int\nu(\mathrm{d}\varrho|\mu)\mu_{\varrho}.$ (3.84)

The set of such convex combinations is denoted by $\mathcal{S}$. Entropy and
energy per site now turns to

$s(\mu)=\int\nu(\mathrm{d}\varrho|\mu)\mathsfit{S}(\varrho),\qquad
e(\mu)=\int\nu(\mathrm{d}\varrho|\mu)\mathsfit{E}(\varrho).$ (3.85)

While the second identity is obvious, the first one looks wrong and is so at
finite volume. But consider
$\alpha\varrho_{1}\times\ldots\times\varrho_{1}+(1-\alpha)\varrho_{2}\times\ldots\times\varrho_{2}$,
each product having $N$ factors. For large $N$ these densities are supported
on essentially disjoint sets in phase space and the integral defining the
total entropy is in leading approximation $\alpha
N\mathsfit{S}(\varrho_{1})+(1-\alpha)N\mathsfit{S}(\varrho_{2})$.

We now choose a (sub)sequence such that $\mu^{N}\to\mu_{\infty}$. Then the
limit

$\lim_{N\to\infty}\frac{1}{N}E(\mu^{N})=e(\mu_{\infty})$ (3.86)

exists.
(i) lower bound. For $\mu^{N}$ the marginal onto $[1,\ldots,n]$ is denoted by
$\mu^{N}_{n}$ and correspondingly by $\mu_{n}$ for some infinite volume
measure by $\mu$. We consider large $N$ and divide the volume into disjoint
segments each of length $n$. Then, by subadditivity of entropy,

$S(\mu^{N})\leq\Big{\lfloor}\frac{N}{n}\Big{\rfloor}S(\mu^{N}_{n})+\mathcal{O}(1)$
(3.87)

with $\lfloor\cdot\rfloor$ indicating integer part. Therefore

$\lim_{N\to\infty}\frac{1}{N}S(\mu^{N})\leq\lim_{n\to\infty}\frac{1}{n}s(\mu_{\infty,n})=s(\mu_{\infty})$
(3.88)

and, since $f=e-s$,

$\lim_{N\to\infty}\frac{1}{N}F(\mu^{N})\geq
e(\mu_{\infty})-s(\mu_{\infty})=f(\mu_{\infty}).\vskip 3.0pt plus 1.0pt minus
1.0pt$ (3.89)

(ii) upper bound. By the standard finite volume variational principle for the
free energy

$\lim_{N\to\infty}\frac{1}{N}F(\mu^{N})\leq\lim_{N\to\infty}\min_{\mu^{\prime}\in\mathcal{S}}\frac{1}{N}F(\mu^{\prime}_{N})=\min_{\mu^{\prime}\in\mathcal{S}}f(\mu^{\prime})\leq
f(\mu_{\infty}).\vskip 3.0pt plus 1.0pt minus 1.0pt$ (3.90)

Combining both bounds

$\lim_{N\to\infty}\frac{1}{N}F(\mu^{N})=f(\mu_{\infty})=\min_{\mu^{\prime}\in\mathcal{S}}\int\nu(\mathrm{d}\varrho^{\prime}|\mu^{\prime})\mathcal{F}(\varrho^{\prime}),$
(3.91)

where the one-particle free energy functional reads

$\mathcal{F}(\varrho)=\big{\langle}\varrho,V+\log\varrho+W\varrho\big{\rangle}.$
(3.92)

In mean-field models from equilibrium statistical mechanics generically there
are several minimizers, indicating a phase transition in the microscopic
system. For integrable many-body systems uniqueness seems to be the rule.
Assuming such property to hold and denoting the unique minimizer by
$\varrho^{*}$, our formulas simplify drastically,

$\lim_{N\to\infty}\frac{1}{N}F(\mu^{N})=\cal{F}(\varrho^{*}).$ (3.93)

Modulo technical points, this is the limit claimed in (3.29). Furthermore,
subsequence can be replaced by sequence and

$\mu^{N}\to\mu_{\varrho^{*}}.$ (3.94)

The density of states, integrated against some test function $f$, equals
$\langle f\varrho_{N}\rangle$ which is still random. From (3.94) one deduces
the limits

$\lim_{N\to\infty}\big{\langle}\langle
f\varrho_{N}\rangle\big{\rangle}_{\mu^{N}}=\langle
f\varrho^{*}\rangle,\qquad\lim_{N\to\infty}\big{\langle}\big{(}\langle
f\varrho_{N}\rangle\big{)}^{2}\big{\rangle}_{\mu^{N}}=\langle
f\varrho^{*}\rangle^{2},$ (3.95)

which establishes the law of large numbers for $\langle f\varrho_{N}\rangle$,
in other words the almost sure convergence of $\langle f\varrho_{N}\rangle$ to
the deterministic limit $\langle f\varrho^{*}\rangle$ as $N\to\infty$.

In statistical mechanics mean-field usually refers to an approximation where a
fluctuating field is replaced by its average. In contrast, for integrable
many-body systems the respective mean-field problem is exact and arrives
unexpectedly. For the Toda lattice the transformation of Dumitriu-Edelman
accomplishes the deal. As to be discussed, for the Calogero fluid the proof
relies on scattering coordinates and for the Ablowitz-Ladik system on the
transformation linked to CMV matrices. For quantum many-body systems a mean-
field type problem naturally arrives through the Bethe ansatz, more generally
through the asymptotic Bethe ansatz.

Notes and references

Section 3.0

At age twenty-four Ludwig Boltzmann [27] wrote his fundamental contribution on
the microcanonical ensemble. He argued that this ensemble provides the natural
description of the long time behavior of a mechanical system. In addition, he
uncovered the connection to thermodynamics. Recommended is his very readable
letter-style account Boltzmann [28]. The way how we teach equilibrium
statistical mechanics today goes back to the must-read book by J. Williard
Gibbs [132].

Section 3.1

The generalized Gibbs ensemble came naturally into focus when studying the
quench of integrable systems starting from a translation invariant state, see
the reviews Polkovnikov et al. [234], Vidmar and Rigol [294], the special
volume edited by Calabrese et al. [43], and the more recent short review with
many references by Alba and Calabrese [11]. In this context the conserved
charges averaged over the initial state determine the parameters of the GGE
obtained in the long time limit. In some circles the notion “generalized” was
initially not so welcome, since there is a long tradition in the study of
Gibbs measures for a general class of potentials, for a minute subset see
Ruelle [246], Lanford [187], Friedli and Velenik [122]. In our context
“generalized” means that for the given mechanical system one has a high
dimensional set of time-invariant measures all of them having the form
anticipated by Statistical Mechanics.

Equivalence of ensembles in its original meaning refers to the property that
the resulting thermodynamic potentials are related to each other through the
respective Legendre transform. This is a widely studied subject, with close
relations to the theory of large deviations, Lanford [187], Simon [261],
Friedli and Velenik [122]. In our context we use a stronger version: In the
limit of infinite volume, two distinct ensembles constructed from the same
bulk hamiltonian yield identical distributions of strictly local observables.
To put it differently, provided the respective parameters are transformed
according to the rules of thermodynamics, correlations of local functions are
identical. Such a stricter version generically fails at phase transitions, for
a detailed discussion see Georgii [130, 131].

Section 3.2

A general problem of nonequilibrium statistical mechanics is to determine the
slowly varying fields. The obvious candidates are the locally conserved
fields. In higher dimensions the equilibrium phase diagram has more structure.
For example, gas and fluid phase may coexist. The respective phase boundary is
then also a slowly varying field. More generally, the order parameters of the
thermodynamic phases have to be included in the list, see Forster [120].
Generalized hydrodynamics (GHD) follows the same approach. For the Toda
lattice the Lax filter is a natural construction, which is applicable to other
integrable systems whenever there is a Lax matrix available.

Section 3.3

The generalized free energy of the Toda lattice has been obtained in Spohn
[269], see also Doyon [89, 87] for complimentary discussions. Originally,
Dumitriu and Edelman [96] investigated how the fully occupied GUE random
matrix transforms isospectrally to a tridiagonal matrix. For this purpose they
iteratively applied the Householder transform, which is a standard numerical
scheme. Once Dumitriu and Edelman had understood how this scheme works for GUE
random matrices, they realized that their algebra holds for arbitrary values
of the parameter $\beta$ and not only for $\beta=2$. In their context the
quadratic confining potential appears naturally. Since the transformation of
volume elements is established, the extension to general confining potentials
is straightforward. The Dumitriu-Edelman identity has been used to obtain a
fairly detailed information on the edge behavior of the general
$\beta$-ensemble, see Ramirez et al. [242] and Bourgade et al. [34] for
further developments. A standard reference on log gases is the monumental
volume by Forrester [118]. Variational problems with a logarithmic kernel are
investigated by Saff and Totik [251].

In a recent contribution Guionnet and Memin [144] prove the generalized free
energy to be given by (3.31) and the almost sure limit of the DOS as stated in
(3.16) with limit $\nu\rho_{\mathsf{p}}$. Their theorem covers confining
potentials such that $\lim_{|w|\to\infty}w^{-2n}V(w)=a_{\infty,n}>0$ for some
positive integer $n$. In particularly, $V$ is not necessarily given by a
convergent power series. From the hydrodynamic perspective, required would be
also the inverse operation. One chooses $\rho_{\mathsf{n}}$ from a specified
class of functions and concludes that $V$ is determined by Eq. (3.50). A
further line of research, not yet covered at all, concerns the spatial
structure of GGEs. This includes the decay of correlation functions and mixing
properties. For such an analysis of use could be Dobrushin’s statistical
mechanics theory, Dobrushin [80], developed specifically for one-dimensional
systems.

Section 3.4

The particular form of the TBA formalism can be better grasped upon reading
Section 13. The double-valuedness of solutions to the TBA equation was pointed
out to me by Bulchandani et al. [39] with numerical confirmation in Cao et al.
[50]. Numerical plots of the densities as they vary through $P_{\mathrm{c}}$
can be found in Mendl and Spohn [204]. The TBA equation (3.50) with quadratic
confining potential was first obtained by Opper [226], who investigated the
classical limit of the quantum Toda chain. He already obtained the solution
(3.71). Our discussion is based on Allez et al. [13], a study of Dyson
Brownian motion with weak interaction, see Chapter 4 for details. A related
study has been carried out for $\beta$-Wishart ensembles by Allez et al. [14].
The recursion relation (3.72) appears already in Duy and Shirai [98] and Duy
[97].

Section 3.5

Out of the many texts and monographs on the theory of large deviations
recommended is the book by Dembo and Zeitouni [73]. The complete proof based
on permutation symmetry can be found in Messer and Spohn [206]. A more recent
study are the lecture notes by Rougerie [245] who also presents an anologous
discussion of quantum mean-field models. De Finetti’s theorems on exchangeable
measures are discussed in Aldous [12]. The generality required in our context
is covered by Hewitt and Savage [153].

##  4 Dyson Brownian motion

In the mid-sixties Freeman Dyson pioneered the study of random matrices, as
initiated by Eugene Wigner and others as a phenomenological model for the
complex energy spectra of highly excited nuclei. More specifically, Dyson
studied the statistics of eigenvalues for the Gaussian orthogonal, unitary,
and symplectic ensembles. For this purpose he introduced a stochastic dynamics
of eigenvalues with the property that its stationary distribution coincides
with the distribution of eigenvalues of the respective random matrix ensemble.
Such a diffusion process is now called Dyson Brownian motion. In this
approach, one studies a physically more intuitive stochastic particle system,
which in the long time limit approaches the desired probability density of
eigenvalues. Particles move in one dimension, their positions being denoted by
$\\{x_{j}(t),j=1,\ldots,N\\}$, and are governed by the coupled stochastic
differential equations,

$dx_{j}(t)=-V^{\prime}(x_{j}(t))dt+2\beta\sum_{i=1,i\neq
j}^{N}\frac{1}{x_{j}(t)-x_{i}(t)}dt+\sqrt{2}db_{j}(t),$ (4.1)

for $j=1,\ldots,N$ and with $\beta\geq 0$. The particles repel each other with
a $1/x$ force and are confined by the external potential $V$, which in fact
will coincide with the generalized chemical potential (3.5) of the Toda
lattice – in the first place the reason for using the same name and symbol.
$\\{b_{j}(t),j=1,\ldots,N\\}$ is a collection of independent standard Brownian
motions. For later purposes the generator, $\mathcal{L}_{N}$, of the coupled
diffusion processes is defined as the linear operator

$\mathcal{L}_{N}\mathsfit{f}(x)=\sum_{j=1}^{N}\Big{(}\partial_{x_{j}}-V^{\prime}(x_{j})+2\beta\sum_{i=1,i\neq
j}^{N}\frac{1}{x_{j}-x_{i}}\Big{)}\partial_{x_{j}}\mathsfit{f}(x),$ (4.2)

acting on functions on configuration space,
$\mathsfit{f}:\mathbb{R}^{N}\to\mathbb{R}$, $x=(x_{1},\ldots,x_{N})$. For
$t\geq 0$, the kernel of the semigroup

$\mathrm{e}^{\mathcal{L}_{N}t}(x,x^{\prime})\mathrm{d}^{N}x^{\prime}$ (4.3)

is the transition probability, in other words, the probability density at time
$t$ given the initial configuration $x$.

The reader might suspect a luxury detour. But this is not the case. Dyson
Brownian motion is a powerful method to numerically solve the TBA equations of
the Toda lattice. Even more importantly, as will be discussed, the GGE
averaged currents are linked to eigenvalue fluctuations. To study such
properties Dyson Brownian motion is a convenient tool.

Dyson investigated the particular values $\beta=1,2,4$, which possess a large
symmetry group. But recently there has been much progress also for general
$\beta$. In our context, we will have to study the case of very small $\beta$,
specifically $\beta=\alpha/N$. In fact, $\alpha=P$ in the application to the
Toda lattice. In (4.1) the drift is the gradient of a potential. Hence the
diffusion process is reversible and its unique stationary measure is given by

$(Z_{\alpha,V,N})^{-1}\exp\\!\Big{(}-\sum_{j=1}^{N}V(x_{j})+\frac{\alpha}{N}\sum_{i,j=1,i\neq
j}^{N}\log|x_{i}-x_{j}|\Big{)},$ (4.4)

which agrees with the normalized version of the probability density function
appearing in (3.3). Our strategy will be to first study the large $N$ limit of
the dynamics and then deduce information on the stationary measure, which is
our real interest.

###  4.1 Macroscopic equation, law of large numbers

We choose some smooth test function $f$ and introduce the empirical density,