$\varUpsilon_{N}(x,t)$, through

$\varUpsilon_{N}(f,t)=\frac{1}{N}\sum_{j=1}^{N}f(x_{j}(t))=\int_{\mathbb{R}}\mathrm{d}x\varUpsilon_{N}(x,t)f(x).$
(4.5)

Then

$\displaystyle
d\varUpsilon_{N}(f,t)=\varUpsilon_{N}(-V^{\prime}\partial_{x}f+\partial_{x}^{2}f,t)dt$
(4.6) $\displaystyle\hskip
10.0pt+\alpha\int_{\mathbb{R}}\mathrm{d}x\int_{\mathbb{R}}\mathrm{d}y\frac{f^{\prime}(x)-f^{\prime}(y)}{x-y}\varUpsilon_{N}(x,t)\varUpsilon_{N}(y,t)dt+\frac{1}{N}\sum_{j=1}^{N}f^{\prime}(x_{j}(t))\sqrt{2}db_{j}(t).$

When integrated in time, the process averaged square of the noise term becomes

$\frac{2}{N^{2}}\sum_{j=1}^{N}\int_{0}^{t}\mathrm{d}s\mathbb{E}\big{(}f(x_{j}(s))^{2}\big{)}.$
(4.7)

Thus the noise is of order $1/\sqrt{N}$ and vanishes in the limit
$N\to\infty$.

We assume a starting measure such that with probability one the initial
empirical density has the limit $\varrho_{0}$,

$\lim_{N\to\infty}\varUpsilon_{N}(f,0)=\int_{\mathbb{R}}\mathrm{d}x\varrho_{0}(x)f(x).$
(4.8)

In the statistical physics literature this property is often referred to as
“self-averaging”, meaning that the limit is non-random and no averaging is
required. Since in (4.6) only the drift term survives, one concludes the limit

$\lim_{N\to\infty}\varUpsilon_{N}(f,t)=\varUpsilon(f,t)=\int_{\mathbb{R}}\mathrm{d}x\varrho(x,t)f(x)$
(4.9)

for all $t>0$, again with probability one. The limit density $\varrho(x,t)$
then satisfies

$\displaystyle\frac{d}{dt}\int_{\mathbb{R}}\mathrm{d}x\varrho(x,t)f(x)$ (4.10)
$\displaystyle=\int_{\mathbb{R}}\mathrm{d}x\varrho(x,t)(-V^{\prime}\partial_{x}f+\partial_{x}^{2}f)(x)+\alpha\int_{\mathbb{R}}\mathrm{d}x\int_{\mathbb{R}}\mathrm{d}y\frac{f^{\prime}(x)-f^{\prime}(y)}{x-y}\varrho(x,t)\varrho(y,t)$

together with the initial data $\varrho_{0}(x)$. Written pointwise,
$\varrho(x,t)$ is governed by the nonlinear Fokker-Planck equation

$\partial_{t}\varrho(x,t)=\partial_{x}\big{(}V_{\mathrm{eff}}^{\prime}(x,t)\varrho(x,t)+\partial_{x}\varrho(x,t)\big{)}.$
(4.11)

The bare confining potential $V$ is modified to an effective potential given
by

$V_{\mathrm{eff}}(x,t)=V(x)-\alpha(T\varrho)(x,t),$ (4.12)

the integral operator $T$ being defined in (3.55).

Our particular interest is the stationary Fokker-Planck equation

$\partial_{x}\Big{(}V^{\prime}(x)\varrho_{\mathfrak{s}}(x)-2\alpha\int_{\mathbb{R}}\mathrm{d}y\frac{1}{x-y}\varrho_{\mathfrak{s}}(y)\varrho_{\mathfrak{s}}(x)+\partial_{x}\varrho_{\mathfrak{s}}(x)\Big{)}=0,\quad\int_{\mathbb{R}}\mathrm{d}x\varrho_{\mathfrak{s}}(x)=1,$
(4.13)

which under the constraint

$\int_{\mathbb{R}}\mathrm{d}x\varrho_{\mathfrak{s}}(x)=1,$ (4.14)

has a unique strictly positive solution $\varrho_{\mathfrak{s}}$. Since Dyson
Brownian motion is time-reversible, the large round bracket itself has to
vanish. We compare to the mean-field free energy functional
$\mathcal{F}^{\circ}_{\mathrm{de}}(\varrho)$ in Eq. (3.28), setting
$P=\alpha$. Its minimizer $\varrho^{\star}$ satisfies the Euler-Lagrange
equation

$V(x)-\mu-2\alpha\int_{\mathbb{R}}\mathrm{d}x^{\prime}\log|x-x^{\prime}|\varrho^{\star}(x^{\prime})+\log\varrho^{\star}(x)+1=0,$
(4.15)

where $w$ has been substituted by $x$. Differentiating with respect to $x$ and
then multiplying by $\varrho^{\star}$ yield

$V^{\prime}(x)\varrho^{\star}(x)-2\alpha\int_{\mathbb{R}}\mathrm{d}y\frac{1}{x-y}\varrho^{\star}(x)\varrho^{\star}(y)+\partial_{x}\varrho^{\star}(x)=0.$
(4.16)

Thus we conclude $\varrho_{\mathfrak{s}}=\varrho^{\star}$.

The nonlinear Fokker-Planck equation depends smoothly on $\alpha$ and so does
$\varrho_{\mathfrak{s}}$. The issue of two-valuedness appears only in the TBA
equation.

###  4.2 Fluctuation theory

The nonlinear Fokker-Planck equation is the result of a law of large numbers.
Thus one would expect to have a central limit type theorem which captures the
order $1/\sqrt{N}$ fluctuations. Quite generically the fluctuations are
governed by a linear Langevin equation for the conserved field, in our case
the density field. The drift term is determined by the linearized macroscopic
equation. But in addition the dynamics generates an effective noise term,
whose precise structure depends on the model.

Our interest is the stationary dynamics, hence the initial state of the
$N$-particle system is the one in Eq. (4.4). To study the fluctuations of the
density it is convenient to introduce the fluctuation field as

$\phi_{N}(f,t)=\frac{1}{\sqrt{N}}\sum_{j=1}^{N}\big{(}f(x_{j}(t))-\langle\varrho_{\mathfrak{s}}f\rangle\big{)}=\int_{\mathbb{R}}\mathrm{d}xf(x)\phi_{N}(x,t).$
(4.17)

As explained in more detail below, there is a Gaussian random field
$\phi(f,t)$, jointly in $f$ and $t$, such that in distribution

$\lim_{N\to\infty}\phi_{N}(f,t)=\phi(f,t).$ (4.18)

Note that the limit is still random. The limit field $\phi$ is governed by the
linear Langevin equation

$\partial_{t}\phi(x,t)=\partial_{x}\big{(}D\phi(x,t)+\sqrt{2\varrho_{\mathfrak{s}}(x)}\xi(x,t)\big{)}.$
(4.19)

Here $\xi(x,t)$ is normalized spacetime Gaussian white noise,
$\mathbb{E}\big{(}\xi(x,t)\xi(x^{\prime},t^{\prime})\big{)}=\delta(x-x^{\prime})\delta(t-t^{\prime})$
and $D$ the linear operator

$D=V^{\prime}_{\mathrm{eff}}+\partial_{x}-\alpha\varrho_{\mathfrak{s}}T^{\prime}.$
(4.20)

Thus $\partial_{x}D$ is the Fokker-Planck evolution operator linearized at
$\varrho_{\mathfrak{s}}$. The effective potential $V_{\mathrm{eff}}$ is still
defined as in (4.12) upon substituting $\varrho(x,t)$ by
$\varrho_{\mathfrak{s}}(x)$.

The Gaussian process $\phi(x,t)$ is stationary in time, has mean zero and is
uniquely characterized by its covariance
$\mathbb{E}\big{(}\phi(x,t)\phi(x^{\prime},t^{\prime})\big{)}$. For the Toda
lattice of interest is the spatial covariance

$\mathbb{E}\big{(}\phi(x,0)\phi(x^{\prime},0)\big{)}=C^{\sharp}(x,x^{\prime}).$
(4.21)

As a general property of linear Langevin equations with time-independent
coefficients, $C^{\sharp}$ is determined by

$\langle D^{*}\partial_{x}f,C^{\sharp}g\rangle+\langle
f,C^{\sharp}D^{*}\partial_{x}g\rangle=2\langle\partial_{x}f,\varrho_{\mathfrak{s}}\partial_{x}g\rangle$
(4.22)

with $D^{*}$ the adjoint operator to $D$. We claim that, as an operator, the
solution is

$C^{\sharp}=(1-\alpha\varrho_{\mathfrak{s}}T)^{-1}\varrho_{\mathfrak{s}}-\big{\langle}(1-\alpha\varrho_{\mathfrak{s}}T)^{-1}\varrho_{\mathfrak{s}}\big{\rangle}^{-1}\big{|}(1-\alpha\varrho_{\mathfrak{s}}T)^{-1}\varrho_{\mathfrak{s}}\big{\rangle}\big{\langle}(1-\alpha\varrho_{\mathfrak{s}}T)^{-1}\varrho_{\mathfrak{s}}\big{|},$
(4.23)

where for simplicity we use the Dirac notation $|\cdot\rangle\langle\cdot|$
for the one-dimensional projector. The subtracted term ensures that the number
of particles does not fluctuate, i.e. $C^{\sharp}\varsigma_{0}=0$.

To confirm the claim, we consider only the left most term of (4.22), the other
one following by symmetry, and have to show that

$\langle\partial_{x}f,D\varrho_{\mathfrak{s}}(1-\alpha
T\varrho_{\mathfrak{s}})^{-1}g\rangle=\langle\partial_{x}f,\varrho_{\mathfrak{s}}\partial_{x}g\rangle.$
(4.24)

Upon replacing $g$ by $(1-\alpha T\varrho_{\mathfrak{s}})g$, one arrives at

$\langle\partial_{x}f,D\varrho_{\mathfrak{s}}g\rangle=\langle\partial_{x}f,\varrho_{\mathfrak{s}}\partial_{x}(1-\alpha
T\varrho_{\mathfrak{s}})g\rangle.$ (4.25)

$\varrho_{\mathfrak{s}}$ satisfies the stationary nonlinear Fokker-Planck
equation (4.16),

$\big{(}V_{\mathrm{eff}}^{\prime}+\partial_{x}\big{)}\varrho_{\mathfrak{s}}=0.$
(4.26)

When inserted in (4.25) this leads to the condition

$\varrho_{\mathfrak{s}}\partial_{x}g-\alpha\varrho_{\mathfrak{s}}T^{\prime}(\varrho_{\mathfrak{s}}g)=\varrho_{\mathfrak{s}}\partial_{x}(1-\alpha
T\varrho_{\mathfrak{s}})g,$ (4.27)

which is easily checked.

$\blackdiamond\blackdiamond$ Martingales and central limit theorem. In
probability theory tremendous efforts have been invested to develop tools for
proving the central limit theorem in case of dependent random variables. Such
techniques can be applied to our infinite-dimensional setting, which is
required since the test functions form an infinite-dimensional linear space.
Only the basic computational steps are explained here. A full proof would be
too technical. For example, we would have to discuss the existence of
solutions to (4.1). In fact, if $\beta>1$, trajectories never touch and the
solution theory is standard. But for $\beta<1$ the crossing probability is
strictly positive and the existence of solutions becomes more intricate. The
proof of the limit (4.18) is based on compactness. In an appropriate function
space one ensures that along some subsequence the stochastic process of
fluctuations has a limit. The limit process satisfies an equation for which
uniqueness is established, which then implies convergence. We consider only
the stationary process, but the time-dependent case can be handled in a
similar fashion.

Using (4.2), one arrives at

$\mathcal{L}_{N}\phi_{N}(f)=\phi_{N}(-V^{\prime}\partial_{x}f+\partial_{x}^{2}f)+\alpha\int_{\mathbb{R}}\mathrm{d}x\int_{\mathbb{R}}\mathrm{d}y\frac{f^{\prime}(x)-f^{\prime}(y)}{x-y}\phi_{N}(x)\varUpsilon_{N}(y)=\varkappa_{[1],N}(f).$
(4.28)

For the quadratic variation we obtain

$\mathcal{L}_{N}\phi_{N}(f)^{2}-2\phi_{N}(f)\mathcal{L}_{N}\phi_{N}(f)=2N^{-1}\sum_{j=1}^{N}f^{\prime}(x_{j})^{2}=\varkappa_{[2],N}(f).$
(4.29)

As before, $\varkappa_{[1],N}(f,t)$ is the function $\varkappa_{[1],N}(f)$
evaluated at the random configuration
$(x_{1}(t),\ldots,x_{N}(t))$, and the same for $\varkappa_{[2],N}(f,t)$. By
the standard theory of Markov processes,

$M_{[1],N}(f,t)=\phi_{N}(f,t)-\phi_{N}(f,0)-\int_{0}^{t}\mathrm{d}s\varkappa_{[1],N}(f,s)$
(4.30)

and

$M_{[2],N}(f,t)=M_{[1],N}(f,t)^{2}-\int_{0}^{t}\mathrm{d}s\,\varkappa_{[2],N}(f,s)$
(4.31)

are martingales.

By compactness one first ensures that the limit in (4.18) exists. Since the
martingale property is preserved in the limit $N\to\infty$, using the law of
large numbers to handle $\varkappa_{[1],N}(f)$ and $\varkappa_{[2],N}(f)$, one
concludes that

$\displaystyle
M_{[1]}(f,t)=\phi(f,t)-\phi(f,0)-\int_{0}^{t}\mathrm{d}s\,\phi(-V^{\prime}\partial_{x}f+\partial_{x}^{2}f,s)$
$\displaystyle\hskip
20.0pt-\int_{0}^{t}\mathrm{d}s\,\alpha\int_{\mathbb{R}}\mathrm{d}x\int_{\mathbb{R}}\mathrm{d}y\frac{f^{\prime}(x)-f^{\prime}(y)}{x-y}\phi(x,t)\varrho_{\mathfrak{s}}(y)$
(4.32)

and

$M_{[2]}(f,t)=M_{[1]}(f,t)^{2}-2t\int_{\mathbb{R}}\mathrm{d}x\varrho_{\mathfrak{s}}(x)f^{\prime}(x)^{2}$
(4.33)

are still martingales under the limit process $\phi(f,t)$. The unique solution
to the latter martingale problem is the linear Langevin equation of (4.19),
thereby confirming the claim. $\blackdiamond\blackdiamond$

Notes and references

Section 4.0

Freeman Dyson introduced his model in Dyson [99]. An interesting account of
the early history is his preface to The Oxford Handbook on Random Matrix
Theory edited by Akemann et al. [10]. The Handbook serves as a useful source
of information and allows one to capture the vastness of the subject. An
introductory monograph on random matrices is Anderson et al. [16]. Dyson
Brownian motion also serves as a technical tool in the study of universality
for Wigner random matrices, see Erdős et al. [105], and the edge behavior of
the density of states of the beta random matrix ensembles, see Bourgade et al.
[34].

Section 4.1

For general $\beta>0$, Dyson Brownian motion is properly constructed in Cépa
and Lépingle [56], where also the law of large numbers is proved.

Section 4.2

Fluctuation theory at strong coupling is established in Israelsson [165]. The
methods used can be extended to the case discussed here. Fluctuation theory
for a wider class of stochastic particle systems is discussed in Spohn [266]
with more detailed arguments and references on the martingale method. For
thermal equilibrium, eigenvalue fluctuations are proved by Nakano and Trinh
[218] by a moment method. Hardy and Lambert [148] use arguments from optimal
transport to more directly estimate the distance between the true density
function and the approximating Gaussian.

##  5 Hydrodynamics for hard rods

So far our focus has been the generalized free energy. But hydrodynamics
relies on further building blocks. The hard rod system will serve to
illustrate the method as such, thereby building a bridge towards the Toda
lattice.

Figure 6: Displayed are the trajectories, $q_{j}(t)$, of $9$ hard rod
particles with length $0.3$ in a volume of size $5$. Very roughly the density
$\bar{\rho}=1.8$ and inverse temperature $\beta=1$, which corresponds to a
pressure $P=3.9$. The trajectory of the 8-th particle is shown in magenta and
its quasiparticle trajectory in red. The scattering shift is indicated on the
right. Courtesy of A. Kundu.

###  5.1 Hard rod fluid

The fluid consists of hard rods of length $\mathsfit{a}$,
$\mathsfit{a}\in\mathbb{R}$, which interact through elastic collisions. The
lattice version will be discussed in section 5.2. Positions are constrained by
$r_{j}\geq\mathsfit{a}$ and quasiparticles maintain their velocity through
collisions. At a collision, i.e. $r_{j}=\mathsfit{a}$, $q_{j}$ jumps
instantaneously to $q_{j}+\mathsfit{a}$ and $q_{j+1}$ to
$q_{j+1}-\mathsfit{a}$. Besides the number of particles, any sum function of
the velocities is conserved. This field is labelled by some general function,
$\phi$, and has a density given by

$Q^{[\phi]}(x)=\sum_{j\in\mathbb{Z}}\delta(q_{j}-x)\phi(p_{j}).$ (5.1)

The particle number is included and corresponds to the constant function,
$\phi(w)=1$. The rod length can be negative. In fact, if $\mathsfit{a}<0$,
then for sufficiently high pressure the free volume turns negative, a property
resembling the Toda lattice.

For the current one has to differentiate the density with respect to $t$. Then
there is the flow term resulting from the $\delta$-function and a collision
term resulting from $\phi$. For the latter we consider a small time
$\mathrm{d}t>0$ and the change $\phi(p_{j}(\mathrm{d}t))-\phi(p_{j})$. If
$r_{j}-\mathsfit{a}<\delta$ with sufficiently small $\delta$ and if
$p_{j+1}<p_{j}$, then $p_{j}(t)=p_{j+1}$. Similarly, if
$r_{j-1}-\mathsfit{a}<\delta$ and if $p_{j-1}>p_{j}$, then $p_{j}(t)=p_{j-1}$,
which implies

$\displaystyle\phi(p_{j}(\mathrm{d}t))-\phi(p_{j})=\delta(r_{j}-\mathsfit{a})(\phi(p_{j+1})-\phi(p_{j}))\chi(\\{p_{j+1}<p_{j}\\})(p_{j}-p_{j+1})\mathrm{d}t$
$\displaystyle\hskip
66.0pt+\delta(r_{j-1}-\mathsfit{a})(\phi(p_{j-1})-\phi(p_{j}))\chi(\\{p_{j}<p_{j-1}\\})(p_{j-1}-p_{j})\mathrm{d}t,$
(5.2)

where $\chi(\\{\cdot\\})$ is the indicator function of the set $\\{\cdot\\}$.
Thus

$\displaystyle
J^{[\phi]}(x)=\sum_{j\in\mathbb{Z}}\Big{(}\delta(q_{j}-x)p_{j}\phi(p_{j})$
(5.3)
$\displaystyle+\big{(}\theta(q_{j}-x)-\theta(q_{j-1}-x)\big{)}\delta(r_{j-1}-\mathsfit{a})\big{(}\phi(p_{j-1})-\phi(p_{j})\big{)}\chi(\\{p_{{}_{j}}<p_{j-1}\\})(p_{j-1}-p_{j})\Big{)}$

with $\theta$ the step function, $\theta(x)=0$ for $x\leq 0$ and $\theta(x)=1$
for $x>0$. Under the hard rod dynamics the local conservation law has the
standard form

$\partial_{t}Q^{[\phi]}(x,t)+\partial_{x}J^{[\phi]}(x,t)=0.$ (5.4)

