Generalized hydrodynamics of the KdV soliton gas

We establish the explicit correspondence between the theory of soliton gases in classical integrable dispersive hydrodynamics, and generalized hydrodynamics (GHD), the hydrodynamic theory for many-body quantum and classical integrable systems. This is done by constructing the GHD description of the soliton gas for the Korteweg-de Vries (KdV) equation. We further predict the exact form of the free energy density and flux, and of the static correlation matrices of conserved charges and currents, for the soliton gas. For this purpose, we identify the solitons' statistics with that of classical particles, and confirm the resulting GHD static correlation matrices by numerical simulations of the soliton gas. Finally, we express conjectured dynamical correlation functions for the soliton gas by simply borrowing the GHD results. In principle, other conjectures are also immediately available, such as diffusion and large-deviation functions for fluctuations of soliton transport.


Introduction
Long wavelength, hydrodynamic theories play a prominent role in physics, from fluids to optics, condensed matter to quantum mechanics, and beyond. Recently, two new hydrodynamic theories have come to the fore, which present a large amount of conceptual and formal similarities, and which have been extremely successful in their respective fields of study: the theory of soliton gases [45], and generalised hydrodynamics [32].
On the one hand, the theory of soliton gases is concerned with the emergent, large-scale behaviours of integrable partial differential equations of dispersive hydrodynamics, such as the Korteweg-de Vries (KdV) and the nonlinear Schrödinger (NLS) equations. Dispersive hydrodynamics theories are of spectacularly different character than their dissipative counterparts. For instance, dispersive shock waves consist of coherent, rank-ordered, nonlinear oscillations that continually expand. Remarkably, the multi-scale nonlinear dynamics of dispersive shock waves are asymptotically described by a new, emergent hydrodynamic theory, arising as a result of averaging over fast oscillations-the famous Whitham modulation theory [86]. The Whitham equations are a system of first-order quasilinear PDEs describing slow evolution of the wave's parameters such as amplitude, wavelength, mean etc. For integrable equations such as KdV or NLS, the Whitham modulations are expressed in terms of slow deformations of hyperelliptic Riemann surfaces associated with spectral finite-gap solutions that locally approximate solutions to initial-value problems for the original dispersive hydrodynamics [64,49].
This new hydrodynamics is in fact of a much more general character. Many physically relevant wave phenomena exhibit complex spatiotemporal behaviours that cannot be restricted to macroscopically coherent dispersive shock waves. In spite of integrability, the inherent randomness of many real-life systems (due to initial and boundary conditions or to complex GHD equations with external force have been confirmed in experiments on cold atomic gases [77,69,68]. GHD is now an extremely active field of research. The kinetic equation for soliton gases, and the (Euler-scale) GHD equation for many-body integrable models, are strikingly similar, something already pointed out both in the contexts of GHD and dispersive hydrodynamics, see e.g. [12,16,38,32,45]. Much like for soliton gases, in GHD it can be seen as an effective propagation equation for the stable quasi-particles of many-body integrable systems, as advocated in [8]. Semi-classical arguments made in [38] and [9] further argue that the GHD equation can be recovered using Zakharov's 1971 idea that scattering shifts lead to modified velocities. This idea is also at the root of the early derivation in the 1970's of the hydrodynamic equations for the hard rod gas [1], which GHD generalises. But besides these remarks, the precise mathematical relation between the kinetic theory of soliton gases in dispersive hydrodynamics, and GHD, has not been established so far.
The main purpose of this paper is to establish an explicit connection between the spectral kinetic theory of soliton gases in dispersive hydrodynamics, and the generalised hydrodynamics of integrable many-body systems. We achieve this in two steps.
First, we develop the GHD description of the KdV soliton gas. This involves making the explicit relation between quantities defined for soliton gases, and those within the context of GHD. In order to make this relation as clear as possible, the respective conventional notations for soliton gases and GHD are kept troughout, and we provide a dictionary between them.
Second, and perhaps most importantly, we give evidence to the fact that the (generalised) thermodynamics of the soliton gas can be described in terms of the classical Thermodynamic Bethe Ansatz (TBA). That is, the general form expressed in [32] for the thermodynamics of integrable many-body systems holds for soliton gases as well. In this, a piece of information about the gas is required, that was not discussed yet in the literature on soliton gases: the "statistics" of the solitons. We find that solitons are associated to the Maxwell-Boltzmann statistics, like classical indistinguishable particles. As a result, we identify the key thermodynamic quantities associated with the KdV soliton gas (including the entropy, the free energy, the temperature). Using the general results of GHD, we then infer the static covariances, the space-integrated two-point correlation functions involving densities and currents. The general statement that ensembles of solitons should admit the Maxwell-Boltzmann statistics was first proposed in [6], from previous results on the sine-Gordon model obtained by taking appropriate classical limits of the quantum TBA. However, here we provide the first (non-rigorous) derivation directly for the classical soliton gas. We also confirm for the fist time the conjecture of the Maxwell-Boltzmann statistics, by comparing static GHD covariances against numerical simulations of KdV soliton gases.
The dictionary between the KdV soliton gas and GHD, and the uncovering of the thermodynamics of the KdV soliton gas, allows us to use results in one field in order to strengthen the other. For instance, the derivation of the KdV soliton gas kinetic equations, from the spectral theory of the KdV equation, gives a mathematically stronger basis for the GHD equation than the principle of local entropy maximisation used up to now to justify it in many-body integrable systems. Also, the various GHD results can now be applied to soliton gases, such as the exact form of dynamical correlation functions at the Euler scale.
The paper is organised as follows: In Section 2 we provide a brief introduction to the theory of soliton gases from the point of view of dispersive hydrodynamics, putting the emphasis on the KdV equation. Section 3 introduces GHD by way of a paradigmatic case study, through its direct application to the Lieb-Liniger model. Section 4 draws explicit parallels between the two aforementioned theories, making the relations between the relevant quantities and notations in both fields. Section 5 builds on the previous one by straightforwardly exploiting results from GHD to construct the thermodynamics of the KdV soliton gas. Section 6 contains a summary of our results and some concluding remarks. The paper is complemented by three appendices providing details of our numerical simulations.

Soliton gas in dispersive hydrodynamics
We start our discussion by introducing the theory of soliton gases in the context of integrable dispersive hydrodynamics. More specifically, we will focus on how it relates to the KdV equation through its N -soliton solutions and further, the thermodynamic limit of finite-gap potentials.

Basic construction
Let us consider the KdV equation in the form This equation (2.1) belongs to the family of completely integrable equations and, for a broad class of initial conditions, its integrability is realised via the inverse scattering transform (IST) method [52]. The inverse scattering theory associates a soliton of the KdV equation with a point of discrete spectrum λ = λ n , of the Schrödinger operator Assuming u → 0 as x → ±∞, the KdV soliton solution corresponding to λ i = −η 2 i , η i > 0, is given by where we denote by a i = 2η 2 i the soliton amplitude, by s i = 4η 2 i its speed, and by x 0 i its initial position or "phase". Note that solitons have finite width ∼ 1/η i , which affects the notion of interaction range, particularly for small-amplitude solitons. In what follows we will be referring to η as a spectral parameter with the understanding that η = √ −λ. Along with the simplest single-soliton solution, the KdV equation supports N -soliton solutions u N (x, t) characterised by N discrete spectral parameters η 1 < η 2 < · · · < η N and the set of initial positions {x 0 i |i = 1, . . . , N }. The integrable structure of the KdV equation has profound implications for the dynamics of soliton interactions.
2. The collision of two solitons with spectral parameters η i and η j , i = j results in their phase (position) shifts given by (2.4) so that the taller soliton acquires shift forward and the smaller one -shift backwards.
3. Solitons interact pairwise, i.e. the resulting phase shift ∆ i of a given soliton with spectral parameter η i after its interaction with M solitons with parameters η j , j = i, is equal to the sum of the individual phase shifts, It is important to stress that the collision phase shifts are far-field effects. Mathematically they are artefacts of the asymptotic representation of the exact two-soliton solution of the KdV equation in the form of a sum of two individual solitons: u 2 (x, t; η 1 , η 2 ) u s (x + ∆ 12 , t; η 1 ) + u s (x + ∆ 21 , t; η 2 ), which is only valid if solitons are sufficiently separated (the long-time asymptotics). The interaction of solitons is a complex nonlinear process [65] and the resulting wave field in the interaction region cannot be represented as a superposition of the phase-shifted one-soliton solutions.
We first introduce soliton gas phenomenologically, as an infinite random ensemble of solitons distributed on R with some non-zero spatial density β, characterised by certain distributions over the spectral parameter η i ∈ Γ = [η min , η max ] ⊂ R + and the phase x 0 j ∈ R. Without loss of generality one can assume Γ = [0, 1], but we shall keep the general notation for the time being. Assuming macroscopic equilibrium (spatial homogeneity), the spectral density of states (spectral DOS) f (η) of soliton gas is introduced in such a way that f (η 0 )dη gives the number of solitons with the spectral parameter η ∈ [η 0 ; η 0 + dη] contained in the portion of soliton gas over a unit interval of x ∈ R (the individual solitons can be counted by "cutting out" the relevant portion of the gas and letting them separate). The corresponding spectral flux density v(η) represents the temporal counterpart of the spectral DOS, i.e. v(η 0 )dη is the number of solitons with the spectral parameter η ∈ [η 0 ; η 0 + dη] crossing any given point x = x 0 per unit interval of time. These definitions are physically suggestive in the context of rarefied soliton gas where solitons are identifiable as individual localised wave structures. The total spatial density of the soliton gas For a rarefied gas, β 1, the phases x 0 i can be assumed to follow Poisson distribution on R with density β. In such a gas, the total spatial shift of a "tracer" soliton with spectral parameter η = η 1 (we shall call it η 1 -soliton) over the time interval dt, due to the interactions with µsolitons, µ ∈ Γ, is given by ∆ 1 ≈ Γ [∆(η 1 , µ)|s 0 (η 1 ) − s 0 (µ)|f (µ)dµ]dt, where s 0 (η) = 4η 2 is the speed of a free, non-interacting soliton. This simple argument was used by Zakharov in 1971 [89].
It turns out that a straightforward generalisation of Zakharov's construction to the case of a dense gas can be made. This can be formulated as the collision rate ansatz, whereby the total position shift of the η 1 -soliton, due to soliton collisions in a gas with DOS f (µ) over the time interval dt, is given by where s(η) is the effective velocity of the 'tracer' soliton with spectral parameter η. In simple terms (2.7) represents an extrapolation of the rarefied gas properties to a dense gas, realised by replacing s 0 (η) → s(η) in the collision rate expression. It is important to stress that the validity of (2.7) for a dense soliton gas is far from being obvious. Indeed, as we mentioned, the very notion of the phase shift in classical soliton theory is only applicable in the context of the long-time asymptotics, i.e. when solitons have sufficient time to separate from each other after the interaction, which can only happen in a rarefied gas. In other words, the interactions between solitons can be treated as short-range only if the typical distance between solitons in a gas is much greater than soliton's width. The assumption that, in a dense gas, where solitons experience significant overlap and continual interaction, the net effect of soliton collisions on the mean velocity is expressed in the same way as in the rarefied gas requires justification. Such a justification has been provided in [43] in the framework of the finite-gap spectral theory [7].
The collision rate ansatz (2.7) implies the integral equation for the effective velocity of a tracer soliton: where the (x, t)-variations of f and s occur on macroscopic, hydrodynamic scales, much larger than the typical scales associated with variations of the wave field u(x, t) in individual solitons. Now, isospectrality of the KdV evolution within the IST framework implies the conservation equation which, together with the equation of state (2.8), provides a spectral hydrodynamic/kinetic description of the KdV soliton gas. As we will now discuss, those last two equations can be rigorously recovered even in the context of a dense gas, in which the collision rate ansatz should a priori not apply.

Nonlinear dispersion relations for soliton gas and spectral kinetic equation
As shown in [43] (see also [45]), equation (2.8) can be derived in the framework of the multiphase modulation (Whitham) theory [86], [49] without making any assumptions about the gas' density or the collision rate. The derivation in [43] makes an extensive use of the spectral properties of the so-called finite-gap potentials [7], the quasi-periodic generalisations of N -soliton solutions of the KdV equation exhibiting the band Lax spectrum, where the spectral bands γ j ⊂ R, j = 1, 2, . . . , N are finite intervals and γ N +1 -a semi-infinite interval. The N -soliton limit of an N -gap solution is achieved by collapsing all the finite bands γ j into double points corresponding to the solitonic spectral values λ j .
It was proposed in [43] that the soliton gas can be described by the thermodynamic limit of the spectral N -gap KdV solutions, achieved by assuming a special band-gap distribution (scaling) of the spectral set S N for N → ∞ on a fixed interval (say λ ∈ [−1, 0]) whereby the finite spectral bands are required to be exponentially narrow compared to the gaps ∼ 1/N between them. Such spectral scaling enables the soliton limit of finite-gap potentials as N → ∞ while preserving finite spectral DOS.
Within the spectral thermodynamic limit approach, the spectral DOS f (η) and the spectral flux density v(η) of an equilibrium (uniform) soliton gas are defined via the nonlinear dispersion relations. The nonlinear dispersion relations for soliton gas are derived by applying the thermodynamic spectral limit to the wavenumber-frequency relations for the finite-gap KdV solutions [49] (see [43,45] for details). Let η = √ −λ, then the spectral DOS f (η) and the corresponding spectral flux density v(η) satisfy Here σ(η) ≥ 0 is the spectral scaling function characterising the Lax spectrum of the soliton gas while Γ is the support of f (η) and v(η). Eliminating σ one obtains the equation of state (2.8) for the effective velocity of a "tracer" soliton in the gas, s(η) = v(η)/f (η). The derivation of the equation (2.8) via the thermodynamic limit thus provides the justification of the collision rate assumption (2.7) for the KdV soliton gas. The continuity (kinetic) equation (2.9) for slowly varying spectral DOS f (η; x, t) in a non-equilibrium (weakly non-unifiorm) soliton gas is obtained as a thermodynamic limit of kinematic modulation equations expressing the "conservation of waves" [43,45]. A similar derivation has been recently performed in [44] for the focusing NLS equation where the spectral parameter λ is complex. We note that the spectral support Γ should not necessarily be a fixed simply connected set: the general case when Γ is given by a union of disjoint intervals whose endpoints are allowed to vary in space-time, has been considered for the KdV gas in [20]. The equation for the evolution of σ(η; x, t) in a non-equilibrium soliton gas follows from (2.10), (2.8) and (2.9). A direct computation yields for a fixed Γ ∂ t σ + s∂ x σ = 0, (2.12) which is equivalent to the spectral kinetic equation. Thus the function σ(η, x, t) plays the role of the Riemann invariant in the soliton gas kinetics/hydrodynamics. Furthermore, it can be shown that the following conditions are satisfied [15]: (2.14) Relations (2.13) and (2.14) are the continuum limit counterparts of the semi-Hamiltonian (integrability) and linear degeneracy properties respectively of the hydrodynamic reductions of the kinetic equation (2.9), (2.8), established in [42], see also [12,48] for further developments on integrability of hydrodynamic reductions. Finally we note that the limits σ → ∞ and σ → 0 correspond to the special cases of soliton gas: the ideal (noninteracting) soliton gas (σ → ∞) and the soliton condensate (σ → 0). In the latter case, equations (2.10) and (2.11) can be integrated explicitly yielding the critical DOS f c (η) and the corresponding spectral flux v c (η) [20]-see Appendix A.

Conserved quantities
One of the fundamental properties of integrable dispersive hydrodynamics is the availability of an infinite set of local conservation laws ∂ t P n + ∂ x Q n = 0, i = 0, 1, 2, . . . , (2.15) where the P i and Q i are functions of the field variable u and its derivatives (the so-called local polynomial functionals). Of particular interest are the first three conserved KdV densities typically associated with conservation of "mass", "momentum" and "energy". Let u(x, t) describe the random KdV wave field in a soliton gas -the "integrable turbulence" [90]. In the spirit of the Whitham modulation theory [86] we invoke scale separation and, assuming ergodicity of the soliton gas, apply the ensemble averaging to the KdV conservation laws (2.15) to obtain the modulation system (note however that the original Whitham theory uses spatial averaging) where, for the first three conserved densities, we have We should remark that the specific coefficients C n in (2.19) are relevant in the context of the particular normalisation (2.16) of the polynomial KdV integrals. However, the resulting conservation system The second constraint, of purely spectral nature, arises due to the requirement of nonnegativity of the spectral scaling function σ(η) in the nonlinear dispersion relations (2.10), (2.11), hence

Thermodynamics and large-scale correlation functions
The main results of this paper, besides making the dictionary between the solition gas and generalised hydrodynamics, is to provide new conjectures for exact results in the soliton gas. These results are concerned with the thermodynamics and large-scale correlation functions, which have not been studied yet in soliton gases. We therefore conclude this section by providing an early heads up for the dispersive hydrodynamics readership to these main results, without the need for a full understanding of GHD. Namely, we present expressions for the temperature, free energy and correlation functions in the KdV solilton gas in terms of the DOS or, equivalently, the spectral scaling function σ(η). The full discussion detailing the GHD conjectures necessary to explicitly derive the expressions below can be found in Section 5. The relevant GHD tools are introduced in Sections 3 and 4.
For the sake of simplicity we shall focus here on expressing the spectral density for thermal states, and on expressing correlations of the KdV field u(x, t) in generic states; generalised Gibbs ensembles and correlations of other densities P n can also be accessed.
The first result is the spectral scaling function for a thermal state. The definition of a thermal state for the KdV equation, or for a soliton gas, is a subtle question. In this paper we take a natural definition based on the soliton gas viewed as the limit of large soliton numbers, and large lengths, of ensembles of solitons. An ensemble of N soliton on an interval of length L is a ensemble of configurations of the KdV field, supported on the interval (and quickly tending to zero beyond it), that decomposes asymptotically at large times into a train of N separate solitons. We distribute these solitons according to a Boltzmann weight, with the temperature being the parameter associated to the Hamiltonian, P 2 (x)dx. See Subsection 5.1. With this definition, the spectral scaling function for a thermal state in terms of the inverse temperature β is In fact for a general spectral function (not necessarily that of a thermal state), the inverse temperature is defined more generally as β = 1 768 At this level of generality, one may also evaluate the entropy (this is the analogue of the Yang-Yang entropy of quantum integrability [87]), which takes the simple form reaching its minimal value in the condensate limit σ → 0 alluded to in Section 2.2.
Finally we can evaluate correlations and fluctuations of the field: the static covariance is given by 27) and the dynamical correlation function of the KdV field at the Euler scale of large position and time separation is where η * (ξ) = {η : s(η) = ξ}.

