Soliton gas in integrable dispersive hydrodynamics

We review spectral theory of soliton gases in integrable dispersive hydrodynamic systems. We first present a phenomenological approach based on the consideration of phase shifts in pairwise soliton collisions and leading to the kinetic equation for a non-equilibrium soliton gas. Then a more detailed theory is presented in which soliton gas dynamics are modelled by a thermodynamic type limit of modulated finite-gap spectral solutions of the Korteweg-de Vries and the focusing nonlinear Schr\"odinger equations. For the focusing nonlinear Schr\"odinger equation the notions of soliton condensate and breather gas are introduced that are related to the phenomena of spontaneous modulational instability and the rogue wave formation. Integrability properties of the kinetic equation for soliton gas are discussed and some physically relevant solutions are presented and compared with direct numerical simulations of dispersive hydrodynamic systems.


Integrable turbulence and soliton gas
Random nonlinear dispersive waves have been the subject of an active research in nonlinear physics for more than five decades, most notably in the contexts of water wave dynamics. A significant portion of the work in this direction has been centred around weak wave turbulence [1]. The wave turbulence theory deals with out of equilibrium statistics of incoherent, weakly nonlinear dispersive waves in non-integrable systems. One of the early and most significant results of the wave turbulence theory was the analytical determination by Zakharov [2] of the analogs of the Kolmogorov spectra describing energy flux through scales in dissipative hydrodynamic turbulence. These spectra, called Kolmogorov-Zakharov spectra, were obtained as solutions of the kinetic equations for the evolution of the Fourier spectra of random weakly nonlinear dispersive waves in multidimensional non-integrable systems.
More recently, a new theme in turbulence theory has emerged in connection with the dynamics of strongly noninear random waves described by one-dimensional integrable systems such as the Korteweg-de Vries (KdV) and 1D nonlinear Schrödinger (NLS) equations. This kind of random wave motion in nonlinear conservative systems, dubbed 'integrable turbulence' [3], has attracted significant attention from both fundamental and applied perspectives. The interest in integrable turbulence is motivated by the complexity of many natural or experimentally observed nonlinear wave phenomena often requiring a statistical description even though the underlying physical model is, in principle, amenable to the well-established mathematical techniques of integrable systems theory such as the inverse scattering transform (IST) or finite-gap theory [4], [5], [6]. Indeed, integrable systems are known to capture essential properties of many wave processes occurring in real-world systems [7]. The integrable turbulence framework is particularly pertinent to the description of modulationally unstable systems whose solutions, under the effect of random noise, can exhibit highly complex nonlinear behaviour that can be adequately described in terms of the turbulence theory concepts, such as probability distribution functions, ensemble averages, power spectra etc. [8], [9], [10]. We stress that the term 'turbulence' in this context is understood as complex spatiotemporal dynamics that requires probabilistic description and is not related to the energy cascades through scales, the prime feature of hydrodynamic and weak turbulence.
Localised nonlinear solitary waves are a ubiquitous feature of nonlinear dispersive wave propagation whose discovery dates back to shallow water wave observations by John Scott Russell in 1845 [11]. If the wave dynamics are described by one of the completely integrable equations the solitary waves exhibit particle-like properties such as elastic, pairwise interactions accompanied by certain phase/position shifts. Such solitary waves are called solitons [12] and have been extensively studied both theoretically [13], [4], [14] and experimentally [15]. The main tool for the analysis of integrable nonlinear dispersive PDEs is the IST [16] based on the reformulation of a nonlinear PDE as a compatibility condition of two linear problems: a stationary spectral (scattering) problem and the evolution problem for the same auxiliary function. Within the scattering problem solitons are associated with discrete values of the spectrum, while the integrable evolution preserves these spectral values in time.
Solitons can form ordered macroscopic coherent structures such as modulated soliton trains and dispersive shock waves [17], [18]. Furthermore, solitons can form irregular, statistical ensembles that can be interpreted as soliton gases. The nonlinear wave field in such gases represents a particular case of integrable turbulence, often called soliton turbulence (with the caveat that the latter term has also been used in the context of nonintegrable wave dynamics, see e.g. [19], [20]). Generally, soliton gas and soliton turbulence represent two complementary aspects of the same physical object, the natural counterparts of the wave-particle duality of a single soliton. In the soliton-gas description the focus is on the collective dynamics/kinetics of solitons as interacting (quasi)particles characterised by certain amplitude (or velocity) distribution function, while the soliton turbulence description emphasises the characterics of the random nonlinear wave field associated with the soliton gas, such as probability density function, power spectrum etc. The observations and analysis of irregular soliton complexes in the ocean have been reported in [21], [22]. Recent laboratory experiments on the generation of shallow-water and deep water soliton gases were reported in [23] and [24] respectively. It has also been demonstrated that soliton gas dynamics in the focusing NLS equation enables a remarkably accurate description of the statistical properties of the nonlinear stage of noise-induced modulational instability [25] as well as provides important insights into the dynamical and statistical mechanisms of the spontaneous formation of rogue waves [26], [27].  Analytical description of soliton gases in nonlinear dispersive wave systems was initiated in the Zakharov's 1971 paper [28], where a spectral kinetic equation for KdV solitons was introduced using an IST-based phenomenological 'flea gas' reasoning enabling the evaluation of an effective adjustment to the soliton's velocity in a rarefied gas due to the interactions (collisions) between individual solitons, accompanied by the well-defined phase-shifts. The kinetic equation for a rarefied soliton gas describes the evolution of the distribution function of the solitons with respect to the spectral parameter and the positions of soliton centres, i.e. the density of states (DOS). In a dense gas, however, the solitons exhibit significant overlap and, as a result, are continuously involved in a strong nonlinear interaction with each other. This can be seen in Fig. 1 displaying a laboratory realisation of a dense soliton gas in a viscous fluid conduit [29], a versatile fluid dynamics platform enabling high precision experiments on the generation and interaction of solitary waves that exhibit nearly elastic collisions [30], [31]. One can appreciate that, for a dense gas the particle interpretation of individual solitons becomes less transparent and the wave aspect of the collective soliton dynamics comes to the foreground. Indeed, a consistent generalisation of Zakharov's kinetic equation for KdV solitons to the case of a dense soliton gas has been achieved in [32] in the framework of the nonlinear wave modulation (Whitham) theory [33]. It was proposed in [34], [32] that the soliton gas can be modelled by the thermodynamic type solitonic limit of the spectral finite-gap KdV solutions and their modulations [35]. The resulting kinetic equation has the form of a nonlinear integro-differential equation for the DOS in the IST spectral phase space. The structure of the kinetic equation obtained in [32] suggested that, remarkably, in a dense gas the net effect of soliton interactions can be formally evaluated using the same phase-shift argument that was used in the rarefied gas theory [28]. This observation, termed the collision rate assumption, has enabled an effective phenomenological theory of a dense soliton gas for the focusing NLS equation [36] and more recently, for the defocusing and resonant NLS equations [37]. The phenomenological theory of soliton gas for the focusing NLS equation proposed in [36] has been confirmed and substantially extended in [38] within the framework of the thermodynamic limit of spectral finite-gap solutions of the focusing NLS equation and their modulations. This latter work has revealed a number of new soliton gas phenomena due to a very different structure of the spectral phase space of the focusing NLS equation compared to the KdV equation. In particular, the generalisation of soliton gas, termed breather gas, was introduced by considering a special family of focusing NLS solitonic solutions called breathers. Such a breather gas represents an intriguing type of integrable turbulence observed in the ocean [39] and recently realised numerically [40].
Mathematical properties of the kinetic equation for soliton gas were studied in [41], [42] where it was proved that it admits an infinite series of integrable linearly degenerate hydrodynamic reductions obtained by a multi-component delta-function ansatz for the DOS. A further study of the classical integrability properties of the soliton gas kinetic equation was undertaken in [43]. A very recent paper [44] presents rigorous results for the uniqueness, existence and non-negativity of the solutions to the integral equations for the DOS in the spectral kinetic theory for KdV and focusing NLS soliton gases [32], [38].
Apart from the above line of research on soliton gases inspired by the Zakharov 1971 work there have been many other developments-analytical, numerical and experimental -exploring various aspects of soliton gas/soliton turbulence dynamics in both integrable and nonintegrable classical wave systems (see e.g. [45], [46], [47], [48], [49], [50], [51], [52], [53], [54]). Additionally, there has been a recent surge of related activity in generalised hydrodynamics (see [55], [56], [57] and references therein), where the equations analogous to those arising in the soliton gas theory became pivotal for the understanding of large-scale, hydrodynamic properties of integrable quantum many-body systems.

Dispersive hydrodynamics
This review considers classical soliton gases from the perspective of dispersive hydrodynamics modelled by hyperbolic or elliptic conservation laws regularized by small conservative, dispersive corrections [58], the KdV equation u t + 6uu x + u xxx = 0 being the simplest example. Smallness of the dispersive term is understood in the sense that the typical coherence length of the medium-determined by a balance between nonlinear and dispersive effects-is much less than the scale of initial, boundary, or intermediate flow conditions where long wavelength, nonlinear, hydrodynamic effects dominate. This distinguishes dispersive hydrodynamic problems from typical formulations in classical soliton theory which usually involve variations of the wave field on the scale of the coherence length. The dispersive hydrodynamic framework arises naturally in the description of problems related to the large-scale dynamics of shallow water or internal waves, but also proves to be extremely useful in the modelling of various phenomena in nonlinear optics including the "atom optics" of quantum fluids such as Bose-Einstein condensates.
In the scalar case general dispersive hydrodynamics are described by the equations of the form where F (u) is the nonlinear hyperbolic flux and D[u] is a differential (generally integro-differential) operator, possibly nonlinear, that gives rise to a real-valued dispersion relation ω = ω 0 (k) for linearised waves. Scalar integrable dispersive hydrodynamics of the form (1.1), such as the KdV, modified KdV, Camassa-Holm or Benjamin-Ono equations often arise as small-amplitude, 'unidirectional' approximations of more general Eulerian bidirectional systems (see [33]) (ρu) t + ρu 2 + σP (ρ) x = (D 2 [ρ, u]) x , σ = ±1, (1.2) where D 1,2 [ρ, u] are conservative, dispersive operators, ρ and u are interpreted as a mass density and fluid velocity, respectively and σP (ρ), where P (ρ) > 0 is the pressure law. For σ = +1 this class of equations generalise the shallow water and isentropic gas dynamics equations while encompassing many of the integrable dispersive hydrodynamic models such as the Kaup-Boussinesq system [59], the hydrodynamic form of the defocusing NLS equation [60], the Calogero-Sutherland system describing the dispersive hydrodynamics of quantum many-body systems [61] and many others. For σ = −1 equations (1.2) describe 'elliptic' dispersive hydrodynamics, the most prominent example being the focusing NLS equation describing, in particular, the evolution of weakly nonlinear narrow-band wave packets on deep water and the propagation of light in optical media in the self-focusing, cubic nonlinearity regime. Integrable bidirectional dispersive hydrodynamics admit Lax representation similar to (1.3) with ψ(x, t, λ) being a vector function.
Integrability of dispersive hydrodynamics (1.1) is understood as the existence of the Lax pair, a system of linear differential equations for an auxiliary (generally vector) function ψ(x, t; λ): where L and A are linear differential operators depending on u(x, t) and its derivatives, and λ is a spectral parameter so that the equation (1.1) can be written in the operator form L t + (LA − AL) = 0 under the additional requirement of isospectrality, λ t = 0. For sufficiently rapidly decaying fields u(x, t) → 0 as |x| → ±∞ the spectrum of the operator L, often called the Lax spectrum, consists of two components: discrete and continuous. The points of discrete spectrum correspond to solitons in the asymptotic solution at t 1, while the continuous spectrum corresponds to dispersive radiation. The solutions whose spectrum consists of N discrete points only are called N -soliton solutions.
A prominent example of a multiscale nonlinear wave structure supported by dispersive hydrodynamics is dispersive shock wave exhibiting coherence at both microscopic (soliton) and macroscopic (wave modulation) scales [18]. Contrastingly, soliton gas as a dispersive hydrodynamic structure exhibits coherence at the microscopic scale while being macroscopically incoherent, in the sense that the values of the wave field at two points separated by a distance much larger than the intrinsic dispersive length of the system (the soliton width), are not dynamically related. The source of randomness in dispersive hydrodynamics is typically related to some sort of stochastic large-scale initial or boundary conditions although one can envisage dynamical mechanisms of the effective randomization of the wave field and soliton gas generation [62], [63]. In particular, statistical soliton ensembles can be naturally generated from both non-vanishing deterministic (e.g. quasiperiodic) and random initial conditions via the processes of soliton fission [64], [65] or modulation instability [25].
The inherent scale separation in dispersive hydrodynamic flows suggests the use of asymptotic methods for their description. One such method, the Whitham modulation theory [33] proved particularly effective and has been extensively developed in the contexts of both integrable and nonintegrable systems. Within the Whitham theory the dispersive hydrodynamic system is asymptotically reduced to a system of quasilinear equations, which describe large-scale variations of the wave's modulation parameters such as the amplitude, the wavenumber, the mean etc. The Whitham hydrodynamic equations can be derived by applying a multiple-scale singular perturbation theory or, equivalently, by averaging the dispersive hydrodynamic conservation laws over rapid oscillations. Although the Whitham method does not rely on integrability, the presence of integrable structure in the original dispersive hydrodynamic system greatly enhances the structure of the modulation system as well. The inherited integrability of the Whitham modulation equations is realised via the generalised hodograph transform [66], [67]. Importantly, the Whitham theory for integrable dispersive PDEs admits spectral formulation within the extension of IST called the finite-gap theory [4]. The finite-gap theory describes an important class of periodic or quasiperiodic (multiphase) solutions to integrable systems, whose Lax spectrum (the set of admissible values of the spectral parameter λ in (1.3) corresponding to L 2 eigenfunctions ψ) lies in the finite union of disjoint bands. In 1980 Flaschka, Forest and McLaughlin showed that the Whitham equations obtained via multiphase averaging of the KdV equation describe slow evolution of the band spectrum of the finite-gap Lax operator L [35]. It has then been shown by Lax and Levermore [68] that the spectral Whitham equations derived in [35] also describe the weak, distribution limits of dispersive conservation laws obtained in the limit of small dispersion.
It has been noted by P. Lax in [69] that the Whitham modulation equations can be viewed as certain analogs of the moment equations in classical statistical fluid mechanics [70]. This observation, combined with the spectral formulation of the Whitham theory [35] and the ideas of Venakides on the continuum limit of theta-functions [71] has enabled the development of the spectral theory of soliton gas [34], [32], [38] providing the justification of the phenomenological approach of [36] and a general mathematical framework for the study of soliton gases/soliton turbulence in classical integrable dispersive hydrodynamic systems.
The structure of the review is as follows. In Section 2 we present a phenomenological theory of soliton gas in unidirectional and bidirectional systems of integrable dispersive hydrodynamics using the KdV and NLS equations as prototype examples. The main result is the construction of the kinetic equation describing the evolution of the density of states in a soliton gas. In the bi-directional case we identify two different types of soliton gases: isotropic and anisotropic, differing in the properties of overtaking and head-on solitons. In Section 3 we develop a detailed spectral theory of soliton gas for the KdV equation and of soliton and breather gases for the focusing NLS equations. This is done by applying a special thermodynamic type limit to the spectral finite-gap solutions and their modulations described by the appropriate Whitham equations. For the focusing NLS equation we introduce the notion of soliton condensate and also consider three types of 'rogue wave' breather gases. In Section 4 hydrodynamic reductions of the spectral kinetic equation are derived by employing a multi-component delta-function ansatz for the density of states and their integrability is investigated. In Section 5 we apply the hydrodynamic reductions to consider a canonical Riemann problem describing collision of two multi-component soliton gases. Finally, the Conclusion provides a brief summary and outlines some directions of future research.