A GGE is characterized by the one-particle distribution function
$f(v)=\bar{\rho}h(v)$. Here $h\geq 0$ and $\langle h\rangle=1$. The velocities
$\\{p_{j},j\in\mathbb{Z}\\}$ are a family of i.i.d. random variables with
common probability density $h$. The positions $\\{q_{j},j\in\mathbb{Z}\\}$ are
statistically invariant under spatial shifts and $\\{r_{j}=q_{j+1}-q_{j}\\}$
are i.i.d., such that $r_{j}$ is exponentially distributed with density
$P\exp\big{(}-P(x-\mathsfit{a})\big{)}\chi(\\{x\geq\mathsfit{a}\\})$, $P>0$.
Then its first moment $\langle x\rangle=\bar{\rho}=P/(1+\mathsfit{a}P)$ and
$\bar{\rho}\mathsfit{a}<1$. Close packing is reached in the limit $P\to\infty$
with $\bar{\rho}=1/\mathsfit{a}$. In the sense of thermodynamics,
$\bar{\rho},h$ are extensive parameters, while $P,V$ are intensive parameters
upon setting $h=\mathrm{e}^{-V}/\langle\mathrm{e}^{-V}\rangle$. Slightly
deviating from our standard convention the GGE average is denoted by
$\langle\cdot\rangle_{\bar{\rho},h}$. The required GGE averages are easily
computed. For the fields one obtains

$\langle Q^{[\phi]}(0)\rangle_{\bar{\rho},h}=\bar{\rho}\langle h\phi\rangle.$
(5.5)

For the currents one uses that

$\int_{\mathbb{R}}\mathrm{d}v\int_{\mathbb{R}}\mathrm{d}wh(v)h(w)\big{(}\phi(w)-\phi(v)\big{)}\chi(\\{v<w\\})(w-v)=\int_{\mathbb{R}}\mathrm{d}v\phi(v)(v-u)h(v)$
(5.6)

with the mean velocity $u=\langle vh(v)\rangle$. Then

$\langle J^{[\phi]}(0)\rangle_{\bar{\rho},h}=\bar{\rho}\,\langle
hv^{\mathrm{eff}}\phi\rangle$ (5.7)

with the effective velocity

$v^{\mathrm{eff}}(v)=\frac{1}{1-\mathsfit{a}\bar{\rho}}(v-\mathsfit{a}\bar{\rho}u).$
(5.8)

The hydrodynamic evolution equations are built from the average conserved
fields and currents, i.e. the input (5.5) and (5.7). For the initial state one
assumes to have locally the statistics of some GGE. In our case, one example
would be to maintain independence but let $f=\bar{\rho}h$ change very slowly
on the scale set by the maximal correlation length. On a more formal level one
introduces the dimensionless parameter $\epsilon$, $\epsilon\ll 1$, of a
typical scale for the ratio “microscopic length / macroscopic length”.
Microscopically the physically relevant length is the maximal correlation
length of the initial GGE. Translated to the initial $f$, the scale separation
is ensured through the functional form $f(\epsilon x;v)$. Given such random
initial conditions, the hard rod fluid is evolving through free motion and
collisions. Independence of momenta is lost immediately. Yet on the
macroscopic scale nothing is moving, at least not for short times. Due to the
conservation laws, changes can be observed only on microscopic times of order
$\epsilon^{-1}t$. Hydrodynamic scaling is ballistic, since changes in space
and in time are both $\mathcal{O}(\epsilon^{-1})$ in microscopic units. GHD
then postulates that even for such long times local GGE remains approximately
valid. If propagation of local GGE holds, one can write down the equations
governing the evolution of the parameters characterizing local GGEs.

The starting point is the exact microscopic conservation law (5.4). Since this
law holds for any dynamical trajectory, one can average over an arbitrary
ensemble, in particular the initial slowly varying state. Instead of averaging
the time $t$ fields over initial conditions, we average the physical fields
over the statistical state at time $t$, average denoted by
$\langle\cdot\rangle_{t,\epsilon}$ with the subscript $\epsilon$ indicating
the parameter of initial slow variation,

$\partial_{t}\langle Q^{[\phi]}(x)\rangle_{t,\epsilon}+\partial_{x}\langle
J^{[\phi]}(x)\rangle_{t,\epsilon}=0,$ (5.9)

which is still exact. Now, switching to macroscopic spacetime
$(\epsilon^{-1}x,\epsilon^{-1}t)$ propagation of local GGE amounts to the
approximation

$\langle
Q^{[\phi]}(\epsilon^{-1}x)\rangle_{\epsilon^{-1}t,\epsilon}\simeq\langle
Q^{[\phi]}(0)\rangle_{f(x,t)},\qquad\langle
J^{[\phi]}(\epsilon^{-1}x)\rangle_{\epsilon^{-1}t,\epsilon}\simeq\langle
J^{[\phi]}(0)\rangle_{f(x,t)},$ (5.10)

where on the right hand side of each equation $Q^{[\phi]}(x)$ is replaced by
$Q^{[\phi]}(0)$, resp. $J^{[\phi]}(x)$ by $J^{[\phi]}(0)$, because of the
spatial translation invariance of the GGE average
$\langle\cdot\rangle_{f(x,t)}$. Using the effective velocity (5.7) evaluated
at $(x,t)$, we arrive at a closed equation for $f(t)=\bar{\rho}(t)h(t)$,

$\partial_{t}f(x,t;v)+\partial_{x}\big{(}(1-\mathsfit{a}\bar{\rho}(x,t))^{-1}(v-\mathsfit{a}\bar{\rho}(x,t)u(x,t))f(x,t;v)\big{)}=0.$
(5.11)

Here the argument of $f$ indicates that $(x,t)$ is the spacetime patch under
consideration and $v$ is the label of the conserved field in that patch. Eq.
(5.11) is the hydrodynamic equation of the hard rod fluid. Through
nonlinearity the hydrodynamic fields are coupled and decouple only in the
ideal gas limit $\mathsfit{a}\to 0$. For $\mathsfit{a}\geq 0$ the function
$f(t)$ counts the number of particles in the macroscopic volume element
$\mathrm{d}x\mathrm{d}v$ at time $t$. But for $\mathsfit{a}<0$, the density
$\bar{\rho}(t)$ can have negative values, and so does $f(t)$.

On purpose we used the symbol $f$ to emphasize that this quantity is very
close to the Boltzmann $f$-function in the sense of counting the number of
particles in one-particle phase space. However the mathematical structure of
the dynamical evolution is distinct from kinetic theory. On the Euler scale
the Boltzmann collision term is not yet visible and Eq. (5.11) has some
similarity with the Vlasov equation. In that case, the current at $(x,t)$
results from the force when averaged over particles, which is linear in $f$
and nonlocal in space. On the other hand for the hydrodynamic equation (5.11),
the current has a more complicated nonlinearity and is strictly local in
$(x,t)$.

As an apparently general feature of generalized hydrodynamics, the
hydrodynamic equation (5.11) can be written in a quasilinear form. For hard
rods the transformation reads

$\rho_{\mathsf{n}}(v)=\frac{1}{\nu-\mathsfit{a}}f(v),\quad\nu\bar{\rho}=1,$
(5.12)

and

$\partial_{t}\rho_{\mathsf{n}}(x,t;v)+v^{\mathrm{eff}}(x,t;v)\partial_{x}\rho_{\mathsf{n}}(x,t;v)=0.$
(5.13)

The fields $\rho_{\mathsf{n}}(v)$ are the normal modes of the hard rod fluid
equations and the convective equation (5.13) identifies their local
propagation velocities as $v^{\mathrm{eff}}(v)$, which depend nonlinearly on
$\rho_{\mathsf{n}}$ through (5.8).

Hydrodynamic equations are expected to become exact in the limit $\epsilon\to
0$. But in practice, hydrodynamics becomes applicable much before the actual
limit. The precise value of $\epsilon$ depends on the particular physical
system and its initial conditions. The range of validity for hydrodynamics is
mostly a qualitative rule of thumb, rather than a sharp error

estimate.
$\blackdiamond\blackdiamond$ Hydrodynamics of simple fluids, the issue of
mathematical rigor. In 1757 Leonhard Euler published his monumental treatise
on fluids, which in particular included the compressible Euler equations. In
this case the locally conserved fields are number, momentum, and total energy.
Euler argued for the appropriate form of the currents on a phenomenological
basis. It was a triumph of early Statistical Mechanics to understand that
these currents can also be obtained from averaging the microscopic currents
over a Gibbs ensemble. The required technical tool is an identity known as
virial theorem. Since the resulting equations agreed with what had been known
already for a long time, the discussion was mostly confined to the textbook
level. Still, the derivation was a further compelling evidence for the
universal validity of Statistical Mechanics. For integrable many-body systems
one cannot rely on phenomenological evidence. The only available method is to
figure out the GGE averaged currents as based on the microscopic model.

A much harder question is the propagation of local GGEs. For the hard rod
fluid Roland Dobrushin and collaborators proved such a result in the mid
80ies. But for more complicated integrable systems, in particular the Toda
lattice, no results are available yet. Due to their physical importance, for
simple fluids there have been many attempts to establish the compressible
Euler equations in a mathematically convincing way. From that perspective the
best result is still the work of Olla, Varadhan, and Yau. They investigated
the kind of dynamical mixing which would be needed for propagation of local
equilibrium. In fact they modified the mechanical model by assuming that the
deterministic collisions between particles are substituted by random ones,
still respecting the mechanical conservation laws. Thus the macroscopic
equations are the compressible version of Euler equations. For such stochastic
dynamics, Olla, Varadhan, and Yau prove the required mixing. But for a
mechanical system mixing would have to come from deterministic chaos. The
resulting gap is so huge that progress is unlikely, even allowing for a long
term perspective. $\blackdiamond\blackdiamond$

###  5.2 Hard rod lattice

Physically the fluid picture can be more easily visualized. But as in the case
of the Toda chain, one can switch to the lattice version. Then the dynamics is
formulated in terms of the stretches $r_{j}$, where $r_{j}(t)\geq\mathsfit{a}$
because of the hard core. If $p_{j}(t)-p_{j+1}(t)<0$, then $r_{j}(t)$ is
decreasing and reaches collision at $r_{j}(t_{\mathrm{c}})=\mathsfit{a}$,
$t_{\mathrm{c}}$ the collision time. Instantaneously the momenta
$(p_{j},p_{j+1})$ are exchanged and $r_{j}(t)$ is increasing for
$t>t_{\mathrm{c}}$. Such local collisions happen for every spatial patch
throughout time. Obviously, the densities of the conserved fields are

$Q^{[0]}_{j}=r_{j},\qquad Q^{[\phi]}_{j}=\phi(p_{j}).$ (5.14)

Note that $\phi=1$ is just the constant function, hence not admissible. Thus,
in contrast to the fluid, the index $0$ will have to be treated differently
from $n\geq 1$. The currents are

$J^{[0]}_{j}=-p_{j}$ (5.15)

and

$J^{[\phi]}_{j}=\delta(r_{j-1}-\mathsfit{a})(\phi(p_{j-1})-\phi(p_{j}))\chi(\\{p_{j}<p_{j-1}\\})(p_{j-1}-p_{j}).$
(5.16)

Compared to (5.3), there is no flow term and the factor
$(\theta(q_{j}-x)-\theta(q_{j-1}-x))$ is missing.

In a hard rod lattice GGE, the velocities are i.i.d. with probability density
$h(v)$, $h\geq 0$, $\langle h\rangle=1$ and the $r_{j}$’s are also i.i.d. with
probability density
$P\exp\big{(}-P(x-\mathsfit{a})\big{)}\chi(\\{x\geq\mathsfit{a}\\})$, $P>0$.
The average stretch is then $\nu=P^{-1}+\mathsfit{a}$. In particular,
$\nu-\mathsfit{a}>0$. Now $\nu,h$, are the extensive parameters, while $P,V$
are still intensive parameters. The GGE average is denoted by
$\langle\cdot\rangle_{\nu,h}$. The averaged conserved fields are given by

$\langle Q^{[0]}_{0}\rangle_{\nu,h}=\nu,\qquad\langle
Q^{[\phi]}_{0}\rangle_{\nu,h}=\langle h\phi\rangle,$ (5.17)

and the averaged currents by

$\displaystyle\langle J^{[0]}_{0}\rangle_{\nu,h}=-\langle hv\rangle=-u,$
$\displaystyle\langle
J^{[\phi]}_{0}\rangle_{\nu,h}=(\nu-\mathsfit{a})^{-1}\big{(}\langle
hv\phi\rangle-u\langle h\phi\rangle\big{)}=\nu^{-1}\langle
h(v^{\mathrm{eff}}-u)\phi\rangle,$ (5.18)

where we used the same computation as in (5.6) and in Eq. (5.8) rewritten as
$v^{\mathrm{eff}}(v)=(\nu-\mathsfit{a})^{-1}(\nu v-\mathsfit{a}u)$. As a
result the hydrodynamic equations of the hard rod lattice take the form

$\partial_{t}\nu(x,t)-\partial_{x}u(x,t)=0,\quad\partial_{t}h(x,t;v)+\partial_{x}\big{(}\nu(x,t)^{-1}(v^{\mathrm{eff}}(x,t;v)-u(x,t))h(x,t;v)\big{)}=0.$
(5.19)

Compared to (5.11), there are now two equations rather than one. But switching
to normal modes through

$\rho_{\mathsf{n}}(v)=\frac{1}{\nu-\mathsfit{a}}h,$ (5.20)

one arrives again at the quasilinear equation (5.13).

###  5.3 TBA, collision rate ansatz

A central model-dependent feature of TBA is the two-particle scattering shift
as appearing in the definition of the $T$ operator. Now for hard rods the two-
particle scattering shift equals $-\mathsfit{a}$. To have the same structure
as the Toda lattice would mean to replace $T$ in (3.50) by the projection
operator

$T=-\mathsfit{a}|1\rangle\langle 1|,$ (5.21)

in other words $Tf(w)=-\mathsfit{a}\,\langle f\rangle$. With this input the
TBA equation for the hard rod lattice reads

$V-\mu+\mathsfit{a}\,\langle\rho_{\mathsf{n}}\rangle+\log\rho_{\mathsf{n}}=0,\quad\langle\rho_{\mathsf{n}}\rangle=P.$
(5.22)

Since $T$ is such a simple operator, the solution is still explicit with

$\mu=\log P+\mathsfit{a}P-\log\langle\mathrm{e}^{-V}\rangle$ (5.23)

and

$h=\mathrm{e}^{-V}/\langle\mathrm{e}^{-V}\rangle,\quad\rho_{\mathsf{n}}=(\nu-\mathsfit{a})^{-1}h,\quad\rho_{\mathsf{p}}=\nu^{-1}h.$
(5.24)

The average of the conserved fields is determined by

$\langle Q_{0}^{[\phi]}\rangle_{\nu,h}=\langle
h\phi\rangle=\nu\langle\rho_{\mathsf{p}}\phi\rangle,$ (5.25)

in agreement with (5.17). Indeed, with adjusted $T$, the TBA formalism holds
also for hard rods.

Qualitatively, the Toda lattice is best approximated by choosing
$\mathsfit{a}<0$. The identities in (5.25) then provide a simple illustration
for the uniqueness Insert below Eq. (3.62). The stretch
$\nu(P)=(1/P)+\mathsfit{a}$ is monotone decreasing with
$P_{\mathrm{c}}=-1/\mathsfit{a}$. Since $\nu-\mathsfit{a}>0$, in addition
$\rho_{\mathsf{n}}\geq 0$. However $\rho_{\mathsf{p}}$ globally flips sign at
$P_{\mathrm{c}}$, while $\nu\rho_{\mathsf{p}}\geq 0$.

We now invert the logic and try to guess a formula for the GGE averaged
currents of hard rods, through using a conventional GHD argument. Employed is
the notion of quasiparticles, most easily visualized in the fluid picture. A
quasiparticle retains its velocity when undergoing a collision. While a
physical particle rattles back and forth between its two neighbors, the
quasiparticle moves with constant velocity interrupted by jumps of size
$\pm\mathsfit{a}$ at collisions. We now prepare a GGE characterized by
$\bar{\rho},h$ and initially start a tracer quasiparticle at the origin with
velocity $v$. Viewed on a somewhat coarser scale, when time-averaging over
many collisions, the quasiparticle travels with a yet to be determined
effective velocity $v^{\mathrm{eff}}(v)$, which turns out to be increasing in
$v$. The tracer quasiparticle collides with a fluid quasiparticle of velocity
$w$. Such fluid particles have density $\rho_{\mathsf{p}}(w)$. If the
collision is from the left, $v<w$, then the tracer jumps by $-\mathsfit{a}$.
Hence the collision rate is
$\rho_{\mathsf{p}}(w)(v^{\mathrm{eff}}(w)-v^{\mathrm{eff}}(v))$. On the other
hand, if the collision is from the right, $v>w$, then the tracer jumps by
$\mathsfit{a}$ and the collision rate is
$\rho_{\mathsf{p}}(w)(v^{\mathrm{eff}}(v)-v^{\mathrm{eff}}(w))$. When
integrating over all fluid quasiparticles, the bare tracer quasiparticle
velocity is modified through collisions according to

$v^{\mathrm{eff}}(v)=v-\mathsfit{a}\int_{v}^{\infty}\mathrm{d}w\rho_{\mathsf{p}}(w)\big{(}v^{\mathrm{eff}}(w)-v^{\mathrm{eff}}(v)\big{)}+\mathsfit{a}\int_{-\infty}^{v}\mathrm{d}w\rho_{\mathsf{p}}(w)\big{(}v^{\mathrm{eff}}(v)-v^{\mathrm{eff}}(w)\big{)}$
(5.26)

and thus

$v^{\mathrm{eff}}(v)=v-\mathsfit{a}\int_{\mathbb{R}}\mathrm{d}w\rho_{\mathsf{p}}(w)\big{(}v^{\mathrm{eff}}(w)-v^{\mathrm{eff}}(v)\big{)}.$
(5.27)

Equation (5.27) is the collision rate ansatz. More concisely in operator form,

$v^{\mathrm{eff}}(v)=v+T(\rho_{\mathsf{p}}v^{\mathrm{eff}})(v)-(T\rho_{\mathsf{p}})(v)v^{\mathrm{eff}}(v).$
(5.28)

For hard rods the rate equation is easily solved with the result

$v^{\mathrm{eff}}(v)=v+\frac{\mathsfit{a}\bar{\rho}}{1-\mathsfit{a}\bar{\rho}}\int_{\mathbb{R}}\mathrm{d}wh(w)(v-w)=(1-\mathsfit{a}\bar{\rho})^{-1}(v-\mathsfit{a}\bar{\rho}u).$
(5.29)

Indeed, the average current equals $\bar{\rho}hv^{\mathrm{eff}}$, in agreement
with the microscopic computation (5.8).

With hindsight, the expression (5.27) for the effective velocity carries
already a recipe of how to extend to other integrable models. The two-particle
scattering shift for the Toda fluid equals $2\log|v-v^{\prime}|$. Thus one
might hope to obtain the correct average currents for the Toda fluid by
considering the motion of a tracer quasiparticle. Analytically this amounts to
substituting the operator $T$ in (5.28) by the operator $T$ from (3.55).

Notes and references

Section 5.1

The hydrodynamic limit for hard rods, i.e. the asymptotic validity of (5.11),
is proved by Boldrighini et al. [25] under the assumption of a sufficiently
random initial state. The limit holds with probability one. A short account is
provided by Spohn [266]. A different proof is constructed by Ferrari et al.
[112]. From the perspective of generalized hydrodynamics, hard rods are
discussed by Doyon and Spohn [93] and by Doyon [89]. In experiments on
integrable many-body systems, particles are moving in a one-dimensional tube
and are confined by an external potential along the tube. Generically this
breaks integrability. To elucidate more precisely the underlying mechanism,
hard rods serve as an instructive model in the contribution by Cao et al.
[49].