Introduction
In the theory of Generalised Hydrodynamics (GHD), a rather different viewpoint is taken from that described in the previous section. The idea is to directly apply the general principles of hydrodynamics (we concentrate on Euler hydrodynamics here) to many-body integrable systems, in order to predict their dynamics away from equilibrium. In the simplest expression of GHD, we consider a system, composed of many particles in interaction, which has the property of integrability. The system can be either quantum or classical, but the theory was initially developed to discuss the quantum case in the original works [18] and [8]. This includes the Lieb-Liniger model, considered in [18], which will be our main example. Additionally, there is no need to restrict to systems of particles, and the Heisenberg spin chain was considered in [8].
In GHD, we assume that the system is initially in a fluctuating state which is, at every point in space, very well described by a state of Gibbs form, where entropy is maximised with respect to all available conservation laws. In integrable models these are the generalised Gibbs ensembles (GGEs), widely studied [40,46,60,62]. For instance, the initial state could be a profile of temperature, where at each point the system is thermal; but it could be something more general, where many conserved charges are involved. The system is taken to be infinite, although GHD with external force fields, which may constrain the system to a finite length, has also been developed. The fluctuations of the full system may be with respect to any ensemble (micro-canonical, canonical, grand-canonical, and more), as this does not affect the description of local states.
Then, the evolution of this state is obtained by applying the hydrodynamic principles: the state is assumed to be, throughout time, well described by local Gibbs-like states, and the dynamics is determined by the continuity equations for all conservation laws admitted by the model.
As usual, the above program requires knowledge of the thermodynamics -all the "equilibrium" states to which the system locally relaxes, in integrable models these are the GGEs -and of the equations of state -how the average "currents" are related to the average "densities", these being the objects involved in the local conservation laws. Solving the latter question was the main technical achievement of the original works [18,8].
In this program, no consideration is made of the explicit microscopic dynamics: one only needs thermodynamic quantities, and one uses hydrodynamic principles to conjecture an emergent evolution equation. The result is supposed to be valid whenever changes of all local probes occur on large enough scales in space and time.
In order to illustrate these ideas, we will consider the Lieb-Liniger model [67] (which can be viewed as a quantum nonlinear Schrödinger equation [76]). This is a system of N bosons with repulsive delta interactions described by the Hamiltonian where c is a positive constant parametrising the strength of interactions. This system is integrable in the sense that its Hamiltonian is part of a family of conserved quantities (also called charges) in involutionĤ This includes the number of particles and momentum, Higher conserved charges have a more complicated form [24,25]. Their explicit expressions are not important, except for noting that they are deformations of the natural conserved charges of free-particle models, the higher-power of momenta: HereV n (x i − x 1 , . . . , x i − x N ) is a differential operator with derivatives up to a maximal order n − 1, and with coefficients that are functions of x i − x 1 , . . . , x i − x N . It does not depend on x j whenever |x i − x j | is large enough, and tends to zero when |x i − x j | tend to infinity for all j = i. The normalisation in (3.4) is purely conventional. The existence of this set of conserved charges ensures that many-particle scattering is elastic (it preserves all momenta), and factorizes into separate two-body scattering events [2]. A crucial notion is that of extensivity of the chargesQ n : locality or quasi-locality of their densities [61]. That is, we may writeQ where the densityq n (x, 0) is an operator "supported" at, or on some finite region around, the point x. A full and accurate notion of support would need more care, but, as a part of this notion,q n (x, 0) should commute withq 0 (x , 0) = i δ(x − x i ) for |x − x | large enough. By conservation of the charges ∂ tQn = 0, these satisfy the microscopic conservation laws for some local or quasi-local currentsĵ n (x, t).
We will start by discussing the thermodynamic limit of such a model and how integrability impacts the properties of its stationary states. To that end we will rely on the celebrated Thermodynamic Bethe Ansatz (TBA). This was initially developed to specifically deal with the Lieb-Liniger model [87], but it turns out that its applications are significantly more far-reaching [50]. We will then use these results to construct a hydrodynamic theory of integrable systems according to [32]. As we shall see in the next section, parallels are to be drawn between quantum particles and solitons.
Note that, in the following, we consider systems featuring only one type of particle but all of this can easily be extended as discussed for instance in [30].