Kinetic equation for soliton gas: phenomenological construction
We first introduce soliton gas phenomenologically, as an infinite ensemble of interacting solitons randomly distributed on the line with non-zero density. This intuitive definition lacks precision but it is sufficient for the purposes of this section. A more elaborate mathematical model of a soliton gas will be described in Section 3.

Unidirectional soliton gas
It is convenient to introduce the basic ideas of soliton gas theory using the KdV equation as a prototype model of nonlinear dispersive wave propagation in a broad variety of physical systems. We consider the KdV equation in the 'physical' form The KdV equation (2.1) belongs to the family of completely integrable equations. For a broad class of initial conditions its integrability is realised via the IST method [4]. The inverse scattering theory associates soliton of the KdV equation (2.1) with a point of discrete spectrum λ n , of the Schrödinger operator, which is the Lax operator for the KdV equation (cf. (1.3)), Assuming u → 0 as x → ±∞ the KdV soliton solution corresponds to λ i = −η 2 i , η i > 0, is given by where the soliton amplitude a i = 2η 2 i , the speed s i = 4η 2 i and x 0 i is the initial position or 'phase'. 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 }. Since the discrete spectrum of the (self-adjoint) Schrödinger operator (2.2) is non-degenerate (i.e. η i = η j ⇐⇒ i = j) all solitons in the N -soliton solution have different velocities and the long-time asymptotic solution assumes the form of a rank-ordered soliton train, where the phase (position) shifts ∆ i occur due to the interaction of individual solitons at the initial stage of the evolution [72], [4]. These phase shifts play important role in our consideration as detailed below. The integrable structure of the KdV equation has profound implications for the soliton interaction dynamics: (i) the KdV evolution preserves the IST spectrum, ∂ t η j = 0, implying that solitons retain their 'identity' (amplitude, speed) upon interactions; (ii) the collision of two solitons with spectral parameters η i and η j , i = j results in their phase (position) shifts given by so that the taller soliton acquires shift forward and the smaller one -shift backwards; (iii) solitons interact pairwise, i.e. the resulting, accumulated 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 the far-field effects and mathematically, the 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 [73] and the resulting wave field in the interaction region cannot be represented as a superposition of the phase-shifted one-soliton solutions. Motivated by the above properties of N -soliton solutions we now introduce a rarefied soliton gas by generalising the asymptotic soliton train solution (2.4) in the following way. We let N → ∞ in (2.4) and introduce the density of states (DOS) f (η, x, t) such that f (η 0 , x 0 , t 0 )dηdx is the number of solitons found at t = t 0 in the element [η 0 , η 0 + dη] × [x 0 , x 0 + dx] of the phase space S = Γ × R, where Γ = [η min , η max ] ⊂ R + is the spectral support of DOS (it is assumed in the above definition of DOS that the interval [x 0 , x 0 + dx] contains a sufficiently large number of solitons). Also, without loss of generality one can assume Γ = [0, 1]. The parameter controlling the total spatial density of the soliton gas is then For a rarefied gas β 1 (this criterion is understood in the asymptotic sense since the actual, numerical, value of β depends on the definition of the unit interval of x. We further assume (to be validated later) the Poisson distribution with density β for the soliton centres x j = x 0 j + 4η 2 t ∈ R so that the distances d k = x k − x k−1 between solitons are independent random values distributed with probability density P(d) = β exp (−βd). We thus arrive at a 'stochastic soliton lattice' that can be viewed as an approximate model of a rarefied soliton gas. Each realisation of the random process u ∞ (x, t) (2.8) satisfies the KdV equation (2.1) almost everywhere in x ∈ R. Due to the small spatial density β, the individual solitons in a soliton gas overlap only in the regions of their exponential tails, except for the rare events of their collisions where a complex nonlinear interaction occurs [73] affecting the statistical characteristics of the random field (2.8) [74].
We first look at the equilibrium, or uniform, soliton gas for which ∂f /∂x = 0. In view of isospectrality of the KdV evolution, the spatially independent density of states f (η) at t = 0 also implies that ∂f /∂t = ∂f /∂x = 0 ∀t > 0.
Consider propagation of a 'tracer' soliton with spectral parameter η = η 1 ∈ [0, 1] in a uniform soliton gas with a given DOS f (η), η ∈ [0, 1] and assume that the gas is rarefied, 1. Due to the collisions of the tracer 'η 1 -soliton' with the 'µ-solitons' (µ = η 1 ) in the gas, each collision leading to the phase shift (2.5), the effective (mean) velocity of the trial soliton is approximately evaluated as The integral correction term in (2.9) follows from the continuum limit of eq. (2.6) assuming that in a rarefied gas s(η 1 ) = 4η 2 1 to leading order. Then the total spatial shift of η 1 -soliton over the time interval dt due to the interactions with µ-solitons, µ ∈ [0, 1] is given by dµ is the average collision rate of η 1 -soliton with µ-solitons having their spectral parameter µ ∈ [µ 0 , µ 0 + dµ]. Then using the expression (2.5) for the KdV phase shift we obtain (2.9). The effective velocity s(η 1 ) has a natural interpretation of the transport velocity in soliton gas.
We note that the spectral parameter η 1 of the 'special' soliton should not necessarily belong to the support Γ = [0, 1] of f (η). If η 1 / ∈ Γ we shall call such a soliton a 'trial' soliton. The important distinction between tracer and trial solitons will become more transparent later. The modification of the 'free' velocity of a trial soliton due its propagation through a uniform soliton gas is illustrated in Fig. 2.
We now consider a non-equilibrium (non-uniform) soliton gas with f (η) ≡ f (η, x, t), s(η) ≡ s(η, x, t) so that spatiotemporal variations of f and s occur on the scales much larger than those associated with variations of the nonlinear wave field u(x, t) in (2.8). Isospectrality of the KdV dynamics implies the conservation equation for the density in the phase space S (i.e. DOS), which, together with the expression for the effective transport velocity s(η) given by (2.9), can be viewed as kinetic equation for rarefied soliton gas first introduced by Zakharov in 1971 [28].  We shall see later that this analogy has a more precise meaning, with a hierarchy of spatiotemporal scales involved. Zakharov's approximate equation (2.9) for the effective transport velocity in a soliton gas was generalised in [32] to the case of dense (β = O(1)) gas. This was done by evaluating the thermodynamic limit of the nonlinear dispersion relations associated with the spectral finite-gap solutions of the KdV equation (2.1), see Section 3.2 below. The result has the form of a linear integral equation which essentially provides the relation between the spectral flux density v = f s and the DOS f in the transport equation (2.10) and so can be viewed as the equation of state of the soliton gas. The expression (2.9) for the effective soliton velocity in a rarefied soliton gas represents an approximate first order solution of the equation of state (2.11) obtained by assuming that the integral term is a small correction to the free soliton velocity 4η 2 . One can notice that equation (2.11) can be formally written down by utilising the same phase-shift argument used in the derivation of the approximate expression (2.9), i.e. by a formal replacement in (2.9) of the leading order value 4η 2 for the soliton velocity with its effective value s(η). The validity of this replacement in the KdV equation case suggests the general collision rate assumption whereby the total position shift of the soliton with η = η 1 due to soliton collisions in a gas with DOS f (µ) over the time interval dt is given by (2.12) where Γ is a spectral support of the DOS f (η). This assumption was used in [36] for the construction of the kinetic equation for the dense soliton gas for the focusing NLS equation, see Section 2.2.2 below. We note that in a different context, the collision rate assumption is at heart of the generalised hydrodynamics (GHD), the theory for the large-scale dynamics of quantum many-body integrable systems. In GHD the pointlike quasiparticles are subject to the instantaneous velocity-dependent spatial shifts upon colliding (see [55], [56], [76], [77], [57], [78] and references therein). It is important to stress, however, that the validity of (2.12) in the context of classical soliton gases, although intuitively suggestive, is far from being obvious. Indeed, as we have already mentioned, the very notion of the phase shift in 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. The fact 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 is quite remarkable. Along with the notion of the phase shift, the phenomenological definition of the DOS f (η) introduced above for a rarefied gas also requires a more careful treatment as the procedure of identifying individual solitons within a dense soliton gas is not obvious (as a matter of fact, the dense soliton gas is no longer described by an approximate 'soliton lattice' expression (2.8)). This issue can be partially addressed by assuming that for any sufficiently broad interval x ∈ [x 0 , x 0 + L], L 1 the soliton gas can be approximated by an exact N -soliton solution and then using the 'windowing' procedure introduced in [37] to 'release' the solitons contained within this interval and count them when they get sufficiently separated. This procedure will be described in Section 2.3. Later, in Section 3 the DOS for dense soliton gases for the KdV and FNLS equations will be introduced in a more mathematically satisfactory way via the thermodynamic limit of the finite-gap spectra of nonlinear multiphase solutions.
With all the above caveats, the equation of state for a general unidirectional soliton gas can be formulated. Let the soliton solution of integrable dispersive hydrodynamics be characterised by the value λ = λ i of the discrete spectrum of the Lax operator, and the position shift in the overtaking two-soliton collisions be ∆(λ, µ) = sgn(λ − µ)G(λ, µ), where G(λ, µ) > 0. Then we obtain upon using (2.12), under the additional assumption s (λ) > 0 (to be verified in concrete cases). Here s 0 (λ) is the velocity of a free soliton and Γ is the support of the DOS f (λ). We use the general notation λ for the spectral parameter in (2.13) with the understanding that it could be some function of the eigenvalue of the Lax operator (as is the case for the KdV equation, cf. (2.2), (2.3)). By re-arranging integral equation (2.13) a useful representation for s(λ) can be obtained (2.14) Although the equation of state (2.13), as written, implies that s(λ) is defined for λ ∈ Γ, its consequence (2.14) can be used to find the speed s(λ 1 ) of the 'trial' soliton with λ = λ 1 / ∈ Γ propagating through the soliton gas with a given DOS f (λ) and transport velocity s(λ), λ ∈ Γ. This extension can be readily justified by formally replacing f (λ) → f (λ) + wδ(λ − λ 1 ) in (2.13) and subsequently letting w → 0.
It is suggested by (2.14) that the case of soliton gas with f (λ), s(λ) satisfying is special. Indeed, this case corresponds to the peculiar type of soliton gas termed soliton condensate [38]. Soliton condensate will be considered in Section 3.3.4 in the context of the focusing NLS equation.

Bidirectional soliton gas
In this section, following [37], we will derive the kinetic equation for integrable bidirectional Eulerian dispersive hydrodynamics (1.2) using the phenomenological construction outlined in the previous section. For convenience, we reproduce system (1.2) here in a slightly simplified form covering majority of dispersive hydrodynamic systems arising in applications: where D[ρ, u]) is dispersion operator that gives rise to a real-valued dispersion relation ω = ω 0 (k) for linearised solutions. Bidirectionality of (2.16) implies that the linear dispersion relation has two branches: ω − 0 (k) and ω + 0 (k)-corresponding to slow and fast waves respectively, i.e. ω − 0 (k)/k < ω + 0 (k)/k in the long wavelength limit k → 0. Suppose that system (2.16) supports a family of bidirectional soliton solutions that bifurcate from the two branches of the linear wave spectrum ω = ω ± 0 (k). We denote the corresponding soliton families (ρ − s , u − s ) and (ρ + s , u + s ). Let these soliton solutions be parametrised by a realvalued spectral (IST) parameter λ so that λ ∈ Γ + for the "fast" branch and λ ∈ Γ − for the "slow" branch, where Γ ± are simply-connected subsets of R with one intersection point at most. Let the respective soliton velocities be c ± (λ). For convenience we assume that c ± (λ) > 0, and The generalisation to the case c ± (λ) < 0, c − (λ 1 ) > c + (λ 2 ) will be straightforward.
One can distinguish between two types of pairwise collisions in a bidirectional soliton gas: the overtaking collisions between solitons belonging to the same spectral branch and characterised by the position shifts ∆ ++ and ∆ −− respectively, and the head-on collisions between solitons of different branches, characterized by the position shifts ∆ +− and ∆ −+ . Let λ = µ, and ∆ ±± (λ, µ) and ∆ ±∓ (λ, µ) denote the position shifts of a λ-soliton due to its collision with a µ-soliton, with the first and the second signs ± in the subscript indicating the branch correspondence of the λ-soliton and the µ-soliton respectively, e.g. ∆ −+ (λ, µ) is the position shift of a λ-soliton with λ ∈ Γ − in a collision with a µ-soliton with µ ∈ Γ + . The typical wave patterns in the overtaking and head-on soliton interactions in a bidirectional shallow-water soliton gas are shown in Fig Figure 3: Soliton collisions in a bidirectional shallow-water soliton gas: a) overtaking; b) headon.