For a many-particle hamiltonian systems with short range interactions, C.B.
Morrey launched an interesting early attempt, Morrey [213]. The modification
of adding randomness in collisions is studied by Olla et al. [223]. The
hydrodynamic limit on the Euler scale can be viewed as a stability result, in
the sense that on the ballistic scale there is no deviation yet from local
equilibrium, resp. local GGE.

Starting in the early 1980ies it was realized that the hydrodynamic limit is
also a challenging problem for stochastic lattice gases and interacting
diffusions. The dynamics is overdamped in the sense that particles have only a
position and evolve through random update rules. An early account is Spohn
[266].

##  6 Equations of generalized hydrodynamics

###  6.1 Average currents

The microscopic version of the currents has been obtained already, see (2.25)
and (2.27). Now the goal is to compute their GGE average in the limit
$N\to\infty$. As before the lattice size is $N$ and we use the pressure
ensemble (3.6). Then in terms of the spectral resolution of $L_{N}$ the $n$-th
microscopic current reads

$J^{[n],N}=\mathrm{tr}\big{[}(L_{N})^{n}L_{N}^{\scriptscriptstyle\downarrow}\big{]}=N\langle\rho_{\mathrm{J},N}\varsigma_{n}\rangle,$
(6.1)

$n\geq 1$, with

$\rho_{\mathrm{J},N}(w)=\frac{1}{N}\sum_{j=1}^{N}\delta(w-\lambda_{j})\Big{(}\sum_{i=1}^{N}a_{i}\psi_{j}(i)\psi_{j}(i+1)\Big{)}.$
(6.2)

The $\delta$-peaks of the Lax density of states are weighted by coefficients
with arbitrary sign. In analogy, not to duplicate names, we call
$\rho_{\mathrm{J},N}$ the empirical current DOS. For given GGE the current DOS
is self-averaging and has the deterministic limit $\rho_{\mathrm{J}}$. In
particular

$\lim_{N\to\infty}\frac{1}{N}\langle J^{[n],N}\rangle_{P,V,N}=\langle
J^{[n]}_{0}\rangle_{P,V}=\langle\rho_{\mathrm{J}}\varsigma_{n}\rangle.$ (6.3)

Setting $n=0$ in (6.1) right hand side, one concludes
$\langle\rho_{\mathrm{J}}\rangle=0$. Hence, as physically expected,
$\rho_{\mathrm{J}}$ cannot have a definite sign. From the few numerical
simulations available, very qualitatively $\rho_{\mathrm{J}}(w)$ has the shape
of $-\partial_{w}\rho_{\mathrm{Q}}(w)$.

Since the Dumitriu-Edelman identity worked so well for the conserved fields,
one is tempted to use the same strategy for the currents. But now not only the
distribution of eigenvalues is in demand. While we have an explicit form for
the joint distribution of $\\{\lambda_{j},\psi_{j}(1),j=1,\ldots,N\\}$, to
apply the map $\Phi$ of Dumitriu-Edelman to the weights
$\sum_{i=1}^{N}a_{i}\psi_{j}(i)\psi_{j}(i+1)$ seems to be a complicated
enterprise. A new approach is needed. Unexpectedly, the key idea will come
from the susceptibility matrices.

We start with the fields and define the infinite volume field-field correlator

$C_{m,n}(j-i)=\langle Q_{j}^{[m]}Q_{i}^{[n]}\rangle_{P,V}^{\mathrm{c}}$ (6.4)

for $m,n\geq 0$, where the superscript c denotes truncation, respectively
connected correlation, $\langle gf\rangle^{\mathrm{c}}_{P,V}=\langle
gf\rangle_{P,V}-\langle g\rangle_{P,V}\langle f\rangle_{P,V}$. Truncated
correlations decay rapidly to zero and the field-field susceptibility matrix
is given by

$C_{m,n}=\sum_{j\in\mathbb{Z}}C_{m,n}(j)=\langle
Q^{[m]};Q^{[n]}\rangle_{P,V},$ (6.5)

where the second equality is merely a convenient notation. $C_{m,n}$ is the
matrix of second derivatives of the generalized free energy. Correspondingly
we introduce the field-current correlator and the field-current susceptibility
matrix

$B_{m,n}(j-i)=\langle J_{j}^{[m]}Q_{i}^{[n]}\rangle_{P,V}^{\mathrm{c}},\quad
B_{m,n}=\sum_{j\in\mathbb{Z}}B_{m,n}(j).$ (6.6)

Despite its apparently asymmetric definition, $B$ satisfies

$B_{m,n}(j)=B_{n,m}(-j).$ (6.7)

To prove, we employ the conservation law and spacetime stationarity to arrive
at

$\displaystyle\partial_{j}\langle
J_{j}^{[m]}(t)Q_{0}^{[n]}(0)\rangle_{P,V}^{\mathrm{c}}=-\partial_{t}\langle
Q_{j}^{[m]}(t)Q_{0}^{[n]}(0)\rangle_{P,V}^{\mathrm{c}}$ $\displaystyle\hskip
37.0pt=-\partial_{t}\langle
Q_{0}^{[m]}(0)Q_{-j}^{[n]}(-t)\rangle_{P,V}^{\mathrm{c}}=\partial_{j}\langle
Q_{0}^{[m]}(0)J_{-j}^{[n]}(-t)\rangle_{P,V}^{\mathrm{c}},$ (6.8)

denoting the difference operator by $\partial_{j}f(j)=f(j+1)-f(j)$. Setting
$t=0$, the difference $\langle
J_{j}^{[m]}Q_{0}^{[n]}\rangle_{P,V}^{\mathrm{c}}-\langle
J_{-j}^{[n]}Q_{0}^{[m]}\rangle_{P,V}^{\mathrm{c}}$ is constant in $j$. Since
truncated correlations decay to zero, this constant has to vanish, which
yields (6.7). In particular, the field-current susceptibility matrix is
symmetric,

$B_{m,n}=B_{n,m}.$ (6.9)

Using this symmetry and restricting to $n\geq 1$, we consider the
$P$-derivative of the average current

$\partial_{P}\langle(L^{n}L^{\scriptscriptstyle\downarrow})_{0,0}\rangle_{P,V}=-B_{n,0}=-B_{0,n}=\langle
Q^{[1]};Q^{[n]}\rangle_{P,V},$ (6.10)

since $J^{[0]}=-Q^{[1]}$ by (2.27). We easily arrived at a very surprising
identity. The $P$-derivative of the average current equals a particular matrix
element of the field-field susceptibility. Such susceptibility depends only on
the eigenvalues of the Lax matrix. In fact, this quantity has been studied
already in Section 4.1 and one only has to adjust the results from there,
setting $\alpha=P$.

As in the case of the free energy, since the pressure is varying as $1/N$, the
fluctuation covariance is adding up, resulting in

$\langle\varsigma_{1},C^{\sharp}\varsigma_{n}\rangle=\int_{0}^{1}\mathrm{d}u\langle
Q^{[1]};Q^{[n]}\rangle_{uP,V}$ (6.11)

with the operator $C^{\sharp}$ defined in (4.23). Inserting from (6.10), one
concludes that

$\partial_{P}\big{(}\langle
J^{[n]}_{0}\rangle_{P,V}-P\langle\varsigma_{1},C^{\sharp}\varsigma_{n}\rangle\big{)}=0.$
(6.12)

The term in the round brackets has to be independent of $P$, in particular
equal to its value at $P=0$. Since the covariance is bounded, the second
summand vanishes at $P=0$. For $\langle J^{[n]}_{0}\rangle_{P,V}$, we use the
random walk expansion (2.18). Each walk contains at least one factor of
$a_{j}$ and thus vanishes because under the a priori measure $a_{j}\to 0$ in
the limit $P\to 0$. Therefore

$\langle
J^{[n]}_{0}\rangle_{P,V}=P\langle\varsigma_{1},C^{\sharp}\varsigma_{n}\rangle.$
(6.13)

Inserting from (4.23) and substituting as
$\alpha\rho_{\mathfrak{s}}=\rho_{\mathsf{n}}$,

$\langle
J^{[n]}_{0}\rangle_{P,V}=\langle\varsigma_{1},(1-\rho_{\mathsf{n}}T)^{-1}\rho_{\mathsf{n}}\varsigma_{n}\rangle-\nu\langle\varsigma_{1}(1-\rho_{\mathsf{n}}T)^{-1}\rho_{\mathsf{n}}\rangle\langle\varsigma_{n}(1-\rho_{\mathsf{n}}T)^{-1}\rho_{\mathsf{n}}\rangle.$
(6.14)

Since $(1-\rho_{\mathsf{n}}T)^{-1}\rho_{\mathsf{n}}$ is a symmetric operator,
one arrives at

$\langle
J^{[n]}_{0}\rangle_{P,V}=\langle\rho_{\mathsf{n}}\varsigma_{1}^{\mathrm{dr}}\varsigma_{n}\rangle-
q_{1}\langle\rho_{\mathsf{p}}\varsigma_{n}\rangle,\qquad
q_{1}=\nu\langle\rho_{\mathsf{p}}\varsigma_{1}\rangle.$ (6.15)

Using (6.3), one finally concludes

$\langle\rho_{\mathrm{J}}\varsigma_{n}\rangle=\langle(\rho_{\mathsf{n}}\varsigma_{1}^{\mathrm{dr}}-q_{1}\rho_{\mathsf{p}})\varsigma_{n}\rangle$
(6.16)

valid for $n\geq 1$ and, in addition,

$\langle\rho_{\mathrm{J}}\varsigma_{0}\rangle=-q_{1}.$ (6.17)

In the previous chapter we entertained the possibility that, by analogy, the
Toda current can be written as

$\rho_{\mathrm{J}}=\rho_{\mathsf{p}}(v^{\mathrm{eff}}-q_{1}),$ (6.18)

compare with (5.2) and (5.24). Using the scattering shift of the Toda lattice,
this proposal amounts to

$v^{\mathrm{eff}}(v)=v+2\int_{\mathbb{R}}\mathrm{d}w\log|v-w|\rho_{\mathsf{p}}(w)\big{(}v^{\mathrm{eff}}(w)-v^{\mathrm{eff}}(v)\big{)}.$
(6.19)

If such an analogy holds, then
$\rho_{\mathsf{p}}v^{\mathrm{eff}}=\rho_{\mathsf{n}}\varsigma_{1}^{\mathrm{dr}}$
and hence

$v^{\mathrm{eff}}=\frac{\varsigma_{1}^{\mathrm{dr}}}{\varsigma_{0}^{\mathrm{dr}}}\,.$
(6.20)

To confirm, we start from

$\displaystyle
T(\rho_{\mathsf{p}}v^{\mathrm{eff}})-(T\rho_{\mathsf{p}})v^{\mathrm{eff}}=T(\rho_{\mathsf{n}}\varsigma_{1}^{\mathrm{dr}})-(T\rho_{\mathsf{p}})v^{\mathrm{eff}}=(T\rho_{\mathsf{n}}-1+1)(1-T\rho_{\mathsf{n}})^{-1}\varsigma_{1}-(T\rho_{\mathsf{p}})v^{\mathrm{eff}}$
$\displaystyle\hskip
68.0pt=-\varsigma_{1}+\varsigma_{1}^{\mathrm{dr}}-\big{(}\rho_{\mathsf{n}}^{-1}\rho_{\mathsf{p}}-1\big{)}v^{\mathrm{eff}}=v^{\mathrm{eff}}-\varsigma_{1},$
(6.21)

as claimed.

At low pressure the Toda particles are far apart and interact mostly through
isolated two-particle collisions. Quasiparticles are then defined as for hard
rods, namely through maintaining their velocity upon excluding the time span
for collisions. The tracer quasiparticle jumps by a distance regulated by the
Toda two-particle scattering shift. The validity of the collision rate ansatz
(6.19) can be argued as we did already for a hard rod fluid. However, for a
dense Toda fluid with $P$ not too close to $0$ the notion of a quasiparticle
is fuzzy, even more so the meaning of a tracer quasiparticle. The validity of
(6.19) for arbitrary GGE parameters will be concluded only indirectly.

Conceptually the collision rate ansatz can be read in a somewhat different
way, as inspired by the time-of-flight method employed in recent experiments
on the rapidity distribution for the $\delta$-Bose gas. One prepares the
$N$-particle Toda system in a GGE and runs the dynamics over the time span
$[0,t]$ with periodic boundary conditions, $t$ of the same order as $N$. To
identify the tracer quasiparticle, a past and future diagnostic step is added,
to say in the time interval $[-\tau,0]$ backwards and $[t,t+\tau]$ forwards
the dynamics is executed with open boundary conditions, $\tau>0$ and of order
1. Then at times $-\tau,t+\tau$ the particles are well separated and have
definite velocities. At time $-\tau$, we pick one particle with velocity $w$
and call it tracer quasiparticle. Then at time $t+\tau$, there is exactly one
particle with the same velocity $w$. For $\tau\ll t$ the displacement of the
tracer quasiparticle should be $v^{\mathrm{eff}}(w)t$ in good approximation.

###  6.2 Hydrodynamic equations

On the hydrodynamic scale the local GGE is characterized by the stretch $\nu$
and the Lax DOS $\nu\rho_{\mathsf{p}}$, both of which now become spacetime
dependent. Merely inserting the average currents, one arrives at the Euler
type hydrodynamic evolution equations,

$\displaystyle\partial_{t}\nu(x,t)-\partial_{x}q_{1}(x,t)=0,$
$\displaystyle\partial_{t}\langle\nu(x,t)\rho_{\mathsf{p}}(x,t)\varsigma_{n}\rangle+\partial_{x}\langle(v^{\mathrm{eff}}(x,t)-q_{1}(x,t))\rho_{\mathsf{p}}(x,t)\varsigma_{n}\rangle=0$
(6.22)

for $n\geq 1$. The latter equation extends to $n=0$, since the bracket under
$\partial_{t}$ equals $1$ and, as can be concluded from (6.1), under
$\partial_{x}$ equals $0$. Therefore one can switch to the pointwise version
of the hydrodynamic equations,

$\displaystyle\partial_{t}\nu(x,t)-\partial_{x}q_{1}(x,t)=0,$
$\displaystyle\partial_{t}\big{(}\nu(x,t)\rho_{\mathsf{p}}(x,t;v)\big{)}+\partial_{x}\big{(}(v^{\mathrm{eff}}(x,t;v)-q_{1}(x,t))\rho_{\mathsf{p}}(x,t;v)\big{)}=0.$
(6.23)

As a most remarkable feature of generalized hydrodynamics, these equations can
be transformed explicitly to a quasilinear form. For this purpose we recall
the identity (3.61), which expresses $\rho_{\mathsf{n}}$ in terms of
$\rho_{\mathsf{p}}$,

$\rho_{\mathsf{n}}=\rho_{\mathsf{p}}(1+(T\rho_{\mathsf{p}}))^{-1}.$ (6.24)

Then Eq. (6.2) acquires the normal form

$\partial_{t}\rho_{\mathsf{n}}+\nu^{-1}(v^{\mathrm{eff}}-q_{1})\partial_{x}\rho_{\mathsf{n}}=0.$
(6.25)

Thus the linearization operator is in fact merely a multiplication by
$\nu^{-1}(v^{\mathrm{eff}}-q_{1})$, in other words the operator is diagonal.
For nonintegrable chains the corresponding operator is a $3\times 3$ matrix,
which in general is not diagonal. Its diagonalization would yield the normal
modes. The solution of the respective hyperbolic conservation laws can develop
shock discontinuities, since locally the solution may have to jump to another
eigenvalue. In our case $\nu^{-1}(v^{\mathrm{eff}}-q_{1})$ is expected to be a
smooth function. Thus the spectrum of the linearization operator is continuous
and has no eigenvalues. Hence for smooth initial data solutions of (6.25) are
expected to stay smooth. No mathematical analysis has been attempted to verify
this conjecture. In fact, the situation could be more subtle since $\nu$ can
take either sign and $\nu^{-1}$ might be singular.

To verify (6.25) we start from

$\rho_{\mathsf{p}}\partial_{t}\nu+\nu\partial_{t}\rho_{\mathsf{p}}+\partial_{x}\big{(}(v^{\mathrm{eff}}-q_{1})\rho_{\mathsf{p}}\big{)}=0,$
(6.26)

which together with the continuity equation yields

$\nu\partial_{t}\rho_{\mathsf{p}}+\partial_{x}(v^{\mathrm{eff}}\rho_{\mathsf{p}})-q_{1}\partial_{x}\rho_{\mathsf{p}}=0.$
(6.27)

Using this identity and (6.24), one obtains

$\displaystyle\nu\partial_{t}\rho_{\mathsf{n}}=\frac{\rho_{\mathsf{n}}}{\rho_{\mathsf{p}}}\big{(}-\partial_{x}(v^{\mathrm{eff}}\rho_{\mathsf{p}})+q_{1}\partial_{x}\rho_{\mathsf{p}}\big{)}-\frac{\rho_{\mathsf{n}}^{2}}{\rho_{\mathsf{p}}}\big{(}T(-\partial_{x}(v^{\mathrm{eff}}\rho_{\mathsf{p}})+q_{1}\partial_{x}\rho_{\mathsf{p}})\big{)}$
$\displaystyle\hskip
12.0pt=\frac{\rho_{\mathsf{n}}}{\rho_{\mathsf{p}}}\big{(}-\partial_{x}(v^{\mathrm{eff}}\rho_{\mathsf{p}})+\rho_{\mathsf{n}}\partial_{x}T(v^{\mathrm{eff}}\rho_{\mathsf{p}})\big{)}+q_{1}\partial_{x}\rho_{\mathsf{n}}.$
(6.28)

The effective velocity solves the integral equation

$v^{\mathrm{eff}}=v+T(v^{\mathrm{eff}}\rho_{\mathsf{p}})-(T\rho_{\mathsf{p}})v^{\mathrm{eff}}$
(6.29)

and, by inserting (6.24),

$\partial_{x}\big{(}\rho_{\mathsf{n}}^{-1}v^{\mathrm{eff}}\rho_{\mathsf{p}}\big{)}=\partial_{x}T(v^{\mathrm{eff}}\rho_{\mathsf{p}}).$
(6.30)

Hence the first summand in (6.2) reads

$\frac{\rho_{\mathsf{n}}}{\rho_{\mathsf{p}}}\big{(}-\partial_{x}(v^{\mathrm{eff}}\rho_{\mathsf{p}})+\rho_{\mathsf{n}}\partial_{x}(\rho_{\mathsf{n}}^{-1}v^{\mathrm{eff}}\rho_{\mathsf{p}})\big{)}=-v^{\mathrm{eff}}\partial_{x}\rho_{\mathsf{n}},$
(6.31)

thereby confirming (6.25).

The normal form (6.25) of the hydrodynamic equations is a crucial insight. A
priori, the coupled conservation laws might be so complicated that it becomes
an impossible task to extract information of interest. Here the normal form is
of considerable help. It allows to obtain some partially analytic solutions,
as the domain wall problem and the linearized version to be discussed in the
next Chapters. The normal form suggests also how to devise a numerical scheme
for solving the hydrodynamic equations. In its most basic version, at the
current time, $t$, one keeps $\nu^{-1}(v^{\mathrm{eff}}-q_{1})$ fixed and
solves the resulting linear equation for one further time step $\mathrm{d}t$.
With the so obtained $\rho_{\mathsf{n}}(t+\mathrm{d}t)$ one updates
$\nu^{-1}(v^{\mathrm{eff}}-q_{1})$ according to (6.20) for $v^{\mathrm{eff}}$,
$\nu^{-1}=\langle\rho_{\mathsf{p}}\rangle$,
$\nu^{-1}q_{1}=\langle\rho_{\mathsf{p}}\varsigma_{1}\rangle$, and
$\rho_{\mathsf{p}}=\rho_{\mathsf{n}}\varsigma_{0}^{\mathrm{dr}}$. No explicit
use of TBA is involved in this step. The iteration is then repeated many times
until the desired final time is reached.