Stationary states
One of the most important principles of statistical mechanics is that the set of all states to which a many-body system may relax, as seen by local probes, is expected to be the set of all stationary, clustering states. A clustering state is one where local observables have vanishing correlations at large distances; clustering state are "ergodic" (in some sense). In typical, nonintegrable gases, a so-called "ergodic principle" says that these are the Gibbs states and their Galilean (or relativistic) boosts, with density matrix ρ ∝ e −β(Ĥ−µN +νP ) , Trρ = 1 . (3.7) Thus, in the inertial frame, these are equilibrium states. In integrable models, they may involve any combination of the local chargesρ ∝ e −Ŵ , (3.8) for any extensive conserved quantityŴ , typically of the form These are the GGEs. This in fact implies that stationary, clustering states include those which carry nonzero currents which cannot be made to vanish by simple boosts, and thus are nonequilibrium steady states. Extensivity ofŴ guarantees that the resulting state is clustering. As the Hamiltonian is homogeneous, the stationary, clustering states are invariant under space-time translations. An equivalent description of the stationary, clustering states (3.8) is as those which maximise entropy with respect to the conserved charges. The β n 's are the generalised "inverse temperatures" associated to the conserved charges Q n (the Lagrange multipliers for the constraints of conserved charges).
The GGEs are most efficiently described via the TBA.