Isotropic and anisotropic soliton gases
We call the bidirectional soliton gas isotropic if the position shifts for the overtaking and head-on collisions between λ-and µ-solitons satisfy the following sign conditions: i.e. the λ-soliton experiences a shift of a certain sign, say shift forward (and the µ-solitonthe shift of an opposite sign) irrespectively of the type of the collision-overtaking or head-on. If conditions (2.17) are not satisfied, i.e. the sign of the phase shift depends on the type of the collision, we shall call the corresponding soliton gas anisotropic. The difference between the phase shifts in isotropic and anisotropic collisions is illustrated in Fig. 5 using concrete examples. Following the construction for unidirectional soliton gas outlined in Section 2.1, we now consider bidirectional soliton gases for integrable Eulerian equations (1.2). We introduce two separate DOS's f − (λ, x, t) and f + (λ, x, t) for the populations of solitons whose spectral parameters belong to the slow (Γ − ) and fast (Γ + ) branches of the spectral set Γ respectively. The isospectrality of integrable evolution implies now two separate conservation laws: where s − (λ, x, t) and s + (λ, x, t) are the transport velocities associated with slow and fast spectral branches Γ − and Γ + respectively. We derive the equations of state for s ± by extending the phenomenological approach of Section 2.1 based on the collision rate assumption. Consider a λ-soliton from the slow branch, λ ∈ Γ − , and compute its displacement in a gas over a 'mesoscopic' time interval dt, sufficiently large to incorporate a large number of collisions, but sufficiently small to ensure that the spatiotemporal field f ± (λ, x, t) is stationary over dt and homogeneous on a typical spatial scale c ± (λ)dt. Having this in mind, we drop the (x, t)-dependence for convenience. Each overtaking collision with a µ-soliton of the same branch, µ ∈ Γ − , shifts the λ-soliton by the distance ∆ −− (λ, µ). Then, invoking the collision rate assumption (2.12) the displacement of the λ-soliton over the time dt due to the overtaking collisions with µ-solitons, where µ ∈ [µ 0 , µ 0 + dµ], is given by Additionally, each head-on collision with a fast soliton, µ ∈ Γ + , shifts the slow λ-soliton with λ ∈ Γ − by ∆ −+ (λ, µ), and the resulting displacement after a time dt is given Γ + ∆ −+ (λ, µ)f + (µ)|s + (λ) − s − (µ)|dt dµ. A similar consideration is applied to the fast soliton branch, λ ∈ Γ + , in the gas. Equating the total displacements of the slow and fast λ-solitons to s − (λ)dt and s + (λ)dt respectively, we obtain the equation of state of a bidirectional gas in the form of two coupled linear integral equations: where λ ∈ Γ − for the first equation and λ ∈ Γ + for the second equation.
If the spectral support Γ = Γ − ∪ Γ + ⊂ R is a simply connected set and the gas is isotropic, the distinction between the fast and slow branches becomes unnecessary and the kinetic equation (2.18), (2.19) for bidirectional soliton gas is naturally reduced to the unidirectional gas equation (2.10),(2.13) for a single DOS f (λ) defined on the entire set Γ. It was shown in [37] that the dynamics governed by the kinetic equations (2.10),(2.13) and (2.18), (2.19) are in a very good agreement with the results of direct numerical simulations of isotropic and anisotropic bidirectional soliton gases respectively.

Kinetic equations for soliton gas in NLS dispersive hydrodynamics
As a representative example of bidirectional integrable dispersive hydrodynamics (2.16), we consider the system For µ = 1, system (2.20) is equivalent to the NLS equation: with the mapping between the two representations being realised by the so-called Madelung transform: The case σ = +1 in (2.21) corresponds to the defocusing NLS equation describing the propagation of light beams through optical fibres in the regime of normal dispersion, as well as nonlinear matter waves in quasi-1D repulsive Bose-Einstein condensates, see for instance [79]. If σ = −1, equation (2.21) is the focusing NLS equation, which is a canonical model for the description of modulationally unstable wave systems such as deep water waves or the propagation of light in optical fibres in the regime of anomalous dispersion.
For σ = +1, µ = −1 system (2.21) is equivalent to the NLS equation in a 'quantum potential' also known as the resonant NLS equation [80], [81], This equation, in particular, describes long magneto-acoustic waves in a cold plasma propagating across the magnetic field [82]. It is also directly related to the Kaup-Boussinesq system, an integrable model for bidirectional shallow water waves [59] (see [37] for the transformation between the resonant NLS and the Kaup-Boussinesq system).
We now look at the soliton gas descriptions for each of the above NLS systems The inverse scattering theory of the defocusing NLS equation was constructed in [83]. It was shown that the defocusing NLS equation supports a family of dark (or grey) soliton solutions on the finite background, which, up to the initial position and phase, can be most conveniently represented in terms of the hydrodynamic variables ρ, u as where λ ∈ R is the discrete spectral parameter in the linear scattering problem associated with the defocusing NLS equation [83] , Γ − = (−1, 0] for the slow soliton branch and Γ + = [0, +1) for the fast soliton branch; note that solutions (ρ + s , u + s ) and (ρ − s , u − s ) have the same analytical expression. Also, despite the same analytical expression for the soliton velocities c + and c − in the defocusing NLS case we keep the the formal notational distinction to retain the connection with general equation of state (2.19) for bidirectional gas, where the expressions for c + (λ) and c − (λ) can be, in principle, different. Additionally, without loss of generality we assumed in (2.24) the unit density background.
Typical dark soliton solutions for ρ and u are displayed in Fig. 4a. The position shifts in the defocusing NLS overtaking and head-on soliton collisions are given by the same analytical for all λ, µ ∈ (−1, 1). Expressions (2.25) were obtained in [83]. Note that, unlike in the KdV equation, one needs to distinguish between the position and phase shifts for soliton collisions in the NLS dispersive hydrodynamics. The expressions for the phase shifts can also be found in [83] but we do not consider them here. One can readily verify that the soliton position shifts given by (2.25) satisfy the isotropy conditions (2.17). The variation of ∆(λ, µ) with respect to µ for a fixed λ is displayed in Fig. 5a. One can see that the position shifts for the head-on and overtaking collisions lie on the same curve with ∆(λ, µ) being continuous at λ = 0, the point of intersection of Γ − and Γ + . Due to the isotropic nature of the defocusing NLS soliton interactions the coupled kinetic equation (2.18), (2.19) for the bidirectional defocusing NLS gas reduces to the single spectral transport equation (2.10) with the equation of state One can verify by direct computation that the condition s (λ) > 0 necessary for the validity of (2.26) is satisfied.
The resonant NLS equation (2.23) is reducible to the well-known integrable Kaup-Boussinesq system for shallow water waves for which the inverse spectral theory was constructed in [59]. It is not difficult to show that the resonant NLS equation has a family of spectral anti-dark soliton solutions given by [81] However, now one can verify that, unlike in the defocusing NLS case, the isotropy condition , that is in a head-on collision between a λ-soliton and a µ-soliton with λ > µ, the λ-soliton's position is now shifted backwards. The variation of ∆ ±± (λ, µ) for the resonant NLS equation is shown in Fig. 5b. One can see that it is qualitatively different from the variation of ∆ ±∓ (λ, µ) for the defocusing NLS equation (Fig. 5a). The kinetic equation for the anisotropic resonant NLS soliton gas then takes the form of two continuity equations (2.18) complemented by the coupled equations of state assuming that s ± (λ) > 0 and s + > s − (verified by direct computation).

(iii) Focusing NLS equation
The case of the focusing NLS equation (equation (2.21) with σ = −1) is special as it represents an example of the focusing dispersive hydrodynamics, where the long-wave, 'hydrodynamic' motion is described by an elliptic system. The focusing NLS equation is a canonical model for the description of modulationally unstable systems in fluid dynamics and nonlinear optics. Nevertheless, the focusing NLS equation supports stable soliton solutions that can propagate in both directions so the focusing NLS soliton gas should be classified as bidirectional.
We consider the focusing NLS equation in the following normalisation: which is standard in the context of the IST analysis. The IST for the focusing NLS equation was introduced in the celebrated paper of Zakharov and Shabat [84], where it was shown that the spectral problem associated with (2.30) has the form where λ is the spectral parameter and Y = Y (x, t, λ) is a vector; ψ denotes complex conjugate. The Lax operator L in (2.31) is often called the Zakharov-Shabat operator. The main difference of the Zakharov-Shabat operator from the Lax operators for the KdV, defocusing and resonant NLS equations is that it is not self-adjoint, implying that the spectral parameter λ is complex, λ ∈ C.
A single-soliton solution of equation (2.30) is characterised by a discrete complex eigenvalue, λ 1 = a + ib and c.c., of the Zakharov-Shabat operator and is given by where x 0 is the initial position of the soliton and φ 0 is the initial phase. One can see that the focusing NLS soliton represents a localised wavepacket with the envelope propagating with the group velocity c g = −4a = −4Reλ 1 and the carrier wave having the phase velocity c p = (b 2 − a 2 )/a = −2Re(λ 2 1 )/Reλ 1 . Also, similar to the defocusing and resonant NLS equations, and in contrast with KdV equation, the amplitude and velocity of the focusing NLS soliton are two independent parameters. We identify the two families (±) of focusing NLS solitons according to the sign of their group velocity, sgn(c g ) = −sgn(Reλ).
Similar to other integrable NLS models, the solitons of the focusing NLS equation interact pairwise and experience both position and phase shifts upon the interaction. The position shifts in the focusing NLS overtaking and head-on soliton collisions are given by the same expression, One can see that the position shifts (2.33) satisfy conditions (2.17) so we classify the focusing NLS soliton collisions as isotropic.
The DOS f (λ) of the focusing NLS soliton gas is generally supported on some compact Schwarz symmetric 2D set Λ ⊂ C so it is sufficient to consider only the upper half plane part Λ + . Here Schwarz symmetry means that if λ ∈ C is a point of the spectrum then so is the c.c. pointλ. The kinetic equation for the focusing NLS soliton gas then assumes the form [36] f t + (f s) x = 0, The case Λ + ⊂ iR + requires a separate consideration. Solitons of the focusing NLS equation can form stationary complexes described by the special N -soliton solutions, called bound states, for which all discrete spectrum points are located on the imaginary axis [84]. Since for the corresponding bound state soliton gas Reλ = 0 the equation of state in (2.34) immediately yields s(λ) = 0 resulting in the equilibrium DOS, f t = 0. It was shown in [25] that the turbulent wave field ψ(x, t) associated with the DOS f (λ) in a dense bound state soliton gas exhibits very peculiar statistical properties, shedding light on the fundamental phenomenon of spontaneous modulational instability. This special soliton gas will be considered in Section 3.3.4.