The collision rate equation suggests a distinct approach, called flea gas
algorithm. The basic idea is to replace the true motion by an instantaneous
exchange of velocities $v_{1},v_{2}$ whenever the two particles first reach
the distance $2\log|v_{1}-v_{2}|$. Thereby the simulation of the particle
dynamics is massively speeded up and it is ensured that on large scales the
Toda generalized hydrodynamics is realized. In fact, the true dynamics is more
complicated, since a third particle or even more particles might interfere
with the collision between particles 1 and 2. To properly set up the particle
model requires a separation into isolated particle clusters combined with an
instantaneous many-body collision. For the Toda fluid, and also integrable
quantum many-body systems, one can use such a scheme as an alternative to
direct numerical simulations of the hydrodynamic equations.

Physically, the initial conditions have to satisfy the constraints resulting
from GGE. This amounts to $\nu\in\mathbb{R}$, $\nu\rho_{\mathsf{p}}\geq 0$,
and $\rho_{\mathsf{n}}\geq 0$. Such conditions should be propagated by the
Euler equations. For $\rho_{\mathsf{n}}$, this property follows from (6.25).
The prefactor of $\partial_{x}$ looks singular at $\nu=0$ and the propagation
of $\rho_{\mathsf{p}}$ properties is less obvious. A specific example is
discussed in Chapter 8.

Comparing Eqs. (5.19) and (6.2) and identifying $h=\nu\rho_{\mathsf{p}}$,
$u=q_{1}$, Toda and hard rods have structurally identical hydrodynamic
evolution equations. On the hard rod side the hydrodynamic limit has been
proved, thereby supporting the validity of the corresponding property for the
Toda lattice.

Notes and references

Section 6.1

The computation leading to (6.1) is taken from De Nardis et al. [75], see also
Karevski and Schütz [169]. The symmetry of the $B$-matrix is discussed in
Grisi and Schütz [142] and Spohn [267]. Rather likely there is earlier work.
The connection between the symmetry of the $B$-matrix and the average current
was first noted by Spohn [271] and extended to quantum systems by Yoshimura
and Spohn [298]. For integrable quantum systems, the appropriately adjusted
rate equation (6.19) has been discovered by Castro-Alvaredo et al. [53]. Very
quickly it was understood that the scheme originally developed in Bertini et
al. [22] and Castro-Alvaredo et al. [53] for specific models has a much wider
applicability, see Doyon et al. [95]. There are further attempts to
microscopically justify the form of the average currents, to mention the ones
by Vu and Yoshimura [295], Borsi et al. [31], Pristyák and Pozsgay [239], and
the review Cubero et al. [65].

Section 6.2

The transformation to normal form is established by Castro-Alvaredo et al.
[53]. The flea gas algorithm is discussed by Doyon et al. [95] and Mestyán and
Alba [207] as an alternative approach to more standard PDE discretization
schemes. A numerically oriented contribution is Møller et al. [212]. In
Bulchandani [38] it is argued that GHD equations are continuum integrable
systems, which is related to the geometric viewpoint developed by Doyon et al.
[94].

##  7 Linearized hydrodynamics and GGE dynamical correlations

In equilibrium statistical mechanics a central theme is to understand the
structure of static correlations, with particular focus on critical points in
thermodynamic parameter space in the neighborhood of which the correlation
length is large on microscopic scales. A natural extension are time-dependent
correlations, which can be viewed as the propagation of initially small
perturbations in the equilibrium state. Now conservation laws and broken
symmetries will play a crucial role. The most elementary approach is the
Landau-Lifshitz theory which uses the link to macroscopic equations linearized
at thermal equilibrium, as will be discussed below. Near critical points the
more sophisticated techniques of critical dynamics would come into play.

###  7.1 Equilibrium spacetime correlations for nonintegrable chains

The general structure behind the Landau-Lifshitz theory can be explained
already in the context of nonintegrable systems with a few conservation laws.
For concreteness we consider nonintegrable anharmonic chains. In essence, the
extension to integrable chains consists in substituting the respective
$3\times 3$ matrices by operators on Hilbert spaces with infinite dimension.
To distinguish from the infinite dimensional case, we use $\vec{u}$ for $3$
-vectors and $\mathsfit{A}$ for $3\times 3$ matrices.

The infinitely extended chain is governed by the hamiltonian

$H_{\mathrm{ch}}=\sum_{j\in\mathbb{Z}}\big{(}\tfrac{1}{2}p^{2}_{j}+V_{\mathrm{ch}}(r_{j})\big{)}$
(7.1)

with equations of motion

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

The chain potential, $V_{\mathrm{ch}}$, is bounded from below and increases at
least one-sided linearly at infinity so to have a finite partition function.
Stretch, momentum, and energy are defined by the local fields

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

and the respective current densities are

$\vec{J}_{j}=\big{(}-p_{j},-V_{\mathrm{ch}}^{\prime}(r_{j-1}),-(p_{j-1}+p_{j})V_{\mathrm{ch}}^{\prime}(r_{j-1})\big{)}.$
(7.4)

It is assumed that there are no further local conservation laws, which is the
essence of non-integrability.

The canonical equilibrium state factorizes with one factor given by

$Z_{0}(P,\bar{u},\beta)^{-1}\exp\\!\big{(}-\beta\big{(}\tfrac{1}{2}(p_{0}-\bar{u})^{2}+V_{\mathrm{ch}}(r_{0})\big{)}-Pr_{0}\big{)},$
(7.5)

where we introduced the intensive dual parameters $P,\bar{u},\beta$.
$\mathfrak{p}=P/\beta$ is the physical pressure, $\bar{u}$ the mean velocity,
and $\beta$ the inverse temperature. The chain free energy per site equals

$F_{\mathrm{ch}}=-\log Z_{0}(P,\bar{u},\beta).$ (7.6)

Of course, more explicit expressions could be provided, but this is not needed
at the moment. In accordance with conventions for the Toda chain, the indices
run over $m,n=0,1,2$, throughout. The static field-field correlator in
infinite volume is defined by

$\mathsfit{C}_{m,n}(j)=\langle
Q_{j}^{[m]}Q_{0}^{[n]}\rangle_{P,\bar{u},\beta}^{\mathrm{c}}=\delta_{0j}\langle
Q_{0}^{[m]}Q_{0}^{[n]}\rangle_{P,\bar{u},\beta}^{\mathrm{c}}$ (7.7)

and the static field-field susceptibility matrix by

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

$\mathsfit{C}$ is the matrix of second derivatives of the chain free energy
$F_{\mathrm{ch}}$. In the same fashion the field-current correlator is given
by

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

and the static field-current susceptibility matrix by

$\mathsfit{B}_{m,n}=\sum_{j\in\mathbb{Z}}\mathsfit{B}_{m,n}(j).$ (7.10)

The field-field spacetime correlator is the matrix

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

which for $t\neq 0$ is no longer $\delta$-correlated in $j$. It is convenient
to also introduce its Fourier transform

$\hat{\mathsfit{S}}_{m,n}(k,t)=\sum_{j\in\mathbb{Z}}\mathrm{e}^{\mathrm{i}kj}\mathsfit{S}_{m,n}(j,t)$
(7.12)

with $k\in[-\pi,\pi]$.

Switching to the macroscopic continuum scale, the average fields and currents
are denoted by
$\vec{\mathsfit{u}}=\langle\vec{Q}_{0}\rangle_{P,\bar{u},\beta}$,
$\vec{\mathsfit{j}}=\langle\vec{J}_{0}\rangle_{P,\bar{u},\beta}$. According to
the rules of thermodynamics there is a one-to-one map between
$\vec{\mathsfit{u}}$ and $(P,\bar{u},\beta)$. Therefore the hydrodynamic
equations for the chain can be written as

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

with quasilinear form

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

The matrix $\mathsfit{A}$ is known as flux Jacobian and given through

$\mathsfit{A}_{m,n}(\vec{\mathsfit{u}})=\partial_{\mathsfit{u}_{n}}\langle
J^{[m]}_{0}\rangle_{\vec{\mathsfit{u}}},$ (7.15)

which by the chain rule yields

$\mathsfit{A}=\mathsfit{B}\mathsfit{C}^{-1},$ (7.16)

viewed as an identity in $(P,\bar{u},\beta)$. The matrix $\mathsfit{A}$ is not
symmetric in general. As established in (6.1), essentially using only
spacetime stationarity, the matrix $\mathsfit{B}$ turns out to be symmetric.
By construction $\mathsfit{C}$ is symmetric and also $\mathsfit{C}>0$. Thus

$\mathsfit{A}=\mathsfit{C}^{\frac{1}{2}}\mathsfit{C}^{-\frac{1}{2}}\mathsfit{B}^{-1}\mathsfit{C}^{-\frac{1}{2}}\mathsfit{C}^{-\frac{1}{2}},$
(7.17)

which implies that $\mathsfit{A}$ has real eigenvalues, $c_{\alpha}$,
$\alpha=0,1,2$, and a complete system of right, $|\psi_{\alpha}\rangle$, and
left eigenvectors, $\langle\tilde{\psi}_{\alpha}|$.

Following Landau and Lifshitz the average (7.11) is viewed as a small initial
perturbation of equilibrium. Then to leading order the correlator may be
approximated by Eq. (7.14) linearized around a constant background with
parameters $P,\bar{u},\beta$. Denoting this perturbation by $\vec{u}$ the
nonlinear hydrodynamic equation turns to its linearized version

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

where the flux Jacobian $\mathsfit{A}$ is now evaluated at the background,
thus constant in spacetime. However the initial conditions, $\vec{u}(x)$, are
random, their average being denoted by $\mathbb{E}$. According to thermal
equilibrium, because of independence, on the macroscopic scale the random
field $\vec{u}(x)$ is Gaussian white noise with mean zero and covariance

$\mathbb{E}\big{(}u_{m}(x)u_{n}(0)\big{)}=\mathsfit{C}_{m,n}\delta(x)=\mathsfit{S}^{\sharp}_{m,n}(x,0),$
(7.19)

which defines the continuum static correlator $\mathsfit{S}^{\sharp}(x,0)$. To
obtain its spacetime version one has to solve (7.18) with initial conditions
(7.19), which yields

$\mathbb{E}\big{(}u_{m}(x,t)u_{n}(0,0)\big{)}=\mathsfit{S}^{\sharp}_{m,n}(x,t),$
(7.20)

thereby defining the right hand side. In Fourier space

$\hat{\mathsfit{S}}^{\sharp}_{m,n}(k,t)=(\mathrm{e}^{-\mathrm{i}k\mathsfit{A}t}\mathsfit{C})_{m,n}.$
(7.21)

Using the eigenvectors of $\mathsfit{A}$, in position space the correlator
(7.20) has the explicit form

$\mathsfit{S}^{\sharp}_{m,n}(x,t)=\sum_{\alpha=0}^{2}\delta(x-c_{\alpha}t)(|\psi_{\alpha}\rangle\langle\tilde{\psi}_{\alpha}|\mathsfit{C})_{m,n}.$
(7.22)

To establish the relation to the microscopic definition (7.11), one introduces
the dimensionless scale parameter $\epsilon$. Then under ballistic scaling it
is expected that for small $\epsilon$,

$\epsilon^{-1}\mathsfit{S}_{m,n}(\lfloor\epsilon^{-1}x\rfloor,\epsilon^{-1}t)\simeq\mathsfit{S}^{\sharp}_{m,n}(x,t).$
(7.23)

Here $\lfloor\cdot\rfloor$ denotes the integer part and the prefactor is
chosen such that the sum rule

$\sum_{j\in\mathbb{Z}}\mathsfit{S}(j,t)=\sum_{j\in\mathbb{Z}}\mathsfit{S}(j,0)$
(7.24)

holds. In Fourier space, more compactly,

$\lim_{\epsilon\to 0}\hat{\mathsfit{S}}_{m,n}(\epsilon
k,\epsilon^{-1}t)=\hat{\mathsfit{S}}^{\sharp}_{m,n}(k,t).$ (7.25)

The eigenvectors of $\mathsfit{A}$ are called normal modes. The $\alpha$-th
mode travels with velocity $c_{\alpha}$ and the initial condition determines
the particular linear combination of normal modes as encoded by the
susceptibility matrix $\mathsfit{C}$. The spacetime correlator has three delta
peaks moving ballistically. If the background has zero average momentum, then
the eigenvalues are $\vec{c}=(-c,0,c)$ with $c$ the isentropic speed of sound.
There are two sound peaks with equal speed moving in opposite directions and a
heat peak standing still. For an integrable system, there are so to speak
infinitely many peaks, each moving with its own velocity. We thus expect the
corresponding time correlator to have a broad spectrum which expands
ballistically.

In Figure 7 shown is molecular dynamics of a fluid with shoulder interaction
potential. Clearly visible is the anomalous broadening of the delta-peaks with
power laws $t^{3/5}$ and $t^{2/3}$.

Figure 7: Heat and sound peak of the equilibrium spacetime correlations for a
fluid interacting through a shoulder potential. Shoulder means
$V_{\mathrm{sh}}(x)=\infty$ for $|x|<\tfrac{1}{2}$, $V_{\mathrm{sh}}(x)=1$ for
$\tfrac{1}{2}<|x|<1$, and $V_{\mathrm{sh}}(x)=0$ for $1<|x|$. Heat peak is on
top and sound peak, displaced by $ct$, below. System size is $N=4096$, time
$t=1024$, $\beta=2$, and $P=1.2$ for heat respectively $P=1$ for sound. The
sound speed is $c=1.74$. There is a single scale-parameter in comparison with
theoretical predictions. The right panels display the corresponding
logarithmic plot. Numerically, the shoulder potential has the simplification
that particle trajectories are piecewise linear. From [201].

###  7.2 GGE spacetime correlations for the Toda lattice

For the Toda chain we follow step by step the road map provided by the finite
mode case, where we recall that $\varsigma_{n}(w)=w^{n}$, including $n=0$. The
field-field and field-current correlator have been introduced before in (6.5),
resp. (6.6). Since we have computed already the GGE average of fields in
(3.62) and of currents in (6.16), the matrix elements of $C,B$ can be
determined by differentiating these averages with respect to $P$ and
$\mu_{n}$. The computation is somewhat lengthy and stated is merely the final
result. For the field-field correlator one obtains

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

$m,n\geq 1$, and for the field-current correlator

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

Both matrices have a two-block structure. Rather than using the discrete
index, $n=1,2,\ldots$, we introduce some real valued function, $\phi$. Our
linear space consists now of two-vectors $(r,\phi)$ with $r$ the coefficient
for the index $0$ and $\phi$ for all other modes. More formally, the
appropriate linear space is
$\mathbb{C}\oplus\big{(}L^{2}(\mathbb{R},\mathrm{d}w)\ominus\\{\varsigma_{0}\\}\big{)}$,
the second summand consisting of all square-integrable functions orthogonal to
$\varsigma_{0}$. We also introduce the linear operator, $F$, in Dirac
notation,

$F=(1-T\rho_{\mathsf{n}})^{-1}-\nu|\varsigma_{0}^{\mathrm{dr}}\rangle\langle\rho_{\mathsf{p}}|.$
(7.28)

Note that $F\varsigma_{0}=0$. Then

$C=\begin{pmatrix}\nu^{3}\langle\rho_{\mathsf{p}}\varsigma_{0}^{\mathrm{dr}}\varsigma_{0}^{\mathrm{dr}}\rangle&-\nu^{2}\big{\langle}\rho_{\mathsf{p}}\varsigma_{0}^{\mathrm{dr}}\big{|}F\\\\[2.15277pt]
-\nu^{2}F^{*}\big{|}\rho_{\mathsf{p}}\varsigma_{0}^{\mathrm{dr}}\big{\rangle}&\nu
F^{*}\rho_{\mathsf{p}}F\end{pmatrix}.$ (7.29)

For some general operator $A$, as introduced before, $A^{*}$ stands for the
adjoint operator of $A$ with respect to the standard $L^{2}$ inner product
$\langle\cdot,\cdot\rangle$. In particular,

$F^{*}=(1-\rho_{\mathsf{n}}T)^{-1}-\nu|\rho_{\mathsf{p}}\rangle\langle\varsigma_{0}^{\mathrm{dr}}|$
(7.30)

With the same notation, the matrix $B$ is given by

$B=\begin{pmatrix}\nu^{2}\langle\rho_{\mathsf{p}}(v^{\mathrm{eff}}-q_{1})\varsigma_{0}^{\mathrm{dr}}\varsigma_{0}^{\mathrm{dr}}\rangle&-\nu\big{\langle}\rho_{\mathsf{p}}(v^{\mathrm{eff}}-q_{1})\varsigma_{0}^{\mathrm{dr}}\big{|}F\\\\[2.15277pt]
-\nu
F^{*}\big{|}\rho_{\mathsf{p}}(v^{\mathrm{eff}}-q_{1})\varsigma_{0}^{\mathrm{dr}}\big{\rangle}&F^{*}\rho_{\mathsf{p}}(v^{\mathrm{eff}}-q_{1})F\end{pmatrix}.$
(7.31)

Following the road map we are supposed to compute $\mathrm{e}^{At}C$ with
$A=BC^{-1}$, which does not seem to be completely straightforward. But instead
one might guess the entire solution $S(x,t)$ by noting that
$S(x,0)=\delta(x)C$ and
$\partial_{t}S(x,t)\big{|}_{t=0}=-\delta^{\prime}(x)B$. Indeed, the two
conditions can be satisfied by setting

$S(x,t)=\begin{pmatrix}\nu^{3}\langle\rho_{\mathsf{p}}\,\delta(x-t\nu^{-1}(v^{\mathrm{eff}}-q_{1}))\varsigma_{0}^{\mathrm{dr}}\varsigma_{0}^{\mathrm{dr}}\rangle&-\nu^{2}\big{\langle}\rho_{\mathsf{p}}\,\delta(x-t\nu^{-1}(v^{\mathrm{eff}}-q_{1}))\varsigma_{0}^{\mathrm{dr}}\big{|}F\\\\[2.15277pt]
-\nu^{2}F^{*}\big{|}\rho_{\mathsf{p}}\,\delta(x-t\nu^{-1}(v^{\mathrm{eff}}-q_{1}))\varsigma_{0}^{\mathrm{dr}}\big{\rangle}&\nu
F^{*}\rho_{\mathsf{p}}\delta(x-t\nu^{-1}(v^{\mathrm{eff}}-q_{1}))F\end{pmatrix}.$
(7.32)

As a control check, the correlator scales indeed ballistically,

$S(x,t)=t^{-1}S(x/t,1).$ (7.33)

Our computation misses that matching simply the value and first derivative at
$t=0$ does not determine the full solution. Leaving for the moment a more
convincing argument aside, let us reflect on the resulting predictions. For
example, within the stated approximations, the spacetime stretch-stretch
correlation function is predicted as

$S_{0,0}(x,t)=\nu^{2}\int_{\mathbb{R}}\mathrm{d}w\delta(x-t\nu^{-1}(v^{\mathrm{eff}}(w)-q_{1}))\nu\rho_{\mathsf{p}}(w)\varsigma_{0}^{\mathrm{dr}}(w)^{2},$
(7.34)

and similarly for the momentum-momentum correlation,

$S_{1,1}(x,t)=\int_{\mathbb{R}}\mathrm{d}w\delta(x-t\nu^{-1}(v^{\mathrm{eff}}(w)-q_{1}))\nu\rho_{\mathsf{p}}(w)(\varsigma_{1}-q_{1}\varsigma_{0})^{\mathrm{dr}}(w)^{2}.$
(7.35)