Bethe Ansatz
Consider a N -particle eigenstate of the Hamiltonian (3.1) on a circle of length L. As the Hamiltonian is free whenever coordinates are disjoint, in each region of a given ordering, the wave-function is a solution of the free Schrödinger equation. Thus it can simply be written as a superposition of plane waves where the set of "quasi-momenta" p i , no two of which are identical [87], will be determined by the periodicity conditions. The sum is over the group of all permutations P of N elements; P represents permutations of the quasi-momenta. Q is also a permutation of N elements; it is not summed over, but rather is fixed by the positions x 1 , . . . , x N (thus it is a function of the positions) and represents the "particle order". By definition, it is the permutation satisfying Amplitudes are permutation invariant, in the sense that A P (Q) = A P•P (Q • P ) for all permutations P ; this in fact guarantees that the resulting wave-function is indeed bosonic, i.e. invariant under permutations of its arguments (the particles' positions). Amplitudes are fixed by the delta-function interaction, and their ratios give rise to scattering phases, as the scattering problem can be obtained from the limit L → ∞. With P differing from P only through permutation of quasi-particles i and j, the ratio involves the two-particle phase shift φ(p i , p j ), The above form of the wave function implies that three-body and higher processes factorise into two-body processes [2]. This is the blueprint for every Bethe-ansatz integrable model. If we have N particles on a circle, the wave function must be periodic. Assuming the particles did not interact, the wave function would simply acquire a factor e ip i L if we were to take one particle of momentum p i around the circle and back to its position, leading to the usual quantization conditions Taking into account the interactions, however, the same particle will have to scatter through all the others, accumulating phase-shifts and yielding the Bethe equations [66] where I i ∈ Z. Indeed, if φ = 0, without interactions, equation (3.13) becomes equivalent to equation (3.12), of which it can be considered an extension. As mentioned, the Lieb-Liniger model admits an infinite number of conserved charges. Although these are difficult to construct in generality, as they are local deformations of higher powers of momenta, their eigenvalues can be evaluated by acting on asymptotic states, where particles are well-separated and do not interact; it turns out that the finite-L effects on the wavefunction (3.10) identically vanish, hence these are the correct eigenvalues in general. Eigenvalues of chargesQ n on eigenstates |{p i } take the form where h n (p i ) is the amount of charge Q n carried by particle i. In the case of the Lieb-Liniger model, we have the total number of particlesQ 0 = N 1, the total momentumQ 1 = −i i ∂ ∂x i , and the total energyQ 2 =Ĥ LL , with More generally eigenvalues take the form Q n ∝ i p n i . As the momentum operatorQ 1 =P and the energyQ 2 =Ĥ LL play an important role, we give their one-particle eigenvalues special names: TBA can be written in universally applicable forms using these objects.

Thermodynamics
We are looking to take the thermodynamic limit of the states (3.8), in the sense that both L → ∞ and N → ∞ while N/L is kept constant, In fact it will be simpler to allow N to fluctuate and keep N /L constant; thusQ 0 is included in the exponential in (3.8). We can characteriseρ by its eigenvalues on the eigenstates |{p i } : The function w(p) replaces the Lagrange multipliers, as There is no explicit need for the Lagrange multipliers {β n }, the GGE is fully fixed by the spectral function w(p). This defines a good GGE under the condition that it grows fast enough as |p| → ∞, and is bounded from below. The set {β n }, and the function w(p), form systems of coordinates in the manifold of maximum entropy states. The thermodynamic limit of the Lieb-Liniger model was first investigated in [87]. For a generic state, we can approximate the distribution of momenta by a spectral, or phase space, density ρ p (p), which is such that ρ p (p i ) ≈ (L(p i+1 − p i )) −1 . Defining a density of holes by ρ h (p i ) = (I i+1 − I i )ρ p (p i ), the Bethe equations (3.13) become [84] 2πρ s (p) = P (p) − dp ϕ(p, p )ρ p (p ) , is the total density or "asymptotic space density" (a notion that will become clear later), where P = dP/dp is the derivative of the momentum eigenvalue (= 1 in the Lieb-Liniger model), and ϕ(p, p ) = ∂ p φ(p, p ) is the differential scattering phase. Eigenvalues of conserved charges can now be computed as In a "micro-canonical ensemble", there are several states with almost identical charges {Q n }, i.e. with the same charge density, as averaged densities are not sensitive to the redistribution of particles within infinitesimal regions of quasi-momenta. Equivalently, these are the states that are consistent with a given pair (ρ p , ρ h ). In an interval dp the number of different configurations N yielding the same densities ρ p and ρ h is given by We may then define an entropy from log N (ρ p , ρ h ) which can be evaluated in the thermodynamic limit (L → ∞) using Stirling's formula Following Yang and Yang's argument [87] and its generalization [70], the GGE distribution can be described using a partition function where, in the entropy S, ρ h has been expressed as a functional of ρ p using the thermodynamic Bethe equation (3.20), and w[ρ] = dp w(p)ρ p (p). We evaluate this partition function under the saddle point approximation, with saddle point condition This function is referred to as the "pseudo-energy". This is because it allows us to define an occupation function and recover the form it takes for the Fermi-Dirac statistics In a similar fashion, the free energy contribution of a Fermi particle of energy [91] F , (3.29) yields the free energy density F = −L −1 log Z saddle point of the system, in the form The Fermi-Dirac statistics is natural for the Lieb-Liniger model, even though it is originally a bosonic theory, because of the delta-function interaction, which forbids particles from being at the same point. For a given function w(p), or set of Lagrange multipliers {β n }, equation (3.26) completely determines the stationary state. Although the state is fluctuating, the fluctuations accounted for do not affect the charge densities (they are fixed to the saddle-point value), hence their averages still take the form (3.22), q n = dp h n (p)ρ p (p) , (3.31) where ρ p (p) is evaluated at the saddle point. The value of ρ p (p) can be established by expressing averages as derivatives of the free energy density This gives us the relation between phase-space densities and Lagrange multipliers β n , or spectral weight w(p), which is the main expression of the thermodynamics of the model. We obtain it in the form of integral equations, 2πρ p = n(P ) dr , (3.33) where the "dressing operation" is a linear operation defined by the integral equation Recall that in the Lieb-Liniger model P (p) = 1 is the constant unit function. Conversely, if ρ p is known, then ρ s is evaluated by (3.20) and therefore the pseudo-energy is known, from which the Lagrange multipliers can be computed by differentiating the saddle-point condition (3.26), . (3.35)

Equation of state
In order to construct a hydrodynamic theory we need a way to relate averages of currents to that of charges: we need an equation of state. In statistical mechanics, this information is fully encoded within the "free energy flux", first considered in [18]. This is a function G of the state whose derivatives give rise to average currents It is shown in [18] that such a generating function always exist. In free particle models, one obtains a form similar to (3.30), where the free energy F ( (p)) controls the contribution of each mode, but with a measure determined by the energy E(p) instead of the momentum. It is then natural to expect that this generalises to integrable models as This expression was first derived in the Lieb-Liniger model in [18] using a property of relativistic quantum field theory called crossing symmetry, and taking the non-relativistic limit. Simple calculations then give the rather intuitive expression Here v gr = E /P (= p in the Lieb-Liniger model) is the group velocity (and recall that P (p) = 1 in the Lieb-Liniger model). This effective velocity can be interpreted as the large-scale modification of the group velocity of a particle going through a gas as it accumulates scattering shifts when colliding with other particles. The expression (3.38) for the currents was given an alternative proof in quantum field theory [23], and a full Bethe ansatz derivation in quantum spin chains in [73]. Recently a general proof was obtained, valid in quantum and classical systems, using self-conserved currents in [80,88].

Hydrodynamic approximation
As in regular hydrodynamics, GHD relies heavily on the separation of scales: systems are divided in fluid cells considered thermodynamically large, so that we assume entropy to be locally maximised, but small compared to the laboratory scales over which variations of local averages occur. That is, we assume that the averages of local observables can be well approximated, at large times, by averages evaluated in local GGEs Ifq n andj n are differentiable, macroscopic conservation laws (3.41) can be re-written in their differential form where, this time, ∂ x and ∂ t represent large-scale derivatives encoding variations between fluid cells. Note that hydrodynamics is a derivative expansion; we will here only focus on the Euler scale, giving the lowest order approximation, and neglect for example diffusive, dispersive and higher order terms.

Fundamental equations
With the above, the fundamental equations of GHD are easy to derive. To that end we shall now characterize the states of locally maximised entropy by using results from section (3.2.2). Now that the Lagrange multipliers depend on (x, t), we must introduce the space-time dependence The space of pseudo-local charges being complete [29] we can finally obtain the fundamental equation of (Euler scale) generalised hydrodynamics It is important to note that we can diagonalise this continuity equation using the occupation function defined in (3.27) Thus, the TBA occupation function plays the role of a continuous Riemann invariant for the GHD equation.

Relation between KdV soliton gas and GHD
One can see some striking parallels between the spectral soliton gas and the GHD theories. In particular, the equation of state (3.39) and the continuity equation (3.45) in GHD have the same structure as the spectral kinetic equations (2.8), (2.9) in the soliton gas theory (up to the form of the phase shift kernel). Furthermore, the thermodynamic Bethe equation (3.20) and the equation (3.46) for the occupation function have their counterpart equations (2.10) and (2.12) respectively in the spectral soliton gas theory.
We shall now make explicit connections between the formalisms of soliton gases and GHD. We will re-contextualise the relevant quantities of one theory in terms of the other, methodically identifying the various notations, to make apparent what each approach can bring to the table. As direct by-product of this consideration, in the next section we shall then develop the thermodynamic characterisation of the KdV soliton gas in terms of entropy, free energy, temperature etc., and from this obtain and verify new conjectures for static covariances in the soliton gas.

"Classical Bethe Ansatz" and KdV
We start this discussion by a heuristic, TBA-inspired, approach to soliton gases. More explicitly, we follow the classical TBA approach, already discussed in the context of Toda gas in [31], as a way to heuristically recover the dispersion relation (2.10), analogous to the thermodynamic Bethe ansatz equation (3.20), albeit in the context of a soliton gas. The aim of this section is to build intuition and further reinforce the aforementioned parallels between soliton gas theory and GHD. It will also serve as a vector to naturally introduce notions such as the "asymptotic space", and a geometric interpretation of the classical TBA which will provide much insight into the hydrodynamics and thermodynamics of the soliton gas.
As we previously discussed, the KdV equation supports multi-soliton solutions which, at sufficiently large time, can be asymptotically represented as a superposition of well-separated single-soliton solutions (2.3), Here x − i and x + i refer to the phases of soliton i as t → −∞ and t → +∞, respectively. The "phase" of soliton i is the shift of the soliton's center with respect to the straight trajectory of velocity 4η 2 i passing by the origin. Then, x − i is the phase before the soliton interacted with any other soliton in the gas, and x + i is the phase after it accumulated shifts through interactions with every other soliton. This distinction is important as will be made clear below. These solitons will serve as the quasi-particles of our "gas", making for a rather literal interpretation of the "soliton gas" appellation.
In order to describe this gas, we consider two different sets of dynamical variables that can be extracted from the multi-soliton solution; their relations will allow us to derive the TBA.
One set is simply composed of the {η i }'s and the {x − i }'s. That is, we consider the positions at time t = 0 obtained by extending the linear trajectories of asymptotic solitons, and their spectral parameters. These asymptotic coordinates can be seen as dynamical variables under the KdV equation: under the KdV evolution we have But the {x − i }'s do not help in measuring averages of physical observables in the KdV gas, or in establishing its thermodynamics. The second way is to assume that, even at intermediate times t, solitons can be considered as quasi-particles indexed by i, of position x t i and spectral parameter η i . The KdV soliton physics is implemented, on these quasi-particles, by assuming that quasi-particles undergo 2-soliton scattering shifts (2.4), at all scattering events. Note that, by convention, in the GHD literature δ(η, µ) > 0 means that if the object η is to the left (right) of the object µ before the collision, then it will be to the right (left) after collision, and its trajectory has been shifted with respect to the non-interacting trajectory towards the left (right). In other words, for δ(η, µ) > 0 the shift is "backward" with respect to the direction of the collision. For the KdV soliton gas, δ(η, µ) < 0. With this prescription, the KdV evolution now induces a nontrivial dynamics U t x 0 i = x t i . This dynamics has an explicit algorithmic implementation as the "flea gas" [38], but its microscopic realisation is not important here. It satisfies the asymptotic condition , thus the quasi-particle describes correctly the trajectory of the soliton i at large negative times. Using this, we may relate the coordinates x 0 i of the quasi-particles at time 0, to the asymptotic coordinates x − i , taking into account all 2-soliton scattering events that occurred at negative times. As such, the quasi-particle gas, with coordinates x 0 i , is akin to the Bethe wave function (3.10).
The exact total scattering shift x + i − x − i of the i-th soliton, after full scattering has occurred, is given by the sum of individual phase shifts from the collisions with other solitons in the gas, as per the well known factorised scattering theory of KdV. Therefore, by the definition of the dynamics of quasi-particles, x t i ∼ x + i + 4η 2 i t (t → +∞), hence the quasi-particle also describes correctly the trajectory of the soliton i /at large positive times. The main assumption is that the quasi-particle description is in fact good enough also for intermediate times, even if actual solitons are not well separated. In other words, the effect of the i th soliton at time t is restricted to an area around x t i , even if the KdV field is of a rather complicated form and does not look at all like separate solitons; from the viewpoint of the large-scale physics and thermodynamics, this knowledge is sufficient. See Figure 1 for an illustration of the meaning of x ± i and x 0 i Now, consider x 0 i , the position of quasi-particle i at t = 0. In a gas of finite density lying on the interval [0, L], we take N solitons and require for the gas at time 0 to be confined within this interval, i.e. {∀i = 1 . . . N, x 0 i ∈ [0, L]}, for large N and L along with fixed N/L. We would like to "measure" the resulting spatial extent on the asymptotic, free-particle coordinates x − i given the constraint x 0 i ∈ [0, L]. The dynamics being 2-body reducible, to obtain x 0 i from x − i we need to add all spatial shifts incurred by quasi-particle i as it crosses other quasi-particles in the gas from time t → −∞ to reach x 0 i at time t = 0. Now, consider a choice of value of x − i such that quasi-particle i, at time 0, is all the way to the left, x 0 i = 0. We denote x − i = x left i . Since for t → −∞ quasi-particles are ordered from fastest to slowest, and since at t = 0 they all lie within [0, L], collisions have almost surely all come from the left and therefore involve all faster quasi-particles. Denoting by T left i the set of all quasi-particles that are faster than i, we find , which is such that the same quasi-particle i is all the way to the right at time 0, x 0 i = L, collisions have come from the right and involve all slower quasi-particles, from the set T right i : Writing both equations (4.6) and (4.7) allows us to compute the length of the asymptotic space i (the "free length" as perceived by quasi-particle i). Noting that the set T left i ∪ T right i contains all quasi-particles except for i (recall from Section 2.1 that all solitons have different velocities since η i = η j for i = j), we have This is illustrated in Figure 2. Let L N (η) be a smooth function interpolating L i 's for i = 1, . . . N and define an "asymptotic space density" K N as Taking the thermodynamic limit of (4.9) we obtain in which the DOS f (η) is the density of the spectral points η i per unit distance and per unit spectral interval (see Section 2.1). Here and elsewhere in this section the spectral integration is performed over the support Γ of the DOS f (µ).
(recall the definition (4.5) for δ(η, µ)). The first two equations follow by identifying Bethe roots with quasi-particle (soliton) spectral parameters; the last two follow from comparing (3.20) with (4.10). Recall that for the Lieb-Liniger model P (η) = 1, however we keep it arbitrary here as this is a more universal expression -we provide an in-depth discussion of P (η) in the context of the KdV soliton gas in Section 4.5. The classical TBA equation is also equivalent to one half of the dispersion relations, equation (2.10), from the soliton gas theory, upon the identification In particular, with the occupation function n = ρ p /ρ s , this leads to . (4.13) A few important comments are in order.
1. We assumed that spatial density of solitons does not depend on their spectral parameter. This is appropriate for the thermodynamics, and is what happens in the spectral theory of "periodic" soliton gases [85]. In the more general, quasi-periodic setting of [43], [44] the DOS is given by f (η) = φ(η)κ(η), where φ(η) is the spectral measure on Γ (the distribution of solitons over the spectrum) and κ(η) is the function responsible for the variable spatial density of solitons depending on their spectral parameter (in the periodic setting κ(η) = α -the inverse period).
2. The formula (4.10), as derived, is, strictly speaking, valid only for a rarefied gas. However, as mentioned, its validity in the context of a dense gas may be argued for by the idea that the effect of the i th soliton is indeed restricted to where the i th quasi-particle is located. Full justification is obtained by the thermodynamic finite-gap derivation of [43] (see also [45] and the discussion in Section 2.1).
3. The classical interpretation of K(η) in (4.11) as a density of asymptotic space was first proposed in [36] in the context of a geometric, semi-classical interpretation of the GHD equations. According to the definition (4.9) taken in the thermodynamic limit (L, N → ∞, L/N fixed), the quantity K(η) is the finite ratio of the perceived length L N (which grows to infinity) to the real length L (which also grows to infinity). As we assume the gas to be homogeneous, this translates into the ratio of a length element dx − (η) in the asymptotic space of a soliton with spectral parameter η, to the length element dx in real space (in the thermodynamic limit). Thus, there is a change of coordinates dx − (η) = K(η)dx as viewed by quasi-particle η, and we may identify the metric as g(η) = K(η) 2 . Because the KdV equation translates into the free evolution in asymptotic coordinates, the fluid equation in terms of asymptotic coordinates is simple, and the change of coordinates to real space gives rise to the GHD equation. This is reviewed in Subsection 4.3.2. Interestingly, it relates the classical concept of a metric in asymptotic space, with the quantum Bethe ansatz concept of the total density ρ s (η). In the soliton gas theory, it gives a geometric interpretation of the spectral scaling function σ(η).

Conserved quantities in KdV
It is important to identify the charges of GHD with the conserved quantities of KdV soliton gases (c.f. section 2.3). It is often convenient to work with their local densities. Here, the notations in the KdV and GHD contexts collide. In the GHD notation, for instance for the quantum Lieb-Liniger model (see Section 3), the densities are denoted q i (x, t), the total charges Q i , and the index usually starts at i = 0 with an "ultra local" charge (if the model admit such a charge), whose Hamiltonian flow is trivial. With the GHD convention, the first few conserved densities are defined, according to (2.16), as q 0 = u, q 1 = u 2 , q 2 = −u 3 + 1 2 u 2 x , . . . . (4.14) The latter two, when integrated, reproduce the momentum and energy (Hamiltonian in one canonical formulation), so we have While Q 1 and Q 2 respectively generate space and time translations using the canonical Poisson bracket, Q 0 generates a trivial Hamiltonian flow: it is, in the GHD language, an "ultra-local charge", much likeQ 0 for the Lieb-Liniger model We may express those charges in terms of the spectral parameter of the elementary solution that is the soliton u η (·, ·) = u s (·, ·; η), given by equation ( Given any uniform density of solitons ρ p (η) per unit distance and per unit spectral parameter -recall that this is equivalent to f (η) within the soliton gas notations -then the average of conserved densities can be expressed as q n = dη ρ p (η)h n (η) , (4.18) just as it was in the quantum case (3.31); the h-functions are model dependent but the equations take the same form. This follows from the invariance of Q n under the KdV time evolution, and the asymptotically separated form (4.1) of the N -soliton solution. This is perfectly equivalent to equation (2.18).
As the densities are conserved, they have associated local currents, which can be evaluated from the KdV equation since each q n only depends on u and its derivatives, hence This is just equation (2.15), the classical version of (3.6). In particular, the KdV equation itself directly leads to the relation (4.20)

GHD equation for the gas of solitons from classical TBA
Following the logic of the quantum TBA, we would now construct the free energy. However, this has no direct equivalent in the theory of soliton gas as developed until now. But within the classical TBA, it is in fact possible to directly justify the GHD equations, without the need for the full thermodynamics as was done in the quantum context. Hence, we keep the development of the KdV thermodynamics for Section 5, and here complete the parallels between TBA and the KdV soliton gas theory of Section 2. We would like to derive a hydrodynamic equation for the gas of quasi-particles of the classical TBA developed in Section 4.1. We present two separate arguments.

Collision rate ansatz
First, the soliton collision rate ansatz (2.7) is immediately applicable to the quasi-particles of the classical TBA. Within the gas every quasi-particle has "bare" group velocity (4.4), v gr (η) = 4η 2 , while the two-body scattering is characterised by the scattering shift (2.4), or equivalently (up to a sign convention) δ(η, µ) (4.5). The shift is then used to calculate the effective velocity of a quasi-particle with given spectral parameter as a modification to its group velocity to account for collisions with other quasi-particles within the gas,  The collision rate ansatz is then used to express the average currents within the gas, which, given equation (4.18), has a very intuitive interpretation: the average current is obtained by multiplying the amount h n (η) of charge carried by a soliton of parameter η by the density of solitons of the same parameter and by the velocity at which they move within the gas. Note that from (4.20) we have j 0 = 3 q 1 , (4.24) such that KdV admits a self-conserved current [80,88]. Finally, following from Section 3.3.2, in a long-wavelength inhomogeneous gas, the GHD Euler equation agrees with the KdV soliton gas kinetic equation (2.9). The relation (4.13) between the solitongas quantity σ(η) and the TBA quantity n(η) makes it clear that the "continuum" of Riemann invariants found in soliton gas, equation (2.12), and in GHD, equation (3.46), are essentially just inverse of each other.

Change of metric
The second way does not involve the collision rate ansatz, but rather the metric interpretation of K that arose from the classical TBA. For a weakly inhomogeneous gas confined to a large volume [0, L], we may repeat the classical TBA arguments of Section 4.1, and obtain a relation between the asymptotic coordinates of quasi-particles x − i of spectral parameter η i , and the real coordinates x 0 i . Recall from (4.3) that under KdV evolution, and x 0 i evolves to x t i . Denoting by x − (η; x, t) the asymptotic coordinate for spectral parameter η corresponding to real coordinates x at time t, we have, according to the results of Section 4.1, Integrating, and far enough on the left, x → −∞, where the gas is absent, K(η; y, t) = 1, which gives a full characterisation of x − (η; , x, t) up to an overall constant shift. As the evolution of the asymptotic coordinates is trivial, the density ρ − p (η; x − , t) of quasiparticles in the asymptotic space satisfies the Liouville equation, Clearly, taking into account the transformation property of a density, we have under the relation (4.27). Using the results established above, this gives . (4.31) The combination of (4.27), (4.29) and (4.31) leads to equation (2.12), equivalently equation (3.46); the derivation is provided in [36]. Therefore, the soliton gas kinetic equation, equivalently the GHD equation, is nothing else but the Liouville equation for asymptotic coordinates, under the transformation of metric that accounts for the relation between the asymptotic and real spaces induced by the two-body scattering shifts.

Nonlinear dispersion relations
Let us recall the relations obtained until now: from (4.11) and (4.13) and from (4.22), The left-hand sides are quantities from soliton gas theory: the spectral parameter η, the DOS f (η), the spectral scaling function σ(η), the two-body soliton scattering shift δ(η, µ) = sgn(α − η)δ(η, µ) (up a sign convention), and the tracer velocity s(η). The right-hand sides are those from TBA and GHD: the Bethe root, or quasi-momentum, p, the quasi-particle phase-space density ρ p (η), the occupation function n(η), the differential scattering phase ϕ(η, µ), and the effective velocity v eff (η). The right-hand side also involves a momentum function P (η), which has no equivalent in the soliton gas theory and which we discuss below. The last missing element of the soliton gas theory that needs comparison to GHD are the nonlinear dispersion relations, equations (2.10) and (2.11). We have seen that (2.10) is in fact the TBA equation (4.10) for the density of asymptotic space in the classical TBA, or (3.20) for the total density of quasi-particles and holes in the quantum TBA. Using the relation (3.27), ρ p = nρ s , the total density takes the form, in terms of soliton gas objects,   It is, using v gr (η) = 4η 2 , the dressing equation This justifies the introduction of the energy function E(η), whose derivative is In this sense, (4.38) and (4.40) really are "dressed" dispersion relations. They lead to the well-known GHD relation for the effective velocity which is a dressed version of the usual definition of the group velocity in terms of the dispersion relation.
In the Lieb-Liniger model, the energy function as defined here agrees with the energy oneparticle eigenvalue E(p) = p 2 /2 in (3.16). In the soliton gas theory, it is an object that depends on our choice of the momentum function P (η), which we now discuss.
Note that the dressing operation has no explicit equivalent in soliton gas theory, although the y α functions introduced in [57, Eq 49] can be intrepreted as the dressed quantities h dr α .

Momentum function
As seen above, in many relations between objects from soliton gas theory and from GHD, the momentum function P (η) (4.43) is involved. This function does not enter any of the physical equations of soliton gas theory: the kinetic equation (2.9), the dispersion relations (2.10), (2.11) or the effective velocity relation (2.8). Therefore, within the soliton gas theory, it is superfluous: the choice of momentum function does not alter any of the physical quantities evaluated in the soliton gas. Nevertheless, it is a part of its relation to GHD. This is because GHD was developed for quantum systems. In quantum systems, the momentum function P (η) has an unambiguous physical meaning. Quantum scattering is described by a phase S = e iφ , and in order to obtain a scattering shift in space, we need to know "what space is". This information is given by the one-particle eigenvalue P (η) of the physical momentum operatorP . In continuous quantum systems with translation-invariant dynamics such as the Lieb-Liniger model -the systems of interest here -the momentum operatorP always exist and defines the notion of space, as operatorsÔ are displaced viaÔ(x) = e −iP x Oe iP x (with = 1). Here one-particle states are parametrised by some "rapidity" η; the precise choice of parametrisation is not important, but the momentum eigenvalue is unambiguous. From it, the semi-classical scattering shift is The choice of parametrisation η is determined by the momentum function P (η). Once a parametrisation is chosen, the conventional definition of the differential scattering shift is in terms of this parametrisation, ϕ(η, µ) = −i ∂ log S(η, µ) ∂η . (4.45) This gives the relation between differential scattering phase and scattering shift, The parametrisation of quasi-particles in GHD in terms of a "rapidity" η is in fact an arbitrary choice. Indeed, all GHD equations are invariant under re-parametrisation, when written in universal ways in terms of the fundamental quantities; this justifies the writing in "universal" forms used above. These quantities themselves transform, under re-parametrisation, either as vector fields (like derivatives ∂/∂η, hence getting an additional P factor) or as scalar fields (without getting any additional factors). Under a change of parametrisation, all equations written in their universal form stay invariant if the following quantities transform as follows: ρ p (η), ρ s (η), P (η), E (η) as vector fields (4.47) and v gr (η), v eff (η), n(η), K(η), δ(η, µ) , h n (η) as scalar fields (in all variables); (4.48) recall that h n (η) are the one-particle eigenvalues used in evaluating average conserved densities. As a consequence of (4.46), ϕ(η, µ) is a vector field in the first variable, and scalar field in the second. Reparametrisation invariance was discussed in the context of GHD in [30]. The relation to soliton gas, however, makes clear that the universal formulation of GHD is in fact redundant, something which had not been pointed out before. Indeed, once a parametrisation is chosen, all physical quantities (other than the momentum itself) in fact do not involve the momentum function, and can be recast in terms of momentum-independent, "fundamental" quantities. Note how the following combinations only depend on the fundamental quantities of soliton gases ρ p (η), δ(η, µ), v gr (η), and not on P (η): , n(η)P (η), σ(η), 2πρ s P (η) , K(η), v eff (η) = (E ) dr (η) (P ) dr (η) are independent of P (η). (4.49) This means that, when re-writing the soliton gas theory in terms of GHD, we may in fact choose the momentum function that we wish; there is no "more meaningful" choice coming from physical scattering phases, as there are no scattering phases, just scattering shifts. As this is a structural observation, it is true for the GHD of any quantum system as well: it can first be written in terms of fundamental soliton-gas quantities, and then re-written in a universal form using any chosen momentum function. For quantum systems, however, choosing the incorrect momentum, although still giving correct hydrodynamic equations, would give the incorrect thermodynamics; this is not so in classical systems, as discussed in the next section. The freedom in choosing the momentum function in the soliton gas theory can be used fruitfully, in order to simplify the relation with GHD. A "physical" choice is the momentum carried by the soliton, from (4.17) With this choice, the energy function defined as E = v gr P is then the actual physical energy of the soliton, E phys (η) = h 2 (η) = 32 5 η 5 (4.51) and the differential scattering phase is In fact, this is not the most convenient. Instead, as in the classical TBA, we can see solitons as classical quasi-particles: it is natural to simply choose the momentum corresponding to unitmass particles of velocities v gr (η) (4.53) In this way, the energy function takes its natural form These, in the soliton gas theory, do not correspond to any conserved charge in the Kruskal series, expressed in terms of extensive functionals of the KdV field. Nevertheless, they are natural in the quasi-particle picture of soliton gases. With this choice, the differential scattering phase becomes This now satisfies the symmetry relation ϕ(η, µ) = ϕ(µ, η) , (4.56) which simplifies many of the equations. This choice also streamlines the relation between two fundamental quantities of both theories: σ and n σ(η) ≡ π 4n(η) , (4.57) making one simply the inverse of the other, up to constant factors. In the TBA literature, one typically chooses a parametrisation that guarantees symmetry of the differential scattering phase. As discussed in [30], the dressing operation in general is different for quantities h n (η) that are scalar fields, and those that are vector fields such as P (η): the differential scattering phase appears with opposite order of its arguments. This guarantees that the dressing operation preserves the transformation property. But with the symmetric choice, the dressing has the same form in both case.
Finally, we note that average densities and currents may be re-written using the dressing operation in various ways: (4.59) The last two equalities in equations (4.58) and (4.59) are "symmetry relations" consequences of the dressing definition. One can deduce from this symmetry that, in average, the current of solitonic charge j , with h-function h(η) = P (η), is equal to the density of charge q , with h-functionh(η) = E (η). Equality of the current of a certain conserved charge with the density of another conserved charge is in fact "boost-invariance relation". For instance, in Galilean systems, the particle current equals the momentum density, and in relativistic systems, the energy current equals the momentum density. Many of the thermodynamic consequences of boost invariance can be deduced from such equality. Given our current conventions P (η) = 2h 0 (η) and E (η) = 6h 1 (η), we conclude from this symmetry that j 0 = 3 q 1 (4.60)

Quantity
Soliton gas theory Generalized Hydrodynamics Conservation laws recovering (4.24), a direct consequence of the KdV equation itself. We note that relation (4.60) can be used, in combination with the exact expression of the free energy (5.10) derived in the next section, in order to obtain the exact expression of the currents (4.59) using the methods of [80,88], instead of directly using a collision rate ansatz. However, here this does not increase the rigour of the derivation, as the classical TBA leading to the free energy of the soliton gases is not on fully rigorous grounds.

Overview of the soliton gas -GHD relation
The parallels drawn in this section are summarised in Table 1 under the conventional choice of momentum P (η) = v gr (η). With this formulation, we can now argue as in [31] in order to derive the thermodynamics of the soliton gas.

Thermodynamics of the KdV soliton gas: conjectures from GHD
In this section we shall use the previously established correspondence in order to construct the thermodynamics of soliton gases by directly applying results from GHD and TBA in this new context. As such, we will exclusively use notations associated with soliton gas theory, otherwise referring the reader to the dictionary Table 1.

Generalised Gibbs Ensemble for KdV soliton gases
The thermodynamics for a soliton gas can be constructed by writing the generalised Gibbs measure for an ensemble of N -soliton solutions (4.1), using the asymptotic, action-angle coordinates (η n , x − n ) discussed in Section 4.1. As these evolve as free particles, the a priori invariant measure for a single particle, with conventional normalisation coming from the choice = 1, is (2π) −1 dP (η n )dx − n , where the momentum function is the quasi-particle velocity, P (η) = v gr (η) = 4η 2 . Therefore, the invariant measures -the GGEs -for a gas on a finite length L are of the form w(η i ) χ |u N (x, 0)| < ε |x|−L/2 for |x| > L/2 , (5.1) for functions w(η). Here ε y → 0 + quickly enough as y → ∞, say as ε y = e −y/a for some a > 0. This measure is to be taken on η n > 0, x − n ∈ R. The measure factorises into the scattering data except for a single non-trivial factor, which confines the gas (at time t = 0) to lie on the interval [−L/2, L/2], up to exponentially small corrections, and plays the role of an entropy term.
The measure may be seen as coming from the conserved charges in the GGE way, as was done in the quantum context in Section 3.2. It involves some conserved charge W such that As in the quantum case, there is no necessity for w(η) to be expressible as a Taylor series for the measure to be a valid GGE. A valid GGE is a space-time translation invariant measure which is clustering at large distances (which is ergodic), and a large space of functions w(η), or equivalently conserved charges W , will work, going beyond those which possess a Taylor expansion.
We are interested in taking the limit L → ∞; note that for "good" w(η), the average number of solitons will increase in proportion to L, to give a finite density, in a way that is determined by the precise form of w(η).

Free energy and entropy of KdV soliton gases
The GGEs now characterised, the next step in discussing the thermodynamics of our soliton gas would naturally be to construct its free energy and entropy.
For this purpose, we need to adapt (3.26) and (3.30), derived in the quantum context, to the classical TBA. As stressed in [32], even if the TBA was originally developed to tackle quantum integrable systems, its general structure is not limited to those: classical limits have been investigated in [26] or [83], and [82] explored TBA approaches to soliton-bearing systems. In all cases, the TBA structure remains valid, with F ( ) being replaced by the free energy function representing the correct statistics of the asymptotic particles.
The universality of the TBA has, in fact, a geometric underpinning which follows from the geometric interpretion [36] of the classical TBA developed in Section 4.1. Recall that in the geometric viewpoint, the interaction can be implemented as a density-dependent change of metric on a system of free particles. This is because the scattering shifts δ(η, µ) effectively change the perceived space. In this interpretation, each "mode" µ induces an extra amount of physical space given by δ(η, µ) as perceived by the mode η. This change of metric is exactly what we need to take into account when recasting the entropy factor in (5.1) as a constraint on the asymptotoic coordinates x − n . The full calculation from this geometric idea was explained in [31]. It is easy however, from this perspective, to re-interpret the formulae obtained in the quantum context, in order to explain the classical result. In the quantum context, equation (3.30) is the free energy density for a system of free particles where each quasi-momentum p is independently distributed in accordance with the Boltzmann weight determined by the "energy" (p). But the Yang-Yang equation (3.26) says that the interaction "dresses" this energy: it modifies the probability for the level to be occupied by taking into account the free energy within the additional space introduced by mode p.
With this interpretation, for the soliton gas, the free energy density (per unit length) is where, because of the measure (5.1), we must use the Maxwell-Boltzamann statistics of classical particles, and (η) free energy in additional space introduced by quasi-particle η . (5.7) Note that dµ ϕ(µ, η) = dP (µ)δ(µ, η) and, as such, the right-hand side actually involves the scattering shift. Hence the extra length perceived by quasi-particle (or soliton) µ, as it should.
In thermodynamics the occupation function is fundamentally defined by (3.29), and thus as per the Maxwell-Boltzmann statistics. Recalling the equivalence (4.57) we can straightforwardly obtain our conjectured expressions for the thermodynamics of the soliton gas. We have (η) = log 4σ(η) π , F ( (η)) = − π 4σ(η) , (5.9) which is directly related to the free energy density of the soliton gas The classical Yang-Yang equation (5.7) then allows for re-contextualisation of the notions of inverse temperature and chemical potential in terms of the soliton gas characteristic quantity σ(η), Recalling the definition (5.4) of w(η), and that h k (η) = C k η 2k+1 for KdV -given the normalisation (2.18) -we may write 12) and, in particular, the inverse temperature of a soliton gas (i.e. the Lagrange multiplier associated with its energy) is given by Finally, the free energy density also provides a way to compute the entropy density S F = W − S , (5.14) hence we infer by successively making use of Yang-Yang equation (5.7) and then of the first dispersion relation (2.10), after commuting integrations on η and µ. The soliton gas spectral entropy density S(η) = f (η)[1 + log 4σ(η) π ] is then consistent with both the classical entropy, that typically takes the from S cl (η) = ρ p (η)[1 − log n(η)], and the notion of condensate (examined in more details in Appendix A), since it reaches its minimal value for σ = 0.
Above we have used the Maxwell-Boltzmann statistics: we are explicitly working under the assumption that solitons behave as classical particles. The validity of this assumption is examined in Appendix A. A further justification is that the Maxwell-Boltzmann statistics (5.6) is the only choice which has the property that, when (5.7) is written in terms of the fundamental soliton gas quantity δ(α, η), a change of momentum function dP (η) → r(η)dP (η) is equivalent to a change of bare weight w(η) → w(η) − log r(η), in agreement with the GGE measure (5.1). This property is desirable, as the choice of the momentum function is irrelevant to the physics of the soliton gas, besides its effect on the chosen a priori measure.

Correlation functions
With the above results, we can now borrow the GHD results from [35,30]. If we take (5.1) as describing the correct soliton gas measure, assuming solitons behave as classical particles, then we are able to conjecture fluctuations (thermodynamics) and correlations (hydrodynamic linear response) in the soliton gas. These directly follow from the expression of the free energy (5.10) as per standard statistical mechanics. Recalling from definition (4.58) and using the KdV notations P n and Q n for the charges and currents space-integrated connected correlation functions take the form where θ(η) characterises the statistics, details of the computations being found in [30]. In particular, for the Maxwell-Boltzmann statistics, Correlation functions, by way of the statistical factor θ(η), thus provide an avenue for investigation into the statistics of solitons, which we perform in Appendix A. Similarly, we can construct the field-current correlator by differentiating the free energy flux G such that and hence The spatial integrals are taken over any region that is large enough surrounding the point 0.
We expect correlation functions to decay exponentially, so the integration should go over "a few correlation lengths". Note that general principles of statistical mechanics imply [18,28] B ab = B ba . (5.21) The above correlations represent the thermodynamics, but more interesting are maybe the quantities that involve the dynamics. The Drude weight, for instance, is an important quantity for characterising ballistic transport. It can be expressed in terms of a time integral of correlation functions using the Kubo formula and takes the form For this result to hold, the integration time t typically needs to be larger than the relaxation time to pure ballistic effects, although it may be difficult to assess how large that is. Going further, we may evaluate the actual two-point functions of conserved densities and currents. The usual results from hydrodynamics apply at the Euler scale 23) where η * (ξ) = {η : s(η) = ξ}. By convention, evaluation of two-point functions (5.23) often involves an averaging (usually either over space, time or space-time along the ray) to eliminate fluctuactions. We therefore predict a decaying behaviour in 1/t.
To conclude with an example, we shall consider the static covariance C 00 and the two-point auto-correlation function for the KdV wave field P 0 (η) = u(η), with h 0 (η) = 4η. Explicitly, we start by evaluating the dressed identity in terms of the soliton gas notations obtained from definition (3.34). Comparison with the dispersion relation (2.10) allows for direct identification id dr (η) = σ(η)f (η). In the more general case one typically needs to evaluate dressed quantities numerically (c.f. Appendix B). As such the simplest correlations readily take the form (5.25)