Ensemble averages and modulation equations for soliton turbulence
Within the spectral kinetic description of soliton gas the solitons are viewed as quasiparticles moving with the speeds determined by the nonlocal equation of state (2.13). The DOS f (λ) in this description represents a comprehensive spectral characteristics of soliton gas. However, of ultimate interest in dispersive hydrodynamics is the description of the turbulent nonlinear wave field u(x, t) associated with the underlying spectral soliton gas dynamics. A natural question arises then: how to solve the inverse scattering problem i.e. how to 'translate' the spectral characterisation of soliton gas (i.e. the DOS) into the statistical characterisation (ensemble averages, probability density function, correlations, etc.) of the associated nonlinear random wave field u(x, t) satisfying an integrable dispersive hydrodynamic equation, e.g. the KdV equation?
Within the classical, deterministic IST setting the inverse spectral problem is solved via the Gelfand-Levitan-Marchenko equation (see e.g. [72], [85], [4]) or utilising the Riemann-Hilbert problem approach (see e.g. [86]). While the construction of a comprehensive extension of the IST for random potentials seems to be far away at present (see e.g. [87], [88] for some of the important developments), some valuable quantitative results can be obtained by elementary means. We shall describe some of these results using the KdV equation as the simplest accessible example, the generalisation to other integrable unidirectional and bidirectional dispersive hydrodynamic equations being straightforward.
As is well known (see e.g. [70], [1]) a turbulent wave field is characterised by the moments u n over the statistical ensemble, which, assuming ergodicity, can be computed as spatial averages u n = 1 0 u n (x, t)dx, over a sufficiently large 'mesoscopic' interval l L, where l is the typical scale for spatial variations of u(x, t) in a soliton gas, i.e. the typical soliton width, and L is the typical scale for spatial variations of the density of states f (λ, x, t) in the non-uniform soliton gas. One of the fundamental properties of integrable dispersive hydrodynamics is the availability of an infinite set of local conservation laws (c) Variation of the truncated distributionũ(y, τ ) at time τ t * . Figure 6: Schematic for the evaluation of the integral (2.37) in soliton gas using the 'windowing' procedure.
where the P i and Q i are functions of the field variable u and its derivatives. For decaying initial data, the integrals ∞ −∞ P i dx are conserved in time. For the KdV equation the existence of an infinite series of local polynomial conservation laws was established in [89]. Their conserved densities sometimes called Kruskal integrals can be deduced from the Lax pair via appropriate expansions for large values of the spectral parameter, see e.g. [4], [90]. E.g. for KdV (2.1) Here we describe a simple heuristic approach proposed in [37] that enables one to link the spectral DOS f (λ) of a soliton gas with the ensemble averages of P i of the integrable dispersive hydrodynamic system (1.1). We describe the general principle using the KdV equation as a prototype example.
Consider a uniform KdV soliton gas, i.e. a gas whose statistical properties, particularly the DOS f (η), do not depend on x, t (we return to the conventional use of the spectral variable η = √ −λ in the context of the KdV solitons, see Section 2.1). We now make a natural assumption that the nonlinear wave field in a homogeneous soliton gas represents an ergodic random process, both in x and t (we note in passing that ergodicity is inherent in the spectral model of soliton gas that will be presented in Section 3). The ergodicity property implies that the ensemble-averages P i [u] in the soliton gas can be replaced by the corresponding spatial averages. Generally, for any functional H[u(x, t)] we have for a single representative realization of soliton gas. We define where P j [u] is the conserved density and L 1.
. Let u(y, t) be a representative realization of a soliton gas and letũ(y, t) be defined in such a way that for some t = t * one has (ũ(y, t * ) = u(y, t * ) for y ∈ (x−L, x+L) and (ũ(y, t * ) = 0 outside of this interval so that the conserved densities P i [ũ(y, t * )] and respective fluxes Q i [ũ(y, t * )] also vanish outside (x − L, x + L). To avoid complications we assume that the transition between the two behaviors is smooth but sufficiently rapid so that such a 'windowed' portion of a soliton gas (see Fig. 6a,6b) can be approximated by the spectral N -soliton solution for some N 1, with the discrete IST spectrum points η i , i = 1, . . . N being distributed on Γ = [0, 1] with density ϕ(η) ≈ 2Lf (η) (recall the definition of DOS in Sec. 1).
The integrals (2.37) then can be re-written as Since the integrals (2.38) do not depend on time they can be computed at t = τ t * when the 'windowed' solutionũ(x, t) gets resolved into a train of N well-separated solitons u s (x, t; η i ), i = 1, . . . , N , see Fig. 6c. Assuming that the overlap between (exponentially decaying) solitons in the train becomes negligible as t → ∞ the integral (2.38) can be represented as a sum of integrals over the individual soliton solutions, (2.39) E.g. for the KdV equation the two first integrals in (2.39) are readily evaluated using the one-soliton solution (2.3) to give: These are the 'spectral mass' and 'spectral momentum' of the η-soliton respectively. Now, using (2.3) and taking the continuous limit as N → ∞ in (2.39) according to i F (η i ) → Γ F (η)ϕ(η)dη , we obtain the expressions for the first two statistical moments in the KdV soliton turbulence defined according to (2.36): One fundamental restriction imposed on the distribution function f (η) follows from nonnegativity of the variance Some consequences of this restriction have been explored in [91]. In particular, considering the 'cold' soliton gas characterised by the Dirac delta-function DOS, f (η) = wδ(η − η 0 ) we obtain from (2.41) u = 4η 0 w, u 2 = 16 3 (η 0 ) 3 w. Then the condition (2.42) yields the constraint We note that the method presented here only requires one to integrate the single-soliton solution and thus can be readily applied to any integrable dispersive hydrodynamic system supporting the soliton resolution scenario. In particular, the ensemble averages for the density ρ , velocity u and momentum ρu in the bidirectional soliton gases for the defocusing and resonant NLS equations were obtained in [37]. Next, we note that the above simple derivation implies that the expressions (2.41) apply to both dense and rarefied gas. Indeed, the same expressions were obtained in [50] for a rarefied KdV gas (see also [92] for the similar modified KdV equation results). In Section 3 we will show how these expressions follow from the general spectral construction of soliton gas, not invoking the heuristic 'windowing' procedure.
In the above consideration of homogeneous soliton gases the ensemble averages (2.36) are constant. For a non-uniform (non-equilibrium) gas the DOS is a slowly varying function of x, t and so are the ensemble averages that now need to be interpreted as "local averages" in the spirit of modulation theory [33]. Essentially, one invokes the "mesoscopic" scale : l L -so that the DOS is approximately constant on any interval (x − , x + ). Then the constant ensemble averages (2.36) are replaced by the slowly varying quantities: The 'local' averages H[u] do not depend on at leading order, and their spatiotemporal variations occur on x, t-scales that correspond to the scales associated with variations of f (η) and are much larger than typical (fast) variation scales of the nonlinear wave field u(x, t). Using the scale separation and applying the ensemble averaging to the infinite set of conservation laws for integrable dispersive hydrodynamics we obtain the modulation system for soliton gas Equations (2.45) can be viewed as a stochastic version of the Whitham modulation equations [33] first discussed in [69]. Equations (2.45) can also be viewed as the integrable turbulence counterpart of the moment equations in the classical hydrodynamic turbulence theory [70]. As a matter of fact the infinite-component stochastic modulation system (2.45) is consistent with the integro-differential kinetic equation (2.10), (2.11). We note that the local ensemble averages of the infinite set of conserved densities P i [u] (x, t) satisfying system (2.45) can be evaluated in terms of the spectral moments of the DOS:

Nonlinear spectral theory of soliton gas
The phenomenological kinetic theory of soliton gas described in the previous section is essentially based on the interpretation of solitons as quasi-particles experiencing short-range pairwise interactions accompanied by the well-defined phase/position shifts. As was already stressed, although this theoretical framework is justifiable in the case of a rarefied gas, it is not quite satisfactory for a dense gas where solitons experience significant overlap and continual nonlinear interactions so that they could become indistinguishable as separate entities. Clearly it would be desirable to have a more mathematically robust approach to the description of a general, dense soliton gas that would, in particular, provide a formal justification of the collision rate assumption (2.12) and, consequently, of the equation of state (2.13), at least in some concrete cases of integrable dispersive hydrodynamics.
The above discussion strongly suggests that we need to look at the wave side of the soliton's 'dual identity'. As we already mentioned, there are insightful parallels between kinetic theory of soliton gas and the nonlinear wave modulation theory introduced by G.B. Whitham in 1965 [93], [33], just before the discovery of the IST method. The Whitham theory describes slow evolution (modulation) of the parameters characterising nonlinear periodic and quasiperiodic waves, such as amplitude, wavenumber, frequency, mean etc. For integrable equations such as the KdV and NLS equations there is an elegant and powerful spectral approach to the derivation of the modulation equations first developed by Flaschka, Forest and McLaughlin for the KdV equation [35] and then extended to other integrable equations such as the NLS, Benjamin-Ono, sine-Gordon equations etc. The spectral modulation theory is based on the extension of the IST method called the finite-gap theory [4], [5]. In this section we will show how the finite-gap theory and modulation equations can be used to construct a mathematical model of soliton gas and justify the kinetic equation (2.10), (2.11) for the soliton gas of the KdV equation and further, equation (2.34) for the soliton gas of the focusing NLS equation.

The Big Picture
We start with a high-level description of the spectral theory of soliton gas based on properties of the so-called finite-gap modulated solutions of integrable dispersive hydrodynamic equations. For simplicity we shall refer to scalar dispersive hydrodynamics (1.1) although the general principles are equally applicable to the vector, bidirectional case. Our qualitative description will be substantiated in the following sections by considering the examples of the KdV and the focusing NLS equations. As we have seen in Section 2, the description of soliton gas includes two aspects: (i) the 'microscopic' structure of an equilibrium (uniform) gas, characterised by the equation of state (2.13); (ii) the slow evolution of the non-equilibrium (non-uniform) soliton gas' parameters described by the spectral transport equation (2.10). This setting bears a strong resemblance to the Whitham modulation theory in which the local microscopic structure is given by an exact periodic or quasiperiodic solution of a dispersive hydrodynamic equation while the slow evolution is described by the modulation equations obtained via some period-averaging procedure or formal multiple scales analysis.

Spectral modulation theory of multiphase waves
The simplest yet instructive example of modulation theory arises when considering the evolution of linear modulated waves in dispersive media, the so-called kinematic wave theory [33]. Dispersive hydrodynamics (1.1) upon linearisation about an equilibrium u = u 0 admits harmonic multiperiodic solutions in the form of a finite Fourier series whereũ = u − u 0 , a j are the amplitudes of the Fourier modes, k = (k 1 , k 2 , . . . , k N ), ω = (ω 1 , ω 2 , . . . , ω N ) are the wavenumber and frequency vectors respectively, and θ 0 = (θ 0 1 , . . . , θ 0 N ) is the initial phase vector. The frequency components in (3.1) satisfy the linear dispersion relation of (1.1), i.e. ω j = ω 0 (k j ). Assuming non-commensurability of the components k j and ω j of the respective wavenumber and frequency vectors in (3.1), this solution is defined on N -dimensional 2π-torus, For a modulated linear wave all parameters in the Fourier series (3.1) become slow functions of x, t. To capture the modulations one introduces a small parameter ε and writes an approximate solution asũ(θ, , where X = εx, T = εt are the slow x and t variables respectively, and θ j (x, t) = ε −1 S j (X, T ) are the fast generalised phases. We now define the local wavenumbers and frequencies by so that in the absence of modulations the approximate solution reduces to (3.1). The phase consistency conditions ∂ xt θ j = ∂ tx θ j then yield N wave conservation laws, System (3.4) represents the simplest example of modulation system. We note that the second part of (3.4) implies that, despite the slow x, t-variations, the components of the frequency vector remain to be locally (on the typical scale of the fast oscillations in θ) related to the wavenumber vector components by the same linear dispersion relations as in the non-modulated wave (3.1). This result is non-trivial and is established by a multiple scales analysis of the linearised dispersive hydrodynamics with the approximate solutionũ(θ, X, T ) as the leading order term in the perturbation series in ε. This analysis also yields the modulation equations for the amplitudes a j (X, T ), which we do not consider here. We also note that the initial phases θ 0 j can be viewed as independent random values, each uniformly distributed on [0, 2π). Then solution (3.1) becomes a random process, whose slow modulations are described by eq. (3.4) complemented by the appropriate amplitude modulation equations.
We now turn to the nonlinear analogue of the multiphase solution (3.1). The one-phase nonlinear periodic solutions are quite common in dispersive hydrodynamics and describe travelling waves. However, the existence of nonlinear multiphase, quasiperiodic solutions is an exceptional feature of integrable PDEs. Let the integrable dispersive hydrodynamics (1.1) admit multiphase solutions 20 i.e. θ = (θ 1 , θ 2 , . . . , θ N ) ∈ T N . Quasiperiodic solution (3.5) can be viewed as a nonlinear analogue of the linearised Fourier solution (3.1). The existence and properties of nonlinear multiphase solutions to the KdV equation were established in 1970-s in the series of pioneering works [94], [95], [96], [97] where the finitegap theory, a nontrivial extension of the IST to periodic and quasiperiodic potentials has been developed (see also the monographs [4], [5] and a historical review [98]). It was shown that the most natural parametrisation of the multiphase solutions (3.5) is achieved in terms of the spectrum of the corresponding Lax operator. For the preliminary discussion of this section it is convenient to assume that the Lax operator is self-adjoint so that its spectrum is real valued. This is the case for the (unidirectional) KdV and (bidirectional) defocusing NLS equations. The case of complex band spectrum arises for the focusing NLS equation, and this case will be considered separately in Section 3.3. The fundamental result of the finite gap theory is that the Lax spectrum S N of the N -phase solution (3.5) lies in the union of N + 1 disjoint bands γ j = [λ 2j−1 , λ 2j ] (one of which could be semi-infinite, then γ N +1 = [λ 2N +1 , +∞)) separated by N finite gaps c j = (λ 2j , λ 2j+1 ), For that reason multiphase solutions of integrable equations are often called finite-gap potentials (we recall that the Lax spectrum of a general periodic potential consists of an infinite number of bands and has even more complicated structure for almost periodic potentials [99]). Thus the real spectrum of a finite-gap potential is fully parametrised by the state vector λ = (λ 1 , λ 2 , . . . , λ D ), where D = 2N + 1 or D = 2N + 2 depending on the presence or absence of a semi-infinite band.
One of the outcomes of the finite-gap theory are the nonlinear dispersion relations linking the physical parameters of the multiphase solution (3.5) such as the wavenumbers, frequencies, mean etc. with the components of the D-dimensional spectral state vector λ. In particular, for the N -component wavenumber and frequency vectors k and ω in (3.5) we have k j = K j (λ), ω j = Ω j (λ), j = 1, . . . , N. (i) Harmonic (linear wave) limit is achieved by collapsing spectral gaps, |c j | = λ 2j+1 − λ 2j → 0, j = 1, . . . , N . In this limit the N -gap solution converts into the linear quasiperiodic solution (3.1), and the nonlinear dispersion relations (3.8) into the dispersion relation ω j = ω 0 (k j ) for linearised wave modes. This is not surprising as the inverse scattering theory, including its finite-gap extension, essentially represents a nonlinear analogue of the Fourier method [85].
(ii) Solitonic limit. In this limit the N -gap quasiperiodic solution (3.5) transforms into the Nsoliton solution exponentially decaying as |x| → ∞ [94] (see also monograph [4]). Spectrally, the solitonic limit is realised by collapsing all finite bands into double points, |γ j | = λ 2j − λ 2j−1 → 0, while keeping the gaps open, i.e. one requires vanishing of the band/gap ratio, r j = |γ j |/|c j | → 0. The collapsed bands correspond to the discrete spectrum points in the traditional IST for decaying potentials. On the other hand, the spectral limit r j → 0 implies vanishing of the corresponding wavenumber, k j → 0, which is another, physically suggestive way to describe the solitonic limit of finite-gap potentials, particularly useful in the context of soliton gases.
As a simple example illustrating the harmonic and solitonic limits in finite-gap potentials we consider the one-gap (single-phase) KdV solution. For N = 1 the KdV Lax spectrum S 1 = [λ 1 , λ 2 ] ∪ [λ 3 , ∞) and the corresponding KdV solution (3.5) assumes the form of a 'cnoidal wave' (see e.g. [72]); without loss of generality we set λ 3 = 0, where cn[·] is the Jacobi elliptic function, m = λ 2 /λ 1 is the modulus, K(m) is the complete elliptic integral of the first kind, and θ = kx − ωt + θ 0 is the phase with the wavenumber k and the frequency ω given by the nonlinear dispersion relations (cf. (3.8)) , ω = −2k(λ 1 + λ 2 ) . (3.10) The band/gap ratio in the one-gap solution is given by r = 1 − m. The solitonic limit r → 0 (λ 2 , λ 1 → −η 2 ) of (3.9), (3.10) is then evaluated using the asymptotic behaviour of elliptic functions [100], so that solution (3.9) takes the form of a soliton (2.3), (3.12) The opposite, harmonic limit of the cnoidal wave solution (3.9) is realised by closing the spectral gap (i.e. letting m → 0) but we do not consider this limit here as it does not play a role in the soliton gas construction.
The above consideration provides a good intuition for what happens in general case N > 1. Indeed, a N -gap KdV solution can be represented as a nonlinear superposition of N cnoidal waves [101], [6] with the interaction between the nonlinear modes described by the so-called Riemann period matrix B (see eq. (3.25) below). In the N -soliton limit (r j → 0 ∀j = 1, 2, . . . , N ), one has (see e.g. [98], [6]) so that the off-diagonal elements of the interaction matrix transform into the normalised twosoliton phase shifts (cf. (2.5)).
Modulation theory of nonlinear multiphase waves F N (θ; λ) describes slow evolution of the endpoints of the band spectrum, λ = λ(X, T ), where X = εx, T = εt, ε 1 [35], [67]. Similar to the linear modulation theory, slow modulations necessitate the introduction of the generalised phase vector θ = S(X, T )/ε so that the local wavenumber and local frequency vectors are defined by (cf. 14) The consistency condition θ xt = θ tx then leads to the system of N wave conservation equations, an analogue of the kinematic equation (3.4) for linear waves with an important difference that, instead of the linear dispersion relation ω = ω 0 (k) linking the components of the frequency and wavenumber vectors it now includes the nonlinear dispersion relations (3.8) for finite-gap potentials, Importantly, the kinematic modulation system (3.15) for a general case of multiphase nonlinear waves is not closed since dim(λ) = D > N . Indeed, the full modulation system contains D equations for λ j (X, T ), and the wave conservation equations are always consistent with (but not equivalent to) the full system (see [35] for the complete description of the KdV spectral modulation theory). However, in the harmonic and soliton limits corresponding to the collapsed spectral gaps or bands the dimension D of the state vector λ decreases, enabling the necessary closure for the system of wave conservation laws in (3.15) under the additional constraint of constant background (see [17], [102] for the relevant theory of the dynamic wave-mean flow interaction showing how the effects of nonconstant background can be included). In particular, in the harmonic limit system (3.15) transforms into the kinematic system (3.4). The solitonic limit is a singular one and requires a more delicate treatment.

Thermodynamic limit of finite-gap spectral solutions
The main idea of the spectral construction of soliton gas is to take simultaneously the solitonic limit and the limit N → ∞ of N -gap potential in such a way that with a similar behaviour for the frequency components ω j . The limit (3.16) is the thermodynamic type limit for nonlinear multiphase waves (as we shall see, β in (3.1.2) agrees with the soliton gas density introduced earlier in (2.7)). This limit suggests the following scaling for the wavenumbers and frequencies: Analysis of the nonlinear dispersion relations (3.8) for the KdV and NLS equations yields the asymptotic structure of the spectral set S N compatible with the thermodynamic scaling (3.17). It turns out that the large N behaviour (3.17) of k j , ω j can be achieved by introducing the special distribution of spectral bands γ j and gaps c j on a fixed spectral interval Γ = [λ 1 , λ D ] (an arc in the complex plane for the focusing NLS) such that i.e. by making the bands exponentially narrow compared to the gaps. We note that the band-gap scaling (3.18) is inspired by the spectrum of the periodic Lax operator (2.2) in the semiclassical limit [103], [71]. We shall call the band-gap scaling (3.18) the exponential thermodynamic spectral scaling and denote the limit as N → ∞ of a function F (λ) considered on this scaling as T-lim F (λ). For the thermodynamic scaling (3.18) we have r j = |γ j |/|c j | → 0 ∀j, yielding essentially an infinitesoliton limit, in full agreement with our intention to describe soliton gas. Other meaningful scalings (sub-exponential, super-exponential) compatible with (3.17) and leading to special cases of soliton gases are possible and will be discussed later.
We will show that the application of the thermodynamic limit to the kinematic modulation system (3.15) results in the kinetic equation for soliton gas fT + (sf )X = 0, s(λ) = S[f (λ)], (3.19) where S[. . . ] is a functional and are the density of states and the transport velocity respectively, both depending on the continuous parameter λ ∈ Γ and on the 'superslow' spacetime variablesX = δx,T = δt, where δ ε 1 is a small parameter that scales the typical spatiotemporal modulations of soliton gas, which are much slower than those associated with the typical Whitham modulations of finite-gap potentials. In practice we will not be introducing the small parameters ε and δ explicitly, assuming that the typical x, t-variations of u and f (λ) occur on disparate micro-and macroscopic scales respectively.
The outlined construction was concerned with the spectral characterisation of soliton gas defined as a thermodynamic limit of finite-gap potentials and can be symbolically represented as S N −→ f (λ). However, this description is incomplete as the Lax spectrum S N determines the finite-gap potential (3.5) only up to N initial phases, θ 0 ∈ T N . For a general finite-gap potential the respective incommensurability of the components of the wavenumber and frequency vectors in (3.5) implies dense winding on the phase torus. Then the natural assumption for the construction of soliton gas would be to let the components of θ 0 be independent random values, each distributed uniformly on [0, 2π). This is the so-called Random Phase Approximation, the standard assumption in the wave turbulence theory [1], which was also used for the construction of random finite-gap solutions of the KdV equation in [104] and of the focusing NLS equation in [105], [106]. Following the classical construction of the configuration space of the ideal onedimensional gas of non-interacting particles (see e.g. [107]) it can be shown that, upon the thermodynamic limit the uniform distribution of the initial phase vector θ 0 over the invariant torus T N transforms into the Poisson distribution on R with the density β = 1 0 f (λ)dλ for the position phases x j = θ 0 j /k j [34], which is consistent with the distribution of the soliton centres in the phenomenologically introduced rarefied soliton gas (2.8). The derivation of the Poisson distribution for the position phases in the KdV soliton gas will be presented in Section 3.2.3.

Thermodynamic limit and nonlinear dispersion relations for soliton gas
We now realise the thermodynamic spectral limit construction of soliton gas for the KdV equation (2.1) following [34], [32]. The Lax spectrum (3.7) of the N -phase KdV solution (3.5) lies in the union of bands, The state vector λ = (λ 1 , λ 2 , . . . , λ 2N +1 ) parametrises the N -gap KdV solution up to N initial phases θ 0 j , which we assume to be independent random values, each uniformly distributed on [0, 2π). In what follows, we take advantage of some known results from the KdV finite-gap theory [4], [5], [35] and apply them to the description of the KdV soliton gas.

(i) Nonlinear dispersion relations for finite-gap potentials
To formulate the nonlinear dispersion relations (3.8) for the multiphase KdV solutions we need to introduce several fundamental objects underlying the algebraic structure of finite-gap potentials expressible in terms of multidimensional Riemann theta-functions [108]. In the finitegap theory the vector λ = (λ 1 , λ 2 , . . . , λ 2N +1 ) of the endpoints of spectral bands defines the two-sheeted hyperelliptic Riemann surface R of genus N via with R ∼ z N +1/2 as z → ∞. We make the branch cuts of R(z) along the spectral bands [λ 1 , λ 2 ], . . . [λ 2j−1 , λ 2j ], . . . , [λ 2N +1 , ∞) and introduce a canonical homology basis on R as follows (see Fig. 7): the α j -cycle surrounds the j-th band clockwise on the upper sheet, and the β j -cycle is canonically conjugated to α j such that the closed contour β j starts at λ 2j , goes to +∞ on the upper sheet and returns to λ 2j on the lower sheet. We introduce a basis of holomorphic differentials on R: where the coefficients κ jk (λ) are determined by the normalisation over the α-cycles β k w j = δ jk , (3.24) while the integrals over the β-cycles give the entries of the symmetric N × N Riemann period matrix B, with positive definite imaginary part. The components k j , ω j of the wavenumber and frequency vectors are expressed in terms of the branch points λ j of the spectral Riemann surface (3.22) by the relations [35]  (ii) Thermodynamic spectral scaling We now introduce the exponential thermodynamic scaling (3.18) for the spectrum (3.21) of the KdV finite-gap solutions. We fix λ 1 = −1 and λ 2N +1 = 0 and, following [71], consider the lattice of points 0 are centres of spectral bands. We then define two positive smooth functions on [0, 1]: (a) The normalised density ϕ(η) > 0 of the lattice points η j , introduced such that i.e. ϕ(η)dη is the probability measure on [0, 1]. (b) The normalised logarithmic band width distribution τ (η) > 0 defined by The functions ϕ(η) and τ (η) asymptotically define the Riemann surface R (3.22) for N 1. The asymptotics (3.28), (3.29) imply the exponential spectral scaling (cf. (3.18)) Other spectral scalings of interest are the sub-exponential scaling: e −aN |γ j | 1 N , j = 1, . . . , N for any a > 0, and super-exponential scaling: e −aN |γ j |, j = 1, . . . , N for any a > 0. The latter corresponds to the case of an 'ideal' soliton gas consisting of noninteracting solitons, and the former to the special kind of soliton gas termed soliton condensate. These scalings were introduced in [38] in the context of the focusing NLS soliton gas and will be considered in Section 3.3.
Consider a partial sum K M = 1 2π M j=1 k j over the spectral lattice {η 1 , η 2 , . . . , η N }, where 1 ≤ M ≤ N . Given that 2πk j is the spatial density of waves (which we treat as quasiparticles) associated with the j-th spectral band, the quantity K M has the meaning of the integrated density of states [109], [99]. Invoking the scaling (3.34) and passing to the continuum limit N → ∞, η M → η, we obtain Then the DOS is given by Similarly, for the temporal counterpart of the integrated DOS-the spectral flux-we have (3.42) so that the spectral flux density in a soliton gas is given by

Equation of state and spectral kinetic equation
Eliminating σ from the nonlinear dispersion relations (3.37), (3.38) we obtain for s(η) = v(η)/f (η) which is exactly the equation of state (2.11) obtained in Section 2 under the collision rate assumption (2.12). Hence this assumption is now justified. As suggested by the phenomenological derivation in Section 2 the function s(η) in (3.44) has the meaning of the effective soliton velocity (i.e. the transport velocity) in a soliton gas. We will now show how this interpretation is justified within the thermodynamic limit framework. To this end we consider non-equilibrium soliton gas with f (η) ≡ f (η, x, t), s(η) ≡ s(η, x, t) and derive the evolution equation for the DOS. We go back to the original, discrete wavenumber and frequency components k j (λ), ω j (λ) of the finite-gap potential, defined in terms of the fixed branch points λ of the Riemann surface R of (3.22). Let us now consider a slowly modulated finite-gap potential with λ = λ(x, t). The modulation system describing the evolution of 2N +1 parameters λ 1 (x, t), λ 2 (x, t), . . . λ 2N +1 (x, t) has been derived in [35]. This system admits an infinite number of hyperbolic conservation laws, that include a finite subset of N wave conservation laws (3.15), which can be manipulated into the equivalent system In particular, for n = 1, 2 the Kruskal integrals coincide with the respective statistical moments of u-the observables of the KdV soliton gas field-and have the form in full agreement with the result (2.41) obtained by the heuristic 'windowing' procedure in Section 2.3. Concluding this section we note that the presented construction of soliton gas it was assumed that soliton propagate on a fixed (zero) background, which was achieved by fixing the endpoint λ 2N +1 = 0 of the spectrum (3.21). A generalisation to a slowly varying background is possible following the modulation construction of solitonic dispersive hydrodynamics in [17], [110]. Such a generalisation could provide interesting insights into new soliton gas phenomena.

Poisson distribution for position phases
Having defined the thermodynamic limit for the spectrum of finite-gap potentials F N (θ), θ ∈ T N , we now need to determine what happens with the phases θ j = k j x − ω j t + θ (0) j , j = 1, . . . , N in this limit.
As discussed in Section 3.1.2 we adopt the Random Phase Approximation in which the initial phases θ (0) j are assumed to be N independent random values uniformly distributed on [0, 2π). Under the random phase approximation the finite-gap potential u(x, t) transforms into an ergodic random process [99], both as function of x and t. We now fix t = t * and represent the cyclic phases θ j ( mod 2π) in the form θ j = k j (x − x j ), where the position phases x j = x 0 j + c j t * with c j = ω j /k j being the phase velocities. The initial position phases x 0 j = −θ (0) j /k j are independent random values, each distributed uniformly on the respective period interval [0, 2π/k j ). The deterministic shifts c j t * can be absorbed into the random initial positions so that we shall replace x (0) j → x j . The thermodynamic spectral scaling (3.30) implies that in the limit N → ∞ all wavenumbers vanish, k j → 0, i.e. all spatial periods of the finite-gap solution become infinite. We now demonstrate that, upon the thermodynamic limit (3.16) the uniform distribution of θ (0) on T N transforms into the Poisson distribution for the position phases x j ∈ R.
Since ξ j are independent random variables the probability-generating function for ξ (N ) is where Hence the probability of having n points x i in the interval ∆ is given by the Poisson distribution with parameter β, the total integrated DOS in the soliton gas. The position phases x j coincide with the soliton centres in the stochastic soliton lattice model (2.8) of a rarefied gas. The random process (2.8) can then be associated with a compound Poisson process [111], a stochastic process with jumps. The jumps, associated with solitons, are distributed on the line randomly according to a Poisson distribution with probability (3.55) and the size of the jumps (the soliton amplitudes a i = 2η 2 i ) is also random, where the spectral parameter η has the probability distribution ϕ(η) = β −1 f (η).

Focusing NLS equation: soliton and breather gases
The spectral theory of soliton gas for the KdV equation was generalised in [38] to the case of the focusing NLS equation (2.30). As was mentioned in Section 2.2.2 the essential difference between the KdV and focusing NLS equations at the level of the IST is that the Lax (Zakharov-Shabat) operator (2.31) for the focusing NLS has a complex-valued spectrum satisfying the Schwartz symmetry property. This results in the availability of the qualitatively new, compared to KdV, families of localised wave solutions and, consequently, to new types of soliton gases. Along with the individual fundamental solitons propagating on zero background the focusing NLS equation supports a family of the so-called bound state N -soliton solutions, the complexes of strongly interacting, co-propagating (i.e. having the same, possibly zero, velocity) solitons. Another important family of localised solitonic structures supported by the focusing NLS equation are breathers, that can be viewed as solitons on finite background. All these localised solutions can serve as quasiparticles in the respective soliton or breather gases, which we consider in this section.

Solitons and breathers
The fundamental soliton solution of the focusing NLS equation is given by the formula (2.32) (see Fig. 8(a) for the profile of |ψ(x)| in such soliton). The inherent feature of the fundamental solution is that it can only exist on the zero background, which makes it very different from KdV solitons and solitons in the defocusing NLS/resonant NLS dispersive hydrodynamics. The IST spectral portrait of the fundamental soliton consists of two complex conjugate points λ ∈ {λ 1 ,λ 1 } (see the inset in Fig. 8 (a)). Apart from the fundamental solitons (2.32) the focusing NLS equation supports a more general family of localised solutions, called breathers, that can be viewed as solitons on finite background. The Lax spectrum of a breather consists of a vertical band γ 0 = [−iq, iq], q > 0 complemented by two complex conjugate solitonic discrete points λ 1 ,λ 1 . In physical space such a breather, sometimes called the Tajiri-Watanabe (TW) breather, represents a localised nonlinear wavepacket propagating on a uniform nonzero background described at x → ±∞ by the plane wave solution ψ = qe i2q 2 t (see Fig. 8 (b) for the typical behaviour of |ψ| in the TW breather along with the spectral portrait in the inset). The analytical expression for the TW breather solution is available elsewhere (see e.g. [112], [113], [114]), here we only present its group (envelope) and phase (carrier wave) velocities: where λ 1 is the IST solitonic spectral point characterising the TW breather and R 0 (λ) = λ 2 + q 2 . The transition from the TW breather solution to the fundamental soliton (2.32) is achieved by vanishing the background, q → 0.  The TW breather has three special reductions, spectrally realised by placing the soliton discrete eigenvalues λ 1 ,λ 1 on the imaginary axis. These reductions are the Akhmediev breather (AB), the Kuznetsov-Ma (KM) breather and the Peregrine soliton (PS) and they, particularly the PS, are often considered as 'analytical prototypes' for rogue, or anomalous waves (see e.g. [116], [117], [118], [119] and references therein). All these three special families of breathers possess specific localisation properties: the AB is localised in time and is periodic in space, the KM breather is localised in space and is periodic in time, and finally, the PS exhibits both temporal and spatial localisation. Their IST spectral characterisation is as follows (see e.g. [115]).
Let the spectral band corresponding to the plane wave background of the breather be γ 0 = [−iq; iq] for some q > 0, and the solitonic discrete spectrum points λ 1 = ip,λ 1 = −ip, p > 0. Then p < q corresponds to AB, p > q to KM breather and p = q to PS. The behaviours of |ψ(x, t)| in the AB, KM and PS breathers, along with their spectral portraits, are displayed in Fig. 9.

Finite-gap solutions
Similar to the KdV equation, the focusing NLS equation (2.30) supports finite-gap solutions which transform into the solitonic solutions when the spectral bands collapse. An n-gap solution ψ = ψ n (x, t) of the focusing NLS equation (2.30) is defined by a fixed set of 2n + 2 endpoints {λ j ,λ j , j = 0, 1, 2, . . . , n} of spectral bands γ j , j = 0, 1, . . . , n + 1, and depends on n real phases θ(x, t) = kx − ωt + θ 0 with the initial phase vector θ 0 ∈ T n , so that |ψ n (x, t)| = F n (θ(x, t)), where F n is a multi-phase (quasiperiodic) function in both x and t, that can be expressed in terms of the Riemann theta-functions [120]. The n-component wavenumber k and the frequency ω vectors depend on the endpoints {λ j , j = 0, 1, 2, . . . , n} of the spectral bands, which define a hyperelliptic Riemann surface R of genus n given by (cf. equation (3.22) for the KdV equation): (3.57) z ∈ C being a complex spectral parameter in the Zakharov Shabat scattering problem; R n (z) ∼ z n+1 as z → ∞. The branch cuts of R(z) are made along spectral bands which will be specified below. An example of the behaviour of |ψ(x, t)| in a genus 4 solution is shown in Fig. 10. There are two qualitatively different types of spectral bands characterising finite-gap NLS solutions: (i) the bands that cross the real axis in the complex spectral plane z ∈ C and connect the complex conjugate spectral points λ j andλ j (these are often called the Stokes bands, see e.g. [6]); (ii) the Schwartz symmetric bands that do not cross the real axis (we shall call them the solitonic bands). By manipulating the spectral bands one can achieve various wave configurations for the field ψ n (x, t). E.g. collapsing a pair of solitonic bands into two complex conjugate double points of the spectrum gives rise to a localised mode in the solution, a soliton on the finite-gap background, a generalised breather (see e.g. [121]). Collapsing a Stokes band into a double point on the real axis implies vanishing of one of the background plane wave modes. By collapsing all Stokes and solitonic bands into double points of the spectrum a finite-gap solution is transformed into a multi-soliton solution. Keeping one of the Stokes bands finite, but collapsing all other bands implies a multi-breather solution.
The finite-gap analogues of the TW breather solutions and their reductions (AB, KM and PS) have genus n = 2 and contain one Stokes band in their Lax spectrum. We shall consider a generalisation of these solutions for an arbitrary genus and, having in mind the construction of soliton and breather gas, make two simplifying assumptions: (i) assume an even genus n = 2N of the spectral Riemann surface; (ii) assume that the Lax spectrum of the finite-gap potential is located on a Schwarz symmetric, simply connected 1D curve Γ ⊂ C, see Fig. 11.
Thus we have one Stokes band γ 0 between λ 0 and λ 0 and 2N solitonic bands: γ j , j = 1, . . . , N between λ 2j−1 and λ 2j , and their Schwarz symmetric counterparts γ −j , j = 1, . . . , N between λ 2j−1 and λ 2j in the lower half-plane. We note that the described spectral geometry also covers the case of an odd genus (achieved by collapsing the Stokes band) and a 'bound state' configuration, when Γ lies on a vertical line so that all the solitons or breathers corresponding to the collapsed bands have the same velocity. Due to the symmetry of the curve Γ it is sufficient to consider only the upper complex half-plane (C + ) part of it, which we denote Γ + (so that Γ + = Γ ∩ C + ) While appearing quite restrictive, the described 1D spectral geometry provides a major insight into the properties of breather and soliton gases and admits a straightforward generalisation in the more physically realistic case, where the spectral bands γ j are located in some (Schwarz symmetric) 2D region Λ ⊂ C. It also admits generalisation to the case of more than one Stokes bands. 3 e. if z ∈ γ j is the point of the spectrum in lf plane C + of the complex plane then the elongs to the spectrum. olution ψ = ψ n (x, t) of (2) is defined by a + 1 spectral bands γ j , j = 1, . . . , n + 1, on n real phases η(x, t) = kx + ωt + θ 0 al phase vector θ 0 ∈ T n , so that |ψ n (x, t)| = here F n , expressed in terms of the Riemann ns, is a multi-phase (quasiperiodic) funcand t [32]. The n-component wavenumber quency ω vectors depend on the endpoints 2, . . . , n} of the spectral bands, which define c Riemann surface R of genus n given by: (3) complex spectral parameter in the ZS scatm. The branchcuts of R(z) will be specified tailed description of finite-gap solutions of NLS equation can be found in [32][33][34]. We t the finite-gap theory provides a natural r the construction of random solutions to tion by assuming a uniform distribution of ase vector θ (0) ∈ T n [21,35].
assume that all spectral bands lie on a fiartz-symmetric curve Γ (see Fig. 2); later n will be removed to allow the bands to a 2D compact subset of C. The curve Γ e to be a connected curve. We also assume s n is even, n = 2N , N ∈ N; the transid genus case will be described below. We e branch points according to the following the spectral bands. The spectral band γ j , , is defined as the segment of Γ in the upe C + connecting the branch points α 2j and ectral band γ −j is Schwartz symmetric to ects, naturally, the branch pointsᾱ 2j and for convenience, we will denote as α −2j and ctively. Finally, there is an exceptional band ing the real axis and connecting the branch α −1 =ᾱ 1 . An odd genus case can then forsidered by collapsing the exceptional band on the real axis. Note that generally one e than one exceptional band, We notice that the contours A j , |j| = 1, · · · , N , and −B j , B −j , j = 1, . . . , N , form a homology basis of the Riemann surface R with the branchcuts along γ j , |j| = 0, 1, . . . , N .
Collapsing a single pair of bands γ ±j , j > 0 into a pair of double points (α ±(2j+1) → α ±2j ) implies the appearance of a soliton on a (n − 1)-gap solution background. Shrinking all bands γ j , |j| = 0, 1, . . . , N to points corresponds to the transformation of a 2N -gap solution into a 2N -soliton solution with b j corresponding to soliton amplitudes and a j -to their velocities at t → ∞ (cf. Eq. (3) for the definitions of a j , b j ). If all a j are equal (without loss of generality one can assume that a j = 0 ∀j, i.e. each γ j is an interval on the imaginary axis) the limiting N -soliton solution is a bound state [25,36], in which the solitons do not separate as t → ∞. We note that in the soliton limit the special "reference" band γ 0 collapses to the origin (which can be associates with a zero-amplitude soliton) and so does not contribute to the limiting Nsoliton solution. If the exceptional band γ 0 remains finite, then collapsing all the other bands into double points correspond to a (multi)-breather limit of a 2N -gap solution with the finite reference band γ 0 being "responsible" for the background (indeed, the genus zero solution with a single band γ 0 = [−iq, iq], q > 0, in the spectrum is the so-called plane wave, or the condensate given by ψ = qe i2q 2 t ). The generic "elementary" breather corresponding to a degenerate genus 2 solution is the so-called Tajiri-Watanabe (TW) breather [27] with the typical behavior of the amplitude |ψ T W (x, t)| shown in Fig. 3. The spectral portrait of the TW breather (shown in the inset of Fig. 3) consists of the vertical branchcut connecting the points ±i, and two double points: α 2 = i cos(µ + iΨ)

Wavenumbers and Frequencies
The wavenumber and frequency vectors, k and ω respectively, associated with a given finitegap solution ψ 2N (x, t) of the focusing NLS equation, are not uniquely defined since any linear combination of the wavenumber (frequency) vector components with integer coefficients is also a wavenumber (frequency). The construction of the thermodynamic limit outlined in Section 3.1.2 and realised in Section 3.2 for the KdV equation relies on the availability of the wavenumbers and frequencies that vanish in the limit when the relevant bands collapse into double points. Such wavenumbers and frequencies automatically appear in the finite-gap KdV theory [35]. The situation with the focusing NLS equation is different since the sets of the wavenumbers and frequencies that arise 'naturally' in the focusing NLS finite-gap theory (see e.g. [122], [123], [124]) do not necessarily possess the requisite property.
The special, fundamental wavenumber and frequency vectors k = (k 1 , . . . , k N ,k 1 , . . . ,k N ) and ω = (ω 1 , . . . , ω N ,ω 1 , . . . ,ω N ), (3.58) that possess the properties necessary for the application of the thermodynamic limit have been identified in [38]. Following the KdV construction in Section 3.2 we introduce two new quantities characterising the finite-gap potentials instead of the endpoints of spectral bands (we consider only the upper half-plane with posistive indices of spectral bands γ j and gaps c j for convenience): (3.59) where j = 1, . . . , N . We shall call the point η j the centre of the j-th band γ j and 2|δ j | the j-th band width. Also for the Stokes band we have δ 0 = iImλ 0 . Note that the notations and the band numeration we are using here are slightly different from those used in [38]. It then follows that the fundamental wavenumbers and frequencies defined by (A.1) and (A.2) have drastically different asymptotic properties in the soliton/breather limit, when one of the solitonic spectral bands collapses into a double point, λ 2i−1 , λ 2j → η j . Namely, In particular, for N = 1 (genus 2), the limit (3.60) (k 1 → 0, ω 1 → 0) with non-zero Stokes band γ 0 (i.e. δ 0 = 0) corresponds to the TW breather limit of the two-phase nonlinear wave solution. The remaining wavenumber and frequencyk 1 = O(1),ω 1 = O(1) correspond to the "carrier" wave of the TW breather (see Fig. 8). If we further tend δ 0 → 0, thenk 1 → 0,ω → 0 and the TW breather transforms into a fundamental soliton. Motivated by these properties for N = 1 we shall call the components k j , ω j of the wavenumber and the frequency vectors k and ω the solitonic components and the componentsk j ,ω j -the carrier components. The solitonic components k j , ω j are the focusing NLS counterparts of the conventional KdV wavenumbers and frequencies specified by (3.26) while the carrier componentsk j ,ω j do not have analogues in the KdV theory. Generally the limit (3.60) for finite N corresponds to N -breather solution if δ 0 = 0 and to N -soliton solution if δ 0 = 0.

Nonlinear dispersion relations
It was shown in [38] that the solitonic components k m , ω m satisfy the following nonlinear dispersion relations (again, considering only the upper half plane): Re λ k + κ j,2 ), where and κ i,j are the coefficients of the normalized holomorphic differentials w j defined by: The contours α j , β j on the Riemann surface R are defined as follows (see Fig. 11): the α jcontours surround the bands γ j clockwise, and the β j -contours start at the Stokes band γ 0 , go to the band γ j on the upper sheet clockwise and return to the band γ 0 on the lower sheet of the Riemann surface. The derivation of (3.61) is outlined in Appendix.
Relations (3.61) are the analogs of the nonlinear dispersion relations (3.26) for the finite-gap potentials of the KdV equation. Similar relations are available for the carrier components [38] but we do not present them here as they play a secondary role in the thermodynamic limit construction.
As for the scaling of the band widths, we consider the following options (cf. eq. (3.30) for KdV): (i) exponential spectral scaling: the band widths |δ j | are exponentially small in N : where τ (λ) is a smooth positive function on Γ + having the meaning of the normalised logarithmic band width (τ (η j ) ∼ − ln |δ j |/N ).
(ii) sub-exponential spectral scaling: for any a > 0 It is clear that in this limit τ (λ) → 0.
We also note that the exponential and sub-exponential spectral scalings have the 'thermodynamic' property in the sense that they preserve finiteness of the total density of waves K N = N j=1 k j in the limit N → ∞ so that lim N →∞ K N = β, where 0 < β < ∞. Note that for the super-exponential scaling β → 0.

Nonlinear dispersion relations and kinetic equation
We first present the results for the general case of breather gas by evaluating the thermodynamic limit of the nonlinear dispersion relations (3.61) for the exponential spectral scaling (3.64). Without much loss of generality we assume that the Stokes band lies on the imaginary axis, γ 0 = [−iq, iq], q > 0 (so that δ 0 = iq) and that η i ∈ arc(a, b) ⊂ Γ + , and introduce the simultaneous scaling for the solitonic wavenumbers and frequencies (cf. (3.34)) so that κ j = κ(η j ) and ν j = ν(η j ), where the functions κ(λ) ≥ 0 and ν(λ) are smooth interpolations of κ j , ν j . Similar to the KdV case, the scaling (3.67) provides a balance of terms in relations (3.61) for N 1. The resulting nonlinear dispersion relations for breather gas have the form (we refer the reader to [38] for details of the derivation): where R 0 (z) = z 2 + q 2 (with the branch cut [−iq, iq], and the branch of the radical defined by R 0 (z) → z as z → ∞), and with a slight abuse of notation, we denoted b a . . . |dµ| ≡ Γ + . . . |dµ|. Further, Equations (3.68), (3.69) are the focusing NLS counterparts of the nonlinear dispersion relations (3.37), (3.38) for the KdV soliton gas with f (λ) and v(λ) having the meanings of the DOS and the spectral flux density respectively, and the function σ(λ) encoding the Zakharov-Shabat spectrum of the finite gap potentials in the thermodynamic limit. The soliton gas limit in (3.68), (3.69) is achieved by vanishing the Stokes band, q → 0, resulting in is the effective velocity of the tracer soliton (breather) in the gas, and s 0 (λ), G(λ, µ) are defined as follows. For breather gas: and for soliton gas: One can see that s 0 (λ) in (3.75) coincides with the group velocity c g (λ) of an isolated TW breather (3.56) and s 0 (λ) in (3.76) coincides with the group velocity of the fundamental soliton of the focusing NLS equation. Furthermore, the integral kernel G(λ, µ) in (3.76) coincides with the absolute value of the position shift (2.33) in the two-soliton collision for the focusing NLS equation. As one may expect, the expression for G(λ, µ) in eq. (3.75) is identified with the position shift in two-breather collisions [125], [114] (see the proof in [40]). Thus the outlined derivation provides the justification of the collision rate assumption (2.12) for dense soliton and breather gases of the focusing NLS equation.
The integral equations (3.72), (3.73) have been recently studied in a rather general spectral geometry [44]. The existence and uniqueness of their solutions was investigated and the nonnegativity of the solutions for the DOS f (λ) was proved.

Spectral kinetic equation for nonequilibrium gas
The derivation of the transport equation for non-equilibrium soliton and breather gases for the focusing NLS equation follows the general modulation construction outlined in Section 3.1 and realised for the KdV equation in Section 3.2.2, albeit with some important differences.

Generalisation to 2D Case
The above derivation of the spectral kinetic equation for the focusing NLS soliton and breather gases involves the basic assumption that the DOS f (λ) is supported on a 1D symmetric curve Γ in the complex plane. This restriction can be readily removed as we show below.
In the case when the shrinking bands γ j , j > 0 fill a compact 2D region Λ + of the upper complex half-plane, the counterpart of the exponential scaling (3.64) is where τ (λ) is a positive smooth function on Λ + . The scaling of the gaps remains O(1/N ), where by the gap width we understand the closest distance between the bands. In this case ϕ(λ) > 0 is the 2D density of bands (and we also distinguish the cases of exponential, sub-exponential and super-exponential scalings of bands, similarly to the 1D case). We now assume one of the 2D spectral thermodynamic spectral scalings (exponential, subexponential, and super-exponential) when the shrinking bands γ j fill a 2D region Λ of the complex plane, see Section III B. For the wave numbers and frequencies instead of (3.67) we introduce where κ j = κ(η j ) and ν j = ν(η j ), and the interpolating functions κ(λ) ≥ 0, ν(λ) are assumed to be smooth on Λ + . Then the 2D thermodynamic limit of the nonlinear dispersion relations (3.61) leads to the same integral equations (3.68)-(3.74) but with the line integration along Γ + replaced by the integration over a 2D compact domain Λ + : where µ = ξ + iζ.
For convenience of the exposition we shall be using the notation Γ + . . . |dµ| in both 1D and 2D cases keeping in mind that in the 2D case the meaning of the integral is given by Eq. (3.85).

Rarefied soliton gas and soliton condensate
The dispersion relations and the equations of state for soliton and breather gases were derived under the assumption of the exponential spectral scaling (3.64) (or (3.83) in 2D case). We now consider the two other scalings of interest: the superexponential scaling (3.66) and subexponential scaling (3.65). We will restrict our exposition to soliton gas but will indicate when the results can be extended to breather gas.
It is convenient to characterise the spectral scalings and the corresponding soliton gases in terms of the function σ(λ) parametrising the nonlinear dispersion relations (3.68) and (3.72). From Eq. (3.68) we have: (there is a similar expression for breather gas which we do not present here). For the exponential scaling σ(λ) = O(1), while the limiting cases σ → ∞ and σ → 0 correspond to the super-and sub-exponential spectral scalings respectively.

Rarefied soliton gas
Rarefied soliton gas represents an infinite random ensemble of weakly interacting solitons characterized by a small DOS, f 1, and therefore, σ 1 by (3.86). We shall refer to the limit f → 0, σ → ∞, f σ = O(1) as the ideal gas limit as it corresponds to the gas of non-interacting breathers (solitons). Spectrally this limit corresponds to the super-exponential spectral scaling (3.66).
For a rarefied gas the interaction (integral) term in the equation of state (3.74) is subdominant so the leading order term s(λ) = s 0 (λ) describes the group velocity distribution in an ideal breather (soliton) gas. Then the first correction to the ideal gas velocity s 0 (λ) is readily computed to give: (3.87) Eq. (3.87) represents the focusing NLS counterpart of the equation (2.9) for the effective soliton velocity in a rarefied KdV soliton gas introduced by Zakharov [28]. All the above results apply to breather gas as well.

Soliton condensate
The inequality σ(λ) ≥ 0 in (3.86) imposes a fundamental constraint on the DOS f (λ). As follows from the discussion in Sec. 3.3.2, the critical value σ(η) = 0 corresponds to the subexponential spectral scaling (3.65). One can see from the nonlinear dispersion relations (3.72), (3.73) that in this case the gas properties are fully determined by the interaction (integral) terms, while the information about the individual quasi-particles (described by the non-integral, secular, terms) is completely lost. By analogy with Bose-Einstein condensation we shall call the soliton gas at σ = 0 the soliton condensate. From (3.72), (3.73) we obtain the dispersion relations for the focusing NLS soliton condensate If the spectral support Γ of the DOS belongs to a vertical line, Γ + ⊂ iR + , the second equation (3.88) implies v(λ) = 0 and, hence, s(λ) = 0. We shall generally call such a non-propagating gas a bound state soliton gas, and it will be the bound state soliton condensate in the present context. Let Γ = [−iq, iq] for some q > 0 and assume that f (μ) = −f (µ) for all µ ∈ [0, iq] (an odd extension of f (µ) onto [−iq, 0]). Then, as was shown in [38] the first equation (3.88) for the DOS reduces to One can observe that the DOS (3.90) coincides with the appropriately normalised semi-classical distribution of discrete Zakharov-Shabat spectrum for the potential ψ ∈ R in the form of a broad rectangular barrier. Eq. (3.90) (up to a norming factor) is obtained as a derivative of the corresponding Weyl's law following from the Bohr-Sommerfeld quantisation rule for the Zakharov-Shabat operator [84]. A numerical realisation of the bound state soliton condensate is shown in Fig. 12. It is achieved by building a dense (as dense as possible) ensemble of 100 strongly interacting solitons with the discrete spectrum eigenvalues distributed according to the Bohr-Sommerfeld rule and random phases of the norming constants (see [26] for the details of the effective numerical implementation of dense N -soliton ensembles with N large). It was shown in [25] that the bound state soliton condensate dynamics reproduces with remarkable accuracy the statistical characteristics of stationary integrable turbulence generated at the nonlinear stage of the development of the noise-induced modulational instability [8], [10], a fundamental and ubiquitous physical phenomenon in focusing nonlinear media [128], [129], [6]. Furthermore, it was demonstrated in [27] that the nonlinear wave field in the bound state soliton condensate exhibits a spontaneous emergence of rogue waves, offering a new perspective on the dynamical and statistical mechanisms underlying this phenomenon observed in water waves, nonlinear optics and matter waves (see [118] and references therein).
The bound state soliton condensate DOS (3.90) is a particular solution of the dispersion relations (3.88). Another explicit solution of the soliton condensate dispersion relations (3.88) is available in the special case when the spectral support Γ of the DOS represents a circle or a circular arc in the complex plane. This solution with nonzero spectral flux density v(λ) obtained in [38] corresponds to a dynamic (non-bound state) soliton condensate.
In conclusion we note that the soliton condensate regime is also available for the KdV soliton gas. Indeed, substitution of σ = 0 in the KdV soliton gas dispersion relations (3.37), (3.38) yields

Special breather gases
In Section 3.3.1 we presented some properties of breather solutions to the focusing NLS equation and indicated three special cases: the Akhmediev breather (AB), the Kuznetsov-Ma (KM) breather and the Peregrine soliton (PS), which are often considered as mathematical models for rogue waves in fluids and nonlinear optical media. These breathers can be viewed as quasiparticles of the special, 'rogue wave breather gases'. Assuming the Stokes band γ 0 = [−iq, iq], q > 0 the AB and KM breather gases are characterised by the DOS f (λ) that is supported on the appropriate intervals of the imaginary axis, λ ∈ Γ + : Γ + = [0, ip], p < q for AB gas and Γ + = [ip, ib], b > p > q for KM gas. For PS gas the DOS is given by f (λ) = wδ(λ − iq), where w > 0 is the PS gas density. The three types of the rogue wave breather gases were realised numerically in [40] by building a random ensemble of N ∼ 50 breathers via the Darboux transform recursive scheme in high precision arithmetic, see Fig. 13. The interaction properties of the constructed breather gases were investigated by propagating through them a trial generic breather (TW) and comparing the effective (mean) propagation velocity with the predictions following from the equation of state (3.74). One of these predictions is that the propagation through the PS breather gas is ballistic, i.e. the PS breather gas does not alter the velocity s 0 (λ) (3.75) of the trial breather due to vanishing of the logarithmic interaction kernel G(λ, µ) (3.75) evaluated at µ = iq [38]. Other predictions, confirmed by detailed quantitative comparisons with numerical simulations in [40], are that the effective velocity of propagation s(λ) of a trial TW breather through a rogue wave gas is greater than s 0 for the KM gas and less than s 0 for the AB gas.
(2.10), (2.13) for unidirectional/bidirectional isotropic soliton gas: where f (λ, x, t), s(λ, x, t), s 0 (λ) and G(λ, µ) are real functions and the parameters λ, µ can be real or complex. For unidirectional gas the spectral support Γ of the DOS f (λ, x, t) is an interval in R. For bidirectional isotropic gas Γ is either a symmetric interval [−a, a] ⊂ R (cf. (2.26)) or a Schwarz symmetric compact domain in C (1D or 2D, see (3.85) for the formal correspondence between the two cases). Consider a soliton gas composed of a finite (M ∈ N) number of distinct components, termed monochromatic, or cold, components modelled by the DOS in the form of a linear combination of the Dirac delta-functions centred at distinct spectral points ζ j ∈ Γ, where w j (x, t) > 0 are the components' weights, and {ζ j } M j=1 ⊂ Γ, (ζ j = ζ k ⇐⇒ j = k). Substitution of (4.2) into the kinetic equation (4.1) reduces it to a system of hydrodynamic conservation laws where the component densities w i (x, t) and the transport velocities s j (x, t) ≡ s(ζ j , x, t) are related algebraically: Here we used the notation We note that, although we derived the hydrodynamic reduction (4.3), (4.4) in the context of unidirectional/isotropic soliton gas, it can be readily generalised to the bidirectional anisotropic case, see [37]. In the latter case the ansatz (4.2) is introduced separately for the slow and fast soliton branches so that the pair of distributions For M = 2 system (4.4) can be solved to give explicit expressions for s 1,2 (w 1 , w 2 ): .
An important remark is in order on the meaning of the delta-function ansatz (4.2) for the DOS f (λ). As a matter of fact, the representation (4.2) is a mathematical idealisation, which has a formal sense in the context of the integral equation of state (4.1), but cannot be applied to the original dispersion relations where it appears in both the integral and the secular terms (cf. (3.37), (3.38) for the KdV equation). In a physically realistic description the delta-functions in (4.2) should be replaced by some narrow distributions around the spectral points ζ j , i.e. we first take the thermodynamic limit N → ∞ and then allow the distributions to become sharply peaked. As a result, the non-condensate condition σ > 0 in the dispersion relation for the DOS (equation (3.37) for KdV or equation (3.68) for focusing NLS) would impose a constraint Γ G(λ, µ)f (λ)dµ < 1 on f (λ) which, among other things, implies that the denominators in (4.7) must be positive. It is interesting to note that a variant of hydrodynamic system (4.3), (4.7) has been recently derived in [130] in the context of TT -deformed conformal field theories out of equilibrium.
We now assume that the integral kernel G(λ, µ) in the equation of state (4.1) can be represented as G(λ, µ) = a(µ) (λ, µ), where (λ, µ) = (µ, λ) > 0, λ = µ (4.8) for some non-singular real function a(µ) = 0. This structure of the soliton interaction kernel is typical for integrable dispersive hydrodynamics and is indeed shared by the equations of state of the KdV, defocusing and focusing NLS unidirectional/isotropic soliton gases (see (3.44), (2.26), (2.34)). E.g. for KdV we have a(µ) = µ, (λ, µ) = 1 λµ ln λ+µ λ−µ , see (2.11). We then have G jk = a k jk , jk = kj , j = k, (4.9) where a k = a(ζ k ), jk = (ζ j , ζ k ). Solving equations (4.7) for w 1,2 gives (4.10) Substituting (4.10) in (4.3) for M = 2 we obtain a diagonal system [36]: with s 1,2 being the Riemann invariants. One can see that system (4.11) is linearly degenerate because its characteristic velocities do not depend on the corresponding Riemann invariants. System (4.11) represents the diagonal form of the so-called Chaplygin gas equations, the system of isentropic gas dynamics with the equation of state p = −A/ρ, where p is the pressure, ρ is the gas density and A > 0 is a constant. The Chaplygin gas equations occur in certain theories of cosmology (see e.g. [131]) and is also equivalent to the 1D Born-Infeld equation [132], [33] arising in nonlinear electromagnetic field theory. As any two-component quasilinear hyperbolic system, system (4.11) is integrable (linearisable) via the classical hodograph transform. However, for any M ≥ 3 integrability of the original system (4.3), (4.4) is no longer obvious. As a matter of fact, a M -component hydrodynamic type system is generally not integrable for M ≥ 3. Below we present the theorem, proved in [41], which summarises the mathematical properties of the hydrodynamic reductions (4.3), (4.4) subject to the structure (4.9) of the interaction matrix G jk .
M -component reductions (4.3), (4.4) of the generalized kinetic equation (4.1) with the interaction kernel satisfying (4.8) are hyperbolic, linearly degenerate integrable hydrodynamic type systems for any M ∈ N.
Specifically, it was shown in [41] that the following fundamental properties are satisfied for the system (4.3), (4.4): (i) Riemann invariants.
The characteristic velocities V j of the diagonal system (4.13) satisfy so that system (4.3), (4.4) is linearly degenerate in the Lax sense [133].
In addition to (4.15) the characteristic velocities V i (r) satisfy the overdetermined system for each three distinct characteristic velocities (∂ k ≡ ∂/∂r k ). Diagonal hydrodynamic type systems satisfying (4.16) are called semi-Hamiltonian [66].
A semi-Hamiltonian hydrodynamic type system can be locally integrated via the generalised hodograph transform due to Tsarev [66], [134]. Its general local solution for ∂ x r i = 0, i = 1, . . . , M is given by the formula where functions W i (r) are found from the linear system of PDEs: Properties of linearly degenerate semi-Hamiltonian (integrable) hydrodynamic type systems have been studied in [135], [136]. Using the results of [135], [136] it was shown in [41] that the general local solution of the hydrodynamic reductions where φ 1 (ξ), . . . , φ M (ξ) are arbitrary functions. Some interesting particular solutions (selfsimilar, quasiperiodic) following from (4.19) can be found in [41]. The Hamiltonian structure of the hydrodynamic reductions (4.3), (4.4) was investigated in [43]. The generalised hydrodynamic reductions arising if one allows the discrete spectral points ζ j in (4.2) be functions of (x, t) were studied in [42].
In conclusion of this section we note that, somewhat counter-intuitevely, the above results on the hyperbolic structure and integrability of hydrodynamic reductions of the kinetic equation for soliton gas apply to the multi-component soliton gas of the focusing NLS equation, whose Whitham modulation system is known to be elliptic for a generic set of modulation parameters and for any genus [137], [122]. This apparent paradox is resolved by noticing that the Whitham system for the focusing NLS equation exhibits real eigenvalues (characteristic speeds) in the soliton limit. E.g. for the genus 1 case, the two pairs of complex conjugate eigenvalues of the modulation matrix collapse into a single, quadruply degenerate real eigenvalue, corresponding to the velocity of the fundamental soliton, see e.g. [63]. 43 Having obtained the hydrodynamic description of multi-component soliton gas it is natural to consider the Riemann problem that plays the fundamental role in classical gas and fluid dynamics [138]. The classical Riemann problem consists in finding solution to a system of hyperbolic conservation laws subject to piecewise constant initial conditions exhibiting discontinuity at x = 0. The distribution solution of the Riemann problem generally consists of a combination of constant states, simple (rarefaction) waves and strong discontinuities (shocks or contact discontinuities) [133]. The discontinuous solutions satisfy the Rankine-Hugoniot jump conditions. If a hyperbolic system is linearly degenerate it does not support simple waves and classical shocks, and the solution of Riemann problem contains only constant states and contact discontinuities [139]. Here, following [75], [37], we present such solutions for different types of soliton gases.

General solution
We consider the hyperbolic system (4.3), subject to the initial data corresponding to two n-component unidirectional or bidirectional soliton gases prepared in the respective uniform states w L ∈ R n and w R ∈ R n , that are initially separated: We emphasize here that the Riemann problem (5.1), (5.2) is formulated for the kinetic equation (via the substitution (4.2)) but not for the original dispersive hydrodynamic system (1.1) or (1.2). For this reason we shall call it the spectral Riemann problem as it essentially describes the spatiotemporal evolution of the spectral DOS of the soliton gas. The spectral Riemann problem for the KdV and focusing NLS two-component soliton gases (n = 2) was investigated in [36], [75], [38] and for rather general isotropic and anisotropic dispersive hydrodynamic soliton gases with an arbitrary number of spectral components in [37]. The Riemann problem of this type has also been considered in the context of generalised hydrodynamics [55], [56], [140], [141], [142], [143]. The spectral Riemann problem (5.1), (5.2) corresponds to the dispersive hydrodynamic Riemann problem that can be viewed as a soliton gas shock tube problem, an analog of the shock tube problem of classical gas dynamics [138], [33]. The shock tube problem represents a good benchmark for the spectral kinetic theory of soliton gas where one can investigate both overtaking and head-on collisions by choosing the appropriate components of solitonic spectra and then comparing the predictions of the kinetic theory with direct numerical simulations of dispersive hydrodynamics.
Following [75], [37] we consider the spectral Riemann problem (analytically) and the corresponding soliton gas shock tube problem (numerically) for n-component unidirectional and bidirectional soliton gases. Due to the scale invariance of the problem (the kinetic equation (5.1) and the initial condition (5.2) are both invariant with respect to the scaling transformation x → Cx, t → Ct), the spectral solution is a self-similar distribution w(x/t). Because of the linear degeneracy (4.15) of the hydrodynamic system (5.1) the only admissible similarity solutions are constant states separated by propagating contact discontinuities, cf. for instance [139]. Discontinuous, weak, solutions are physically acceptable here since the kinetic equation describes the conservation of the number of solitons within any given spectral interval, and the Rankine-Hugoniot type conditions can be imposed to ensure the conservation of the number of solitons across the discontinuities. As a result, the solution of the Riemann problem (4.3), (5.2) for each component w i (x, t) is composed of n + 1 constant states, or plateaus, separated by M discontinuities (see e.g. [133]): where the lower index i indicates the i-th component of the vector w, and the upper index j is the index of the plateau. For clarity we labeled the superscripts j = 1 as "L" (left boundary condition) and j = n + 1 as "R" (right boundary condition). The contact discontinuities propagate at the characteristic velocities [133]: where the plateaus' values w j i are found from the Rankine-Hugoniot jump conditions [33]: where i, j = 1 . . . n. The Rankine-Hugoniot conditions with i = j are trivially satisfied by the definition of contact discontinuity (5.4). We also note that conditions (5.4) are compatible with (5.5) for linearly degenerate systems of conservation laws. Note that, if the solitons were not interacting, the initial step distribution w i (x, 0) for the component λ = ζ i would have propagated at the free soliton velocity s 0i : which is very different compared to the solution (5.3).

Riemann problem: Examples
KdV soliton gas We now consider the simplest concrete Riemann problem for soliton gas. Namely, we shall study the collision of two single-component KdV soliton gases [75]. To be consistent with the notations of Sections 2.1 and 3.2 we shall be using the spectral variable η = √ −λ. Also we simplify the general notation of the previous section to make the specific results more transparent.
We consider the two-component ansatz (4.2) for the DOS, Substitution of (5.7) into the kinetic equation (2.10), (2.11) for KdV soliton gas yields the system of two conservation laws where the transport velocities s 1,2 ≡ s(η 1,2 , x, t) are given by (cf. (4.7)) where It is assumed that w 1 α 1 + w 2 α 2 < 1 (see the remark after eq. (4.7)). We now consider the system (5.9), (5.10) subject to the Riemann data (5.2) with w L = (w 10 , 0), w R = (0, w 20 ), i.e. w 1 (x, 0) = w 10 , w 2 (x, 0) = 0 , x < 0, where w 10 , w 20 > 0 are some constants satisfying w i0 η i /3 (see (2.43)). Note that the initial velocities of the colliding soliton gases are fully determined, via Eq. (5.10), by the density distribution (5.12). An example of one realisation of direct numerical simulation of the collision of two single-component KdV soliton gases modelled by the spectral Riemann problem (5.9), (5.12) is shown in Fig. 14. Details of the numerical implementation of the KdV soliton gas in the simulations can be found in [75]. The required solution for each component w i represents a combination of three constant states (plateaus) separated by two contact discontinuities given by equations (5.3) -(5.5) for n = 2. In the notations of this section, where the velocities s 1c = s 1 (w 1c , w 2c ) and s 2c = s 2 (w 1c , w 2c ) are determined by relations (5.10).  The piece-wise constant solution for β(x, t) is shown in Fig. 14 above the soliton gas plot. It is not difficult to show that the density of the two-component soliton gas in the interaction region, β c = w 1c + w 2c > w 10 , w 20 but β c < w 10 + w 20 . We emphasize that, although the soliton gas is initially prepared in a rarefied regime where solitons are spatially well-separated, the total density of solitons increases in the interaction region, and a more dense soliton gas is formed for which solitons exhibit significant overlap. The expanding interaction region in the soliton gas Riemann problem can be viewed as a stochastic counterpart of the traditional, coherent dispersive shock wave forming due to a dispersive regularisation of the Riemann initial data in the KdV equation [144], [18].
Comparisons of the theoretical results (5.16), (5.17), (5.19) for the soliton gas physical 'observables' c ± , β, u , u 2 with direct numerical simulations of the shock tube problem for the KdV soliton gas are presented in Figs. 15, 16.

Solition gases for the defocusing and resonant NLS equations
The spectral Riemann problem (5.1) (5.2) for bidirectional soliton gases for the defocusing and resonant NLS equations (see Section 2.2) was considered in [37], where the spectral Riemann problem solution (5.3) -(5.5) was constructed for the collisions of two-and three-component gas and comparisons of some analytically obtained 'observables' (the speeds of contact discontinuities and the spatial profiles of the mean density ρ ) with direct numerical simulations of the soliton gas collision have been made.  Figure 17: (Adapted from [37]). Spatiotemporal plots of the field ρ(x, t) for one realization of the anisotropic soliton gas collision. Trajectories of the solitons appear in solid lines. Red dash-dotted lines correspond to the trajectories of the contact discontinuities: x = z 1 t, x = z 2 t, cf. (5.4), and blue dashed lines to the free soliton trajectories: x = s 01 t, x = s 02 t. (a) Overtaking collisions (ζ 1 , ζ 2 ) = (1.05, 1.10), cf. initial condition (i) in Table 1. (b) Head-on collisions (ζ 1 , ζ 2 ) = (−1.05, 1.10), cf. initial condition (ii) in Table 1. Fig. 17 displays the numerically obtained x − t diagrams extracted from a single realisation of the soliton gas collision and illustrating the difference between overtaking and head-on interactions in the anisotropic soliton gas of the resonant NLS equation.
In Fig. 18 the comparisons are shown for the profiles of ensemble average ρ of the wave field ρ(x, t) associated with the spectral Riemann problem solution for the collision of threecomponent soliton gases for the defocusing and resonant NLS equations. The spectral data used for the numerical realisation of the soliton gases are collected in Table 1. The numerical evaluation of the ensemble average ρ involves averaging over 100 realisations. As expected, the solution is composed of 4 plateaus, where regions II and III contain at least two distinct  Figure 18: (Adapted from [37]). Comparison between the ensemble average ρ(x, t) obtained by averaging over 100 realisations of the numerical solution of the three-component soliton gas shock tube problem (solid line) for bidirectional dispersive hydrodynamics and the analytical solution obtained via the corresponding spectral Riemann problem (dash-dotted line). The dashed line corresponds to the average field density in the artificial soliton gas composed of noninteracting solitons with the spectral distribution given by (5.6). (a) Defocusing NLS soliton gas at t = 2000, case (iii) Table 1. (b) Resonant NLS soliton gas at t = 5000, case (iv) in Table 1.  soliton components and are region of interactions. The comparison shows an excellent agreement confirming the efficacy of the developed spectral kinetic theory. Note that in Fig. 18 the discontinuities in ρ(x, t) have a finite slope which is an artefact of the statistical averaging procedure detailed in [37].

Summary and Outlook
In this article we have reviewed the results on the spectral theory of soliton gases obtained by the author and his collaborators over the last two decades, although many of the developments are relatively recent. The theory is inspired by Zakharov's 1971 paper [28] where the ideas of the spectral kinetic description of a weakly interacting, rarefied soliton gas for the KdV equation were introduced for the first time. It has transpired later that the perspective of dispersive hydrodynamics, particularly the Whitham theory of modulated spectral finite-gap solutions provides an appealing mathematical framework for the generalisation of Zakharov's kinetic equation to the case of strongly interacting, dense soliton gas [32], [36], [38]. It has also been realised that soliton gas represents a ubiquitous physical phenomenon, a kind of strongly nonlinear turbulent wave motion, that can observed in the environmental conditions [21], [39] and realised in laboratory experiments [23], [24]. The recently discovered connections of the soliton gas dynamics and statistics with the topical areas of modulational instability and rogue wave formation [25], [26], [10] have provided the soliton gas research with further relevance. One can envisage many important and interesting mathematical and physical problems arising in connection with the developed spectral theory of soliton gas. Probably the most pertinent one is the determination of the statistical parameters of the integrable turbulence associated with a given spectral DOS of a soliton gas, in particular, obtaining the probability density function, the power spectrum and the autocorrelation function of the nonlinear wave field. This would provide the theoretical explanation of the intriguing statistical properties of integrable turbulence observed in the numerical simulations and physical experiments [8], [10], [25], [27], [145]. Further, the spectral hydrodynamics/kinetics of non-equilibrium soliton gas could provide the means for the manipulation of the integrable turbulence statistics, and therefore, control of the rogue wave emergence in physical systems described at leading order by integrable equations.
Another pertinent question is the generation of soliton gas out of 'non-gas' initial or boundary conditions. This problem is important from both the theoretical and experimental points of view. One possibility is to use the spontaneous 'soliton fission' mechanism [64], [31] inspired by the original Zabusky-Kruskal work [12] and employed in the experiments on the creation of a shallow water bidirectional anisotropic soliton gas in [23]. In modulationally unstable media a semiclassical scenario of the transition to integrable soliton turbulence via a chain of topological bifurcations of the finite-gap spectra was theoretically proposed in [63] and to some extent realised in an optics experiment reported in [54]. Of particular relevance is the IST based approach to the controlled (as opposed to the spontaneous) experimental realisation of soliton gas developed in [24]. This approach enables the generation of a soliton gas with a prescribed DOS, which is crucial for the verification and utilisation of the spectral kinetic theory.
Yet another prospective area of the further development of the soliton gas theory is the interaction of soliton gas with external potentials, either 'static' as in the soliton gas in a trapped BECs [48], [146] or dynamic, as in the recently proposed hydrodynamic soliton tunnelling framework [17], [147], [110]. Related to this, the development of the theory of soliton gas in perturbed integrable systems is of particular importance for applications where the higher order, nonintegrable corrections to the core integrable dynamics are always present and generally result in the slow evolution of the DOS, that is distinct from the spectral kinetic transport by the continuity equation, see [24].
Finally, the recently emerged intriguing parallels between the spectral theory of soliton gas and generalised hydrodynamics [57], [78] provide yet another avenue for the novel crossdisciplinary research that could be of significant benefit for both fields of research. near z = ∞ on the main sheet respectively, and the normalisation conditions requiring that all the periods of dp, dq are real (real normalised differentials). The homology basis α j , β j is defined in Section 3.3.1, see Fig. 11. The signs of integrals in (A.1), (A.2) will be opposite if we replace j by −j. The signs of the integrals in (A.1), (A.2) will be opposite if we replace j by −j.
The wavenumbers and frequencies are symmetrically extended to negative indices by k −j = k j , ω −j = ω j , j = 1, . . . , N, and similar equations for the 'tilded' quantities. They also satisfy the corresponding equations (A.1), (A.2), but the signs of integrals in (A.1), (A.2) will be opposite when we replace j by −j.
The proof that k j ,k j , ω j ,ω j defined by (A.1), (A.2) are indeed the wavenumbers and frequencies of the finite-gap focusing NLS solution associated with the spectral surface R of (3.57) can be found in [38]. It was then shown in [38] that for the wavenumbers and frequencies (A.1), (A.2) associated with the finite-gap focusing NLS solution the nonlinear dispersion relations are given byk whereγ is a large clockwise oriented contour containing Γ, and P j (z) and w j are defined in (3.62), (3.63) Taking real and imaginary parts of (A.4) and using the residues in the right hand side, we obtain separate systems for the solitonic components k m , ω m and the carrier components k m ,ω m .