Both identities can be deduced most easily by going back to (7.26) upon
inserting the $\delta$-function. Our result is in perfect analogy to the
finite mode case. The modes are labelled by the spectral parameter $w$ and
propagate with velocity $\nu^{-1}(v^{\mathrm{eff}}(w)-q_{1})$. The last factor
provides the weights and depends on the particular choice for the matrix
elements of $S(x,t)$.

To confirm the guess, one has to figure out the spectral representation of the
operator $A=BC^{-1}$. As key observation, in (6.24) we already introduced a
nonlinear map which transforms the system of conservation laws into its
quasilinear version in such a way that the operator corresponding to $A$ is
manifestly diagonal. This property is expected to persist when linearizing the
map (6.24). In the $\nu,\rho_{\mathsf{p}}$ variables the map is given by

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

Both sides are linearized as $\rho_{\mathsf{n}}+\epsilon g$,
$\nu\rho_{\mathsf{p}}+\epsilon r$, $\nu+\epsilon\phi$, $\langle\phi\rangle=0$.
To first order in $\epsilon$ this yields the linear map $R:g\mapsto(r,\phi)$
given by

$Rg=\nu\begin{pmatrix}-\nu\langle
g(\varsigma_{0}^{\mathrm{dr}})^{2}\rangle\\\\[2.15277pt]
F^{*}(g\varsigma_{0}^{\mathrm{dr}})\end{pmatrix}$ (7.37)

with $F^{*}$ as in (7.30). Note that indeed $\langle F^{*}\phi\rangle=0$. Eq.
(7.36) can be inverted as

$\rho_{\mathsf{p}}=(1-\rho_{\mathsf{n}}T)^{-1}\rho_{\mathsf{n}},$ (7.38)

thereby deriving, by a similar argument as before,

$R^{-1}\begin{pmatrix}r\\\
\phi\end{pmatrix}=(\nu\varsigma_{0}^{\mathrm{dr}})^{-1}(-\rho_{\mathsf{n}}r+(1-\rho_{\mathsf{n}}T)\phi).$
(7.39)

Indeed one checks that

$RR^{-1}=1,\quad R^{-1}R=1,$ (7.40)

the first $``1"$ standing for the identity operator as a $2\times 2$ block
matrix and the second $``1"$ for the identity operator in the space of scalar
functions.

The operators $C,B$ can be written in the new basis with the result

$R^{-1}CR=\nu\big{|}\rho_{\mathsf{p}}\big{\rangle}\big{\langle}(\varsigma_{0}^{\mathrm{dr}})^{2}\big{|}+\nu(\varsigma_{0}^{\mathrm{dr}})^{-1}\rho_{\mathsf{p}}FF^{*}\varsigma_{0}^{\mathrm{dr}}$
(7.41)

and

$\displaystyle
R^{-1}BR=\big{|}(v^{\mathrm{eff}}-q_{1})\rho_{\mathsf{p}}\big{\rangle}\big{\langle}(\varsigma_{0}^{\mathrm{dr}})^{2}\big{|}+(v^{\mathrm{eff}}-q_{1})(\varsigma_{0}^{\mathrm{dr}})^{-1}\rho_{\mathsf{p}}FF^{*}\varsigma_{0}^{\mathrm{dr}}\big{)}$
$\displaystyle\hskip 44.0pt=\nu^{-1}(v^{\mathrm{eff}}-q_{1})R^{-1}CR.$ (7.42)

As anticipated, in the new basis $R^{-1}AR$ is simply multiplication by
$\nu^{-1}(v^{\mathrm{eff}}-q_{1})$. We conclude that

$\mathrm{e}^{At}C=R\exp\\!\big{(}\nu^{-1}(v^{\mathrm{eff}}-q_{1})t\big{)}R^{-1}C.$
(7.43)

Working out the algebra, one arrives at

$\mathrm{e}^{At}C=\begin{pmatrix}\nu^{3}\langle\rho_{\mathsf{p}}\exp\\!\big{(}\nu^{-1}(v^{\mathrm{eff}}-q_{1})t\big{)}(\varsigma_{0}^{\mathrm{dr}})^{2}\rangle&-\big{\langle}\rho_{\mathsf{p}}\exp\\!\big{(}\nu^{-1}(v^{\mathrm{eff}}-q_{1})t\big{)}\varsigma_{0}^{\mathrm{dr}}\big{|}F\\\\[2.15277pt]
-\nu^{2}F^{*}\big{|}\rho_{\mathsf{p}}\exp\\!\big{(}\nu^{-1}(v^{\mathrm{eff}}-q_{1})t\big{)}\varsigma_{0}^{\mathrm{dr}}\big{\rangle}&\nu
F^{*}\rho_{\mathsf{p}}\exp\\!\big{(}\nu^{-1}(v^{\mathrm{eff}}-q_{1})t\big{)}F\end{pmatrix}.$
(7.44)

Adding the spatial dependence, first in Fourier space, the propagator
$\mathrm{e}^{At}C$ is modified to $\mathrm{e}^{\mathrm{i}kAt}C$, which in
position space yields (7.32). For completeness the matrix $A$ is recorded as

$A=\nu^{-1}R(v^{\mathrm{eff}}-q_{1})R^{-1}=\begin{pmatrix}\langle\rho_{\mathsf{n}}(v^{\mathrm{eff}}-q_{1})\varsigma_{0}^{\mathrm{dr}}\rangle&-\big{\langle}(v^{\mathrm{eff}}-q_{1})\varsigma_{0}^{\mathrm{dr}}(1-\rho_{\mathsf{n}}T)\big{|}\\\\[2.15277pt]
-\nu^{-1}F^{*}\big{|}\rho_{\mathsf{n}}(v^{\mathrm{eff}}-q_{1})\big{\rangle}&\nu^{-1}F^{*}(v^{\mathrm{eff}}-q_{1})(1-\rho_{\mathsf{n}}T)\end{pmatrix}.$
(7.45)

In Figure 8 we report on molecular dynamics simulations of the Toda lattice in
thermal equilibrium for $3$, out of $9$ published, distinct parameter values.
While there is a large body of work in providing microscopic confirmation of
GHD in the context of a variety of models, the results below presumably
constitute the most extensive check, so far. There are several points in
favor. The initial statistical state is sharply defined. Modulo numerically
computing the quantities appearing in the TBA formalism, there is the concrete
prediction (7.32) based on linearized GHD. The comparison uses no adjustable
parameter. Errors are smaller than $3.5\%$, which is of the order of
statistical noise. For the shorter time $t=150$ (data not displayed) one still
observes ballistic scaling for our range of parameters. However, as systematic
deviation, the peaks are slightly lower than predicted by the theory.

At the highest pressure the peaks are shifted towards the boundaries, which is
even more pronounced for $\beta=2.0,P=3.53$ (data not displayed). This
behavior reflects that for $\beta,P\to\infty$, at fixed ratio, the Toda DOS
converges to the arcsine law which exhibits an inverse square root singularity
at either boundary, compare with (3.75). To be emphasized is the middle frame
at which $\nu=-0.03$, hence $|\rho_{\mathsf{f}}|=33$. Closer to $\nu=0$ the
TBA simulation becomes difficult. Of course for molecular dynamics one would
merely have to set $q_{1}=q_{N}$ at which point the typical distance of
particles is order $1/\sqrt{N}$, which equals $0.02$ in our case. At such high
densities the collision rate ansatz becomes dubious. But the simulation
confirms that GHD remains valid even under such extreme conditions.

Figure 8: Toda spacetime correlation functions in thermal equilibrium. The
dashed line is the GHD prediction (7.32) for $m,n=0,1,2$ and the colored lines
are molecular dynamics simulations. The system size is $N=3000$ and time
$t=600$. The parameters $\beta,P$ are indicated and $\nu=2.58,-0.03,-0.42$.
Thermal average is based on $3\times 10^{6}$ samples. From [200]. Toda
spacetime correlation functions in thermal equilibrium. The dashed line is the
GHD prediction (7.32) for $m,n=0,1,2$ and the colored lines are molecular
dynamics simulations. The system size is $N=3000$ and time $t=600$. The
parameters $\beta,P$ are indicated and $\nu=2.58,-0.03,-0.42$. The thermal
average is over samples is $3\times 10^{6}$. From [200].

Actually, the peak structure on the Euler scale is only the starting step of
the Landau-Lifshitz theory. Its main focus is the broadening of the peaks due
to dissipation and molecular noise. For a simple fluid in three dimensions the
broadening is of order $\sqrt{t}$ with a Gaussian shape function. In one
dimension the standard Landau-Lifshitz theory would make the same prediction.
But the linear fluctuation theory has to be corrected by expanding the
nonlinear Euler term to second order, which then results in a well-confirmed
superdiffusive spreading as $t^{2/3}$. More details can be found in Section
15.2. But even before there is a generic obstacle. For a well-isolated peak
the broadening is easily observed. But for the Toda lattice the Euler term
yields already a smeared out background. So one has to specifically tune to
physical situations for which the anticipated dissipative corrections can be
detected.

$\blackdiamond\blackdiamond$ Equilibrium spacetime correlations for hard rods.
For hard rods the $T$ operator can be inverted and identities as (7.32) become
more explicit. We consider thermal equilibrium, hence
$h(w)=h_{\beta}(w)=\sqrt{\beta/2\pi}\exp(-\beta w^{2}/2)$, $q_{1}=0$, and

$v^{\mathrm{eff}}(v)=\frac{\nu}{\nu-\mathsfit{a}}v.$ (7.46)

The dressing operator reads

$(1-T\rho_{\mathsf{n}})^{-1}=1+\frac{\mathsfit{a}}{\nu}|\varsigma_{0}\rangle\langle
h_{\beta}|$ (7.47)

and hence

$\varsigma_{0}^{\mathrm{dr}}=1+\frac{\mathsfit{a}}{\nu},\quad\varsigma_{1}^{\mathrm{dr}}=v,\quad\varsigma_{2}^{\mathrm{dr}}=v^{2}+\frac{\mathsfit{a}}{\beta\nu}.$
(7.48)

Rescaling the spatial coordinate to $\mathsfit{x}=(\nu-\mathsfit{a})x/t$ one
concludes

$\\{S_{00}^{\mathrm{hr}}(x,t),S_{11}^{\mathrm{hr}}(x,t),S_{22}^{\mathrm{hr}}(x,t)\\}=(\nu-\mathsfit{a})t^{-1}h_{\beta}(\mathsfit{x})\\{1+(\mathsfit{a}/\nu))^{2},\mathsfit{x}^{2},(\mathsfit{x}^{2}+(\mathsfit{a}/\beta\nu))^{2}\\}.$
(7.49)

The interesting term is the energy-energy spacetime correlation
$S_{22}^{\mathrm{hr}}(x,t)$. For the ideal gas, $\mathsfit{a}=0$, the scaling
function $\mathsfit{x}^{4}h_{\beta}(\mathsfit{x})$ vanishes at
$\mathsfit{x}=0$ and has a double peak. Turning on the interaction the minimum
at $\mathsfit{x}=0$ is shifted to $h_{\beta}(0)(\mathsfit{a}/\beta\nu)^{4}$,
resulting in three peaks. $\blackdiamond\blackdiamond$

Notes and references

Section 7.0

The book by Forster [120] is still a very readable exposition of the
fluctuation theory developed by Landau and Lifshitz [186]. Critical dynamics
is covered in the classic review of Hohenberg and Halperin [155].

Section 7.1

The existence of the infinite volume dynamics is proved by Lanford et al.
[188]. Fluctuation theory for simple fluids is discussed in Spohn [266]. As a
crucial distinction, the hydrodynamic equations for a few conservation laws
may develop shocks. An instructive account is the most readable lecture course
by Bressan [36]. In Das et al. [67] a spin chain is considered, for which at
low temperatures phase differences emerge as an almost conserved field, on top
of the two strictly conserved fields.

Section 7.2

A more detailed account of GGE spacetime correlations can be found in Spohn
[270]. Early molecular dynamics simulations of the Toda lattice and related
theoretical investigations have been carried out by Schneider and Stoll [254],
Schneider [253], Diederich [77], and Cuccoli et al. [66]. Molecular dynamics
simulations of thermal spacetime correlations are reported by Kundu and Dhar
[183]. Amongst other features, investigated are the low density regime and the
low temperature harmonic approximation. A wider range of parameters is covered
by Mazzuca et al. [200] including self- and cross-correlations of stretch,
momentum, and energy. The numerical results are compared with the GHD
prediction in Eq. (7.32). A more general perspective on fluctuations at the
Euler scale and beyond is developed by Doyon and Myers [91].

##  8 Domain wall initial states

Spacetime correlations encode small deviations from a GGE. But the
hydrodynamic scale covers also initial states very far from stationarity,
thereby exhibiting a rich portfolio of time-dependent behavior. From a
theoretical perspective, a very natural and much studied nonequilibrium
initial state is enforced by joining two semi-infinite systems, the left half
line in one GGE and the right half line in another GGE. The fields are
spatially constant except for a single jump at the origin. On the level of
GHD, the solution scales exactly as $x/t$, a simplification which promises to
yield exact solutions. In the mathematical theory of finitely many hyperbolic
conservation laws, such set-up is known as Riemann problem. Its solution
consists of flat pieces, smooth rarefaction waves, and jump discontinuities,
called shocks. Their precise spatial sequence can be complicated and is, in
principle, encoded by the eigenvectors and eigenvalues of the matrix
$\mathsfit{A}$ in their dependence on the thermodynamic parameters, compare
with (7.15). As discussed in Section 6.2 in normal form the flux Jacobian $A$
is multiplication by $\nu^{-1}(v^{\mathrm{eff}}-q_{1})$. This form admits only
rarefaction waves. No spontaneous generation of jump discontinuities is
expected.

Besides its intrinsic interest, we pick this example to also illustrate the
predictive power of GHD. For the domain wall problem, the primary goal is to
figure out macroscopic behavior, i.e. how stretch, momentum, and energy vary
spatially given their particular initial jump discontinuity. Apparently GHD is
the only theoretical scheme through which such a macroscopic behavior can be
computed.

We slightly rewrite Eq. (6.25) as

$\partial_{t}\rho_{\mathsf{n}}(x,t;v)+\tilde{v}^{\mathrm{eff}}(x,t;v)\partial_{x}\rho_{\mathsf{n}}(x,t;v)=0,\quad$
(8.1)

with

$\tilde{v}^{\mathrm{eff}}(v)=\nu^{-1}(v^{\mathrm{eff}}(v)-q_{1}).$ (8.2)

Domain wall initial conditions read

$\rho_{\mathrm{dw}}(x,0;v)=\chi(\\{x<0\\})\rho_{\mathsf{n}-}(v)+\chi(\\{x\geq
0\\})\rho_{\mathsf{n}+}(v).$ (8.3)

Instead of $\rho_{\mathsf{n}\pm}$, physically it might be more natural to
prescribe the DOS of the Lax matrix and the average stretch. But
mathematically the normal form (8.1) together with (8.3) is more accessible.

Since the solution to (8.1), (8.3) scales ballistically, we set
$\rho_{\mathrm{dw}}(x,t;v)=\mathsf{g}(t^{-1}x;v)$ and
$\tilde{v}^{\mathrm{eff}}(x,t;v)=\tilde{\mathsf{v}}^{\mathrm{eff}}(t^{-1}x;v)$.
Without loss of generality one adopts $t=1$ and arrives at

$(x-\tilde{\mathsf{v}}^{\mathrm{eff}}(x;v))\partial_{x}\mathsf{g}(x;v)=0,\qquad\lim_{x\to\pm\infty}\mathsf{g}(x;v)=\rho_{\mathsf{n}\pm}(v).$
(8.4)

Therefore $x\mapsto\mathsf{g}(x;v)$ for fixed $v$ has to be constant except
for jumps at the zeros of $x-\tilde{\mathsf{v}}^{\mathrm{eff}}(x;v)$ as
function of $x$. At this stage, it is not clear which level of generality to
adopt. Physically, one would expect to have a unique solution. Thus in the
$(x,v)$-plane there should be a contact line which divides the plane into two
domains, one characterized to contain the set $\\{-\infty\\}\times\mathbb{R}$
and the other the set $\\{\infty\\}\times\mathbb{R}$. The solution
$\mathsf{g}(x;v)$ is constant in either domain and jumps from
$\rho_{\mathsf{n}-}(v)$ to $\rho_{\mathsf{n}+}(v)$ across the contact line.
Uniqueness means a unique contact point for every $v$. Then the contact line
is represented as the graph of a function defined by
$v\mapsto\tilde{\phi}(v)$. As discussed below, for hard rods, $\tilde{\phi}$
is invertible with inverse denoted by $\phi$. Hence the contact line can be
written also as the graph of $x\mapsto\phi(x)$. We present our argument for
the latter case, since hard rods are reasonably close to the Toda lattice.
More general contact lines could be handled in a similar fashion.

With our assumptions, for every $x$ there is a unique contact point
$v=\phi(x)$, which satisfies

$x-\tilde{\mathsf{v}}^{\mathrm{eff}}(x;v),$ (8.5)

and the solution ansatz reads

$g^{\phi}(v)=\chi(\\{v>\phi\\})\rho_{\mathsf{n}-}(v)+\chi(\\{v\leq\phi\\})\rho_{\mathsf{n}+}(v).$
(8.6)

The superscript $\phi$ will be used to generically indicate that in the TBA
formalism $\rho_{\mathsf{n}}$ is replaced by $g^{\phi}$. For example,

$f^{\mathrm{dr},\phi}=\big{(}1-Tg^{\phi}\big{)}^{-1}f,\qquad\nu^{\phi}\langle\rho_{\mathsf{p}}^{\phi}\rangle=1,$
(8.7)

compare with (3.57) and (3.60). Note that
$g^{\pm\infty}(v)=\rho_{\mathsf{n}\pm}(v)$. In particularly, we define

$\tilde{\mathsf{v}}^{\mathrm{eff},\phi}(v)=(\nu^{\phi})^{-1}\Big{(}\frac{\varsigma_{1}^{\phi}(v)}{\varsigma_{0}^{\phi}(v)}-q_{1}^{\phi}\Big{)}.$
(8.8)

From the solution ansatz it follows that

$\mathsf{g}(x;v)=g^{\phi(x)}(v),\qquad\tilde{\mathsf{v}}^{\mathrm{eff}}(x;v)=\tilde{\mathsf{v}}^{\mathrm{eff},\phi(x)}(v).$
(8.9)

The condition (8.5) turns into

$x=\tilde{\mathsf{v}}^{\mathrm{eff},\phi(x)}(\phi(x)).$ (8.10)

Defining

$G(\phi)=\tilde{\mathsf{v}}^{\mathrm{eff},\phi}(\phi),$ (8.11)

we conclude that

$x=G(\phi(x)),\qquad\phi(x)=G^{-1}(x).$ (8.12)

The contact line is the inverse of the function $G$ in (8.11). To determine
$G$ numerically, one has to compute the dressing depending on $\phi$. This
then yields $\mathsf{g}(x;v)$ from which the observables of interest, as
average stretch, momentum, and energy, at location $x$ are deduced.

Qualitatively $G$ is linear for large $\phi$. But this is a special feature of
fluid-like models with a quadratic kinetic energy. For the discrete linear
wave equation with a general dispersion relation, because of the bounded
momentum space, the concept of contact line as such remains untouched, but the
line is no longer equal to the graph of a function.

The predicted sharp step appears to be oversimplified. Indeed on longer time
scales the step is expected to be smeared roughly as an error function,
compare with Section 15.3.