Conclusions and Outlook
We have reviewed the main aspects of two theories for long-wavelength, hydrodynamic behaviours that were recently developed: that of soliton gases in integrable dispersive hydrodynamic equations (or "integrable turbulence"), and that of generalised hydrodynamics (GHD) for many-body integrable systems. For the former we took the example of the KdV equation, and for the latter, the quantum Lieb-Liniger model, both paradigmatic models in their respective domains.
As noticed in the literature, these theories share a number of similarities, both conceptual and technical. In fact, it is natural to expect that the theory of GHD be able to describe soliton gases, and that the theory of soliton gases be at the basis of GHD. We have shown that this is indeed the case, establishing the precise dictionary between the quantities and concepts studied in soliton gases and in GHD. The crucial observations were the parallels between the soliton gas kinetic equation and the Euler-scale GHD equation, and between the soliton gas dispersion relations and the (classical) thermodynamic Bethe ansatz (TBA) equations. This dictionary has allowed us to clarify aspects of each theories. For instance, we emphasised the geometric interpretation of the soliton gas kinetic equation; and we pointed out that the momentum function, previously seen as a fundamental quantity in GHD, is in fact superfluous and ambiguously defined in the soliton gas interpretation of GHD.
Most importantly, using this dictionary and known results in quantum and classical integrable many-body systems, we went beyond the kinetic equation for the KdV soliton gas, and proposed its exact thermodynamics: free energy, entropy, generalised temperatures. For this purpose, we argued that the classical (and generalised) Yang-Yang equations of TBA correctly describes the soliton gas thermodynamics, if the statistics of the solitons is taken to be that of classical particles, with Maxwell-Boltzmann distribution. The resulting conjectures for susceptibilities of various conserved quantities in the KdV soliton gas were verified numerically, using state-of-the-art numerical methods for generating soliton gases. We compared against other common statistics found in many-body systems (Bose-Eintein, Fermi-Dirac), and confirmed that the Maxwell-Boltzmann statistics gives the most precise predictions.
One of the pertinent questions arising in relation with soliton gas theory is the issue of "completeness" of the soliton gas/GHD description of the KdV hydrodynamics. It is known, by a classification based on the standard inverse scattering method, that the KdV equation admits both solitonic modes (discrete part of the spectrum) and radiative modes (continuous part of the spectrum). Shouldn't a typical KdV gas contain, then, also an ensemble of radiative modes? It is of course possible, naively, to include radiative modes in the hydrodynamic and thermodynamic description of the gas, as the scattering shifts are known.
But we note that a soliton gas lying on the full line does not exhibit decay as |x| → ∞, and so its spectrum cannot be rigorously classified in the framework of the conventional discrete/continuous dichotomy. The intuitive notion of the purely discrete spectrum of soliton gas comes from the physically natural assumption that, on any sufficiently large interval x ∈ [x 0 , x 0 + L] ∈ R, a soliton gas can be approximated by an appropriate N -soliton solution with N 1. This assumption is at the heart of the phenomenological construction of soliton gas outlined in Section 2.1. The GHD interpretation of soliton gas in Section 3 is consistent with this "N -soliton paradigm". The dispersive hydrodynamic spectral theory of soliton gas developed in [43,45] and outlined in Section 2.2 is based on a more general definition of soliton gas as the thermodynamic limit of finite gap, non-decaying KdV solutions exhibiting solitons and radiative modes as special cases. Remarkably, this "wave" approach to soliton gas bridging solitons and radiation leads to the same hydrodynamic description as the "particle" GHD construction, which suggests that the GHD in fact captures the "full" KdV hydrodynamics.
But perhaps the strongest argument in favour of the "completeness" of the soliton gas description of random KdV fields (the integrable turbulence [90]) is the evidence, at the level of spectral data, that in the limit N → ∞ solitons can provide an approximation to the general solution of the integrable model in question, including the continuous spectrum modes (the latter are approximated by a large number of solitons with small amplitudes), see [47], §8. Although the discussion in [47] refers to the NLS equation, it is clear that the same arguments are directly applicable to the KdV equation as well. We also point out the relevant recent paper [54] where a soliton gas solution with non-zero reflection coefficient was constructed as the limit as N → ∞ of the N -soliton solutions of the KdV equation.
Of course, a full understanding of the completeness of the soliton gas set of "action-angle variables" is a significant open problem, not only in terms of a rigorous proof, but also in terms of its very formulation. One question is as to the existence of consistent subsets of hydrodynamic modes. It is well known that in GHD, one may restrict the initial state to a subset of hydrodynamic modes, which is invariant under time evolution. For instance, in the hard rod gas, one may form an ensemble composed of a finite number of velocities (yet an infinite number of rods). Such restricted ensembles are invariant under evolution (in general, this will be true at least at the Euler scale), and thus have a consistent hydrodynamics. Likewise, a soliton gas certainly gives a consistent hydrodynamics. Another question is as to the definition of the "right" maximal-entropy ensemble for the KdV field, for instance the thermal states. What is a Boltzmann weight for the KdV field? What Hamiltonian should one use, and what a priori measure on the field itself? We leave this for future research.
GHD makes a number of far-reaching predictions for the large-scale behaviours of many-body systems, once the statistics of the underlying quasi-particles is known. With Maxwell-Boltzmann statistics for the KdV soliton gas, we made explicit predictions for the matrix of Drude weights, describing ballistic transport of KdV conserved quantity, and for the Euler-scale, dynamical correlation functions.
GHD can be used as the basis of many more predictions, which can immediately be written for soliton gases. As an example, the presence of external potentials coupling to any conserved density in GHD was described in [37], and suggests, for instance, that the kinetic soliton gas equation for evolution under u t + (α(x)(3u 2 + u xx )) x = 0 (6.1) with long-wavelength α(x), takes, instead of (2.9), the form for an effective acceleration a (whose expression can be inferred from [37] using the dictionary described in the present work). More general forms of large-scale space-time dependence in the evolution equation can be introduced and dealt with using GHD [4]. Another example is the exact expression for the diffusion operator, a second-derivative additional term in the kinetic equation. This is known in GHD [27,28], and is conjectured to be applicable to large classes of integrable models, only requiring the knowledge of the scattering shift and the statistics. Thus, the results of the present work can directly be used to make a prediction for diffusion in the KdV soliton gas. A final example is the large-deviation function for fluctuations of soliton transport, which can be obtained from the theory developed in [34,71]. Another aspect is the opportunities that the KdV soliton gases give for a more precise control on its GHD, thanks to the extensive mathematical structure and understanding available. Some of the techniques used in KdV soliton gases may also be transported to other models with a GHD description. For instance, can we adapt Whitham modulation theory to quantum many-body systems?
Finally, it is likely that the arguments and results of the present paper can be extended to soliton gases in other integrable dispersive hydrodynamic systems, particularly the focusing NLS equation, a canonical model for the description of modulationally unstable physical systems. The spectral theory of soliton and breather gases for the focusing NLS equation has been developed in [44] but the counterpart GHD description can provide further insights into the fundamental phenomena of integrable turbulence and rogue wave formation.
A numerical realisation of the genus 0 soliton condensate as a dense ensemble of KdV solitons constructed according to the spectral DOS (A.3) is presented Figure 3. The method for the numerical realisation of KdV soliton gases used here, based on N -soliton solutions with N large, has been developed in [20]; we present a brief description of it in Appendix C. We should stress that, to achieve the required spectral DOS (A.3) in the numerical soliton gas, one has to "place" all solitons in a single point, where the location of a soliton in the gas is defined in terms of the phase of the corresponding norming constant in the N -soliton solution. In other words, this means that, in the context of the TBA approach introduced in section 4.1, for any soliton i of the gas This implies that, given a condensate supported on an interval of size L at t = 0, for any two solitons i and j of the condensate ∆ i = ∆ j = L.
As one can expect, the soliton condensate is not the most interesting case to investigate the statistics of solitons since σ → 0 implies that all the dressed quantities, and by extension the correlation functions, vanish. This agrees with the above statement that u(x, t) = 1 almost everywhere for such a gas, and can also be readily interpreted in GHD terms in the simplest case by looking at equations (2.10)-(5.24) from which we make the identification for any integer n, and whatever the statistics.

A.2 Diluted condensate
A more interesting case for the matter at hand would be the "diluted" condensate [20] constructed from the spectral DOS f (η) = α η where 0 < α < 1. Notably, even for dilution parameters α ≈ 1, the resulting soliton gas exhibits spectacularly different qualitative features compared to the "genuine" condensate, as exemplified in Figure 4. From both the general expression of the dispersion relation (2.10) and the particular case of the condensate (A.1) we infer Note that σ > 0 implies α < 1 which, since α can be directly related to the mean wave field in the gas (2.19) agrees with the notion of the condensate as the densest possible gas, and that the spectral DOS (A.6) actually corresponds to a diluted version of it. In terms of the TBA approach of section 4.1, this dilution corresponds to a widening of the asymptotic space in which the impact parameters (i.e. the phases of the norming constants) are distributed. From here it is possible to compute the simplest correlation function, recalling the identifi-cation (A.4) in the classical case. However, it is typically difficult to obtain explicit analytical expressions, and one usually has to resort to numerical integration (details of which are found in Appendix B). Here we start by solving the second integral equation (2.11) of the dispersion relations, with σ(η) given by (A.7); this is equivalent to computing 4η 3 dr ≡ 3h dr 1 /4. Recalling that in (2.11) v(η) ≡ s(η)f (η), the previous computation directly yields the effective velocity and we recover the canonical GHD expression (4.42) for the choice of momentum function (4.53). From there, using expressions (5.17) and (5.22), we can compute, among other quantities, C 01 , C 11 , B 00 and D 00 for different statistical factors, which we then compare, in both Figure 5 and Table 2, to correlations obtained from a simulated gas. As we can see, even though the data is admittedly noisy, simulations seem to indicate solitons behave as classical particles. Figure 5: Relative difference δ = (X GHD − X simulations )/X GHD , between the GHD predictions and simulations regarding the correlation functions X = C 00 , C 01 , C 11 , B 00 and D 00 . The horizontal axis represents the size of the box over which correlations are computed; for large boxes the gas cannot be considered homogeneous anymore due to finite size effects. Full curves correspond to Maxwell-Boltzmann statistics, dashed to Bose-Einstein and dotted to Fermi-Dirac. Results are obtained from the average over 6000 realisations of a diluted condensate of 150 solitons with parameter α = .91.
Computing correlations for charges of order n ≥ 2 would require solving the integral equation defined by the dressing operation for η 2n+1 dr , which at this point should be rather straightforward.