Figure 9: Domain wall initial state for the Toda lattice at uniform inverse
temperature $\beta=1$, $P_{-}=\frac{1}{2}$, and $P_{+}=2$. Displayed in blue
are $\tilde{\mathsf{v}}^{\mathrm{eff},\phi}(v)$ at $\phi=2$ and the contact
line $G(\phi)$. On the left, the dotted magenta line is
$\tilde{\mathsf{v}}^{\mathrm{eff}}(v)$ computed by using $\rho_{\mathsf{n}+}$
and the olive line using $\rho_{\mathsf{n}-}$. Correspondingly on the right
displayed is $\tilde{\mathsf{v}}^{\mathrm{eff}}(\phi)$ using either
$\rho_{\mathsf{n}+}$ or $\rho_{\mathsf{n}-}$. From [204].

$\blackdiamond\blackdiamond$ Hard rod domain wall. For the hard rod lattice
the function $G$ can still be computed analytically. We recall that $\nu$ is
the average stretch, $\nu>\mathsfit{a}$, and $h(v)$ denotes the normalized
velocity distribution,

$\rho_{\mathsf{p}}(v)=\nu^{-1}h(v),\quad\rho_{\mathsf{n}}(v)=(\nu-\mathsfit{a})^{-1}h(v),\quad
q_{1}=\int_{\mathbb{R}}\mathrm{d}wh(w)w.$ (8.13)

In normal form the hard rod GHD reads

$\partial_{t}\rho_{\mathsf{n}}(x,t;v)+(\rho(x,t)-\mathsfit{a})^{-1}(v-q_{1}(x,t))\partial_{x}\rho_{\mathsf{n}}(x,t;v)=0,$
(8.14)

compare with the structurally identical Eq. (8.1). The initial condition is
the same as in (8.3) and the $G$-function is computed by the method explained
above with the result

$G_{\mathrm{hr}}(\phi)=\Big{(}\int_{\phi}^{\infty}\mathrm{d}v\rho_{\mathsf{n}-}(v)+\int^{\phi}_{-\infty}\mathrm{d}v\rho_{\mathsf{n}+}(v)\Big{)}\phi-\int_{\phi}^{\infty}\mathrm{d}v\,v\rho_{\mathsf{n}-}(v)-\int^{\phi}_{-\infty}\mathrm{d}v\,v\rho_{\mathsf{n}+}(v).$
(8.15)

Clearly, $dG_{\mathrm{hr}}(\phi)/d\phi>0$. Thus $G_{\mathrm{hr}}$ is
invertible. For $\phi\to\infty$ one obtains

$G_{\mathrm{hr}}(\phi)\simeq\frac{1}{\nu_{+}-\mathsfit{a}}(\phi-q_{1,+})$
(8.16)

and correspondingly for $\phi\to-\infty$ with $\nu_{+},q_{1+}$ replaced by
$\nu_{-},q_{1-}$. If the boundary velocity distributions have an exponential
decay, then the error in (8.14) is also exponentially small. The contact line
is asymptotically linear with a well localized monotone interpolation.

If $\mathsfit{a}<0$, one can tune the boundary values such that
$\nu_{-}<0<\nu_{+}$. There is then some intermediate point, $x_{0}$ with
$\nu(x_{0})=0$, a case studied in more detail at the end of Section 9.1.
$\blackdiamond\blackdiamond$

Notes and references

Section 8

This chapter is based on the study of Mendl and Spohn [204], where in
particular numerical simulations of GHD with domain wall initial conditions
are reported. Domain wall for hard rods are studied in Doyon and Spohn [92].
On the Euler scale the step at the contact line is sharp. Earlier work of
Bulchandani et al. [39] considers the case of low pressure,
$P_{-}<P_{+}<P_{\mathrm{c}}$. Molecular dynamics is compared with predictions
of GHD. Also quantum mechanically corrections are obtained to lowest order.
For the Toda lattice, diffusive corrections are discussed in Section 15.3. The
broadening is convincingly observed in the XXZ model, see De Nardis et al.
[74].

Domain wall initial conditions have been studied numerically for the discrete
sinh-Gordon model by Bastianello et al. [18]. The most detailed investigations
are available for the XXZ model, see Piroli et al. [233] and Misguich et al.
[209]. In this case the spectral parameter space contains in addition the type
of string states and the contact line refers to a larger space. Gamayun et al.
[124] study domain wall initial conditions for a classical spin chain and its
quantized version.

Due to momentum conservation, the contact line of the Toda lattice is one-to-
one and covers the full real line. On the other hand, for XXZ, and other
discrete models, the contact line takes values only in a bounded $v$-interval,
which implies left and right edges in $x$-space up to which boundary values
remain constant. The behavior near the edge often shows intricate oscillatory
decay, which is beyond the hydrodynamic scale and has been elucidated in
considerable detail for a variety of models, see Collura et al. [60],
Bulchandani and Karrasch [40], and Grava et al. [138].

##  9 Toda fluid

Toda particles might as well be viewed to move on the real line, which is the
fluid picture, see Fig. 1. The name “fluid” is somewhat misleading, since
particles are distinguishable. Still the dynamical behavior is fluid-like. In
a way, we have to start from the beginning and redo the computation for
average fields and their currents. Since we rely on chain results, the
discussion can be more compressed. To avoid duplicating symbols, the fluid is
distinguished from the chain through the index “f”. $x\in\mathbb{R}$ stands
for the coordinate of the one-dimensional physical space.

###  9.1 Euler equations

For the Toda fluid conserved fields have a density given by

$Q^{[0]}_{\mathsf{f}}(x)=\sum_{j\in\mathbb{Z}}\delta(q_{j}-x),\qquad
Q^{[n]}_{\mathsf{f}}(x)=\sum_{j\in\mathbb{Z}}\delta(q_{j}-x)Q^{[n]}_{j}$ (9.1)

for $\;n\geq 1$. Taking their time derivative,

$\displaystyle\frac{d}{dt}Q^{[0]}_{\mathsf{f}}(x)=\sum_{j\in\mathbb{Z}}\delta^{\prime}(q_{j}-x)p_{j},$
$\displaystyle\frac{d}{dt}Q^{[n]}_{\mathsf{f}}(x)=\sum_{j\in\mathbb{Z}}\big{(}\delta^{\prime}(q_{j}-x)p_{j}Q^{[n]}_{j}+\delta(q_{j}-x)(J_{j}^{[n]}-J_{j+1}^{[n]})\big{)}.$
(9.2)

Hence fluid current densities read

$\displaystyle
J^{[0]}_{\mathsf{f}}(x)=\sum_{j\in\mathbb{Z}}\delta(q_{j}-x)p_{j},$ (9.3)
$\displaystyle
J^{[n]}_{\mathsf{f}}(x)=\sum_{j\in\mathbb{Z}}\big{(}\delta(q_{j}-x)p_{j}Q^{[n]}_{j}+\big{(}\theta(q_{j}-x)-\theta(q_{j-1}-x)\big{)}J_{j}^{[n]}\big{)}$

with $\theta$ the step function, $\theta(x)=0$ for $x\leq 0$ and $\theta(x)=1$
for $x>0$. The index $0$ has been treated separately. But in fact,
$Q^{[n]}_{j}$ naturally extends to $n=0$ by setting $Q^{[0]}_{j}=1$. This
property foreshadows a single hydrodynamic equation.

The next item is GGE, where we start from the canonical ensemble with
partition function

$Z_{\mathrm{can}}(N,\ell)=\int_{\Gamma_{N}}\prod_{j=1}^{N}\mathrm{d}r_{j}\mathrm{d}p_{j}\delta\big{(}Q^{[0],N}-\ell\big{)}\exp\\!\big{(}-\mathrm{tr}[V(L_{N})]\big{)}.$
(9.4)

The ensemble is canonical, since volume, $q_{N+1}-q_{1}$, and number of
particles are fixed. In our discussion the confining potential will play a
passive role and is hence dropped from the notation, except when appearing as
GGE index. For the chain we considered $\ell=\nu N$ and, by the equivalence of
ensembles, switched to $\exp\\!\big{(}-PQ^{[0],N}\big{)}$ with the pressure
$P$ dual to the stretch $\nu$. In the limit $N\to\infty$ the free energy per
site, $F_{\mathrm{to}}(P)$, has been obtained already, compare with (3.40).
Turning to the fluid we want to keep the volume fixed and allow for
fluctuations in $N$. The corresponding chemical potential is denoted by $\mu$.
Then

$Z_{\mathsf{f}}(\mu,\ell)=\sum_{N=1}^{\infty}\mathrm{e}^{\mu
N}\int_{\Gamma_{N}}\prod_{j=1}^{N}\mathrm{d}r_{j}\mathrm{d}p_{j}\delta\big{(}Q^{[0],N}-\ell\big{)}\exp\\!\big{(}-\mathrm{tr}[V(L_{N})]\big{)},$
(9.5)

compare with (3.1). One has to distinguish the cases $\ell>0$ and $\ell<0$,
corresponding to either increasing or decreasing particle labelling, on
average. Thus the particle number $Q_{\mathsf{f}}^{[0],N}>0$, while the fluid
density $\rho_{\mathsf{f}}=N/\ell$ can take either sign. Since $\nu=\ell/N$,
the relation $\nu\rho_{\mathsf{f}}=1$ holds always. Close to $P_{\mathrm{c}}$,
defined by $\nu(P_{\mathrm{c}})=0$, the fluid density $\rho_{\mathsf{f}}$
jumps from $\infty$ to $-\infty$ when increasing $P$. This looks singular but
merely reflects a particular choice of coordinates. $\nu(P)$ is a smooth
function. Of interest is the limit

$\lim_{\ell\to\pm\infty}-\frac{1}{\ell}\log
Z_{\mathsf{f}}(\mu,\ell)=F_{\mathsf{f},\pm}(\mu),$ (9.6)

where $\pm$ refers to the sign of $\ell$. The fluid free energy has two
branches labeled by $\mathrm{sgn}(\ell)$. Correspondingly there are two
distinct infinite volume GGE averages denoted by
$\langle\cdot\rangle_{\mu,\pm,V}$. For a convergent partition function, one
has to require $\mu<\mu_{\mathrm{max}}$ for either branch with
$\mu_{\mathrm{max}}=F_{\mathrm{to}}(P_{\mathrm{c}})$, as displayed in Figure
10.

Figure 10:  For thermal equilibrium at $\beta=1$, to the left is the Toda
lattice free energy $F_{\mathrm{to}}=\mu$ as a function of the pressure $P$.
To the right are the two branches of the fluid free energy
$F_{\mathsf{f},\pm}=-P$ as a function of $\mu$, see (9.10).

Thermodynamically the dual to $P$ is the stretch $\nu$ and to $\mu$ the fluid
density $\rho_{\mathsf{f}}$. Let us denote by $\mathcal{L}F$ the Legendre
transform of some function $F$. By the equivalence of ensembles on the level
of free energies one concludes that

$\lim_{N\to\infty}-N^{-1}\log Z_{\mathsf{can}}(N,\nu
N)=\mathcal{L}F_{\mathrm{to}}(\nu)$ (9.7)

and

$\lim_{\ell\to\pm\infty}-\ell^{-1}\log
Z_{\mathrm{can}}(\rho_{\mathsf{f}}\ell,\ell)=\mathcal{L}F_{\mathsf{f},\pm}(\rho_{\mathsf{f}}).$
(9.8)

Thus the two free energies are related by

$\mathcal{L}F_{\mathrm{to}}(\nu)=\rho_{\mathsf{f}}^{-1}\mathcal{L}F_{\mathsf{f},\pm}(\rho_{\mathsf{f}}),\quad\nu\rho_{\mathsf{f}}=1,$
(9.9)

which implies

$F_{\mathrm{to}}(P)=\mu,\quad-F_{\mathsf{f},\pm}(\mu)=P.$ (9.10)

Stated differently, $F_{\mathrm{to}}$ and $-F_{\mathsf{f},\pm}$ are inverse
functions of each other. The first identity was already noted in (3.40) and
the argument of $F_{\mathsf{f},\pm}$ is in fact the parameter appearing in the
TBA equation. $F_{\mathrm{to}}(P)$ is concave with maximum at
$P_{\mathrm{c}}$. On the other hand, $F_{\mathsf{f},+}$ is concave and
$F_{\mathsf{f},-}$ convex with both graphs being smoothly joined at
$\mu_{\mathrm{max}}$.

The free energy of the Toda fluid can be written in such a way as to suggest
the extension to a larger class of integrable models. Starting from (9.10) and
using (3.41), the free energy equals

$F_{\mathsf{f}}(\mu,V)=-\int_{\mathbb{R}}\mathrm{d}w\mathrm{e}^{-\varepsilon(w)},$
(9.11)

compare with (9.10) and (3.37). Differentiating the TBA equation with respect
to $\mu$ one arrives at

$-\partial_{\mu}F_{\mathsf{f}}(\mu,V)=\int_{\mathbb{R}}\mathrm{d}w\rho_{\mathsf{n}}(w)(1-T\rho{n})^{-1}\varsigma_{0}(w)=\nu^{-1}=\rho_{\mathsf{f}}.$
(9.12)

Similarly by perturbing $V$ as $V+\kappa\varsigma_{n}$, one concludes

$\partial_{\kappa}F_{\mathsf{f}}(\mu,V+\kappa\varsigma_{n})|_{\kappa=0}=\int_{\mathbb{R}}\mathrm{d}w\rho_{\mathsf{n}}(w)(1-T\rho_{\mathsf{n}})^{-1}\varsigma_{n}(w)=\langle\rho_{\mathsf{p}}\varsigma_{n}\rangle=\langle
Q_{\mathsf{f}}^{[n]}(0)\rangle_{\mu,V},$ (9.13)

which confirms the defining derivatives for the fluid free energy
$F_{\mathsf{f}}$.

For the current averages we follow the strategy as developed for the lattice
case by introducing the field-current susceptibility matrix through

$B_{\mathsf{f},m,n}=\int_{\mathbb{R}}\mathrm{d}x\langle
J_{\mathsf{f}}^{[m]}(x)Q_{\mathsf{f}}^{[n]}(0)\rangle^{\mathrm{c}}_{\mu,\pm,V}.$
(9.14)

By the same argument as leading to (6.9), $B_{\mathsf{f}}$ is a symmetric
matrix, $B_{\mathsf{f},m,n}=B_{\mathsf{f},n,m}$. Since
$J_{\mathsf{f}}^{[0]}=Q_{\mathsf{f}}^{[1]}$, one starts from

$\partial_{\mu}\langle
J_{\mathsf{f}}^{[n]}(0)\rangle_{\mu,\pm,V}=B_{\mathsf{f},n,0}=B_{\mathsf{f},0,n}=\int_{\mathbb{R}}\mathrm{d}x\langle
Q_{\mathsf{f}}^{[1]}(x)Q_{\mathsf{f}}^{[n]}(0)\rangle^{\mathrm{c}}_{\mu,\pm,V}.$
(9.15)

The last expression can be written as derivative,

$B_{\mathsf{f},0,n}=-\partial_{\kappa}\langle
Q_{\mathsf{f}}^{[n]}(0)\rangle_{\mu,\pm,V+\kappa
w}\big{|}_{\kappa=0}=-\partial_{\kappa}\langle\rho_{\mathsf{p},\kappa}\varsigma_{n}\rangle\big{|}_{\kappa=0},$
(9.16)

where $\rho_{\mathsf{p},\kappa}$ is determined by the solution of the TBA
equation

$V+\kappa\varsigma_{1}-\mu-T\rho_{\mathsf{n},\kappa}+\log\rho_{\mathsf{n},\kappa}=0.$
(9.17)

In (3.60) we differentiated TBA with respect to $\varsigma_{0}$. This time it
is with respect to $\varsigma_{1}$. Denoting by ′ the derivative at
$\kappa=0$, one obtains

$\rho_{\mathsf{n}}^{\prime}=-\rho_{\mathsf{n}}(1-T\rho_{\mathsf{n}})^{-1}\varsigma_{1}$
(9.18)

and

$\rho_{\mathsf{p}}^{\prime}=\partial_{\mu}\rho_{\mathsf{n}}^{\prime}=-\partial_{\mu}(\rho_{\mathsf{n}}\varsigma_{1}^{\mathrm{dr}})=-\partial_{\mu}\big{(}\rho_{\mathsf{p}}v^{\mathrm{eff}}\big{)},$
(9.19)

where, as before, $v^{\mathrm{eff}}$ is defined in (6.20). Therefore

$\partial_{\mu}\big{(}\langle
J_{\mathsf{f}}^{[n]}(0)\rangle_{\mu,\pm,V}-\langle\rho_{\mathsf{p}}v^{\mathrm{eff}}\varsigma_{n}\rangle\big{)}=0,$
(9.20)

implying that the difference does not depend on $\mu$. To fix the constant,
the current vanishes in the limit of vanishing density of particles, in other
words for $\mu\to-\infty$, and so does $\rho_{\mathsf{p}}$ since
$\nu\langle\rho_{\mathsf{p}}\rangle=1$. We conclude that

$\langle
J_{\mathsf{f}}^{[n]}(0)\rangle_{\mu,\pm,V}=\langle\rho_{\mathsf{p}}v^{\mathrm{eff}}\varsigma_{n}\rangle,$
(9.21)

which holds for all $n$, including $n=0$. Compared to the lattice currents
(6.18), only the shift $\rho_{\mathsf{p}}q_{1}$ has been removed. This can be
understood by starting from the fluid side. Then switching to the lattice
dynamical equations, the lateral motion is subtracted and such a term has to
appear in the average lattice currents.

Since $B_{\mathsf{f}}$ is a symmetric matrix, there has to be a current
potential, denoted by $G_{\mathsf{f}}$, such that

$\partial_{\kappa}G_{\mathsf{f}}(\mu,V+\kappa\varsigma_{n})|_{\kappa=0}=\langle
J_{\mathsf{f}}^{[n]}(0)\rangle_{\mu,V}.$ (9.22)

Indeed, one only has to slightly modify the free energy as

$G_{\mathsf{f}}(\mu,V)=\int_{\mathbb{R}}\mathrm{d}ww\mathrm{e}^{-\varepsilon(w)}.$
(9.23)

Then, as before by differentiating TBA,

$\partial_{\kappa}G_{\mathsf{f}}(\mu,V+\kappa\varsigma_{n})|_{\kappa=0}=\int_{\mathbb{R}}\mathrm{d}ww\rho_{\mathsf{n}}(w)(1-T\rho_{\mathsf{n}})^{-1}\varsigma_{n}(w)=\langle\varsigma_{1}^{\mathrm{dr}}\rho_{\mathsf{n}}\varsigma_{n}\rangle.$
(9.24)

Since in (9.21) the term $n=0$ no longer plays a special role, there is only a
single hydrodynamic equation which reads

$\partial_{t}\rho_{\mathsf{p}}(x,t;v)+\partial_{x}\big{(}v^{\mathrm{eff}}(x,t;v)\rho_{\mathsf{p}}(x,t;v)\big{)}=0.$
(9.25)

Using (6.24) to transform to quasilinear form yields

$\partial_{t}\rho_{\mathsf{n}}(x,t;v)+v^{\mathrm{eff}}(x,t;v)\partial_{x}\rho_{\mathsf{n}}(x,t;v)=0,$
(9.26)

which is identical to the lattice version, except for the term accounting for
the lateral motion. By construction $\rho_{\mathsf{n}}(x,t;v)\geq 0$, a
property which is correctly propagated by (9.26). On the other hand the
overall sign of $\rho_{\mathsf{p}}(x,t;v)$ is governed by $\nu(x,t)$, to say,
if $\nu(x,t)>0$, then $\rho_{\mathsf{p}}(x,t;v)\geq 0$, and if $\nu(x,t)<0$,
then $\rho_{\mathsf{p}}(x,t;v)\leq 0$. An example for $\nu(x,t)$ changing sign
is discussed in the Insert just below.

$\blackdiamond\blackdiamond$ High-low pressure domain wall. As discussed in
Chapter 8, a domain wall initial state consists of two half lines with
homogeneous states $\rho_{\mathsf{n}-}$, resp. $\rho_{\mathsf{n}+}$. If
$\nu_{-}<0<\nu_{+}$, then there is a sign change of the particle density
$\rho_{\mathsf{p}}$. While in that Chapter we used the Toda lattice, it is
instructive to view the high-low pressure domain wall from the perspective of
a fluid. To be concrete, in the initial state the particle with label $0$ is
placed at the origin, $q_{0}=0$, the particles with negative indices
correspond to a thermal state with parameters $(\nu_{-},\beta)$ and those with
positive indices to $(\nu_{+},\beta)$. Thus the momenta are i.i.d. Gaussian
with mean zero and variance $\beta^{-1}$, while $\\{q_{j+1}-q_{j},j\leq-1\\}$
are i.i.d. $\chi_{2P_{-}}$ distributed random variables, and
$\\{q_{j+1}-q_{j},j\geq 0\\}$ are i.i.d. $\chi_{2P_{+}}$, compare with Chapter
8. The case of interest is $P_{-}>P_{\mathrm{c}}>P_{+}$, thus
$\nu_{-}<0<\nu_{+}$. Then, by assumption, $\langle
q_{j+1}-q_{j}\rangle_{P_{-}}=\nu_{-}<0$ for $j\leq-1$, while $\langle
q_{j+1}-q_{j}\rangle_{P_{+}}=\nu_{+}>0$ for $j\geq 0$.

A typical ordering is of the form $\dots>q_{-1}>q_{0}=0<q_{1}<\dots$ . Thus
initially the average particle density equals $|\nu_{-}|^{-1}+\nu_{+}^{-1}$ on
$[0,\infty)$ and decays exponentially on $(-\infty,0]$. Under the dynamics,
the point at which the ordered domain touches the anti-ordered domain is
moving in time, its label being denoted by $\kappa(t)$. Close to
$q_{\kappa(t)}$ particles pile up. The domain boundary acts as bottleneck for
particles. The position of the bottleneck is
$x_{\mathrm{bn}}(t)=q_{\kappa(t)}$. Numerical simulations suggest that
$\kappa(t)$ and $x_{\mathrm{bn}}(t)$ change linearly in time, at least
approximately. Also to the left of the bottleneck the particle density
vanishes rapidly. As observed numerically, the scaling function $G$, see
(8.12), has a nonzero slope at $x_{\mathrm{c}}$. Thus, near the bottleneck the
particles should be distributed according to equilibrium with a linearly
varying pressure. From this property one infers that in physical space the
particle density at the bottleneck diverges as an inverse square root. It
would be of interest to better understand the precise particle statistics
close to the bottleneck. $\blackdiamond\blackdiamond$

###  9.2 Generalized free energy - again

The fluid picture suggests an alternative approach to obtain the generalized
free energy, thereby circumventing entirely the Dumitriu-Edelman change of
volume elements. The same method will be applied to Calogero-Moser models as
to be discussed in Chapter 11. As further spin-off, a closer structural
similarity to the Lieb-Liniger $\delta$-Bose gas will be accomplished.

One starts from the Toda fluid, $(q_{j},p_{j})\in\mathbb{R}^{2}$,
$j=1,\ldots,N$. The open chain is governed by the hamiltonian of Eq. (2.59)
and has the Lax matrix

$(L_{N}^{\diamond})_{j,j}=p_{j},\qquad(L_{N}^{\diamond})_{j,j+1}=(L_{N}^{\diamond})_{j+1,j}=\mathrm{e}^{-((q_{j+1}-q_{j})/2)}$
(9.27)

with $j=1,\ldots,N$ in the former case, $j=1,\ldots,N-1$ in the latter case,
and $(L_{N}^{\diamond})_{i,j}=0$ otherwise. Some aspects of the scattering
theory for the open Toda chain have been discussed already in Section 2.3. The
explicit canonical transformation to scattering coordinates, alias action-
angle variables, is stated in Eqs. (2.66) - (2.68). By construction, under the
transformation $\Phi$, one has the identity

$\mathrm{tr}\big{[}V(L^{\diamond}_{N})\big{]}\circ\Phi=\sum_{j=1}^{N}V(\lambda_{j}).$
(9.28)

Inserting from (2.66), (2.67), in addition one obtains

$\mathrm{e}^{-q_{1}}\circ\Phi=\frac{\sigma_{N-1}}{\sigma_{N}}=\sum_{j=1}^{N}Y_{j}\mathrm{e}^{-\phi_{j}},\qquad\mathrm{e}^{q_{N}}\circ\Phi=\frac{\sigma_{1}}{\sigma_{0}}=\sum_{j=1}^{N}Y_{j}\mathrm{e}^{\phi_{j}}$
(9.29)

with

$Y_{j}=\prod_{i=1,i\neq j}^{N}|\lambda_{i}-\lambda_{j}|^{-1}.$ (9.30)

We first consider the case when the average stretch $\nu>0$. To study GGEs,
the end particles are subject to suitable boundary potentials,
$V^{+}_{\ell,1}(q_{1})$, $V^{+}_{\ell,N}(q_{N})$, such that all particles are
confined to a volume of approximate size $\ell$ with $\ell>0$. Our observation
above suggests to choose the boundary potential as

$V^{+}_{\ell,1}(q_{1})+V^{+}_{\ell,N}(q_{N})=\mathrm{e}^{-\ell/2}\big{(}\mathrm{e}^{-q_{1}}+\mathrm{e}^{q_{N}}\big{)},$
(9.31)

since $\mathrm{e}^{-q_{1}}$ and $\mathrm{e}^{q_{N}}$ are given by rather
simple formulas. Then in approximation $-\ell/2<q_{1}$ and $q_{N}<\ell/2$,
which yields the desired confinement because $\nu>0$. In case $\nu<0$, one has
to choose $\ell<0$ and the boundary potential becomes

$V_{\ell,1}^{-}(q_{1})+V_{\ell,N}^{-}(q_{N})=\mathrm{e}^{\ell/2}\big{(}\mathrm{e}^{q_{1}}+\mathrm{e}^{-q_{N}}\big{)}.$
(9.32)

But now $Y_{j}$ is switched to the denominator, which makes the analysis more
delicate. In the following we restrict ourselves to the case $\nu>0$.

Next considered is a canonical type ensemble with $\ell$ and $N$ as
parameters. Since $\Phi$ is a canonical transformation, its Jacobian
determinant equals $1$ which implies for the volume elements
$\mathrm{d}^{N}p\,\mathrm{d}^{N}q=\mathrm{d}^{N}\phi\,\mathrm{d}^{N}\lambda$.
To the energy-like term $\mathrm{tr}\big{[}V(L^{\diamond}_{N})\big{]}$ one
adds the boundary potential (9.31). For the partition function one thereby
arrives at

$\displaystyle
Z_{\mathrm{to},N}(\ell,V)=\int_{\Gamma_{N}}\mathrm{d}^{N}p\,\mathrm{d}^{N}q\exp\big{(}-\mathrm{tr}\big{[}V(L^{\diamond}_{N})\big{]}-\mathrm{e}^{-\ell/2}\big{(}\mathrm{e}^{-q_{1}}+\mathrm{e}^{q_{N}}\big{)}\big{)}$
$\displaystyle\hskip
20.0pt=\frac{1}{N!}\int_{\mathbb{R}^{N}}\mathrm{d}^{N}\phi\,\int_{\mathbb{R}^{N}}\mathrm{d}^{N}\lambda\exp\Big{(}-\sum_{j=1}^{N}V(\lambda_{j})-2\mathrm{e}^{-\ell/2}Y_{j}\cosh\phi_{j}\Big{)}.$
(9.33)

Here we used that under the map $\Phi^{-1}$ the eigenvalues are ordered. Since
the integrand is symmetric in the action variables, their integration is
extended over $\mathbb{R}^{N}$ at the expense of the factor $1/N!$. The
$\phi_{j}$-integrations can be carried out explicitly and yield

$Z_{\mathrm{to},N}(\ell,V)=\frac{1}{N!}\int_{\mathbb{R}^{N}}\mathrm{d}^{N}\lambda\prod_{j=1}^{N}\mathrm{e}^{-V(\lambda_{j})}\prod_{j=1}^{N}2K_{0}\big{(}2\mathrm{e}^{-\ell/2}Y_{j}\big{)},$
(9.34)

where

$K_{0}(x)=\int_{0}^{\infty}\mathrm{d}t\exp(-x\cosh{t})$ (9.35)

is the zero order modified Bessel function of the second kind.

As a consequence, the probability distribution of eigenvalues of the Lax
matrix $L_{N}^{\diamond}$ under the normalized measure

$Z_{\mathrm{to},N}(\ell,V)^{-1}\exp\big{(}-\mathrm{tr}\big{[}V(L^{\diamond}_{N})\big{]}-\mathrm{e}^{-\ell/2}\big{(}\mathrm{e}^{-q_{1}}+\mathrm{e}^{q_{N}}\big{)}\big{)}\mathrm{d}^{N}p\,\mathrm{d}^{N}q$
(9.36)

is given by

$Z_{\mathrm{to},N}(\ell,V)^{-1}\exp\Big{[}-\sum_{j=1}^{N}V(\lambda_{j})+\sum_{j=1}^{N}\log\Big{(}2K_{0}\big{(}2\exp\big{[}-\tfrac{1}{2}N\big{(}\nu-\frac{2}{N}\sum_{i=1,i\neq
j}^{N}\log Y_{j}\big{)}\big{]}\big{)}\Big{)}\Big{]}$ (9.37)

relative to $(1/N!)\mathrm{d}^{N}\lambda$. Note that

$-\frac{2}{N}\log Y_{j}=\frac{1}{N}\sum_{i=1,i\neq
j}^{N}\phi_{\mathrm{to}}(\lambda_{j}-\lambda_{i})$ (9.38)

which, except for the diagonal term, equals
$\phi_{\mathrm{to}}(\lambda_{j}-\cdot)$ integrated over the empirical DOS.
Only the functional in (9.37) is more complicated than before and has an
explicit $N$-dependence in addition.

The next step is the $\ell\to\infty$ asymptotics at fixed $\nu=\ell/N$. The
Bessel function behaves as $K_{0}(x)=-\log x$ for $x\to 0$ and
$K_{0}(x)=(\pi/2x)^{1/2}\mathrm{e}^{-x}$ for $x\to\infty$. Hence for large
$\ell$

$K_{0}(\mathrm{e}^{-a\ell})\simeq a\ell\hskip 8.0pt\mathrm{for}\hskip
3.0pta>0,\qquad\log K_{0}(\mathrm{e}^{-a\ell})\simeq-\mathrm{e}^{-a\ell}\hskip
8.0pt\mathrm{for}\hskip 3.0pta<0.$ (9.39)

Let us define

$\Upsilon_{j}=\nu+\frac{1}{N}\sum_{i=1,i\neq
j}^{N}\phi_{\mathrm{to}}(\lambda_{j}-\lambda_{i}).$ (9.40)

Then for large $N$

$\log\Big{(}2K_{0}\big{(}2\exp\big{[}-\tfrac{1}{2}N\Upsilon_{j}\big{]}\big{)}\Big{)}\simeq\begin{cases}\log
N+\log\Upsilon_{j},\,\,\mathrm{if}\,\,\Upsilon_{j}>0,\\\\[4.30554pt]
-2\mathrm{e}^{-(N\Upsilon_{j}/2)},\,\,\mathrm{if}\,\,\Upsilon_{j}<0.\end{cases}$
(9.41)

Following the arguments in Section 3.5, the corresponding free energy
functional reads

$\mathcal{F}_{\mathrm{to},\mathsf{f},+}^{\circ}(\varrho)=\nu^{-1}\int_{\mathbb{R}}\mathrm{d}w\varrho(w)\Big{(}V(w)-1+\log\varrho(w)-\log\Big{(}\nu+\int_{\mathbb{R}}\mathrm{d}w^{\prime}\varrho(w^{\prime})\phi_{\mathrm{to}}(w-w^{\prime})\Big{)}\Big{)}$
(9.42)

in case
$\nu+\int_{\mathbb{R}}\mathrm{d}w^{\prime}\varrho(w^{\prime})\phi_{\mathrm{to}}(w-w^{\prime})>0$
and $\mathcal{F}_{\mathrm{to},\mathsf{f},+}^{\circ}(\varrho)=\infty$
otherwise. The variation is over all probability densities with
$\mathcal{F}_{\mathrm{to},\mathsf{f},+}^{\circ}(\varrho)<\infty$. Momentarily
we ignore this infinite-dimensional constraint, but will return to the issue
below. Denoting the minimizer by $\varrho^{*}$, as before, one concludes that

$\lim_{\ell\to\infty}-\frac{1}{\ell}\log
Z_{\mathrm{to},N}(\ell,V)=F_{\mathrm{to},\mathsf{f},+}(\nu,V)=\mathcal{F}_{\mathrm{to},\mathsf{f},+}^{\circ}(\varrho^{*}).$
(9.43)

It is convenient to substitute $\nu^{-1}\varrho=\rho$. Then the free energy
per unit length of the Toda fluid is determined by the free energy functional

$\mathcal{F}_{\mathrm{to},\mathsf{f},+}(\rho)=\int_{\mathbb{R}}\mathrm{d}w\rho(w)\Big{(}V(w)-1+\log\rho(w)-\log\Big{(}1+\int_{\mathbb{R}}\mathrm{d}w^{\prime}\rho(w^{\prime})\phi_{\mathrm{to}}(w-w^{\prime})\Big{)}\Big{)},$
(9.44)

where the minimization is constrained by $\rho\geq 0$ and

$\int_{\mathbb{R}}\mathrm{d}w\rho(w)=\nu^{-1}.$ (9.45)

As before the constraint is lifted by the Lagrange multiplier $\mu$.
Introducing the space density $\rho_{\mathsf{s}}$ by

$\rho_{\mathsf{s}}=1+T\rho,$ (9.46)

the free energy functional can be written as

$\mathcal{F}_{\mathrm{to},\mathsf{f},+}^{\bullet}(\rho)=\int_{\mathbb{R}}\mathrm{d}w\rho(w)\big{(}V(w)-1-\mu\big{)}+\int_{\mathbb{R}}\mathrm{d}w\rho_{\mathsf{s}}(w)\Big{(}\frac{\rho}{\rho_{\mathsf{s}}}\log\frac{\rho}{\rho_{\mathsf{s}}}\Big{)}(w).$
(9.47)

The first summand can be viewed as energy-like and the second summand as the
entropy of $\rho$ relative to $\rho_{\mathsf{s}}$.

Provisionally we denote the minimizer of (9.47) by $\rho_{\mathsf{p}}$ which
satisfies the Euler-Lagrange equation

$V(w)-\mu+\log\frac{\rho_{\mathsf{p}}}{1+T\rho_{\mathsf{p}}}-T\frac{\rho_{\mathsf{p}}}{1+T\rho_{\mathsf{p}}}=0.$
(9.48)

In view of (3.61) and (3.49) this suggests to set

$\rho_{\mathsf{n}}=\frac{\rho_{\mathsf{p}}}{1+T\rho_{\mathsf{p}}},\qquad\rho_{\mathsf{n}}=\mathrm{e}^{-\varepsilon},$
(9.49)

which leads to

$\varepsilon=V-\mu-T\mathrm{e}^{-\varepsilon},$ (9.50)

in agreement with the TBA equation (3.56). Since
$\rho_{\mathsf{s}}=1+T\rho_{\mathsf{p}}$, the space density satisfies

$\rho_{\mathsf{s}}=\varsigma_{0}^{\mathrm{dr}},\quad\rho_{\mathsf{s}}\rho_{\mathsf{n}}=\rho_{\mathsf{p}}.$
(9.51)

In conclusion

$\lim_{\ell\to\infty}-\ell^{-1}\log
Z_{\mathrm{to},\mathsf{f},+,N}(\ell,\mu,V)=\mathcal{F}_{\mathrm{to},\mathsf{f},+}^{\bullet}(\rho_{\mathsf{p}})=-\int_{\mathbb{R}}\mathrm{d}w\mathrm{e}^{-\varepsilon(w)},$
(9.52)

in agreement with (9.11).

As explained in Section 3.5, the convergence of the DOS to a deterministic
limit follows on general grounds, but can also be confirmed by considering the
first order variational derivatives of the free energy. Returning to the
infinite-dimensional constraint below Eq. (9.42), evaluating it at
$\rho_{\mathsf{p}}$ reads

$0<1+\int_{\mathbb{R}}\mathrm{d}w^{\prime}\rho_{\mathsf{p}}(w^{\prime})\phi_{\mathrm{to}}(w-w^{\prime})=\rho_{\mathsf{s}}(w).$
(9.53)

For $P<P_{\mathrm{c}}$ the number density $\rho_{\mathsf{n}}>0$ and also
$\rho_{\mathsf{p}},\rho_{\mathsf{s}}$ stay positive. Hence (9.53) is
equivalent to $\nu>0$, a condition which was required already at the very
beginning of the argument. Thus, indeed, the minimizer remains unaffected by
the constraint.

The free energy for the Toda fluid has also been obtained by using scattering
coordinates in combination with a specifically designed external potential. We
presented two very different proofs for the generalized free energy and
arrived at the same answer, as it should be. Since the proofs are lengthy, one
might worry about missing factors of $2$ or so. Having two independent
derivations vouches for correctness. Furthermore, the availability of a second
tool opens up the possibility to handle the generalized free energy for a
larger class of integrable models.

$\blackdiamond\blackdiamond$ Local macroscopic conservation laws. By locally
conserved we mean here the identity

$\partial_{t}\mathfrak{n}(x,t;w)+\partial_{x}\mathfrak{j}(x,t;w)=0$ (9.54)

pointwise in $w$. If GHD is considered in a periodic interval, say $[0,\ell]$,
this implies the conservation law

$\partial_{t}\int_{0}^{\ell}\mathrm{d}x\,\mathfrak{n}(x,t;w)=0.$ (9.55)

According to (9.47), a physical example is the entropy with density

$s(x,t;w)=-\rho_{\mathsf{p}}(x,t;w)\log\rho_{\mathsf{n}}(x,t;w).$ (9.56)

If the solution to GHD is free of discontinuities, the entropy is locally
conserved.

Surprisingly, as a general property of GHD, beyond entropy a much larger
family of local macroscopic conservation laws is obtained in replacing the
$\log$ by some arbitrary function $g$, hence

$\mathfrak{n}(x,t;w)=(\rho_{\mathsf{p}}g(\rho_{\mathsf{n}}))(x,t;w).$ (9.57)

To verify one starts from

$\displaystyle-\partial_{t}(\rho_{\mathsf{p}}f(\rho_{\mathsf{n}}))=-(\partial_{t}\rho_{\mathsf{p}})f(\rho_{\mathsf{n}})-\rho_{\mathsf{p}}f^{\prime}(\rho_{\mathsf{n}})\partial_{t}\rho_{\mathsf{n}}$
$\displaystyle\hskip
39.0pt=\partial_{x}(v^{\mathrm{eff}}\rho_{\mathsf{p}})f(\rho_{\mathsf{n}})+\rho_{\mathsf{p}}v^{\mathrm{eff}}\partial_{x}f(\rho_{\mathsf{n}})=\partial_{x}(v^{\mathrm{eff}}\rho_{\mathsf{p}}f(\rho_{\mathsf{n}})),$
(9.58)

which implies