A.3 Uniform spectral distribution
Another case we may consider is that of a uniform spectral distribution f (η) = α, for which some analytical results are still easy to obtain. Recalling again that id dr = σf (η), we have (A.11) The cut-off η 0 is here to ensure that σ (as well as id dr ) stays non-negative. By noting that σ(η) reaches its minimal value for η = η min η min = 1 + η 2 0 e α 1 + e α , (A. 12) and imposing that σ(η min ) = 0, we obtain an implicit relation between η 0 and α can be interpreted as a way to ensure the gas remains less dense than the condensate. Table 3 compiles the results shown in Figure 6: GHD predictions for different correlation functions, assuming either Maxwell-Boltzmann, Fermi-Dirac or Bose-Einstein statistics, are compared to results of the simulations. As with the diluted condensate, computing correlation functions require numerical integration and, as before, simulations favour Maxwell-Boltzmann statistics.

A.4 Linear spectral distribution
The last example we shall consider soliton gas realised as result of a long-time KdV evolution of an infinite sequence of broad parabolic pulses randomly distributed on R [57]. The spectral DOS of such a gas grows linearly with the spectral parameter, f (η) = 2αη. As before, we can compute the dressed identity id dr (η) = 2αησ(η) = η + 2α for which a cut-off needs to be defined as well so that σ > 0 and the average density  Figure 7 and Table 4).
While the previous three examples do not constitute a proof that soliton behave as classical particles, we feel they serve as a good indicator. Ultimately, it does not seem outlandish to assume Maxwell-Boltzmann statistics is the correct best fitted to soliton gases.  Table 3: GHD predictions regarding a soliton gas with uniform DOS for various statistics displayed against the simulated correlation functions (averaged over L B ∈ [100; 400]) with standard deviation. Figure 7: Relative difference δ = (X GHD − X simulations )/X GHD , between the GHD predictions and simulations regarding the correlation functions X = C 00 , C 01 , C 11 , B 00 and D 00 . The horizontal axis represents the size of box over which correlations are computed; for large boxes the gas cannot be considered homogeneous anymore due to finite size effects. Full curves correspond to Maxwell-Boltzmann statistics and dashed to Bose-Einstein; relative difference for Fermi-Dirac statistics are too large for their inclusion to be relevant. Results are obtained from the average over 3000 realisations of a simulated gas of 150 soltions with linear spectral distribution α = .26 and cut-off η 0 = .4.

B Computing dressed quantities numerically
We start by recalling the definition of the dressing operation (3.34) a dr (η) = a(η) +  Table 4: GHD predictions regarding a soliton gas with linear DOS for various statistics displayed against the simulated correlation functions (averaged over L B ∈ [100; 300]) with standard deviation.
easily computed analytically. Once written in matrix form it only takes a few seconds to solve system (B.4) on a mid-range laptop using Matlab for n = 4000. Depending on how well-behaved σ is, it may be relevant to go for quadratic interpolation with graded mesh, but it usually seems unnecessary.

C Numerical generation of a soliton gas
In this Appendix we briefly outline the scheme used to numerically generate the soliton gases presented in Appendix A. This procedure is based on [58] and was originally used to generate the numerical synthesis of soliton and breather gases for the focusing NLS equation in [53,75]. The modification of this scheme for the generation of KdV soliton gas adopted here has been developed in [20]; we refer the interested reader to the latter reference for an in depth presentation. Our scheme generates exact N -soliton solutions of KdV by algebraic recursive procedures. The Darboux transformation relates the N -soliton Jost solution J N (η) of the KdV auxiliary IST system to J N −1 (η), so that and D i (η) is the i-th Darboux operator which depends only on the spectral parameter η i , the norming constant b i ≡ (−1) i exp 2η i x − i and on J i−1 (η i ) [58]. As such, the scheme can be broken down into three main steps 1. We start by sampling a set of N spectral parameters {η i } and N impact parameters {x − i } from two given distributions; in the simulations of Appendix A spectral parameters were sampled from the according DOS, and impact parameters from a uniform distribution over the asymptotic space.
2. We compute the norming constants {b i } and then successively construct the Jost solutions {J i } as well as the Darboux operators {D i }.
3. We compute the N -soliton solution in the form where the {K i }'s, as the {D i }'s, depend on the spectral parameter η i , the norming constant b i and on J i−1 (η i ) in a non-trivial manner [58].
We stress once again that N -soliton solutions thus produced are exact solutions of KdV, approximating the soliton gas in the large N limit. This provides precise control over the numerically generated soliton gas, while also being computationally faster than numerically solving KdV through pseudo-spectral methods. However, this scheme necessitates an extremely high precision: though no rigorous stability analysis has been performed a good rule of thumb seems to be that generating N -soliton solutions requires 2N digits of precision, making the scheme unwieldy for N greater than a few hundreds. That is why, in Appendix A, rather than construct a soliton gas featuring thousands of solitons, we make use of ergodicty and average over several realisations of 150-soliton or 200-soliton solutions.