$\mathfrak{j}=v^{\mathrm{eff}}\rho_{\mathsf{p}}f(\rho_{\mathsf{n}}),$ (9.59)

as to be confirmed. $\blackdiamond\blackdiamond$

Notes and references

Section 9.1

A more detailed discussion of the Toda fluid is presented by Doyon [88, 89].
The high-low pressure domain wall is studied by Mendl and Spohn [204].

Section 9.2

The Toda lattice with boundary potential
$g_{1}(\mathrm{e}^{-q_{1}}+\mathrm{e}^{q_{N}})+g_{2}(\mathrm{e}^{-2q_{1}}+\mathrm{e}^{2q_{N}})$
is still integrable, see van Diejen (1995) for a proof. But the respective Lax
matrix is less accessible. Using scattering coordinates for the computation of
the generalized free energy is anticipated in Doyon [88], where the free
energy functional (9.42) can be found already. The identity (9.2) seems to be
a novel result. The large class of local conservation laws for the
hydrodynamic equations has been noted by Caux et al. [54]. I learned about
(9.2) through the master thesis of Deokule [76].

##  10 Hydrodynamics of soliton gases

In his seminal contribution Morikazu Toda realized that a chain with
exponential interactions allows for travelling waves, called soliton in
analogy to continuum nonlinear wave equations. The soliton is a coherent
motion of many particles which is spatially localized and travels at constant
velocity maintaining its shape. Two incoming solitons undergo an intricate
collision process, to eventually move outwards with velocities unchanged and
characteristic scattering shifts. But qualitatively this resembles the
dynamics discussed for Toda quasiparticles already. Apparently the Toda
dynamics possesses two distinct levels with soliton-like structure. To
distinguish, we refer to either soliton-based scheme or particle-based scheme.
The much studied prime example for a soliton-based scheme is the Korteweg-de
Vries (KdV) equation. This topic will be taken up first because of its
similarity with the soliton-based scheme of the Toda dynamics. Allowing for
both hydrodynamic schemes is shared by Toda and Ablowitz-Ladik lattice. For
other models only either one of the two schemes is known.

###  10.1 Soliton gas of the Korteweg-de Vries equation

The most widely studied soliton gas is based on the KdV equation, written in
its conventional form as

$\partial_{t}u+\partial_{x}(3u^{2}+\partial_{x}^{2}u)=0$ (10.1)

with $u$ a real-valued wave field over $\mathbb{R}$. Physically the Korteweg-
de Vries equation describes the motion of surface waves in a shallow channel
of water under gravity and in the approximation of small amplitudes. Solitary
waves have been observed experimentally back in 1834 and the description in
terms of a specific nonlinear wave equation was accomplished in 1895. While
the study of soliton solutions strongly indicated integrability, the infinite
hierarchy of local conservation laws was first obtained by Miura, Gardner, and
Kruskal in 1968.

As the Toda lattice, KdV admits for a Lax pair. Now $L$ is an operator,

$L=-\partial_{x}^{2}-u(t),$ (10.2)

while the partner operator reads

$B=-\big{(}4\partial_{x}^{3}+6u\partial_{x}+4\partial_{x}u\big{)}.$ (10.3)

Under the KdV time evolution

$\partial_{t}L=[B,L],$ (10.4)

a property signalling integrability. A different viewpoint is to consider, for
given $u$, the eigenvalue problem

$Lv=\lambda v$ (10.5)

and the evolution equation

$\partial_{t}v=Bv.$ (10.6)

This system has solutions only if $u(t)$ satisfies the KdV equation, which
thus plays the role of a consistency relation.

The lowest order conserved fields can still be obtained by hand with the
resulting densities

$Q^{[0]}(x)=u(x),\quad Q^{[1]}(x)=u^{2}(x),\quad
Q^{[2]}(x)=u^{3}(x)-\tfrac{1}{2}(\partial_{x}u)^{2}(x).$ (10.7)

Here we adopted that convention that $u^{n}$ is the term with no derivatives.
A more systematic tool comes from the underlying Poisson structure. The formal
Poisson bracket turns out to be

$\\{u(x),u(x^{\prime})\\}_{\mathrm{KV}}=\delta^{\prime}(x-x^{\prime}).$ (10.8)

More conveniently, functions on phase space are replaced by functionals of $u$
and its derivatives,

$F=\int\mathrm{d}xh(u(x),\partial_{x}u(x),\ldots,\partial_{x}^{n}u(x))$ (10.9)

with some polynomial $h$. Of course, more general functionals could be
considered. But for our purposes (10.9) will do. The gradient of $F$, denoted
by $\nabla F(x)$, is the functional derivative of $F$ with respect to $u$ at
$x$. The functionals $F,G$ have then the Poisson bracket

$\\{F,G\\}_{\mathrm{KV}}=\int_{\mathbb{R}}\mathrm{d}x\nabla
F(x)\partial_{x}\nabla G(x).$ (10.10)

In particular this Poisson bracket satisfies the Jacobi identity. Denoting the
total conserved fields by $Q^{[n]}=\int\mathrm{d}xQ^{[n]}(x)$, their
respective gradients are

$\nabla Q^{[0]}=1,\quad\nabla Q^{[1]}=2u,\quad\nabla
Q^{[2]}=3u^{2}+\partial_{x}^{2}u.$ (10.11)

Higher order conservation laws are obtained recursively through

$(2n-1)\partial_{x}\nabla
Q^{[n+1]}=(n+1)\big{(}u\partial_{x}+\partial_{x}u-\tfrac{1}{2}\partial^{3}_{x}\big{)}\nabla
Q^{[n]}$ (10.12)

The prefactors result from our convention in (10.7) upon disregarding the
operator $\partial_{x}^{3}$. In fact, this kind of recursion is familiar from
integrable quantum many-body systems and called boost operator.

For the hamiltonian we set $H_{\mathrm{kv}}=-Q^{[2]}$, which then generates
the time evolution as

$\partial_{t}u(x,t)=\\{u(x,t),H_{\mathrm{kv}}\\}_{\mathrm{KV}}.$ (10.13)

Hence $Q^{[2]}$ is called energy. Analogously, $Q^{[1]}$ generates a spatial
shift to the right with unit speed, hence called momentum. The KdV equation
possesses a two-parameter family of traveling wave solutions, called soliton,
more specifically one-soliton. Their analytic form is

$u_{\mathrm{s},1}(x,t;\eta,\delta)=2\eta^{2}\mathrm{sech}^{2}\big{(}\eta(x-4\eta^{2}t-\delta)\big{)}.$
(10.14)

Here $\eta$ is the spectral parameter, $\eta>0$, and $\delta$ is a shift,
$\delta\in\mathbb{R}$. The center of the soliton moves on the straight line
$x(t)=4\eta^{2}t+\delta$. Hence the soliton speed is $4\eta^{2}$. The taller
and narrower the shape the faster the soliton moves. The conserved fields take
the values $Q^{[0]}(u_{\mathrm{s},1})=4\eta$,
$Q^{[1]}(u_{\mathrm{s},1})=(4^{2}/3)\eta^{3}$, and
$Q^{[2]}(u_{\mathrm{s},1})=(4^{3}/5)\eta^{5}$.

We now start with localized initial data, in the sense that

$\int_{\mathbb{R}}\mathrm{d}x|x||u(x)|<\infty.$ (10.15)

For long times, $t\to\pm\infty$, the KdV solution splits into a soliton part,
denoted by $u^{\mathrm{s}}$, and dispersive part which decays to zero with
some power law. More precisely for any $\epsilon>0$,

$\lim_{|t|\to\infty}\sup_{\\{x||x|>\epsilon|t|\\}}|u(x,t)-u^{\mathrm{s}}(x,t)|=0.$
(10.16)

This behavior is reminiscent of charges coupled to the Maxwell field which
asymptotically acquire straight line motion through emitting radiation. Hence
the dispersive part is also called radiation. The soliton part consists of $N$
freely moving solitons, $N$ depending on the initial condition. In a formula

$u^{\mathrm{s}}(x,t)=\sum_{j=1}^{N}u_{\mathrm{s},1}(x,t;\eta_{j},\phi_{j}^{\pm}),$
(10.17)

where $\pm$ refers to $t\to\pm\infty$. Solitons behave like quasiparticles.
However for KdV there are only overtaking collisions. Thus for $t\to-\infty$
the incoming solitons are ordered as $x_{j}(t)<x_{j+1}(t)$ and the spectral
parameters as $\eta_{1}>\ldots>\eta_{N}>0$, while outgoing solitons travel in
reverse order. The label $N$ is used in analogy to Section 2.3, where $N$
denotes the number of quasiparticles (which equals the number of physical
particles).

The crucial dynamical information comes from the scattering shifts,
$\kappa_{j}$, which are written as

$\phi_{j}^{-}\hskip
6.0pt\mathrm{for}\,\,t\to-\infty,\qquad\phi_{j}^{+}=\phi_{j}^{-}+\kappa_{j}\hskip
6.0pt\mathrm{for}\,\,t\to\infty,$ (10.18)

where

$\kappa_{j}=\frac{1}{\eta_{j}}\sum_{i=1,i\neq
j}^{N}\mathrm{sgn}(\eta_{j}-\eta_{i})\phi_{\mathrm{kv}}(\eta_{j},\eta_{i})$
(10.19)

with the KdV scattering kernel

$\phi_{\mathrm{kv}}(\eta,\eta^{\prime})=\log\Big{|}\frac{\eta-\eta^{\prime}}{\eta+\eta^{\prime}}\Big{|}.$
(10.20)

The latter structure is already familiar from the hard rod gas and Toda
lattice. As hallmark of integrability, the scattering shifts are determined by
the two-soliton scattering shift. But in contrast to hard rods, the latter
shift depends on the incoming velocities. Note that
$\phi_{\mathrm{kv}}(\eta,\eta^{\prime})<0$ as for physical hard rods.

There are concise formulas for $N$-soliton solutions, which are defined by
requiring the radiation part to vanish identically. Considering the special
case $N=2$ we can illustrate how the identity in (10.20) is obtained. The two-
soliton solution can be written as

$u_{\mathrm{s},2}(x,t)=2\partial_{x}^{2}\log\varphi_{2}(x,t)$ (10.21)

with

$\varphi_{2}=1+A_{1}\mathrm{e}^{2\theta_{1}}+A_{2}\mathrm{e}^{2\theta_{2}}+A_{3}\mathrm{e}^{2(\theta_{1}+\theta_{2})},$
(10.22)

where

$\theta_{i}=\eta_{i}x-\omega_{i}t,\quad\omega_{i}=4\eta_{i}^{3},\quad i=1,2.$
(10.23)

The coefficients $A_{1},A_{2}$ are free parameters which can be converted into
$\phi_{i}^{\pm}$. But $A_{3}$ is determined by

$\log A_{3}=\log A_{1}+\log A_{2}+\phi_{\mathrm{kv}}(\eta_{1},\eta_{2}).$
(10.24)

In (10.19) the actual scattering shift carries the prefactor $1/\eta_{j}$.
This dependence can be traced to the solution ansatz (10.23), where in the
definition of $\theta_{i}$ the factor $\eta_{i}$ has not been bracketed out.

Given previous experience, one may boldly jump ahead to the hydrodynamic
equation for the KdV soliton gas. We introduce the counting function
$f(x,t;\eta)\mathrm{d}x\mathrm{d}\eta$ which is the number of solitons with
spectral parameter in the interval $[\eta,\eta+\mathrm{d}\eta]$ and location
in $[x,x+\mathrm{d}x]$ at time $t$. The counting function is governed by

$\partial_{t}f+\partial_{x}(v^{\mathrm{eff}}f)=0.$ (10.25)

For the effective velocity, as before, we assume the collision rate ansatz
which in our case becomes

$v^{\mathrm{eff}}(\eta)=4\eta^{2}+\eta^{-1}\int_{0}^{\infty}\mathrm{d}\eta\phi_{\mathrm{kv}}(\eta,\eta^{\prime})f(\eta^{\prime})\big{(}v^{\mathrm{eff}}(\eta^{\prime})-v^{\mathrm{eff}}(\eta)\big{)}.$
(10.26)

Still missing are some explanations on how the counting function $f$ is
related to local values of the conserved fields. In principle, the counting
function can be measured by the time-of-flight method. One considers some time
$t$ and a sufficiently large spatial box, still small on the scale of
macroscopic variations. The wave field $u(x)$ in this spacetime cell is cut
out and then smoothly interpolated with the zero-field all the way to
infinity. The such constructed field $u_{\mathrm{b}}(x)$ is evolved under the
KdV dynamics. For sufficiently long times the solution consists of many
spatially well separated one-solitons. In approximation, their statistics is
that of an ideal gas with uniform density
$\bar{\varrho}=\int_{\mathbb{R}_{+}}\mathrm{d}\eta f(\eta)$. Each soliton
carries a spectral parameter, $\eta_{j}$, and these form a family of i.i.d.
random variables with single point probability density
$h(\eta)=f(\eta)/\bar{\varrho}$. To determine the average value of the
conserved field $Q^{[n]}$ under such statistics, the positional average yields
the integral over $Q^{[n]}$ evaluated at $u_{\mathrm{s},1}(x,0;\eta,0))$ and
the i.i.d. average results in the integral over $\eta$ weighted with
$h(\eta)$, in formulas

$\int_{0}^{\infty}\mathrm{d}\eta
h(\eta)\bar{\varrho}\int_{\mathbb{R}}\mathrm{d}xQ^{[n]}(u_{\mathrm{s},1}(x,0;\eta,0))=\frac{2^{2n+2}}{2n+1}\int_{0}^{\infty}\mathrm{d}\eta
f(\eta)\eta^{2n+1}.$ (10.27)

Thereby we connected the average of the $n$-th conserved field with the
$(2n+1)$-th moment of the counting function. The numerical coefficients are
known as Kruskal series. Since $f$ vanishes on $\mathbb{R}_{-}$, it is
uniquely determined by its moments (10.27).

A natural question is to have a more intrinsic characterization of the
statistics of the wave field $u_{\mathrm{b}}(x)$ in the prescribed cell. This
is so to speak an inverse scattering problem and not much is known. For
particle-based hydrodynamics the statistics of the positions and momenta of
the particles in the cell should be given by a GGE with properly adjusted
parameters. At least for the Toda fluid we know that, as regards to the
conserved fields, the data agree with those from the time-of-flight method.
Also the Lax filter could be employed, compare with Section 3.2. For KdV the
corresponding construction would involve bound state eigenvalues of the
Schrödinger operator (10.2) with potential given by the box field
$u_{\mathrm{b}}(x)$.

$\blackdiamond\blackdiamond$ TBA revisited. The Toda lattice consists of
nonrelativistic particles with energy-momentum relation $E(P)=P^{2}/2$. For
relativistic models the energy-momentum relation would have to be modified.
Also for solitons there is no particular reason to have a quadratic energy-
momentum relation. Therefore TBA has to be extended to cover such cases.

In a hamiltonian context, if the energy-momentum relation $E(P)$ is given,
then the velocity is determined by $v=\partial E/\partial P$. In case of the
KdV equation the one-soliton momentum and energy are given in terms of the
spectral parameter $\eta$. To stay in context we still use $\eta$, which
however is now a parameter depending on the model under study. The functions
$E(\eta)$ and $P(\eta)$ are prescribed. Then, with prime referring to
parameter differentiation, the velocity becomes $v=E^{\prime}/P^{\prime}$.
Through the interaction between quasiparticles $E^{\prime},P^{\prime}$ become
dressed and, appealing to the fundamental relation (6.20), one argues that

$v^{\mathrm{eff}}=\frac{(E^{\prime}){{}^{\mathrm{dr}}}}{(P^{\prime}){{}^{\mathrm{dr}}}}\,,$
(10.28)

the dressing transformation being defined as before. Of course, for
$E(P)=P^{2}/2$ the original relation is recovered.

Starting from the collision rate ansatz, it is now reformulated as

$v^{\mathrm{eff}}=\frac{1}{P^{\prime}}\big{(}E^{\prime}+T(\rho_{\mathsf{p}}v^{\mathrm{eff}})-(T\rho_{\mathsf{p}})v^{\mathrm{eff}}\big{)},$
(10.29)

where $T(\eta,\eta^{\prime})$ is some symmetric kernel. As before,
$\rho_{\mathsf{p}}$ is the particle density. The number density is defined by

$\rho_{\mathsf{n}}=\frac{\rho_{\mathsf{p}}}{P^{\prime}+T\rho_{\mathsf{p}}}.$
(10.30)

Then the dressing operation for some function $f$ is still

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

with $\rho_{\mathsf{n}}$ regarded as multiplication operator. The inverse to
(10.30) reads

$\rho_{\mathsf{p}}=\rho_{\mathsf{n}}P^{\prime\,\mathrm{dr}},$ (10.32)

which identifies the space density as

$\rho_{\mathsf{s}}=P^{\prime\,\mathrm{dr}}.$ (10.33)

By a computation similar to the one in (6.1), the effective velocity is indeed
given by (10.28).

Applying our scheme to KdV solitons, one finds $P(\eta)=\eta^{2}/2$ and
$E(\eta)=\eta^{4}$, which differ from the one-soliton momentum and energy.
$\blackdiamond\blackdiamond$

###  10.2 Soliton gas of the Toda lattice

We have already studied the hydrodynamic scale of the Toda lattice within the
particle-based scheme. Since the Toda lattice has also solitons in the more
conventional sense, one would expect to have also a hydrodynamic scale for the
soliton-based scheme. This is the topic to be explored now. Of course, which
scheme to use depends on the choice of initial conditions. For sure, there
will be also initial conditions not covered by either scheme.

The dynamical evolution takes place on the infinitely extended lattice with
equations of motion

$\frac{\mathrm{d}^{2}}{\mathrm{d}t^{2}}r_{j}(t)=2\mathrm{e}^{-r_{j}(t)}-\mathrm{e}^{-r_{j-1}(t)}-\mathrm{e}^{-r_{j+1}(t)}$
(10.34)

where $j\in\mathbb{Z}$. We consider solutions for which
$|p_{j}(t)|,|r_{j}(t)|\to 0$ as $j\to\pm\infty$. In Flaschka variables, the
two-parameter family of one-soliton solutions is given by

$a_{j}^{2}(t;\eta,\delta)-1=\omega^{2}\mathrm{sech}^{-2}(|\eta|j-\omega
t-\delta)$ (10.35)

and

$p_{j}(t;\eta,\delta)=|\sinh\eta|\big{(}\tanh(|\eta|j-\omega
t-\delta)-\tanh(|\eta|(j-1)-\omega t-\delta)\big{)}$ (10.36)

with

$\omega(\eta)=\sinh\eta.$ (10.37)

Deviating from the previous convention the soliton travels along the
trajectory $j(t)=|\eta|^{-1}\omega t+|\eta|^{-1}\delta$. The spectral
parameter is $\eta$, $\eta\in\mathbb{R}$, the shift $\delta\in\mathbb{R}$, and
the soliton velocity

$v_{\mathrm{s}}=|\eta|^{-1}\omega.$ (10.38)

For $\eta>0$ the soliton is moving to the right and for $\eta<0$ to the left.
Note that the velocity is odd and monotone increasing in $\eta$, however with
a jump of size $2$ at the origin. Necessarily the soliton speed is supersonic.

The construction of the two-soliton relies on the discrete analogue of (10.21)
and the same ansatz as in (10.22), using $\theta_{i}=|\eta_{i}|-\omega_{i}t$.
This results in the Toda scattering kernel denoted by
$\phi_{\mathrm{tos}}(\eta_{1},\eta_{2})$. One has to distinguish overtaking
collisions, $\eta_{1}\eta_{2}>0$ and head-on collisions, $\eta_{1}\eta_{2}<0$.
Explicitly,