# Dispersive hydrodynamics of soliton condensates for the Korteweg-de Vries
equation

T. Congy, G.A. El, G. Roberti, and A. Tovbis

###### Abstract

We consider large-scale dynamics of non-equilibrium dense soliton gas for the
Korteweg-de Vries (KdV) equation in the special “condensate” limit. We prove
that in this limit the integro-differential kinetic equation for the spectral
density of states reduces to the $N$-phase KdV-Whitham modulation equations
derived by Flaschka, Forest and McLaughlin (1980) and Lax and Levermore
(1983). We consider Riemann problems for soliton condensates and construct
explicit solutions of the kinetic equation describing generalized rarefaction
and dispersive shock waves. We then present numerical results for “diluted”
soliton condensates exhibiting rich incoherent behaviours associated with
integrable turbulence.

##  1 Introduction

Solitons represent the fundamental localised solutions of integrable nonlinear
dispersive equations such as the Korteweg-de Vries (KdV), nonlinear
Schrödinger (NLS), sine-Gordon, Benjamin-Ono and other equations. Along with
the remarkable localisation properties solitons exhibit particle-like elastic
pairwise collisions accompanied by definite phase/position shifts. A
comprehensive description of solitons and their interactions is achieved
within the inverse scattering transform (IST) method framework, where each
soliton is characterised by a certain spectral parameter related to the
soliton’s amplitude, and the phase related to its position (for the sake of
definiteness we refer here to the properties of KdV solitons). Generally,
integrable equations support $N$-soliton solutions which can be viewed as
nonlinear superpositions of $N$ solitons. Within the IST framework $N$-soliton
solution is characterised by a finite set of spectral and phase parameters
completely determined by the initial conditions for the integrable PDE.

The particle-like properties of solitons suggest some natural questions
pertaining to the realm of statistical mechanics, e.g. one can consider a
soliton gas as an infinite ensemble of interacting solitons characterised by
random spectral (amplitude) and phase distributions. The key question is to
understand the emergent macroscopic dynamics (i.e. hydrodynamics or kinetics)
of soliton gas given the properties of the elementary, “microscopic” two-
soliton interactions. It is clear that, due to the presence of an infinite
number of conserved quantities and the lack of thermalisation in integrable
systems the properties of soliton gases will be very different compared to the
properties of classical gases whose particle interactions are non-elastic.
Invoking the wave aspect of the soliton’s dual identity, soliton gas can be
viewed as a prominent example of integrable turbulence [1]. The pertinent
questions arising in this connection are related to the determination of the
parameters of the random nonlinear wave field in the soliton gas such as
probability density function, autocorrelation function, power spectrum etc.

The IST-based phenomenological construction of a rarefied, or diluted, gas of
KdV solitons was proposed in 1971 by V. Zakharov [2] who has formulated an
approximate spectral kinetic equation for such a gas based on the properties
of soliton collisions: the conservation of the soliton spectrum
(isospectrality) and the accumulation of phase shifts in pairwise collisions
that results in the modification of an effective average soliton’s velocity in
the gas. Zakharov’s spectral kinetic equation was generalised in [3] to the
case of a dense gas using the finite gap theory and the thermodynamic,
infinite-genus, limit of the KdV-Whitham modulation equations [4]. The results
of [3] were used in [5] for the formulation of a phenomenological construction
of kinetic equations for dense soliton gases for integrable systems describing
both unidirectional and bidirectional soliton propagation and including the
focusing, defocusing and resonant NLS equations, as well as the Kaup-
Boussinesq system for shallow-water waves [6]. The detailed spectral theory of
soliton and breather gases for the focusing NLS equation has been developed in
[7].

The spectral kinetic equation for a dense soliton gas represents a nonlinear
integro-differential equation describing the evolution of the density of
states (DOS)—the density function $u(\eta;x,t)$ in the phase space
$(\eta,x)\in\Gamma^{+}\times\mathbb{R}$, where
$\eta\in\Gamma^{+}\subset\mathbb{R}^{+}$ is the spectral parameter in the Lax
pair associated with the nonlinear integrable PDE,

$u_{t}+(us)_{x}=0,\quad
s(\eta,x,t)=s_{0}(\eta)+\int_{\Gamma^{+}}G(\eta,\mu)u(\mu,x,t)[s(\eta,x,t)-s(\mu,x,t)]d\mu\,.$
(1.1)

Here $s_{0}(\eta)$ is the velocity of a “free” soliton, and the integral term
in the second equation represents the effective modification of the soliton
velocity in the gas due to pairwise soliton collisions that are accompanied by
the phase-shifts described by the kernel $G(\eta,\mu)$. Both $s_{0}(\eta)$ and
$G(\eta,\mu)$ are system specific. In particular, for KdV $s_{0}=4\eta^{2}$
and $G(\eta,\mu)=\tfrac{1}{\eta}\ln\left|\tfrac{\mu+\eta}{\mu-\eta}\right|$.
The spectral support $\Gamma^{+}$ of the DOS is determined by initial
conditions. We note that, while $\Gamma^{+}\subset\mathbb{R}^{+}$ for the KdV
equation, one can have $\Gamma^{+}\subset\mathbb{C}^{+}$ for other equatons,
e.g. the focusing NLS equation, see [7] . Equation (1.1) describes the DOS
evolution in a dense soliton gas and represents a broad generalisation of
Zakharov’s kinetic equation for rarefied gas [2]. The existence, uniqueness
and properties of solutions to the “equation of state” (the integral equation
in (1.1) for fixed $x,t$) for the focusing NLS and KdV equations were studied
in [8].

The original spectral theory of the KdV soliton gas [3] has been developed
under the assumption that the spectral support $\Gamma^{+}$ of the DOS is a
fixed, simply-connected interval of $\mathbb{R}^{+}$; without loss of
generality one can assume $\Gamma^{+}=[0,1]$. In [7], [8] this restriction has
been removed by allowing the spectral support $\Gamma^{+}$ to be a union of
$N+1$ disjoint intervals $\gamma_{j}=[\lambda_{2j-1},\lambda_{2j}]$, termed
here s-bands: $\Gamma^{+}=\cup_{j=0}^{N}\gamma_{j}$,
($\gamma_{i}\cap\gamma_{j}=\emptyset$, $i\neq j$). In this paper we introduce
a further generalization of the existing theory by allowing the endpoints
$\lambda_{i}$ of the s-bands be functions of $x,t$. We show that this
generalization has profound implications for soliton gas dynamics, in
particular, the kinetic equation implies certain nonlinear evolution of the
endpoints $\lambda_{j}(x,t)$ of the s-bands. We determine this evolution for a
special type of soliton gases, termed in [7] soliton condensates. Soliton
condensate represents the “densest possible” gas whose DOS is uniquely defined
by a given spectral support $\Gamma^{+}$. The number $N$ of disjoint s-bands
in $\Gamma^{+}$ determines the genus $g=N-1$ of the soliton condensate. We
show that the evolution of $\lambda_{j}$’s in a soliton condensate is governed
by the $g$-phase averaged KdV-Whitham modulation equations [4], also derived
in the context of the semi-classical, zero-dispersion limit of the KdV
equation [9].

We then consider the soliton condensate dynamics arising in the Riemann
problem initiated by a rapid jump in the DOS. Our results suggest that in the
condensate limit the KdV dynamics of soliton gas is almost everywhere
equivalent to the (deterministic) generalised rarefaction waves (RWs) and
generalized dispersive shock waves (DSWs) of the KdV equation. We prove this
statement for the genus zero case and present a strong numerical evidence for
genus one. Our results also suggest direct connection of the “deterministic
KdV soliton gases” constructed in the recent paper [10] with modulated soliton
condensates.

Our work puts classical results of integrable dispersive hydrodynamics
(Flaschka-Forest-McLaughlin [4], Lax-Levermore [9], Gurevich-Pitaevskii [11])
in a broader context of the soliton gas theory. Namely, we show that the KdV-
Whitham modulation equations describe the emergent hydrodynamic motion of a
special soliton gas—a condensate—resulting from the accumulated effect of
“microscopic” two-soliton interactions. This new interpretation of the Whitham
equations is particularly pertinent in the context of generalized
hydrodynamics, the emergent hydrodynamics of quantum and classical many-body
systems [52]. The direct connection between the kinetic theory of KdV soliton
gas and generalized hydrodynamics was established recently in [53] (see also
[54] where the Whitham equations for the defocusing NLS equation were shown to
arise in the semi-classical limit of the generalised hydrodynamics of the
quantum Lieb-Liniger model).

Our work also paves the way to a major extension of the existing dispersive
hydrodynamic theory by including the random aspect of soliton gases. To this
end we consider “diluted” soliton condensates whose DOS has the same spectral
distribution as in genuine condensates but allows for a wider spacing between
solitons giving rise to rich incoherent dynamics associated with “integrable
turbulence” [1]. In particular, we show numerically that evolution of initial
discontinuities in diluted soliton condensates results in the development of
incoherent oscillating rarefaction and dispersive shock waves.

An important aspect of our work is the numerical modelling of soliton
condensates using $n$-soliton KdV solutions with large $n$ configured
according to the condensate density of states. The challenges of the numerical
implementation of standard $n$-soliton formulae for sufficiently large $n$ due
to rapid accumulation of roundoff errors are known very well. Here we use the
efficient algorithm developed in [45], which relies on the Darboux
transformation. We improve this algorithm following the recent methodology
developed in [46] for the focusing NLS equation with the implementation of
high precision arithmetic routine. Our numerical simulations show excellent
agreement with analytical predictions for the solutions of soliton condensate
Riemann problems and provide a strong support to the basic conjecture about
the connection of KdV soliton condensates with finite-gap potentials.

It should be noted that soliton condensates have been recently studied for the
focusing NLS equation, where they represent incoherent wave fields exhibiting
distinct statistical properties. In particular, it was shown in [55] that the
so-called bound state soliton condensate dynamics underies the long-term
behavior of spontaneous modulational instability, the fundamental physical
phenomenon that gives rise to the statistically stationary integrable
turbulence [56, 57].

The paper is organised as follows. In Section 2 we present a brief outline of
the spectral theory of soliton gas for the KdV equation and introduce the
notion of soliton condensate for the simplest genus zero case. In Section 3,
following [8], we generalize the spectral definition of soliton condensate to
an arbitrary genus case and prove the main Theorem 3.2 connecting spectral
dynamics of non-uniform soliton condensates with multiphase Whitham modulation
theory [4] describing slow deformations of the spectrum of periodic and
quasiperiodic KdV solutions. Section 3 is concerned with properties of KdV
solutions corresponding to the condensate spectral DOS, i.e. the soliton
condensate realizations. We formulate Conjecture 4.1 that any realization of
an equilibrium soliton condensate almost surely coincides with a finite-gap
potential defined on the condensate’s hyperelliptic spectral curve. This
proposition is proved for genus zero condensates and a strong numerical
evidence is provided for genus one and two. In Section 5 we construct
solutions to Riemann problems for the soliton gas kinetic equation subject to
discontinuous condensate initial data. These solutions describe evolution of
generalized rarefaction and dispersive shock waves. In Section 6 we present
numerical simulations of the Riemann problem for the KdV soliton condensates
and compare them with analytical solutions from Section 5. Finally, in Section
7 we consider basic properties of “diluted” condensates having a scaled
condensate DOS and exhibiting rich incoherent behaviors. In particular, we
present numerical solutions to Riemann problems for such diluted condensates.
Appendix A contains details of the numerical implementation of dense soliton
gases. In Appendix B we present results of the numerical realization of the
genus 2 soliton condensate and its comparison with two-phase solution of the
KdV equation.

##  2 Spectral theory of KdV soliton gas: summary of results

Here we present an outline of the spectral theory of KdV soliton gas developed
in [3, 12]. We consider the KdV equation in the form

$\varphi_{t}+6\varphi\varphi_{x}+\varphi_{xxx}=0.$ (2.1)

The inverse scattering theory associates soliton of the KdV equation (2.1)
with a point $z=z_{1}=-\eta_{1}^{2}$, $\eta_{1}>0$ of the discrete spectrum of
the Lax operator

$\mathcal{L}=-\partial_{xx}^{2}-\varphi(x,t),$ (2.2)

with sufficiently rapidly decaying potential $\varphi(x,t)$: $\varphi(x,t)\to
0$ as $|x|\to\infty$. The corresponding KdV soliton solution is given by

$\varphi_{\rm
s}(x,t;\eta_{1})=2\eta_{1}^{2}\hbox{sech}^{2}[\eta_{1}(x-4\eta_{1}^{2}t-x_{1}^{0})],$
(2.3)

where the soliton amplitude $a_{1}=2\eta_{1}^{2}$, the speed
$s_{1}=4\eta_{1}^{2}$, and $x_{1}^{0}$ is the initial position or ‘phase’.
Along with the simplest single-soliton solution (2.3) the KdV equation
supports $N$-soliton solutions $\varphi_{n}(x,t)$ characterized by $n$
discrete spectral parameters $0<\eta_{1}<\eta_{2}<\dots<\eta_{n}$ and the set
of initial positions $\\{x_{i}^{0}|i=1,\dots,n\\}$ associated with the phases
of the so-called norming constants [13]. It is also known that $n$-soliton
solutions can be realized as special limits of more general $n$-gap solutions,
whose Lax spectrum $\mathcal{S}_{n}$ consists of $N$ finite and one semi-
infinite bands separated by $n$ gaps [13],

$z\in\mathcal{S}_{n}=[\zeta_{1},\zeta_{2}]\cup[\zeta_{3},\zeta_{4}]\cup\dots\cup[\zeta_{2n+1},\infty).$
(2.4)

The $n$-gap solution of the KdV equation (2.1) represents a multiphase
quasiperiodic function

$\begin{split}&\varphi(x,t)=F_{n}(\theta_{1},\theta_{2},\dots,\theta_{n}),\quad\theta_{j}=k_{j}x-\omega_{j}t+\theta_{j}^{0},\\\
&F_{n}(\dots,\theta_{j}+2\pi,\dots)=F_{n}(\dots,\theta_{j},\dots),\end{split}$
(2.5)

where $k_{j}$ and $\omega_{j}$ are the wavenumber and frequency associated
with the $j$-th phase $\theta_{j}$, and $\theta_{j}^{0}$ are the initial
phases. Details on the explicit representation of the solution (2.5) in terms
of Riemann theta-functions can be found in classical papers and monographs on
finite-gap theory, see [59] and references therein.

The $n$-phase ($n$-gap) KdV solution (2.5) is parametrized by $2n+1$ spectral
parameters—the endpoints $\\{\zeta_{j}\\}_{j=1}^{2n+1}$ of the spectral bands.
The nonlinear dispersion relations (NDRs) for finite gap potentials can be
represented in the general form, see [4] for the concrete expressions,

$k_{j}=K_{j}(\zeta_{1},\dots,\zeta_{2n+1}),\quad\omega_{j}=\Omega_{j}(\zeta_{1},\dots,\zeta_{2n+1}),\quad
j=1,\dots,n,$ (2.6)

—and connect the wavenumber-frequency set $\\{k_{j},\omega_{j}\\}_{j=1}^{n}$
of (2.5) with the spectral set $\mathcal{S}_{n}$ (2.4). These are complemented
by the relation $\langle{\varphi}\rangle=\Phi(\zeta_{1},\dots,\zeta_{2n+1})$,
where $\langle\varphi\rangle=\int F_{n}d\theta_{1}\dots d\theta_{n}$ is the
mean obtained by averaging of $F_{n}$ over the phase $n$-torus
$\mathbb{T}^{n}=[0,2\pi)\times\dots\times[0,2\pi)$, assuming respective non-
commensurability of $k_{j}$’s and $\omega_{j}$’s and, consequently, ergodicity
of the KdV flow on the torus.

The $n$-soliton limit of an $n$-gap solution is achieved by collapsing all the
finite bands $[\zeta_{2j-1},\zeta_{2j}]$ into double points corresponding to
the soliton discrete spectral values,

$\zeta_{2j}-\zeta_{2j-1}\to 0,\ \ \zeta_{2j},\zeta_{2j-1}\to-\eta_{j}^{2},\ \
j=1,\dots,n.$ (2.7)

It was proposed in [3] that the special infinite-soliton limit of the spectral
$n$-gap KdV solutions, termed the thermodynamic limit, provides spectral
description the KdV soliton gas. The thermodynamic limit is achieved by
assuming a special band-gap distribution (scaling) of the spectral set
$\mathcal{S}_{n}$ for $n\to\infty$ on a fixed interval
$[\zeta_{1},\zeta_{2n+1}]$ (e.g. $[-1,0]$). Specifically, we set the spectral
bands to be exponentially narrow compared to the gaps so that
$\mathcal{S}_{n}$ is asymptotically characterized by two continuous
nonnegative functions on some fixed interval
$\Gamma^{+}\subset\mathbb{R}^{+}$: the density $\phi(\eta)$ of the lattice
points $\eta_{j}\in\Gamma^{+}$ defining the band centers via
$-\eta_{j}^{2}=(\zeta_{2j}+\zeta_{2j-1})/2$, and the logarithmic bandwidth
distribution $\tau(\eta)$ defined for $n\to\infty$ by

$\eta_{j}-\eta_{j+1}\sim\frac{1}{n\phi(\eta_{j})},\quad\tau(\eta_{j})\sim-\frac{1}{n}\ln(\zeta_{2j}-\zeta_{2j-1}).$
(2.8)

The scaling (2.8) was originally introduced by Venakides [14] in the context
of the continuum limit of theta functions.

Complementing the spectral distributions (2.8) with the uniform distribution
of the initial phase vector ${\boldsymbol{\theta}}^{0}$ on the torus
$\mathbb{T}^{n}$ we say that the resulting random finite gap solution
$\varphi(x,t)$ approximates soliton gas as $n\to\infty$. An important
consequence of this definition of soliton gas is ergodicity, implying that
spatial averages of the KdV field in a soliton gas are equivalent to the
ensemble averages, i.e. the averages over $\mathbb{T}^{n}$ in the
thermodynamic limit $n\to\infty$. We shall use the notation $\langle
F[\varphi]\rangle$ for ensemble averages and $\overline{F[\varphi]}$ for
spatial averages.

From now on we shall refer to $\eta$ as the spectral parameter and
$\Gamma^{+}$–the spectral support. The density of states (DOS) $u(\eta)$ of a
spatially homogeneous (equilibrium) soliton gas is phenomenologically
introduced in such a way that $u(\eta_{0})d\eta dx$ gives the number of
solitons with the spectral parameter $\eta\in[\eta_{0};\eta_{0}+d\eta]$
contained in the portion of soliton gas over a macroscopic (i.e. containing
sufficiently many solitons) spatial interval
$x\in[x_{0},x_{0}+dx]\subset\mathbb{R}$ for any $x_{0}$ (the individual
solitons can be counted by cutting out the relevant portion of the gas and
letting them separate with time). The corresponding spectral flux density
$v(\eta)$ represents the temporal counterpart of the DOS i.e.
$v(\eta_{0})d\eta$ is the number of solitons with the spectral parameter
$\eta\in[\eta_{0};\eta_{0}+d\eta]$ crossing any given point $x=x_{0}$ per unit
interval of time. These definitions are physically suggestive in the context
of rarefied soliton gas where solitons are identifiable as individual
localized wave structures. The general mathematical definitions of $u(\eta)$
and $v(\eta)$ applicable to dense soliton gases are introduced by applying the
thermodynamic limit to the finite-gap NDRs (2.6), leading to the integral
equations [3, 12]:

$\displaystyle\int_{\Gamma^{+}}\ln\left|\frac{\mu+\eta}{\mu-\eta}\right|u(\mu)d\mu+u(\eta){\sigma}(\eta)=\eta,$
(2.9)
$\displaystyle\int_{\Gamma^{+}}\ln\left|\frac{\mu+\eta}{\mu-\eta}\right|v(\mu)d\mu+v(\eta){\sigma}(\eta)=4\eta^{3},$
(2.10)

for all $\eta\in\Gamma^{+}$. Here the spectral scaling function
$\sigma:\Gamma^{+}\to[0,\infty)$ is a continuous non-negative function that
encodes the Lax spectrum of soliton gas via
$\sigma(\eta)=\phi(\eta)/\tau(\eta)$. Equations (2.9), (2.10) represent the
soliton gas NDRs.

Eliminating $\sigma(\eta)$ from the NDRs (2.9), (2.10) yields the equation of
state for KdV soliton gas:

$s(\eta)=4\eta^{2}+\frac{1}{\eta}\int_{\Gamma^{+}}\log\left|\frac{\eta+\mu}{\eta-\mu}\right|u(\mu)[s(\eta)-s(\mu)]d\mu,$
(2.11)

where $s(\eta)=v(\eta)/u(\eta)$ can be interpreted as the velocity of a tracer
soliton in the gas. It was shown in [3] that for a weakly non-uniform (non-
equilibrium) soliton gas, for which $u(\eta)\equiv u(\eta;x,t)$,
$s(\eta)\equiv s(\eta;x,t)$, the DOS satisfies the continuity equation

$u_{t}+(us)_{x}=0,$ (2.12)

so that $s(\eta;x,t)$ acquires the natural meaning of the transport velocity
in the soliton gas. Equations (2.12), (2.11) form the spectral kinetic
equation for soliton gas. One should note that the typical scales of spatio-
temporal variations in the kinetic equation (2.12) are much larger than in the
KdV equation (2.1), i.e. the kinetic equation describes macroscopic evolution,
or hydrodynamics, of soliton gases.

Let the spectral support $\Gamma^{+}$ be fixed. Then, differentiating equation
(2.9) with respect to $t$, equation (2.10) with respect to $x$, and using the
continuity equation (2.12) we obtain the evolution equation for the spectral
scaling function

$\sigma_{t}+s\sigma_{x}=0,$ (2.13)

which shows that $\sigma(\eta;x,t)$ plays the role of the Riemann invariant
for the spectral kinetic equation.

Finally, the ensemble averages of the conserved densities of the KdV wave
field in the soliton gas (the Kruskal integrals) are evaluated in terms of the
DOS as $\langle{\cal
P}_{n}[\varphi]\rangle=C_{n}\int_{\Gamma^{+}}\eta^{2n-1}u(\eta)d\eta$, where
${\cal P}_{n}[\varphi]$ are conserved quantities of the KdV equation and
$C_{n}$ constants [3, 12] (see also [15] for rigorous derivation in the NLS
context). In particular, for the two first moments we have, on dropping the
$x,t$-dependence [3, 12],

$\langle\varphi\rangle=4\int_{\Gamma^{+}}\eta\,u(\eta)d\eta,\quad\langle\varphi^{2}\rangle=\frac{16}{3}\int_{\Gamma^{+}}\eta^{3}\,u(\eta)d\eta.$
(2.14)

We note that in the original works on KdV soliton gas it was assumed
(explicitly or implicitly) that the spectral support $\Gamma^{+}$ of the KdV
soliton gas is a fixed, simply connected interval (without loss of generality
one can assume that in this case $\Gamma^{+}=[0,1]$). In what follows we will
be considering a more general configuration where $\Gamma^{+}$ represents a
union of $N+1$ disjoint intervals.

A special kind of soliton gas, termed soliton condensate, is realized
spectrally by letting $\sigma\to 0$ in the NDRs (2.9), (2.10). This limit was
first considered in [7] for the soliton gas in the focusing NLS equation and
then in [8] for KdV. Loosely speaking, soliton condensate can be viewed as the
“densest possible” gas (for a given spectral support $\Gamma^{+}$) whose
properties are fully determined by the interaction (integral) terms in the
NDRs (2.9), (2.10).

For the KdV equation, setting $\sigma=0$ and, considering the simplest case
$\Gamma^{+}=[0,1]$ in (2.9), (2.10), we obtain the soliton condensate NDRs
[12]:

$\int_{0}^{1}\ln\left|\frac{\mu+\eta}{\mu-\eta}\right|u(\mu)d\mu=\eta,\quad\int_{0}^{1}\ln\left|\frac{\mu+\eta}{\mu-\eta}\right|v(\mu)d\mu=4\eta^{3}.$
(2.15)

These are solved by

$u(\eta)=\frac{\eta}{\pi\sqrt{1-\eta^{2}}},\quad
v(\eta)=\frac{6\eta(2\eta^{2}-1)}{\pi\sqrt{1-\eta^{2}}},$ (2.16)

as verified by direct substitution (it is advantageous to first differentiate
equations (2.15) with respect to $\eta$). The formula (2.16) for $u(\eta)$ is
sometimes called the Weyl distribution, following the terminology from the
semiclassical theory of linear differential operators [9].

######  Remark 2.1.

The meaning of the zero $\eta_{0}=1/{\sqrt{2}}$ of $v(\eta)$ is that all the
tracer solitons with the spectral parameter $\eta>\eta_{0}$ move to the right,
whereas all the tracer solitons with $\eta<\eta_{0}$ move in the opposite
direction while the tracer soliton with $\eta=\eta_{0}$ is stationary. The
somewhat counter-intuitive “backflow” phenomenon (we remind that KdV solitons
considered in isolation move to the right) has been observed in the numerical
simulations of the KdV soliton gas [16] and can be readily understood from the
phase shift formula of two interacting solitons, where the the larger soliton
gets a kick forward upon the interaction while the smaller soliton is pushed
back. As a matter of fact, the KdV soliton backflow is general and can be
observed for a broad range of sufficiently dense gases (see Fig. 16 in Section
7.1 for the numerical illustration).

##  3 Soliton condensates and their modulations

We now consider the general case of the soliton gas NDRs (2.9), (2.10) by
letting the support $\Gamma^{+}\subset\mathbb{R}^{+}$ of $u(\eta),v(\eta)$ to
be a union of disjoint intervals $\gamma_{k}\subset{\mathbb{R}}^{+}$ with
endpoints $\lambda_{j}>0$, $j=1,2,\dots,2N+1$, where
$\gamma_{0}=[0,\lambda_{1}]$ and $\gamma_{k}=[\lambda_{2k},\lambda_{2k+1}]$,
$k=1,\dots,N$, i.e.

$\Gamma^{+}=[0,\lambda_{1}]\cup[\lambda_{2},\lambda_{3}]\cup\dots\cup[\lambda_{2N},\lambda_{2N+1}].$
(3.1)

We shall call the intervals $\gamma_{k}$ the s-bands, and the soliton gas
spectrally supported on $\Gamma^{+}$ (3.1)— the genus $N$ soliton gas.
Correspondingly, we refer to the intervals
$c_{j}=(\lambda_{2j-1},\lambda_{2j})$ separating the s-bands as to s-gaps.
Note that the s-bands and s-gaps are different from the original bands and
gaps in the spectrum $\mathcal{S}_{n}$ of finite-gap potential (cf. (2.4)) as
they emerge after the passage to the thermodynamic limit: loosely speaking,
one can view the s-bands as a continuum limit of the “thermodynamic band
clusters”, each representing an isolated dense subset of
$\mathcal{S}_{\infty}$ consisting of the collapsing original bands. The
existence and uniqueness of solutions $u(\eta)$, $v(\eta)$ for (2.9), (2.10)
respectively, as well as the fact that $u(\eta)\geq 0$ on $\Gamma^{+}$ with
some mild constraints, was established in [8]. Our goal here is to find
explicit expressions for $u,v$ for the genus $N$ soliton condensate, that is,
solutions of (2.9), (2.10) for the particular case ${\sigma}\equiv 0$ on
$\Gamma^{+}$.

Denote by $\Gamma^{-}$ the symmetric image of $\Gamma^{+}$ with respect to the
origin, i.e., $\Gamma^{-}=-\Gamma^{+}$. If we take the odd continuation of
$u,v$ to $\Gamma^{-}$ (preserving the same notations), we observe that
equations (2.9), (2.10) become

$\displaystyle-\int_{\Gamma}\ln|\mu-\eta|u(\mu)d\mu+u(\eta){\sigma}(\eta)=\eta,$
(3.2)
$\displaystyle-\int_{\Gamma}\ln|\mu-\eta|v(\mu)d\mu+v(\eta){\sigma}(\eta)=4\eta^{3},$
(3.3)

where $\Gamma:=\Gamma^{+}\cup\Gamma^{-}$, for all $\eta\in\Gamma^{+}$. In
fact, if we symmetrically extend ${\sigma}(\eta)$ from $\Gamma^{+}$ to
$\Gamma$, equations (3.2), (3.3) should be valid on $\Gamma$ since every term
in these equations is odd. The expressions (2.14) for the first two moments
(ensemble averages) of the KdV wave field in the soliton gas become

$\langle\varphi\rangle=2\int_{\Gamma}\eta\,u(\eta)d\eta,\
\qquad\langle\varphi^{2}\rangle=\frac{8}{3}\int_{\Gamma}\eta^{3}\,u(\eta)d\eta.$
(3.4)

We now consider soliton condensate of genus $N$ by setting ${\sigma}\equiv 0$
in (3.2), (3.3). Then, differentiating in $\eta$ we obtain

$H[u]=\frac{1}{\pi},\quad H[v]=\frac{12\eta^{2}}{\pi}\quad\text{on $\Gamma$},$
(3.5)

where $H$ denotes the Finite Hilbert Transform (FHT) on $\Gamma$, see for
example [17], [18],

$H[f](\xi)=\frac{1}{\pi}\int\limits_{\Gamma}\frac{f(y)dy}{y-\xi}.$ (3.6)

Equations (3.5) are the (transformed) NDRs for the KdV soliton condensate.

To find $u,v$ for the soliton condensate, it is sufficient to invert the FHT
$H$ on $\Gamma$. Denote by $\mathcal{R}_{2N}$ the hyperelliptic Riemann
surface of the genus $2N$, defined by the branchcuts (s-bands) $\gamma_{k}$,
$k=0,\pm 1,\dots,\pm N$, where $\gamma_{-k}=-\gamma_{k}$. Define two
meromorphic differentials of second kind, $dp$ and $dq$ on $\mathcal{R}_{2N}$
by

$dp=\frac{iP(\eta)}{2\pi R(\eta)}d\eta,\qquad dq=\frac{2iQ(\eta)}{\pi
R(\eta)}d\eta,$ (3.7)

where

$R(\eta)=\sqrt{(\eta^{2}-\lambda_{1}^{2})(\eta^{2}-\lambda_{2}^{2})\dots(\eta^{2}-\lambda_{2N+1}^{2})},$
(3.8)

and $P,Q$ are odd monic polynomials of degree $2N+1$ and $2N+3$ respectively
that are chosen so that all their s-gap integrals are zero, i.e.

$\int\limits_{\lambda_{2j-1}}^{\lambda_{2j}}dp=\int\limits_{\lambda_{2j-1}}^{\lambda_{2j}}dq=0,\quad
j=1,\dots,N.$ (3.9)

Equivalently, one can say that $dp,dq$ are real normalized differentials. Note
that Equations (3.7), (3.9) uniquely define $dp,dq$.

######  Theorem 3.1.

Functions $u(\eta)=dp/d\eta$ and $v(\eta)=dq/d\eta$ defined by (3.7) and (3.9)
satisfy the respective equations (3.5) and are odd and real valued on
$\Gamma$. Thus $u,v$ are the solutions of NDRs (2.9), (2.10) for $\sigma=0$.
Moreover, $u(\eta)\geq 0$ on $\Gamma^{+}$. Here the value of $R(\eta)$ for
$\eta\in\Gamma$ is taken on the positive (upper) shore of the branchcut.

Theorem 3.1 for $u$ was proven in [8], Section 6, for the so-called bound
state soliton condensate. The proof for KdV is analogous. The proof for $v$
goes along the same lines, except $v(\eta)$ attains different signs.

######  Remark 3.1.

The normalization (3.9) requires that both polynomials $P,Q$ have zeros in
every of the $2N$ gaps on $[-\lambda_{2N+1},\lambda_{2N+1}]$. Note also that
$P(0)=Q(0)=0$. That takes care of all the zeros of $P$. The polynomial $Q$ has
two additional symmetric real zeros $\pm\eta_{0}$ that must be located on some
band $\gamma_{k}$ and its symmetrical image $\gamma_{-k}$, see below. In the
case $N=0$ such zeros are $\eta_{0}=\pm\frac{1}{\sqrt{2}}$, see (2.16). Let us
prove that $\eta_{0}$ belongs to a band. It is easy to see that the zero level
curves $\Im\int_{0}^{\eta}dp=0$ consist of all bands and the imaginary axis,
whereas the zero level curves $\Im\int_{0}^{\eta}dq=0$ consist of that of
$\Im\int_{0}^{\eta}dp=0$ with an extra two curves crossing ${\mathbb{R}}$ at
$\pm\eta_{0}$ and approaching $z=\infty$ with angles $\pm\frac{\pi}{6}$ and
$\pm\frac{5\pi}{6}$ respectively. Note that there must be four zero level
curves passing through $\pm\eta_{0}$ and, therefore, they must be on the
bands.

Thus, for the soliton condensate of genus $N$, we obtain, on using Theorem 3.1
and Equation (3.7),

$u(\eta)\equiv
u^{(N)}(\eta;\lambda_{1},\dots\lambda_{2N+1})=\frac{iP(\eta)}{2\pi
R(\eta)},\quad v(\eta)\equiv
v^{(N)}(\eta;\lambda_{1},\dots,\lambda_{2N+1})=\frac{2iQ(\eta)}{\pi R(\eta)}.$
(3.10)

Figure 1: Spectral distributions (3.10) for genus 2 soliton condensate. a)
Density of states $u^{(2)}(\eta;{\boldsymbol{\lambda}})$. b) Tracer velocity
$s^{(2)}(\eta;{\boldsymbol{\lambda}})$. Here
${\boldsymbol{\lambda}}=(\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4},\lambda_{5})=(0.3,0.5,0.7,0.9,1)$.

The velocity of a tracer soliton with the spectral parameter
$\eta\in\Gamma^{+}$ propagating in the soliton condensate with DOS $u(\eta)$
is then found as

$s(\eta)\equiv
s^{(N)}(\eta;\lambda_{1},\dots,\lambda_{2N+1})=\frac{v(\eta)}{u(\eta)}=\frac{4Q(\eta)}{P(\eta)}.$
(3.11)

As an illustrative example we present in Fig. 1 the plots of the DOS and
tracer velocity for the genus 2 soliton condensate.

We now consider slow modulations of non-equilibrium (non-uniform) soliton
condensates by assuming $u\equiv u(\eta;x,t)$, $v\equiv v(\eta;x,t)$,
$\Gamma\equiv\Gamma(x,t)$. Equations (2.12), (3.10) then yield the kinetic
equation for genus $N$ soliton condensate:

$\left(\frac{P}{R}\right)_{t}+\left(\frac{4Q}{R}\right)_{x}=0,$ (3.12)

that is valid for $\eta\in\Gamma=\cup_{k=-N}^{N}(\gamma_{k})$. The velocity
(3.11) then assumes the meaning of the tracer, or transport, velocity in a
non-uniform genus $N$ soliton condensate.

######  Theorem 3.2.

The kinetic equation (3.12) for soliton condensate implies the evolution of
the endpoints $\lambda_{j}$, $j=1,\dots,2N+1$ according to the Whitham
modulation equations

$\partial_{t}\lambda_{j}+V_{j}({\boldsymbol{\lambda}})\partial_{x}\lambda_{j}=0,\quad
j=1,\dots,2N+1,$ (3.13)

where ${\boldsymbol{\lambda}}=(\lambda_{1},\dots,\lambda_{2N+1})$ and

$V_{j}(\boldsymbol{\lambda})=s^{(N)}(\lambda_{j};\lambda_{1},\dots,\lambda_{2N+1})=\frac{4Q(\lambda_{j})}{P(\lambda_{j})}.$
(3.14)

Proof. (See [19]) Multiplying (3.12) by $(\eta^{2}-\lambda_{j}^{2})^{3/2}$ and
passing to the limit $\eta\to\lambda_{j}$ we obtain equations (3.13), (3.14)
for the evolution of the spectral $s$-bands (i.e. the evolution of
$\Gamma(x,t)$). These are the KdV-Whitham modulation equations [4], [9] (see
also Remark 3.2 below).

∎

######  Corollary 3.1.

The endpoints of the “special” band
$\gamma_{k}=[\lambda_{2k},\lambda_{2k+1}]$, $k\neq 0$, containing the point
$\eta_{0}$ of zero tracer speed, $s(\eta_{0})=0$, are moving in opposite
directions, whereas all the endpoints on the same side from $\eta_{0}$ are
moving in the same direction. See Fig. 1 (right) for $N=2$

######  Remark 3.2.

Modulation equations (3.12), (3.13) were originally derived by Flaschka,
Forest and McLaughlin [4] by averaging the KdV equation over the multiphase
(finite-gap) family of solutions. These equations along with the condensate
NDRs (3.5), also appear in the seminal work of Lax and Levermore [9] in the
context of the semiclassical (zero-dispersion) limit of multi-soliton KdV
ensembles (see Section 5 and, in particular, Equation (5.23) in [9]). A
succinct exposition of the spectral Whitham theory for the KdV equation can be
found in Dubrovin and Novikov [19]).

######  Remark 3.3.

The Whitham modulation equations (3.13), (3.14) are locally integrable for any
$N$ via Tsarev’s generalized hodograph transform [20, 19]. Moreover, by
allowing the genus $N$ to take different values in different regions of
$x,t$-plane, $N=N(x,t)$, global solutions of the KdV-Whitham system can be
constructed for a broad class of initial data (see Section 4.2 for further
details). Invoking the definitive property $\sigma\equiv 0$ of a soliton
condensate, the existence of the solution to an initial value problem for the
Whitham system for all $t>0$ implies that this property will remain invariant
under the $t$-evolution, i.e. soliton condensate will remain a condensate
during the evolution, however its genus can change.

The finite-genus Whitham modulation system (3.13), (3.14) can be viewed as an
exact hydrodynamic reduction of the full kinetic equation (2.12), (2.11) under
the ansatz (3.10), (3.11). Recalling the origin of the soliton gas kinetic
equation as a singular, thermodynamic limit of the Whitham equations [3] the
recovery of the finite-genus Whitham dynamics in the condensate limit might
not look surprising. On the other hand, viewed from the general soliton gas
perspective the condensate reduction notably shows that the highly nontrivial
nonlinear modulation (hydro)dynamics emerges as a collective effect of the
elementary two-soliton scattering events. This understanding is in line with
ideas of generalised hydrodynamics, a powerful theoretical framework for the
description of non-equilibrium macroscopic dynamics in many-body quantum and
classical integrable systems [21]. The connection of the KdV soliton gas
theory with generalised hydrodynamics has been recently established in [22].
Relevant to the above, it was shown in [23] that the semiclassical limit of
the generalised hydrodynamics for the Lieb-Liniger model of Bose gases yields
the Whitham modulation system for the defocusing NLS equation.

A different type of hydrodynamic reductions of the soliton gas kinetic
equation defined by the multi-component delta-function ansatz
$u(\eta,x,t)=\sum_{i=1}^{m}w_{i}(x,t)\delta(\eta-\eta_{j})$ for the DOS has
been studied in [24] for $\eta_{j}=const$ and in [25, 26] for
$\eta_{j}=\eta_{j}(x,t)$. One of the definitive properties of the
multicomponent hydrodynamic reductions of this type is their linear degeneracy
which, in particular, implies the absence of the wavebreaking and the
occurrence of contact discontinuities in the solutions of Riemann problems
[27]. In contrast, the condensate (Whitham) system (3.13), (3.14) obtained
under the condition $\sigma\equiv 0$ is known to be genuinely nonlinear,
$\partial V_{j}/\partial\lambda_{j}\neq 0$, $j=1,\dots 2N+1$ [28] implying the
inevitability of wavebreaking for general initial data, which is in stark
contrast with linear degeneracy of the multicomponent “cold-gas” hydrodynamic
reductions. Reconciling the genuine nonlinearity property of soliton
condensates with linearly degenerate non-condensate multicomponent cold-gas
dynamics is an interesting problem which will be considered in future
publications.

Thus, we have shown that the spectral dynamics of soliton condensates are
equivalent to those of finite gap potentials, which naturally suggests a close
connection (or even equivalence) of these two objects at the level of
realizations, i.e. the corresponding solutions $\varphi(x,t)$ of the KdV
equation. This connection will be explored in the next section using a
combination of analytical results and numerical simulations for genus 0 and
genus 1 soliton condensates.

##  4 Genus 0 and genus 1 soliton condensates

Having developed the spectral description of KdV soliton condensates, we now
look closer at the two simplest representatives: genus 0 and genus 1
condensates. In particular, we shall be interested in the characterization of
the realizations of soliton condensates, i.e. the KdV solutions, denoted
$\varphi_{\rm c}^{(N)}(x,t)$, corresponding to the condensate spectral DOS
$u^{(N)}(\eta)$ for $N=0,1$. We do not attempt here to construct the soliton
gas realizations explicitly via the thermodynamic limit of finite gap
potentials (see Section 2), instead, we infer some of their key properties
from the expressions (3.4) for the ensemble averages as integrals over the
spectral DOS. We then conjecture the exact form of soliton condensate
realizations and support our conjecture by detailed numerical simulations.

###  4.1 Equilibrium properties

####  4.1.1 Genus 0

For $N=0$ equations (3.10) for the DOS and the spectral flux density yield
(cf. (2.16))

$u(\eta)=u^{(0)}(\eta;\lambda_{1})\equiv\frac{\eta}{\pi\sqrt{\lambda_{1}^{2}-\eta^{2}}},\quad
v(\eta)=v^{(0)}(\eta;\lambda_{1})\equiv\frac{6\eta(2\eta^{2}-\lambda_{1}^{2})}{\pi\sqrt{\lambda_{1}^{2}-\eta^{2}}},$
(4.1)

so that the tracer velocity (cf. (3.11))

$s(\eta)=s^{(0)}(\eta;\lambda_{1})=6(2\eta^{2}-\lambda_{1}^{2}).$ (4.2)

Next, substituting (4.1) in (3.4) (where $\Gamma=[-\lambda_{1},\lambda_{1}]$
or equivalently, $\Gamma^{+}=[0,\lambda_{1}]$), we obtain for the ensemble
averages:

$\langle\varphi\rangle=\lambda_{1}^{2},\quad\langle\varphi^{2}\rangle=\lambda_{1}^{4},$
(4.3)

where $\varphi\equiv\varphi_{\rm c}^{(0)}(x,t)$. Thus the variance
$\Delta=\sqrt{\langle\varphi^{2}\rangle-\langle\varphi\rangle^{2}}=\sqrt{\langle(\varphi-\langle\varphi\rangle)^{2}\rangle}=0$,
which implies (see, e.g. [29]) that genus 0 soliton condensate is almost
surely described by a constant solution of the KdV equation , i.e.

$\varphi=\langle\varphi\rangle=\lambda_{1}^{2}$ (4.4)

(note that constant solution is classified as a genus 0 KdV potential).

This result can be intuitively understood by identifying soliton condensate
with the “densest possible” soliton gas for a given spectral support $\Gamma$.
The densest “packing” for genus 0 is achieved by distributing soliton
parameters according to the spectral DOS $u(\eta)$ (4.1) which results in the
individual solitons “merging” into a uniform KdV field of amplitude
$\lambda_{1}^{2}$. The numerical implementation of soliton condensate
realizations, using $n$-soliton KdV solution with $n$ large, shows that the
condensate DOS (4.1) is only achievable within this framework if all $n$
solitons in the solution have the same phase of the respective norming
constants. Invoking the interpretation of the phase of the norming constant as
the soliton position in space [13, 30] one can say that in the condensate all
solitons are placed at the same point, say $x=0$ (cf. Appendix A for a
mathematical justification). Details of the numerical implementation of KdV
soliton gas using $n$-soliton solutions can be found in Appendix A. Fig. 2
displays the realization $\varphi_{\rm c}^{(0)}(x)$ of genus 0 soliton
condensate with $\lambda_{1}=1$ modeled by $n$-soliton solutions
$\varphi_{n}(x)$ with $n=100$ and $n=200$, along with the absolute errors
$\varphi_{n}(x)-1$; in the following we refer to these $n$-soliton solutions
as “numerical realizations” of the soliton gas. One can see that the error at
the center of the numerical domain, where the gas is nearly uniform, is very
small: Fig. 3 displays the variation of this error with $n$ and shows that it
decreases with $1/n^{2}$.

Figure 2: a) Comparison between numerical realizations of genus 0 condensate
generated with $100$ solitons (dashed line), $200$ solitons (black solid
line), and the constant KdV solution $\varphi=1$ (red solid line). b)
Corresponding absolute errors $|\varphi_{n}(x)-1|$ obtained with $50$ solitons
(dashed line), $100$ solitons (solid line) and $200$ (dash-dotted line); the
absolute error is evaluated at the extrema of the oscillations.  Figure 3:
Variation of the absolute error $|\varphi_{n}(x)-1|$ at the center of the
numerical domain $x=0$ (cf. Fig. 2). The markers correspond to the error
obtained numerically and the solid line the corresponding fit $\alpha/n^{2}$
where $\alpha\approx 0.25$.

The numerical approximation used here is similar to the approximation of the
soliton condensate of the focusing NLS equation via a $n$-soliton solution
presented in [50]. In the latter case the uniform wavefield limit as a central
part of the so-called “box potential”, is also reached when the complex phases
of the norming constants are chosen deterministically. The absolute error—the
difference between the $n$-soliton solution and the expected constant value of
the wavefield—measured at the center of the numerical realization—follows a
different scaling law and is proportional to $n^{-1/2}$.

####  4.1.2 Genus 1

We now consider the case of genus 1 soliton condensate. For $N=1$

$R(\eta)=\sqrt{(\eta^{2}-\lambda_{1}^{2})(\eta^{2}-\lambda_{2}^{2})(\eta^{2}-\lambda_{3}^{2})}$
(4.5)

is purely imaginary on
$\Gamma=[-\lambda_{3},-\lambda_{2}]\cup[-\lambda_{1},\lambda_{1}]\cup[\lambda_{2},\lambda_{3}]$.
According to Theorem 3.1

$\displaystyle
u(\eta)=u^{(1)}(\eta;\lambda_{1},\lambda_{2},\lambda_{3})\equiv\frac{i\eta(\eta^{2}-w^{2})}{\pi
R(\eta)},$ (4.6) $\displaystyle
v(\eta)=v^{(1)}(\eta;\lambda_{1},\lambda_{2},\lambda_{3})\equiv\frac{12i\eta(\eta^{4}-h^{2}\eta^{2}-r^{2})}{\pi
R(\eta)},$ (4.7)

where $h^{2}=\frac{\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2}}{2}$
follows from the fact that $\left.-i{\rm
Res}\frac{Q(\zeta)}{(\zeta-\eta)R(\zeta)}\right|_{\zeta=\infty}=-6\eta^{2}$.
The normalization conditions (3.9) imply that

$w^{2}=\frac{\int_{\lambda_{1}}^{\lambda_{2}}\frac{y^{3}dy}{R(y)}}{\int_{\lambda_{1}}^{\lambda_{2}}\frac{ydy}{R(y)}},\quad
r^{2}=\frac{\int_{\lambda_{1}}^{\lambda_{2}}\frac{y^{5}-\frac{\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2}}{2}y^{3}}{R(y)}dy}{\int_{\lambda_{1}}^{\lambda_{2}}\frac{ydy}{R(y)}}.$
(4.8)

Using 3.131.3 and 3.132.2 from [31], we calculate

$w^{2}=\lambda_{3}^{2}-(\lambda_{3}^{2}-\lambda_{1}^{2})\mu(m),\quad{\rm
where}\quad\mu(m)=\frac{{\rm E}\left(m\right)}{{\rm K}\left(m\right)}\quad{\rm
and}\quad
m=\frac{\lambda_{2}^{2}-\lambda_{1}^{2}}{\lambda_{3}^{2}-\lambda_{1}^{2}}.$
(4.9)

Calculation of $r^{2}$ is a bit more involved as it is based on the
observation

$\begin{split}\int_{\lambda_{1}}^{\lambda_{2}}\frac{y^{5}}{R(y)}dy&=\frac{1}{2}\int_{\lambda_{1}^{2}}^{\lambda_{2}^{2}}\frac{z^{2}dz}{R(z^{\frac{1}{2}})},\\\
&=\frac{\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2}}{3}\int_{\lambda_{1}^{2}}^{\lambda_{2}^{2}}\frac{zdz}{R(z^{\frac{1}{2}})}-\frac{\lambda_{1}^{2}\lambda_{2}^{2}+\lambda_{1}^{2}\lambda_{3}^{2}+\lambda_{2}^{2}\lambda_{3}^{2}}{6}\int_{\lambda_{1}^{2}}^{\lambda_{2}^{2}}\frac{dz}{R(z^{\frac{1}{2}})}.\end{split}$
(4.10)

Using (4.8), (4.10), we obtain after some algebra

$r^{2}=\frac{1}{6}\left[\lambda_{3}^{2}(\lambda_{3}^{2}-\lambda_{2}^{2}-\lambda_{1}^{2})-2\lambda_{2}^{2}\lambda_{1}^{2}-(\lambda_{3}^{2}+\lambda_{1}^{2}+\lambda_{2}^{2})(\lambda_{3}^{2}-\lambda_{1}^{2})\mu(m)\right].$
(4.11)

Thus, the velocity of a tracer soliton with spectral parameter
$\eta\in\Gamma^{+}$ in the genus 1 soliton condensate, characterized by DOS
(4.6), is given by

$\begin{split}s(\eta)&\equiv
s^{(1)}(\eta;\lambda_{1},\lambda_{2},\lambda_{3})=12\frac{\eta^{4}-\frac{\lambda_{2}^{2}+\lambda_{3}^{2}+\lambda_{1}^{2}}{2}\eta^{2}-r^{2}}{\eta^{2}-w^{2}}\\\
&=12\frac{\eta^{4}-\frac{\lambda_{2}^{2}+\lambda_{3}^{2}+\lambda_{1}^{2}}{2}\eta^{2}-\frac{\lambda_{3}^{2}(\lambda_{3}^{2}-\lambda_{2}^{2}-\lambda_{1}^{2})-2\lambda_{2}^{2}\lambda_{1}^{2}-(\lambda_{3}^{2}+\lambda_{1}^{2}+\lambda_{2}^{2})(\lambda_{3}^{2}-\lambda_{1}^{2})\mu(m)}{6}}{\eta^{2}-\lambda_{3}^{2}+(\lambda_{3}^{2}-\lambda_{1}^{2})\mu(m)}.\end{split}$
(4.12)

We note that a similar expression for the tracer velocity in a dense soliton
gas was obtained in [32] in the context of the modified KdV (mKdV) equation.

For $N=1$ the integrals (3.4) for the mean and mean square of the soliton
condensate wave field $\varphi\equiv\varphi_{\rm c}^{(1)}(x,t)$ can be
explicitly evaluated using (253.11) and (256.11) from [33] and 19.7.10 from
[49]:

$\displaystyle\langle\varphi\rangle=\lambda_{1}^{2}+\lambda_{2}^{2}-\lambda_{3}^{2}+2(\lambda_{3}^{2}-\lambda_{1}^{2})\mu(m),$
(4.13)
$\displaystyle\langle\varphi^{2}\rangle=\frac{2(\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2})}{3}\langle\varphi\rangle+\frac{\lambda_{1}^{4}+\lambda_{2}^{4}+\lambda_{3}^{4}-2\lambda_{1}^{2}\lambda_{2}^{2}-2\lambda_{2}^{2}\lambda_{3}^{2}-2\lambda_{1}^{2}\lambda_{3}^{2}}{3},$
(4.14)

with $\mu(m)$ and $m$ given by (4.9). It is not difficult to verify that,
unlike in the case of genus 0 condensates, the variance
$\Delta=\sqrt{\langle\varphi^{2}\rangle-\langle\varphi\rangle^{2}}$ does not
vanish identically implying that all realizations of the genus 1 soliton
condensate are almost surely non-constant.

A key observation is, that formulae (4.13), (4.14) coincide with the period
averages $\overline{\varphi}$ and $\overline{\varphi^{2}}$ of the genus 1 KdV
solution associated with the spectral Riemann surface $\mathcal{R}_{2}$ of
(4.5) (see e.g. [34, 35]):

$\begin{split}&\varphi(x,t)\equiv
F_{1}(\theta;\lambda_{1},\lambda_{2},\lambda_{3})=\lambda_{1}^{2}+\lambda_{2}^{2}-\lambda_{3}^{2}+2(\lambda_{3}^{2}-\lambda_{1}^{2}){\rm
dn}^{2}\left(\frac{\sqrt{\lambda_{3}^{2}-\lambda_{1}^{2}}}{k}\theta;m\right),\\\
&\theta=k(x-Ut)+\theta^{0},\quad
U=2(\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2}),\quad
k=\frac{\pi\sqrt{\lambda_{3}^{2}-\lambda_{1}^{2}}}{{\rm K}(m)},\end{split}$
(4.15)

where $\theta^{0}\in[0,2\pi)$ is an arbitrary initial phase.

Figure 4: Comparison between the numerical realization of genus 1 condensate
generated with $200$ solitons (black solid line), $200$ solitons (black solid
line), and the exact cnoidal wave solution $F_{1}(kx)$ (4.15) (red dashed
line) for $\lambda_{1}=0.5,\lambda_{2}=0.85,\lambda_{3}=1$ ($m=0.63$); the two
plots are visually indistinguishable from one another. b) Corresponding
absolute errors $\varphi_{n}(x)-F_{1}(kx)$ obtained with $50$ solitons (dashed
line), $100$ solitons (solid line) and $200$ (dash-dotted line); the absolute
error is evaluated at the extrema of the oscillations.

The equivalence between the ensemble averages in genus 1 KdV soliton
condensates and the period averages in single-phase KdV solutions, along with
the established in Section 3 equivalence between the respective modulation
dynamics, strongly suggest that realizations of the genus 1 soliton
condensates are described by the periodic solutions $F_{1}(\theta)$ (4.15) of
the KdV equation. This motivates the following

######  Conjecture 4.1.

For any realization $\varphi=\varphi_{\rm c}^{(1)}(x,t)$ of the genus 1 KdV
soliton condensate associated with the spectral curve $\mathcal{R}_{2}$ of
(4.5) one can find the initial phase $\theta^{0}\in[0,2\pi)$ in the periodic
solution $F_{1}(\theta;\lambda_{1},\lambda_{2},\lambda_{3})$ (4.15) such that
almost surely $\varphi_{\rm
c}^{(1)}(x,t)=F_{1}(\theta;\lambda_{1},\lambda_{2},\lambda_{3})$.

We support Conjecture 4.1 by a detailed comparison of a numerical realization
of KdV soliton condensate (as $n$-soliton solution with $n$ large) spectrally
configured according to the DOS (4.6), and the periodic KdV solution (4.15),
defined on the same spectral curve $\mathcal{R}_{2}$, with the appropriately
chosen initial phase $\theta^{0}$ (see Appendix A for the details of the
numerical implementation of soliton condensate). The comparison is presented
in Fig. 4 and reveals a remarkable agreement, which further improves as $n$
increases.

Conjecture 4.1 can be naturally generalized to an arbitrary genus $N$: for any
realization of the KdV soliton condensate of genus $N$ corresponding to the
density of states $u^{(N)}(\eta;{\boldsymbol{\lambda}})$ (3.10) and associated
with the spectral Riemann surface $\mathcal{R}_{2N}$ of (3.8), one can find
$N$-component initial phase vector
${\boldsymbol{\theta}}^{0}\in\mathbb{T}^{N}$ so that $\varphi^{(N)}_{\rm
c}(x,t)$ almost surely coincides with $N$-phase KdV solution
$\varphi=F_{N}({\boldsymbol{\theta}};{\boldsymbol{\lambda}})$ (2.5). To
support this generalization we performed a comparison of a numerical
realization of the genus 2 soliton condensate with the respective two-phase
(two-gap) KdV solution, see Appendix B.

A rigorous mathematical proof of Conjecture 4.1 and its generalization for an
arbitrary genus will be the subject of future work.

In conclusion we note that Conjecture 4.1 correlates with the results of [10]
where a particular “deterministic soliton gas” solution of the KdV equation
has been constructed by considering the $n$-soliton solution with the discrete
spectrum confined within two symmetric intervals—the analogs of s-bands of our
work—and letting $n\to\infty$. This solution was shown in [10] to represent a
primitive potential [36] whose long-time asymptotics is described at leading
order by a modulated genus 1 KdV solution. A similar construction was realized
for the mKdV equation in [32].

###  4.2 Modulation dynamics

The dynamics of DOS in non-equilibrium (weakly non-homogeneous) soliton
condensates is determined by the evolution of the endpoints $\lambda_{j}$ of
the spectral bands of $\Gamma$ (the s-bands). As proven in Section 3, this
evolution is governed by the Whitham modulation equations (3.13). Properties
of the KdV-Whitam modulation systems are well studied: in particular, system
(3.13) is strictly hyperbolic and genuinely nonlinear for any genus $N\geq 1$
[28]. This implies inevitability of wavebreaking for a broad class of initial
conditions. What is the meaning of the wavebreaking in the context of soliton
condensates, and how is the solution of the kinetic equation continued beyond
the wavebreaking time?

We first invoke the definitive property of a soliton condensate—the vanishing
of the spectral scaling function, $\sigma(\eta)\equiv 0$ in the soliton gas
NDRs (2.9). According to Remark 3.3, if $\sigma(\eta;x,0)\equiv 0$ for all
$x\in\mathbb{R}$, then $\sigma(\eta;x,t)\equiv 0$ for all $x\in\mathbb{R}$,
$\forall t>0$ implying that soliton condensate necessarily remains a
condensate during the evolution (at least of some class of initial data). The
only qualitative modification that is permissible during the evolution is the
change of the genus $N$. The description of the evolution of a soliton
condensate is then reduced to the determination of the spectral support
$\Gamma(x,t)$, parametrizing the DOS via the band edges $\lambda_{j}(x,t)$:
$u=u^{(N)}(\eta;\lambda_{1},\dots\lambda_{2N+1})$ (3.10).

In view of the above, the evolution of soliton condensates can be naturally
put in the framework of the problem of hydrodynamic evolution of multivalued
functions originally formulated by Dubrovin and Novikov [19]. Let
$\Lambda_{N}(x,t)=\\{\lambda_{1}(x,t),\dots,\lambda_{2N+1}(x,t)\\}$ be a
smooth multivalued curve whose branches $\lambda_{j}(x,t)$ satisfy the Whitham
modulation equations (3.13). Then, if wavebreaking occurs within one of the
branches it results in a change of the genus $N$ so that
$\Lambda_{N}\to\Lambda_{N+1}$ in some space-time region $[x^{-}(t),x^{+}(t)]$
that includes the wavebreaking point. The curves $\Lambda_{N}$ and
$\Lambda_{N+1}$ are glued together at free boundaries $x^{\pm}(t)$. Details of
the implementation of this procedure can be found in [19, 37, 38, 39]. The
simplest case of the multivalued curve evolution arises when the initial data
for $\Lambda_{N}$ is a piecewise-constant distribution (both for
$\lambda_{j}$’s and for $N$), with a discontinuity at $x=0$ — a Riemann
problem. In this special case the wavebreaking occurs at $t=0$ (subject to
appropriate sign of the initial jump) and smoothness of $\Lambda_{N}$ is not a
prerequisite.

In this paper, we restrict ourselves to Riemann problems involving only genus
0 and genus 1 modulation solutions and show how the resulting spectral
dynamics are interpreted in terms of soliton condensates. For that we will
need explicit expressions for the Whitham characteristic velocities for $N=0$
and $N=1$. These expressions are known very well (see e.g. [11, 34, 35]) but
here we obtain them as transport velocities for the respective soliton
condensates, using the expressions (4.2), and (4.12) respectively.

(i) $N=0$. Consider a non-equilibrium (non-uniform) soliton condensate of
genus 0, characterized by a space-time dependent DOS $u(\eta;x,t)$. To this
end we set $\eta=\lambda_{1}(x,t)$ in (4.2), then the Whitham system (3.13),
(3.14) assumes the form of the Hopf (inviscid Burgers) equation

$(\lambda_{1})_{t}+6\lambda_{1}^{2}(\lambda_{1})_{x}=0.$ (4.16)

Note that this is exactly the result obtained by Lax and Levermore [9] for the
pre-breaking evolution of semi-classical soliton ensembles.

(ii) $N=1$. We obtain on using (4.12),

$(\lambda_{j})_{t}+V_{j}(\lambda_{1},\lambda_{2},\lambda_{3})(\lambda_{j})_{x}=0,\quad
j=1,2,3,$ (4.17)

where

$\begin{split}&V_{1}(\lambda_{1},\lambda_{2},\lambda_{3})\equiv
s^{(1)}(\lambda_{1};\lambda_{1},\lambda_{2},\lambda_{3})=2(\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2})+\frac{4(\lambda_{2}^{2}-\lambda_{1}^{2})}{\mu(m)-1},\\\
&V_{2}(\lambda_{1},\lambda_{2},\lambda_{3})\equiv
s^{(1)}(\lambda_{2};\lambda_{1},\lambda_{2},\lambda_{3})=2(\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2})+\frac{4(\lambda_{3}^{2}-\lambda_{2}^{2})(\lambda_{2}^{2}-\lambda_{1}^{2})}{\lambda_{3}^{2}-\lambda_{2}^{2}-(\lambda_{3}^{2}-\lambda_{1}^{2})\mu(m)},\\\
&V_{3}(\lambda_{1},\lambda_{2},\lambda_{3})\equiv
s^{(1)}(\lambda_{3};\lambda_{1},\lambda_{2},\lambda_{3})=2(\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2})+\frac{4(\lambda_{3}^{2}-\lambda_{2}^{2})}{\mu(m)},\end{split}$
(4.18)

and $\mu(m)$ is defined in (4.9). System (4.17), (4.18) coincides with the
original Whitham modulation equations derived for $r_{j}=6\lambda_{j}^{2}$ in
[40] by averaging KdV conservation laws over the single-phase, cnoidal wave
family of solutions (see also [11, 19, 34, 35]).

##  5 Riemann problem for soliton condensates

The classical Riemann problem consists of finding solution to a system of
hyperbolic conservation laws subject to piecewise-constant initial conditions
exhibiting discontinuity at $x=0$. The distribution solution of such Riemann
problem generally represents a combination of constant states, simple
(rarefaction) waves and strong discontinuities (shocks or contact
discontinuities) [41]. In dispersive hydrodynamics, classical shock waves are
replaced by dispersive shock waves (DSWs) — nonlinear expanding wavetrains
with a certain, well-defined structure [35]. Here we generalize the Riemann
problem formulation to the soliton gas kinetic equation by considering (1.1)
subject to discontinuous initial DOS:

$u(\eta,x,t=0)=\begin{cases}u^{(N_{-})}(\eta;\lambda_{1}^{-},\dots,\lambda_{2N_{-}+1}^{-}),&x<0,\\\
u^{(N_{+})}(\eta;\lambda_{1}^{+},\dots,\lambda_{2N_{+}+1}^{+}),&x>0,\end{cases}$
(5.1)

where $u^{(N)}(\eta;\lambda_{1},\dots,\lambda_{2N+1})$ is the DOS (3.10) of
genus $N$ condensate and $\lambda_{j}^{\pm}>0$.

As discussed in Section 4.2, soliton condensate necessarily retains its
definitive property $\sigma=0$ during the evolution, with the only qualitative
modification permissible being the change of the genus $N$. The evolution of
the soliton condensate is then determined by the motion of the s-band edges
$\lambda_{j}$ according to the Whitham modulation equations (3.13) subject to
discontinuous initial conditions following from (5.1):

$\\{N;{\boldsymbol{\lambda}}\\}(x,t=0)=\begin{cases}\\{N_{-};(\lambda_{1}^{-},\dots,\lambda_{2N_{-}+1}^{-})\\},&x<0,\\\
\\{N_{+};(\lambda_{1}^{+},\dots,\lambda_{2N_{+}+1}^{+})\\},&x>0.\end{cases}$
(5.2)

Thus the Riemann problem for soliton gas kinetic equation is effectively
reduced in the condensate limit to the Riemann problem (5.2) for the Whitham
modulation equations (3.13). Depending on the sign of the jump
$\lambda_{j}^{-}-\lambda_{j}^{+}$ the regularization of the discontinuity in
$\lambda_{j}$ can occur in two ways: (i) if
$(\lambda_{j}^{-}-\lambda_{j}^{+})>0$ then the regularization occurs via the
generation of a rarefaction wave for $\lambda_{j}$ without changing the genus
$N$ of the condensate; (ii) if $(\lambda_{j}^{-}-\lambda_{j}^{+})<0$ (which
implies immediate wavebreaking for $\lambda_{j}$) the regularization occurs
via the generation of a higher genus condensate whose evolution is governed by
the modulation equations.

Below we consider several particular cases of Riemann problems describing some
prototypical features of the soliton condensate dynamics.

###  5.1 $N_{-}=N_{+}=0$

Consider the initial condition for the kinetic equation in the form of a
discontinuous genus 0 condensate DOS,

$u(\eta,x,t=0)=\begin{cases}u^{(0)}(\eta;q_{-}),&x<0,\\\
u^{(0)}(\eta;q_{+}),&x>0,\end{cases}$ (5.3)

where $q_{\pm}=\lambda_{1}^{\pm}$, and $u^{(0)}>0$ is defined in (4.1). The
DOS distribution (5.3) implies the step initial conditions for the Whitham
modulation system (3.13):

$N(x,t=0)=0,\quad\lambda_{1}(x,t=0)=\begin{cases}q_{-},&x<0,\\\
q_{+},&x>0,\end{cases}$ (5.4)

with $q_{-}\neq q_{+}$. Additionally, since the wave field in a genus 0
soliton condensate is almost surely a constant,
$\varphi(x,t)=(\lambda_{1})^{2}$, we conclude that the DOS distribution (5.3)
gives rise to the Riemann step data

$\varphi(x,t=0)=\begin{cases}q_{-}^{2},&x<0,\\\ q_{+}^{2},&x>0,\end{cases}$
(5.5)

for the KdV equation (2.1) itself.

The Riemann problem for the KdV equation was originally studied by Gurevich
and Pitaevskii (GP) [11] in the context of the description of dispersive shock
waves. The key idea of GP construction was to replace the dispersive Riemann
problem (5.5) for the KdV equation by an appropriate boundary value problem
for the hyperbolic KdV-Whitham system (4.17) which is then solved in the class
of $x/t$-self-similar solutions. Here we take advantage of the GP modulation
solutions and their higher genus analogues to describe dynamics of soliton
condensates. The choice of the genus of the Whitham system and,
correspondingly, the genus of the associated soliton condensate, depends on
whether $q_{-}>q_{+}$ or $q_{+}<q_{-}$.

####  5.1.1 Rarefaction wave ($q_{-}<q_{+}$)

The solution of the Riemann problem (1.1), (5.3) is given globally (for $t>0$)
by the genus 0 DOS $u^{(0)}(\eta;\lambda_{1})$ (4.1) modulated by the centered
rarefaction wave solution of the Hopf equation (4.16) subject to the step
initial condition (5.4):

$\lambda_{1}(x,t)=\begin{cases}q_{-},&x<s_{-}t,\\\
\sqrt{\frac{x}{6t}},&s_{-}t<x<s_{+}t,\\\ q_{+},&s_{+}t<x,\end{cases}$ (5.6)

where

$s_{-}=6q_{-}^{2}t,\quad s_{+}=6q_{+}^{2}t.$ (5.7)

Note that the solution (5.6) is admissible since $s_{-}<s_{+}$. Behavior of
$\lambda_{1}$ in the solution (5.6) is shown in Fig. 5a. The evolution of the
soliton condensate’s DOS associated with the spectral rarefaction wave
solution (5.6) is given by

$u(\eta;x,t)=\frac{\eta}{\pi\sqrt{\lambda_{1}^{2}(x,t)-\eta^{2}}}.$ (5.8)

A contour plot of the DOS (5.8) is presented in Fig. 5(a).

Figure 5:  Solutions to the soliton condensate Riemann problem (5.3). a)
Rarefaction wave (genus $0$) solution (5.6), (5.8) for $q_{-}<q_{+}$. Dashed
line: $\lambda_{1}(x,t)$, colors: DOS $u^{(0)}(\eta;\lambda_{1}(x,t))$. b) DSW
(genus $1$) solution (5.9), (5.10) for $q_{-}>q_{+}$. Dashed line:
$\lambda_{1}(x,t)\leq\lambda_{2}(x,t)\leq\lambda_{3}(x,t)$, colors: DOS
$u^{(1)}(\eta;q_{+},\lambda_{2}(x,t),q_{-})$.

####  5.1.2 Dispersive shock wave ($q_{-}>q_{+}$)

The solution (5.6), (5.8) derived previously is not admissible for
$q_{-}<q_{+}$ since $s_{-}>s_{+}$ in that case. In other words, the
compressive discontinuous initial data (5.4) imply immediate wavebreaking and
necessitate the introduction of the higher genus DOS connecting
$u^{(0)}(\eta;q_{-})$ and $u^{(0)}(\eta;q_{+})$. The requisite DOS is given by
equation (4.6), which we reproduce here for convenience,

$u(\eta;x,t)=u^{(1)}(\eta;\lambda_{1},\lambda_{2},\lambda_{3})=\frac{i\eta(\eta^{2}-w^{2})}{\pi
R(\eta)},$ (5.9)

Here $w(\lambda_{1},\lambda_{2},\lambda_{3})$ is given by (4.9) and
$\lambda_{j}=\lambda_{j}(x,t)$, $j=1,2,3$, are slowly modulated according to
the Whitham equations (4.17), (4.18).

The solution of (4.18) is self-similar, $\lambda_{j}(x/t)$, such that
$u^{(1)}(\eta;\lambda_{1},\lambda_{2},\lambda_{3})$ matches with
$u^{(0)}(\eta;q_{-})$ at the left boundary $x=s_{-}t$, and with
$u^{(0)}(\eta;q_{+})$ at the right boundary $x=s_{+}t$, with $s_{-}<s_{+}$.

The requisite solution is the $2$-wave of the Whitham system (4.17) (only
$\lambda_{2}$ is non-constant)

$\lambda_{1}=q_{+},\quad
V_{2}(\lambda_{1},\lambda_{2},\lambda_{3})=x/t,\quad\lambda_{3}=q_{-},\quad\text{for}\quad
s_{-}t<x<s_{+}t,$ (5.10)

where

$s_{-}=V_{2}(q_{+},q_{+},q_{-})=12q_{+}^{2}-6q_{-}^{2},\quad
s_{+}=V_{2}(q_{+},q_{-},q_{-})=2q_{+}^{2}+4q_{-}^{2}.$ (5.11)

This is the famous GP solution describing the DSW modulations in the KdV step
resolution problem [11]. Indeed, we have $s_{-}<s_{+}$ and, interpreting the
GP solution (5.10) in terms of soliton condensates the limiting behaviors at
the DSW edges is given by

$\begin{split}&x\to s_{-}t,\quad\lambda_{2}\to\lambda_{1}=q_{+},\quad
u^{(1)}(\eta;q_{+},\lambda_{2},q_{-})\to u^{(0)}(\eta;q_{-}),\\\ &x\to
s_{+}t,\quad\lambda_{2}\to\lambda_{3}=q_{-},\quad
u^{(1)}(\eta;q_{-},\lambda_{2},q_{+})\to u^{(0)}(\eta;q_{+}).\end{split}$
(5.12)

###  5.2  $N_{-}+N_{+}=1$

Before considering the soliton condensate Riemann problem (1.1), (5.1) for the
case $N_{-}+N_{+}=1$ we list the admissible solutions to the kinetic equation
connecting a genus $0$ distribution $u^{(0)}(\eta;q)$ to a genus $1$
distribution $u^{(1)}(\eta;\lambda_{1},\lambda_{2},\lambda_{3})$. One can
easily verify for the next four solutions that

$\begin{split}&x\to s_{-}t,\quad
u^{(1)}(\eta;\lambda_{1},\lambda_{2},\lambda_{3})\to
u^{(N_{-})}(\eta;\boldsymbol{\lambda}_{-}),\\\ &x\to s_{+}t,\quad
u^{(1)}(\eta;\lambda_{1},\lambda_{2},\lambda_{3})\to
u^{(N_{+})}(\eta;\boldsymbol{\lambda}_{+}),\end{split}$ (5.13)

with $s_{-}<s_{+}$.

We use the following convention to label the fundamental Riemann problem
solutions: we call $j^{\pm}$-wave, where $j$ is the index of the only varying
Riemann invariant $\lambda_{j}$ in the solution, while the remaining
invariants are constant; $+$ indicates that $N_{+}=1$ i.e. the genus $1$
soliton condensate is initially at $x>0$, and $-$ indicates that $N_{-}=1$
i.e. the genus $1$ soliton condensate is initially at $x<0$.

(i) $3^{+}$-wave

Consider the initial condition for the soliton condensate DOS:

$u(\eta,x,t=0)=\begin{cases}u^{(0)}(\eta;q_{-}),&x<0,\\\
u^{(1)}(\eta;\lambda_{1}^{+},\lambda_{2}^{+},\lambda_{3}^{+}),&x>0,\end{cases}\quad\text{with}\quad\lambda_{1}^{+}=q_{-}.$
(5.14)

The resolution of the step (5.14) is described by

$u(\eta,x,t)=\begin{cases}u^{(0)}(\eta;q_{-}),&x<s_{-}t,\\\
u^{(1)}(\eta;q_{-},\lambda_{2}^{+},\lambda_{3}(x/t)),&s_{-}t<x<s_{+}t,\\\
u^{(1)}(\eta;q_{-},\lambda_{2}^{+},\lambda_{3}^{+}),&x>s_{+}t,\end{cases}$
(5.15)

where $\lambda_{3}(x/t)$ is given by the $3^{+}$-wave solution of the
modulation equations (4.17):

$\begin{split}&\lambda_{1}=q_{-},\quad\lambda_{2}=\lambda_{2}^{+},\quad
V_{3}(\lambda_{1},\lambda_{2},\lambda_{3})=x/t,\quad\text{for}\quad
s_{-}t<x<s_{+}t,\\\
&s_{-}=V_{3}(q_{-},\lambda_{2}^{+},\lambda_{2}^{+})=2(q_{-})^{2}+4(\lambda_{2}^{+})^{2},\quad
s_{+}=V_{3}(q_{-},\lambda_{2}^{+},\lambda_{3}^{+}).\end{split}$ (5.16)

The behavior of the Riemann invariants $\lambda_{j}$ in the $3^{+}$-wave is
shown in Fig. 6a. The associated soliton condensate KdV solution
$\varphi(x,t)$ along with the behavior of the mean $\langle\varphi\rangle$ are
shown in Figs. 10 and 11.

(ii) $2^{+}$-wave

Consider the initial condition:

$u(\eta,x,t=0)=\begin{cases}u^{(0)}(\eta;q_{-}),&x<0,\\\
u^{(1)}(\eta;\lambda_{1}^{+},\lambda_{2}^{+},\lambda_{3}^{+}),&x>0,\end{cases}\quad\text{with}\quad\lambda_{3}^{+}=q_{-}.$
(5.17)

The resolution of the step (5.17) is described by

$u(\eta,x,t)=\begin{cases}u^{(0)}(\eta;q_{-}),&x<s_{-}t,\\\
u^{(1)}(\eta;\lambda_{1}^{+},\lambda_{2}(x/t),q_{-}),&s_{-}t<x<s_{+}t,\\\
u^{(1)}(\eta;\lambda_{1}^{+},\lambda_{2}^{+},q_{-}),&x>s_{+}t,\end{cases}$
(5.18)

where $\lambda_{2}(x/t)$ is given by the $2^{+}$-wave solution of the
modulation equations (4.17):

$\begin{split}&\lambda_{1}=\lambda_{1}^{+},\quad
V_{2}(\lambda_{1},\lambda_{2},\lambda_{3})=x/t,\quad\lambda_{3}=q_{-},\quad\text{for}\quad
s_{-}t<x<s_{+}t,\\\
&s_{-}=V_{2}(\lambda_{1}^{+},\lambda_{1}^{+},\lambda_{2}^{+})=12(\lambda_{1}^{+})^{2}-6(q_{-})^{2},\quad
s_{+}=V_{2}(\lambda_{1}^{+},\lambda_{2}^{+},q_{-}).\end{split}$ (5.19)

The behavior of the Riemann invariants $\lambda_{j}$ in the $2^{+}$-wave is
shown in Fig. 6b.

(iii) $1^{-}$-wave

Consider the initial condition:

$u(\eta,x,t=0)=\begin{cases}u^{(1)}(\eta;\lambda_{1}^{-},\lambda_{2}^{-},\lambda_{3}^{-}),&x>0,\\\
u^{(0)}(\eta;q_{+}),&x<0,\end{cases}\quad\text{with}\quad\lambda_{3}^{-}=q_{+}.$
(5.20)

The resolution of the step (5.20) is described by

$u(\eta,x,t)=\begin{cases}u^{(1)}(\eta;\lambda_{1}^{-},\lambda_{2}^{-},q_{+}),&x<s_{-}t,\\\
u^{(1)}(\eta;\lambda_{1}(x/t),\lambda_{2}^{-},q_{+}),&s_{-}t<x<s_{+}t,\\\
u^{(0)}(\eta;q_{+}),&x>s_{+}t,\end{cases}$ (5.21)

where $\lambda_{1}(x/t)$ is given by the $1^{-}$-wave solution of the
modulation equations (4.17):

$\begin{split}&V_{1}(\lambda_{1},\lambda_{2},\lambda_{3})=x/t,\quad\lambda_{2}=\lambda_{2}^{-},\quad\lambda_{3}=q_{+},\quad\text{for}\quad
s_{-}t<x<s_{+}t,\\\ &s_{-}=V_{1}(\lambda_{1}^{-},\lambda_{2}^{-},q_{+}),\quad
s_{+}=V_{1}(\lambda_{2}^{-},\lambda_{2}^{-},q_{+})=12(\lambda_{2}^{-})^{2}-6(q_{+})^{2}.\end{split}$
(5.22)

The behavior of the Riemann invariants $\lambda_{j}$ in the $1^{-}$-wave is
shown in Fig. 6c.

(iv) $2^{-}$-wave

Consider the initial condition:

$u(\eta,x,t=0)=\begin{cases}u^{(1)}(\eta;\lambda_{1}^{-},\lambda_{2}^{-},\lambda_{3}^{-}),&x<0,\\\
u^{(0)}(\eta;q_{+}),&x>0,\end{cases}\quad\text{with}\quad\lambda_{1}^{-}=q_{+}.$
(5.23)

The resolution of the step (5.23) is described by

$u(\eta,x,t)=\begin{cases}u^{(1)}(\eta;q_{+},\lambda_{2}^{-},\lambda_{3}^{-}),&x<s_{-}t,\\\
u^{(1)}(\eta;q_{+},\lambda_{2}(x/t),\lambda_{3}^{-}),&s_{-}t<x<s_{+}t,\\\
u^{(0)}(\eta;q_{+}),&x>s_{+}t,\end{cases}$ (5.24)

where $\lambda_{2}(x/t)$ is given by the $2^{-}$-wave solution of the
modulation equations (4.17):

$\begin{split}&\lambda_{1}=q_{+},\quad
V_{2}(\lambda_{1},\lambda_{2},\lambda_{3})=x/t,\quad\lambda_{3}=\lambda_{3}^{-},\quad\text{for}\quad
s_{-}t<x<s_{+}t,\\\ &s_{-}=V_{2}(q_{+},\lambda_{2}^{-},\lambda_{3}^{-}),\quad
s_{+}=V_{2}(q_{+},\lambda_{3}^{-},\lambda_{3}^{-})=2(q_{+})^{2}+4(\lambda_{3}^{-})^{2}.\end{split}$
(5.25)

The behavior of the Riemann invariants $\lambda_{j}$ in the $2^{-}$-wave is
shown in Fig. 6d. The associated soliton condensate KdV solution
$\varphi(x,t)$ along with the behavior of the mean $\langle\varphi\rangle$ are
shown in Figs. 12 and 13.

Figure 6: Basic modulation configurations in the Riemann problem (1.1), (5.1)
for soliton condensates with $N_{-}+N_{+}=1$. a) $3^{+}$-wave solution
(LABEL:eq:3-wave). b) $2^{+}$-wave solution (5.19). c) $1^{-}$-wave solution
(5.22). d) $2^{-}$-wave solution (5.25). In all cases the dashed lines show
the variation of the spectral edges
$\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$, and the colors visualize the DOS
$u^{(1)}(\eta;{\boldsymbol{\lambda}})$.

##  6 Riemann problem: numerical results

We consider Riemann problems with $N_{-}+N_{+}\leq 1$. Because of the inherent
limitations of the numerical implementation of soliton gas detailed in
Appendix A, we restrict the comparison to the cases $q_{-}=0$ or $q_{+}=0$.

###  6.1 Rarefaction wave

In this first example, we choose

$\\{N;{\boldsymbol{\lambda}}\\}(x,t=0)=\begin{cases}\\{0;q_{-}=0\\},&x<0,\\\
\\{0;q_{+}=1\\},&x>0.\end{cases}$ (6.1)

A numerical realization of the soliton condensate evolution corresponding to
the steplike initial condition (6.1) is displayed in Fig. 7. The same figure
displays the realization at $t=40$. The realization corresponds to a
$n$-soliton solution with parameters distributed according to the initial DOS
of (5.3), (6.1); details are given in Appendix A. As predicted in Sec. 4.1,
the realization of the condensate corresponds to the vacuum $\varphi=0$ at the
left of $x=0$, and a constant $\varphi=1$ at the right of $x=0$. As
highlighted in Appendix A.1, the $n$-soliton solution displays an overshoot at
$x=0$, regardless of the number of solitons $n$, which is reminiscent of
Gibbs’ phenomenon in the theory of Fourier series. This phenomenon has been
originally observed in the numerical approximation of the soliton condensate
of the focusing NLS equation by a $n$-soliton solution in [50]; see for
instance the similarities between Figs. 7a, 8b and Fig. 2a of [50]. Indeed, in
both cases, the IST spectrum of the step distribution contains a non-solitonic
radiative component (cf. [42]), which is not taken into account by the
$n$-soliton solution; the mismatch between the exact step and the $n$-soliton
solution manifests by the occurrence of the spurious oscillations observed
near $x=0$.

Figure 7: Riemann problem with initial condition (6.1) for DOS $u(\eta;x,t)$.
The plots depict the variation of a condensate’s realization $\varphi(x,t)$ at
$t=0$ (a) and $t=40$ (b, solid line). The red dashed line in depicts the
variation of the rarefaction wave $\varphi=\lambda_{1}(x/t)^{2}$ (5.6).

The solution of the Riemann problem with the initial condition (6.1) is given
by $u^{(0)}(\eta;\lambda_{1}(x,t))$ where $\lambda_{1}(x,t)$ is the
rarefaction wave (genus 0) solution (5.6). We have shown in Sec. 4.1 that the
genus 0 soliton condensate is almost surely described by the constant solution
$\varphi=(\lambda_{1})^{2}$. In the context of the evolution of the step (6.1)
$\lambda_{1}$ varies according to (5.6) so $\lambda(x,t)$ should be treated as
a slowly varying (locally constant) condensate solution. In Fig. 7 we compare
the numerical realization of the evolution of genus $0$ condensate with the
analytical solution (5.6).

###  6.2 Dispersive shock wave

We now consider

$\\{N;{\boldsymbol{\lambda}}\\}(x,t=0)=\begin{cases}\\{0;q_{-}=1\\},&x<0,\\\
\\{0;q_{+}=0\\},&x>0.\end{cases}$ (6.2)

A numerical realization of the genus $0$ soliton condensate corresponding to
the step-initial condition (6.2) is presented in Fig. 8 (a): it corresponds to
the vacuum $\varphi=0$ for $x>0$, and a constant $\varphi=1$ for $x<0$. The
realization at $t=40$ is shown in Fig. 8 (b) and it corresponds to a classical
DSW solution for the KdV equation.

Figure 8: Riemann problem with initial condition (6.2) for DOS $u(\eta;x,t)$.
The plots depict the variation of a condensate’s realization $\varphi(x,t)$ at
$t=0$ (a) and $t=40$ (b).

The solution of the condensate Riemann problem with the initial condition
(5.3), (6.2) is given by the genus 1 DOS (5.9) modulated by the $2$-wave
solution (5.10) of the Whitham equations. In order to make a quantitative
comparison of this analytical solution with the numerical evolution of the
soliton gas displayed in Fig. 8, we compute numerically the mean
$\langle\varphi\rangle$ and the variance
$\sqrt{\langle\varphi^{2}\rangle-\langle\varphi\rangle^{2}}$, the latter being
an amplitude type characteristic of the cnoidal wave. We have conjectured in
Sec. 4.1 that any realization of the uniform genus $1$ condensate corresponds
to a cnoidal wave modulo the initial phase $\theta^{0}\in[0;2\pi)$. In that
case, the ensemble average of the soliton condensate reduces to an average
over the phase $\theta^{0}$, or equivalently, over the period of the cnoidal
wave, which can be performed on a single realization. We assume here that the
result generalizes to non-uniform condensates so that the realization computed
numerically and displayed in Fig. 8(b) can be consistently compared with a
slowly modulated cnoidal wave solution. The averages
$\langle\varphi(x,t)\rangle$ and $\langle\varphi(x,t)^{2}\rangle$ can be
determined via a local phase average of one realization of the condensate. The
local period averages are obtained via

$\langle\varphi(x,t)\rangle=\frac{1}{L(x,t)}\int_{x}^{x+L(x,t)}\varphi(y,t)dy,\quad\langle\varphi(x,t)^{2}\rangle=\frac{1}{L(x,t)}\int_{x}^{x+L(x,t)}\varphi(y,t)^{2}dy,$
(6.3)

where $L(x,t)$ is the local wavelength extracted numerically.

The comparison between the analytically determined averages
(4.13),(4.14),(5.10) and the averages (6.3) obtained numerically is presented
in Fig. 11 and shows a very good agreement.

Figure 9: Mean $\langle\varphi\rangle$ (a) and variance
$\sqrt{\langle\varphi^{2}\rangle-\langle\varphi\rangle^{2}}$ (b) of the
solution of the Riemann problem’s solution with the initial condition (6.2).
The markers correspond to averages extracted from the numerical solution using
(6.3), and the solid black lines to the corresponding analytical averages
(4.13),(4.14),(5.10).

###  6.3 Generalized rarefaction wave

$N_{-}+N_{+}=0$ in the two previous examples. In the next examples, we choose
$N_{-}+N_{+}=1$. Let’s start with $N_{+}=1$:

$\\{N;{\boldsymbol{\lambda}}\\}(x,t=0)=\begin{cases}\\{0;q_{-}=0\\},&x<0,\\\
\\{1;(\lambda_{1}^{+}=0,\lambda_{2}^{+}=1/2,\lambda_{2}^{+}=1)\\},&x>0.\end{cases}$
(6.4)

A numerical realization of the step-initial condition is displayed in Fig. 10.
The same figure displays the realization at $t=40$. The realization of the
condensate corresponds to the “vacuum” $\varphi=0$ for $x<0$, and a cnoidal
wave for $x>0$. Note that the KdV equation does not admit heteroclinic
traveling wave solutions, rendering difficult the numerical implementation of
these “generalized” Riemann problems studied for instance in [43, 44].
Remarkably here, the solution depicted in Fig. 10 is an exact, $n$-soliton
solution of the KdV equation. As highlighted previously (see also Appendix
A.1), the $n$-soliton solution exhibits an overshoot at $x=0$, regardless of
the number of solitons $n$.

Figure 10: Riemann problem for soliton condensate with initial condition (6.4)
for DOS $u(\eta;x,t)$. The plots depict the variation of a condensate’s
realization $\varphi(x,t)$ at $t=0$ (a) and $t=40$ (b).

The solution of the Riemann problem for the kinetic equation with the initial
condition (5.14), (6.4) is given by the $3^{+}$-wave (5.15),(LABEL:eq:3-wave).
The comparison between the analytical averages (4.13),(4.14),(LABEL:eq:3-wave)
and the averages obtained numerically is shown in Fig. 11 and shows a very
good agreement. The modulation depicted in Figs. 10b and 11a resembles the
modulation of a cnoidal wave of an almost constant amplitude but with a
varying mean. The variation of the mean $\langle\varphi\rangle$ is similar to
the variation of the field in a classical rarefaction wave, so we call the
corresponding structure shown in Fig. 10b a generalized rarefaction wave. The
variance of the wavefield $\varphi$ in the generalized rarefaction wave is
shown in Fig. 11b.

Figure 11: Mean $\langle\varphi\rangle$ (a) and variance
$\sqrt{\langle\varphi^{2}\rangle-\langle\varphi\rangle^{2}}$ (b) of the
solution of the Riemann problem’s solution with the initial condition (6.4).
The markers correspond to the averages extracted from the numerical solution
using (6.3), and the solid black lines to the corresponding analytical
averages (4.13),(4.14),(LABEL:eq:3-wave).

###  6.4 Generalized dispersive shock wave

We now consider the “complementary” initial condition

$\\{N;{\boldsymbol{\lambda}}\\}(x,t=0)=\begin{cases}\\{1;(\lambda_{1}^{-}=0,\lambda_{2}^{-}=1/2,\lambda_{2}^{-}=1)\\},&x<0,\\\
\\{0;q_{+}=0\\},&x>0.\end{cases}$ (6.5)

An example of the numerical realization of the soliton gas step-initial
condition and its evolution at $t=40$ are displayed in Fig. 12.

Figure 12: Riemann problem with initial condition (6.5) for DOS $u(\eta;x,t)$.
The plots depict the variation of a condensate’s realization $\varphi(x,t)$ at
$t=0$ (a) and $t=40$ (b).

The solution of the Riemann problem with the initial condition (6.5) is given
by the $2^{-}$-wave (5.24), (5.25). The comparison between the analytically
derived averages (4.13),(4.14),(5.25) and the averages obtained numerically is
displayed in Fig. 13, and shows a very good agreement. The modulation observed
in Figs. 12, 13 resembles the modulation of partial dispersive shock wave: the
modulated cnoidal wave reaches the soliton limit $m=1$ for $x\to s_{+}t$ but
terminates at $m\neq 0$ for $x\to s_{-}t$. The solution then continues as a
non-modulated cnoidal wave for $x<s_{-}t$. This structure differs from the
celebrated dispersive shock wave solution of the KdV equation involving the
entire range $0\leq m\leq 1$ [35]. We call the described structure connecting
a constant state (a genus 0 condensate) at $x\to+\infty$ with a periodic
solution (a genus 1 condensate) at $x\to-\infty$ a generalized DSW. We note
that the soliton condensate structure shown in Fig. 12b exhibits strong
similarity to the “deterministic KdV soliton gas” solution constructed in
[10].

Figure 13: Mean $\langle\varphi\rangle$ (a) and variance
$\sqrt{\langle\varphi^{2}\rangle-\langle\varphi\rangle^{2}}$ (b) of the
solution of the Riemann problem’s solution with the initial condition (6.5).
The markers correspond to averages extracted from the numerical solution using
(6.3), and the solid black lines to the corresponding analytical averages
(4.13),(4.14),(LABEL:eq:3-wave).

##  7 Diluted soliton condensates

###  7.1 Equilibrium properties

We now introduce the notion of a “diluted” soliton condensate by considering
DOS $u(\eta)=Cu^{(N)}(\eta)$, where $u^{(N)}(\eta)$ is the condensate DOS of
genus $N$, and $0<C<1$ is the “dilution constant”.

E.g. the diluted soliton condensate of genus 0 is characterized by DOS

$u(\eta)=C\frac{\eta}{\pi\sqrt{\lambda_{1}^{2}-\eta^{2}}},\quad 0<C<1.$ (7.1)

We recover the genus $0$ condensate DOS (4.1) by setting $C=1$. As $C$
decreases, the “averaged spacing” between the solitons

$\kappa^{-1}=\left(\int u(\eta)d\eta\right)^{-1}\propto C^{-1}$ (7.2)

increases and the condensate gets “diluted”. Comparison between the most
probable realization of the condensate ($C=1$) and a typical realization of a
slightly dilute condensate ($C=0.97$) is given in Fig. 14. Remarkably, one can
see that a slight increase of the average spacing between the solitons within
the condensate results in the emergence of significant random oscillations of
the KdV wave field.

As follows from (2.14) we have
$\langle\varphi\rangle=\langle\varphi^{2}\rangle=C$ for the diluted genus 0
condensate so that the variance is given by:

$\Delta=\sqrt{\langle\varphi^{2}\rangle-\langle\varphi\rangle^{2}}=\sqrt{C(1-C)}.$
(7.3)

The comparison between (7.3) and the variance obtained numerically by
averaging over different diluted condensates is presented in Figure 14.
Assuming ergodicity of a generic uniform soliton gas, the ensemble average
$\langle\dots\rangle$ in Fig. 14a (and Fig. 15) is computed here numerically
with a spatial average of one, spatially broad, gas realization.

Figure 14: a) Realizations soliton gas with the DOS (7.1) and $\lambda_{1}=1$:
$C=1$ in dashed line (genus 0 condensate) vs $C=0.97$ in solid line (diluted
genus zero condensate); in both cases the gas is realized numerically with
$N=100$ solitons. b) Variance for diluted condensates $C<1$. Solid line:
formula (7.3); markers: numerically extracted values of the variance; insets:
typical realizations of the KdV wave field $\varphi(x,t)$ in diluted
condensates.

More generally, the diluted soliton condensate of genus $N$ is characterized
by DOS

$u(\eta)=Cu^{(N)}(\eta;\lambda_{1},\dots,\lambda_{2N+1}),\quad 0<C<1.$ (7.4)

We have in the general case

$\langle\varphi\rangle=C\langle\varphi_{\rm
c}^{(N)}\rangle,\quad\langle\varphi^{2}\rangle=C\langle\big{(}\varphi_{\rm
c}^{(N)}\big{)}^{2}\rangle,$ (7.5)

where $\langle\varphi_{\rm c}^{(N)}\rangle$, $\langle\big{(}\varphi_{\rm
c}^{(N)}\big{)}^{2}\rangle$ are the ensemble averages obtained for the genuine
condensate ($C=1$), functions of $\lambda_{1},\dots,\lambda_{2N+1}$ only; for
instance $\langle\varphi_{\rm c}^{(1)}\rangle$, $\langle\big{(}\varphi_{\rm
c}^{(1)}\big{)}^{2}\rangle$ are given by (4.13),(4.14). Since
$\langle\big{(}\varphi_{\rm c}^{(N)}\big{)}^{2}\rangle\neq\langle\varphi_{\rm
c}^{(N)}\rangle^{2}$ for $N\geq 1$ and distinct $\lambda_{i}$’s, the variance
of diluted genus condensates

$\Delta=\sqrt{C\langle\big{(}\varphi_{\rm
c}^{(N)}\big{)}^{2}\rangle-C^{2}\langle\varphi_{\rm c}^{(N)}\rangle^{2}}$
(7.6)

never vanishes if $N\geq 1$, as can be seen in the example $N=1$ shown in Fig.
15. Thus, in contrast with the genus 0 case, the transition from the genus 1
condensate ($C=1$) to diluted genus 1 condensate ($C<1$) does not see a
drastic change in the oscillations’ amplitude. In particular, the oscillations
seem to remain “almost” coherent – i.e. an average period can be identified –
for the dilution factors $C$ close to 1 as depicted in the inset of Fig. 15.

Figure 15: Variance for diluted genus 1 condensates with DOS (7.4) and
$(\lambda_{1},\lambda_{2},\lambda_{3})=(0.5,0.85,1)$; insets: typical
realizations of the KdV wave field $\varphi(x,t)$ in diluted condensates.
Figure 16: Soliton trajectories in a diluted genus 0 soliton condensate with
$C=0.9$. Highlighted is a small-amplitude tracer soliton moving backwards.

Diluted condensates present a convenient framework to verify the prediction
formulated in Remark 2.1 regarding the “backflow” effect (i.e. the existence
of tracer KdV solitons moving in negative direction) in sufficiently dense
soliton gases. A numerical simulation of the diluted genus $0$ condensate with
$C=0.9$ where one can clearly see the soliton trajectory with a negative slope
is presented in Fig. 16.

###  7.2 Riemann problem

We can now consider the soliton condensate Riemann problem for diluted
condensates for which the initial DOS (5.1) is replaced by

$u(\eta,x,t=0)=\begin{cases}C_{-}\,u^{(N_{-})}(\eta;\lambda_{1}^{-},\dots,\lambda_{2N_{-}+1}^{-}),&x<0,\\\
C_{+}\,u^{(N_{+})}(\eta;\lambda_{1}^{+},\dots,\lambda_{2N_{+}+1}^{+}),&x>0,\end{cases}$
(7.7)

where $0<C_{\pm}<1$.

To be specific, we investigate numerically the evolution of the diluted
condensate initial conditions (7.7) with $N_{-}+N_{+}\leq 1$ and $\lambda_{i}$
chosen from the examples presented in Sec. 6. Numerical realizations of the
step-initial condition and their evolution in time are presented in Fig. 17.
One can see that generally, realizations of the diluted soliton condensate do
not exhibit a macroscopically coherent structure as observed in Sec. 6.
However, in the case $N_{-}+N_{+}=1$, the evolution of the diluted condensate
realizations, despite the visible incoherence, still qualitatively resembles
the evolution of the “genuine” condensates depicted in Figs. 10, 12. One can
see that the recognizable patterns of the generalized rarefaction wave (see
Fig. 17f) and the generalized DSW (see Fig. 17h) persist even if $C<1$.
Indeed, as shown in Sec. 7.1, the oscillations in a realization of the diluted
genus 1 condensate appear almost coherent for a small dilution factor. The
persistence of coherence can also be observed in the case $N_{-}+N_{+}=0$ when
$\lambda_{1}^{-}>\lambda_{1}^{+}$ (Fig. 17d): a DSW develops if $C=1$, and
coherent, finite amplitude oscillations still develop for $C\neq 1$ at the
right edge of the structure where the amplitudes of oscillations are large. In
connection with the above, it is important to note that, although the initial
condition (7.7) is given by the discontinuous diluted condensate DOS,
$u(\eta;x,0)=Cu^{(N)}(\eta)$, the kinetic equation evolution does not imply
that the DOS will remain to be of the same form for $t>0$. In other words,
unlike genuine condensates, the diluted condensates do not retain the spectral
“diluted condensate” property during the evolution.

Figure 17: Riemann problem for diluted soliton condensates with initial
condition (7.7) with $C_{\pm}=0.95$. (a)-(d) $N_{-}+N_{+}=0$ and
$\lambda_{1}=1$; (e)-(h) $N_{-}+N_{+}=1$ and
$(\lambda_{1},\lambda_{2},\lambda_{3})=(0,1/2,1)$. The diluted condensates are
realized with exact $n$-soliton solutions ($n=200$) configured spectrally
according to the respective scaled condensate DOS’s. The evolution results in
the generation of incoherent rarefaction and dispersive shock waves.

##  8 Conclusions and Outlook

We have considered a special kind of soliton gases for the KdV equation,
termed soliton condensates, which are defined by the property of vanishing the
spectral scaling function $\sigma(\eta)$ in the soliton gas nonlinear
dispersion relations (2.9), (2.10). As a result, the density of states
$u(\eta)$ in a soliton condensate is uniquely determined by its spectral
support $\Gamma^{+}\in\mathbb{R}^{+}$. By considering $\Gamma^{+}$ to be a
union of $N+1$ disjoint intervals,
$[0,\lambda_{1}]\cup[\lambda_{2},\lambda_{3}]\cup\dots\cup[\lambda_{2N},\lambda_{2N+1}]$,
and allowing the endpoints $\\{\lambda_{j}\\}_{j=1}^{2N+1}$ vary slowly in
space-time we prove that the kinetic equation for soliton gas reduces in the
condensate limit to the genus $N$ KdV-Whitham modulation for
$\lambda_{j}(x,t)$. The KdV-Whitham equations were originally derived via the
wave averaging procedure in [40], [4] and via the semicalssical limit of the
KdV equation in [9]. These equations have been extensively used for the
description of dispersive shock waves [35], particularly in the context of
dispersive Riemann problem originally introduced by Gurevich and Pitaevskii
[11].

Along with the characterization of the large scale, modulation dynamics of
soliton condensates, our work suggests that they represent “coherent” or
“deterministic” soliton gases whose typical realizations are given by finite-
gap potentials. We prove this conjecture for genus zero condensates and
present a strong numerical evidence for $N=1,2$.

By invoking the results from the modulation theory of dispersive shock waves
we have constructed analytical solutions to several Riemann problems for the
soliton gas kinetic equation subject to discontinuous condensate initial data.
These solutions describe the evolution of generalized rarefaction and
dispersive shock waves in soliton condensates. We performed numerical
simulations of the Riemann problem for the KdV soliton condensates by
constructing an exact $n$-soliton solutions with $n$ large and the spectral
parameters distributed according to the condensate densty of states. A
comparison of the numerical simulations with analytical predictions from the
solutions of the kinetic equation showed excellent agreement.

Finally, we considered the basic properties of “diluted” soliton condensates
having a scaled condensate density of states and exhibiting rich incoherent
behaviors.

There are several avenues for future work suggested by our results. One
pertinent problem would be to consider near-condensate soliton gas dynamics by
assuming the spectral scaling function $\sigma$ to be “small”
($\sigma(\eta)=\epsilon\tilde{\sigma}(\eta)$, $\epsilon\ll 1$,
$\tilde{\sigma}=\mathcal{O}(1)$). Another area of major interest is the
extension of the developed KdV soliton condensate theory to other integrable
equations, particularly the focusing NLS equation, where a number of important
theoretical and experimental results on the soliton gas dynamics have been
obtained recently, see [55, 7, 60, 50]. And last but not least, one of the
most intriguing open questions is the possibility of phase transitions in
soliton gases, i.e. the formation of a soliton condensate from non-condensate
initial data. The generalized hydrodynamics approach to the thermodynamics of
soliton gases [22] provides a promising framework to explore this possibility.
At the same time, this direction of research could require some departure from
integrability and the development of the soliton gas theory for perturbed
integrable equations.

## Acknowledgments

The authors would like to thank the Isaac Newton Institute for Mathematical
Sciences for support and hospitality during the programme “Dispersive
hydrodynamics: mathematics, simulation and experiments, with applications in
nonlinear waves” when the work on this paper was undertaken. This work was
supported by EPSRC Grant Number EP/R014604/1. GE’s and GR’s work was also
supported by EPSRC Grant Number EP/W032759/1 and AT’s work was supported by
NSF Grant DMS 2009647. TC and GR thank Simons Foundation for partial support.
All authors thank T. Bonnemain, S. Randoux and P. Suret for numerous useful
discussions.

##  Appendix A Numerical implementation of soliton gas

###  A.1 Riemann problem

The realizations of the soliton gas are approximated numerically by the
$n$-soliton solution

$\varphi\equiv\varphi_{n}(x,t;\eta_{1},\dots,\eta_{n},x^{0}_{1},\dots,x^{0}_{n}),\quad
n\in\mathbb{N},$ (A.1)

where the $\eta_{i}$’s and $x^{0}_{i}$’s correspond respectively to the
spectral parameters and the “spatial phases” of the solitons;
$\eta_{i}<\eta_{i+1}$ by convention. The numerical implementation of (A.1) is
described in Sec. A.3 below. The numerical solutions presented in this work
are all generated with $n=200$ solitons, unless otherwise stated.

Since $n$ is finite, the $n$-soliton solution reduces to a sum of separated
solitons in the limit $|t|\to\infty$. By construction, we have in the limit
$t\to\pm\infty$

$\varphi_{n}(x,t)\sim\sum_{i=1}^{n}2\eta_{i}^{2}{\rm
sech}^{2}\left[\eta_{i}(x-4\eta_{i}^{2}t-x_{i}^{\pm})\right],$ (A.2)

where $x_{i}^{\pm}$ are the spatial phases of the $i$-th soliton at
$t\to\pm\infty$. We then take the spatial phase in (A.1) to be
$x_{i}^{0}=(x_{i}^{-}+x_{i}^{+})/2$.

Consider a uniform soliton gas with the density of states $u(\eta)$. Let the
spectral parameters $\eta_{i}$ be distributed on $\Gamma^{+}$ with density

$\phi(\eta)=\frac{u(\eta)}{\kappa},\quad\kappa=\int_{\Gamma_{+}}u(\eta)d\eta,$
(A.3)

where the normalization by the spatial density of solitons $\kappa$ ensures
that $\phi(\eta)$ is normalized to $1$. It was shown in [47] that the spatial
density $\kappa$ is obtained if the phases $x^{0}_{i}$ are uniformly
distributed on the interval (denoted “$S$-set” in [47]):

$I_{s}=\left[-\frac{n}{2\kappa_{s}},+\frac{n}{2\kappa_{s}}\right],\quad\kappa_{s}=\int_{\Gamma_{+}}\frac{\eta}{\sigma(\eta)}d\eta,$
(A.4)

where $\sigma(\eta)$ is the spectral scaling function in the NDRs (2.9),
(2.10); $y(\eta)$ in [47] is given here by $y(\eta)=u(\eta)\sigma(\eta)/\eta$.
The derivation of (A.4) has been revisited recently in the context of
generalized hydrodynamics [22]: it was shown that $\kappa_{s}$ corresponds to
the density of spatial phases $x_{i}^{0}$, or equivalently $x_{i}^{\pm}$ which
are well defined asymptotically ($t\to\pm\infty$) where the solitons are “non-
interacting” and their position are given by $x_{i}(t)\sim
4\eta_{i}^{2}t+x_{i}^{\pm}$. In the rarefied gas limit the interaction term in
the NDR (2.9) is small and therefore $\sigma(\eta)u(\eta)\approx\eta$ so that
we obtain $\kappa_{s}=\kappa$ as expected. In the general case though the
density $\kappa_{s}$ of non-interacting phases is different from the
“physical” density $\kappa$, as demonstrated with the soliton condensate
examples below. In the thermodynamic limit $n\to\infty$, the soliton solution
(A.1) represents a realization of the uniform soliton gas.

Since the number $n$ of solitons is finite, the $n$-soliton solution has a
finite spatial extent. By distributing the phases $x^{0}_{i}$ uniformly on the
interval $I_{s}$, the $n$-soliton solution $\varphi_{n}(x,t=0)$ approximates a
realization of the uniform soliton gas for $x\in[-\ell/2,\ell/2]$ where
$\ell=n/\kappa$; $\varphi_{n}(x,t=0)\sim 0$ outside of this interval. This
naturally generates the box-like initial condition for the kinetic equation

$u(\eta;x,t=0)\sim\begin{cases}0,&x<-\ell/2,\\\ u(\eta),&-\ell/2<x<\ell/2,\\\
0,&\ell/2<x.\end{cases}$ (A.5)

Note that $u(\eta;x,t=0)=0$ can be seen as the genus $0$ condensate where the
end point of the central s-band $\lambda_{1}\to 0$. This limits the type of
initial condition that can be implemented for the Riemann problem and we
choose in practice $(N_{-}=0,q_{-}=0)$ or $(N_{+}=0,q_{+}=0)$. For
convenience, we shift the $x$-axis by $\pm\ell/2$ to obtain one of the
discontinuities at the position $x=0$.

The evolution in time of the soliton gas realization is obtained by varying
the parameter $t$. Contrarily to a direct resolution of the KdV equation, via
finite difference or spectral method, the time-evolution presented here is
instantaneous and does not accumulate any numerical errors since the
$n$-soliton solution is an exact solution. For the Riemann problem, the
maximal time is bounded by the finite extent of the $n$-soliton solution:
after a sufficently long time, the two hydrodynamic states originating from
the discontinuities at $x=-\ell/2$ and $x=\ell/2$ start interacting. Longer
times can be reached by choosing a larger number of solitons $n$.

We consider now the density of states of interest for this work:

$u(\eta)=Cu^{(N)}(\eta;\lambda_{1},\dots\lambda_{2N+1}),\quad C\leq 1,$ (A.6)

where $u^{(N)}$ is density of states of the condensate defined in Sec. 3.
(A.4) rewrites

$\kappa_{s}=\frac{\kappa(C)}{1-C},\quad\kappa(C)=C\int_{\Gamma_{+}}u^{(N)}(\eta;\lambda_{1},\dots,\lambda_{2N+1})d\eta.$
(A.7)

Fig. 18 shows the comparison between the spatial density of solitons $\kappa$
and the density of phases $\kappa_{s}$ for the genus $0$ case where
$\kappa(C)=C/\pi$. The phases density $\kappa_{s}$ diverges in the condensate
limit $C\to 1$, and $x_{i}^{0}$’s are all equal to the same phase $x^{0}$
($I_{s}\to\\{x^{0}\\}$). This limit is in agreement with the results obtained
in Sec. 4.1 for genus $0$ and genus $1$ condensates: each realisation of the
condensate ($C=1$) is approximated with a coherent $n$-soliton solution where
$x_{i}^{0}=x^{0}={\rm cst},\;\forall i$.

Figure 18: The solid line represents the variation of $\kappa_{s}$ with
respect to $\kappa$ for a diluted genus $0$ condensate, cf. (A.7). The markers
are obtained using the $100$-soliton solution: $\kappa=\ell/n$ where $\ell$
corresponds to the spatial extension of the $n$-soliton solution.

Examples of numerical realizations of soliton condensates and diluted soliton
condensates are given in Sections 4.1, 6, 7 and B. Figs. 2, 4 and 20 shows
that numerical approximations of condensate via the $n$-soliton solution are
not exactly uniform; realizations become more uniform near the center of the
interval $[-\ell,\ell]$ as the number of soliton $n$ increases.

The realization at $t=0$ in Figs. 7, 8, 10 and 12 also displays the “border
effects” observed at the discontinuities of the Riemann problem initial
condition (located at $x=0$). These border effects, manifesting as overshoots
of the realization, persist regardless of the number of solitons $n$ as shown
by the comparison between the $100$-soliton and $200$-soliton solutions in
Fig. 19. However, because of their finite size, the observed border effects
seem to have no effect on the asymptotic dynamics of the condensate as
demonstrated by the very good agreement between the theory and the numerical
solution in Sec. 6.

Figure 19: $n$-soliton solution approximating a realization of the condensate
$(N;\lambda_{1},\lambda_{2},\lambda_{3})=(1;0,0.5,0.85)$. The solid black line
represent the solution $n=100$ and the red dashed line the solution $n=200$.
Both solutions have been shifted such that the maximum of the solution is
located at $x=0$.

###  A.2 Generation of spectral parameters $\eta_{i}$

The spectral parameters of the $n$-soliton solution are distributed with
probability density $\phi(\eta)$, cf. (A.3). This can be achieved by choosing
the solutions of the nonlinear equation

$\int_{0}^{\eta_{i}}\phi(\mu)d\mu=\frac{i}{n},\quad i=1\dots n.$ (A.8)

For genus 0 condensate whose DOS is given by (4.1), this equation reduces to

$1-\sqrt{1-\frac{\eta_{i}^{2}}{\lambda_{1}^{2}}}=\frac{i}{n}.$ (A.9)

For genus 1 condensate with DOS (4.6), this equation reads

$\frac{U(\eta_{i})}{U(\lambda_{3})}=\frac{i}{n},\quad\text{where:}\quad
U(\eta)=\int_{0}^{\eta}u^{(1)}(\mu;\lambda_{1},\lambda_{2},\lambda_{3})d\mu.$
(A.10)

We have

$U(\eta)=\int_{0}^{\eta}\frac{i(\mu^{2}-p^{2})}{2\pi\sqrt{R(\mu)}}d(\mu^{2})=\int_{0}^{\eta^{2}}\frac{i(x-p^{2})}{2\pi\sqrt{(x-\lambda_{1}^{2})(x-\lambda_{2}^{2})(x-\lambda_{3}^{2})}}dx,$
(A.11)

which yields

$U(\eta)=U_{0}+\begin{cases}\displaystyle-\frac{1}{\pi}\left[\sqrt{\frac{(\lambda_{3}^{2}-\eta^{2})(\lambda_{1}^{2}-\eta^{2})}{\lambda_{2}^{2}-\eta^{2}}}-\frac{\lambda_{1}^{2}-p^{2}}{\sqrt{\lambda_{3}^{2}-\lambda_{1}^{2}}}F(\beta,q)-\sqrt{\lambda_{3}^{2}-\lambda_{1}^{2}}E(\beta,q)\right],&0<\eta<\lambda_{1},\\\
0,&\lambda_{1}<\eta<\lambda_{2},\\\
\displaystyle\frac{1}{\pi}\left[\frac{\lambda_{1}^{2}-p^{2}}{\sqrt{\lambda_{3}^{2}-\lambda_{1}^{2}}}F(\kappa,q)+\frac{\lambda_{2}^{2}-\lambda_{1}^{2}}{\sqrt{\lambda_{3}^{2}-\lambda_{1}^{2}}}\Pi(q,\kappa,q)\right],&\lambda_{2}<\eta<\lambda_{3}.\end{cases}$
(A.12)

where:

$\displaystyle\beta=\sin^{-1}\left(\sqrt{\frac{\lambda_{1}^{2}-\eta^{2}}{\lambda_{2}^{2}-\eta^{2}}}\right),\quad\kappa=\sin^{-1}\left(\sqrt{\frac{(\lambda_{3}^{2}-\lambda_{1}^{2})(\eta^{2}-\lambda_{2}^{2})}{(\lambda_{3}^{2}-\lambda_{2}^{2})(\eta^{2}-\lambda_{1}^{2})}}\right),$
(A.13) $\displaystyle
q=\frac{\lambda_{3}^{2}-\lambda_{2}^{2}}{\lambda_{3}^{2}-\lambda_{1}^{2}},$
(A.14) $\displaystyle
U_{0}=\frac{1}{\pi}\left[\sqrt{\frac{\lambda_{3}^{2}\lambda_{1}^{2}}{\lambda_{2}^{2}}}-\frac{\lambda_{1}^{2}-p^{2}}{\sqrt{\lambda_{3}^{2}-\lambda_{1}^{2}}}F(\beta_{0},q)-\sqrt{\lambda_{3}^{2}-\lambda_{1}^{2}}E(\beta_{0},q)\right],\quad\beta_{0}=\sin^{-1}\left(\frac{\lambda_{1}}{\lambda_{2}}\right).$
(A.15)

###  A.3 Algorithm for the $n$-soliton solution

The algorithm generating the exact $n$-soliton (A.1), originally developed in
[45], relies on the Darboux transformation. This scheme is subject to roundoff
errors during summation of exponentially small and large values for a large
number of solitons $n$. We improve it following [46], with the implementation
of high precision arithmetic routine to overcome the numerical accuracy
problems and generate solutions with a number of solitons $n\gtrsim 10$.

In order to simplify the algorithm, it is suggested to consider simultaneously
the KdV equation (2.1) and its equivalent form

$\varphi-6\varphi\varphi_{x}+\varphi_{xxx}=0$ (A.16)

obtained from (2.1) by the reflection $\varphi\to-\varphi$. The Darboux
transformation presented here relates the Jost solution associated with the
$(n-1)$-soliton solution of one equation, to the $n$-soliton solution of the
other equation.

Considering the direct scattering problem for the Lax pair in the matrix form

$\Phi_{x}=\left(\begin{array}[]{cc}\eta&\mp 1\\\
\varphi&-\eta\end{array}\right)\Phi,$ (A.17)

with $-1$ corresponding to (2.1) and $+1$ to (A.16), the Jost solutions
$J,\tilde{J}\in\mathbb{R}^{2\times 2}$ are defined recursively by the Darboux
transformations $D(\eta)$ and $\tilde{D}(\eta)$ such that:

$J_{n}(\eta)=D_{n}(\eta)J_{n-1}(\eta),\quad\textrm{with:}\quad
D_{n}(\eta)=I+\frac{2\eta_{n}}{\eta-\eta_{n}}P_{n},$ (A.18)
$\tilde{J}_{n}(\eta)=\tilde{D}_{n}(\eta)\tilde{J}_{n-1}(\eta),\quad\textrm{with:}\quad\tilde{D}_{n}(\eta)=I-\frac{2\eta_{n}}{\eta+\eta_{n}}\tilde{P}_{n}.$
(A.19)

$P_{n}(x,t)$ and $\tilde{P}_{n}(x,t)$ are independent of $\eta$ and have the
form:

$P_{n}=\sigma_{2}\tilde{P}^{\textrm{T}}_{n}\sigma_{2}=\frac{J_{n-1}\left(-\eta_{n}\right)\left(\begin{array}[]{c}-b_{n}\\\
1\end{array}\right)\left(\begin{array}[]{ll}b_{n}&1\end{array}\right)\tilde{J}_{n-1}^{-1}\left(\eta_{n}\right)}{\left(\begin{array}[]{ll}b_{n}&1\end{array}\right)\tilde{J}_{n-1}^{-1}\left(\eta_{n}\right)J_{n-1}\left(-\eta_{n}\right)\left(\begin{array}[]{c}-b_{n}\\\
1\end{array}\right)}$ (A.20)

with the real constants $b_{n}$ depending on the spatial phases

$b_{n}=\left(-1\right)^{n}\exp\left(2\eta_{n}x^{0}_{n}\right).$ (A.21)

The Jost solutions for the initial seed solution $\varphi_{0}=0$ are given by

$J_{0}(\eta)=\tilde{J}_{0}(\eta)=\left(\begin{array}[]{cc}\exp\left[\eta
x-4\eta^{3}t\right]&-\exp\left[-\eta x+4\eta^{3}t\right]\\\
0&-2\eta\exp\left[-\eta x+4\eta^{3}t\right]\end{array}\right),$ (A.22)

and one can show that at each recursion step

$\varphi_{n}=\varphi_{n-1}+4\eta_{n}\left(P_{n}\right)_{21},$ (A.23)

where $\varphi_{n}$ is the $n$-soliton solution of (2.1) for $n$ even and
solution of (A.16) for $n$ odd. Recently, a more efficient and accurate
algorithm has been proposed in [51] to generate the $n$-soliton KdV solution
employing a $2$-fold Crum transform.

##  Appendix B Genus $2$ condensate

In Fig. 20 a realization of the genus 2 soliton condensate is compared with
the two-phase KdV solution associated with the same spectral surface and
equipped with an appropriately chosen initial phase vector. The two-phase
solution has been computed numerically using the so-called trace formula [13]:

$\varphi(x,t)=\lambda_{1}^{2}-2\sum_{j=1}^{2}\left(\mu_{j}^{2}(x,t)-\frac{\lambda_{2j}^{2}+\lambda_{2j+1}^{2}}{2}\right),$
(B.1)

where the auxiliary spectra $\mu_{j}(x,t)$ satisfy Dubrovin’s ordinary
differential equations:

$\frac{\partial\mu_{j}^{2}}{\partial
x}=\frac{2\sigma_{j}R\left(\mu_{j}\right)}{\prod_{j\neq
k}^{2}\left(\mu_{j}^{2}-\mu_{k}^{2}\right)},\quad j=1,2$ (B.2)

with $\sigma=\pm 1$ and

$R(\mu)=\sqrt{\left(\mu^{2}-\lambda_{1}^{2}\right)\left(\mu^{2}-\lambda_{2}^{2}\right)\left(\mu^{2}-\lambda_{3}^{2}\right)\left(\mu^{2}-\lambda_{4}^{2}\right)\left(\mu^{2}-\lambda_{5}^{2}\right)}.$
(B.3)

Each $\mu_{j}$ oscillates in the corresponding s-gap
$[\lambda_{2j-1},\lambda_{2j}]$ so that the sign of $\sigma_{j}$ changes every
time $\mu_{j}$ changes direction of motion at the gap end point. We observe
that, to compute the KdV finite-gap solution corresponding to a given
realization soliton condensate, all the initial phases $\mu_{j}(x_{0},t)$ must
be placed at the edges of the corresponding s-gaps while the choice of the
gap’s edge (right/left) is determined by the number of discrete eigenvalues
(odd/even) that are located within the s-band $[\lambda_{2j},\lambda_{2j+1}]$
in the numerical construction of the condensate.

Figure 20: a) Comparison between KdV genus 2 soliton condensate realized via
$n=204$-soliton solution (solid line) and the exact 2-phase KdV solution
(dashed line) for
$\lambda_{1}=0.3,\lambda_{2}=0.5,\lambda_{3}=0.7,\lambda_{4}=0.9,\lambda_{5}=1.$;
the two plots are visually indistinguishable from one another; b) Absolute
error for the condensate generated with $50$ solitons (dashed line), $100$
solitons (solid line) and $204$ solitons (dash-dotted line); for readability
the absolute error is evaluated at the extrema of the solutions.

## References

  * [1] V. E. Zakharov, “Turbulence in integrable systems,” Studies in Applied Mathematics, vol. 122, pp. 219–234, 2009.
  * [2] V. Zakharov, “Kinetic equation for solitons,” Journal of Experimental and Theoretical Physics, vol. 33, pp. 538–541, 1971.
  * [3] G. El, “The thermodynamic limit of the Whitham equations,” Physics Letters A, vol. 311, pp. 374–383, 2003.
  * [4] H. Flaschka, M. G. Forest, and D. W. McLaughlin, “Multiphase averaging and the inverse spectral solution of the Korteweg-de Vries equation,” Comm. Pure Appl. Math., vol. 33, pp. 739–784, 1980.
  * [5] G. A. El and A. M. Kamchatnov, “Kinetic equation for a dense soliton gas,” Physical Review Letters, vol. 95, p. 204101, Nov 2005.
  * [6] T. Congy, G. El, and G. Roberti, “Soliton gas in bidirectional dispersive hydrodynamics,” Physical Review E, vol. 103, p. 042201, 2021.
  * [7] G. El and A. Tovbis, “Spectral theory of soliton and breather gases for the focusing nonlinear Schrödinger equation,” Physical Review E, vol. 101, p. 052207, 2020.
  * [8] A. Kuijlaars and A. Tovbis, “On minimal energy solutions to certain classes of integral equations related to soliton gases for integrable systems,” Nonlinearity 34 no. 10 7227 (2021)
  * [9] P. D. Lax and C. D. Levermore, “The small dispersion limit of the Korteweg-de Vries equation: 2,” Comm. Pure Appl. Math., vol. 36, no. 5, pp. 571–593, 1983.
  * [10] M. Girotti, T. Grava, and K. D. T.-R. McLaughlin, “Rigorous asymptotics of a KdV soliton gas,” Comm. Math. Phys., vol. 384, pp. 733–784, 2021.
  * [11] A. V. Gurevich, L. P. Pitaevskii, Nonstationary structure of a collisionless shock wave, Sov. Phys. JETP 38 (2) (1974) 291–297, translation from Russian of A. V. Gurevich and L. P. Pitaevskii, Zh. Eksp. Teor. Fiz. 65, 590-604 (August 1973).
  * [12] G.A. El, “Soliton gas in integrable dispersive hydrodynamics,” J. Stat. Mech.: Theor. Exp. 114001, 2021.
  * [13] S. P. Novikov, S. Manakov, L. P. Pitaevskii, and V. Zakharov, Theory of Solitons: The Inverse Scattering Method.  Monographs in Contemporary Mathematics, Springer, 1984.
  * [14] S. Venakides, “The continuum limit of theta functions,” Communications on Pure and Applied Mathematics, vol. 42, no. 6, pp. 711–728, 1989.
  * [15] A. Tovbis and F. Wang, “ Recent developments in spectral theory of the focusing NLS soliton and breather gases: the thermodynamic limit of average densities, fluxes and certain meromorphic differentials; periodic gases,” arXiv:2203.03566.
  * [16] E. Pelinovsky and E. Shurgalina.  KDV soliton gas: interactions and turbulence.  In Advances in Dynamics, Patterns, Cognition, 295-306,  Springer, Cham., 2017
  * [17] F. Tricomi, “On The Finite Hilbert Transformation,” Q J Math, vol. 2, no. 1, pp. 199–211, 1951.
  * [18] S. Okada and D. Elliot, “The finite Hilbert transform in $\mathcal{L}^{2}$,” Mathematische Nachrichten, vol. 153, no. 1, pp. 43–56, 1991.
  * [19] B. A. Dubrovin and S. P. Novikov, “Hydrodynamics of weakly deformed soliton lattices. Differential geometry and Hamiltonian theory,” Russian Mathematical Surveys, vol. 44, no. 6, pp. 35–124, 1989.
  * [20] S. P. Tsarëv, “The geometry of Hamiltonian systems of hydrodynamic type. the generalized hodograph method,” Mathematics of the USSR-Izvestiya, vol. 37, pp. 397–419, Apr 1991.
  * [21] B. Doyon.  Lecture notes on Generalised Hydrodynamics.  In SciPost Physics Lecture Notes, page 18, 2020.
  * [22] T. Bonnemain, B. Doyon and G. A. El , Generalized hydrodynamics of the KdV soliton gas, Journ. Phys. A: Math. Theor. 2022 (accepted), arXiv:2203.08551
  * [23] E. Bettelheim.  The Whitham approach to the $c\to 0$ limit of the Lieb–Liniger model and generalized hydrodynamics.  Journal of Physics A: Mathematical and General, 53:205204, 2020.
  * [24] G. A. El, A. M. Kamchatnov, M. V. Pavlov, and S. A. Zykov, “Kinetic Equation for a Soliton Gas and Its Hydrodynamic Reductions,” Journal of Nonlinear Science, vol. 21, pp. 151–191, 2011.
  * [25] M. V. Pavlov, V. B. Taranov, and G. A. El, “Generalized hydrodynamic reductions of the kinetic equation for a soliton gas,” Theoretical and Mathematical Physics, vol. 171, pp. 675–682, 2012.
  * [26] E.V. Ferapontov and M.V. Pavlov, Kinetic Equation for Soliton Gas: Integrable Reductions, J. Nonlin. Sci., 32, 26, 2022
  * [27] F. Carbone, D. Dutykh, and G. A. El, “Macroscopic dynamics of incoherent soliton ensembles: Soliton gas kinetics and direct numerical modelling,” EPL (Europhysics Letters), vol. 113, p. 30003, 2016.
  * [28] C.D. Levermore, The hyperbolic nature of the zero dispersion KdV limit, Comm. Part. Diff. Eq., 13, 495-514, 1988.
  * [29] V. K. Rohatgi and A.K.Md. Ehsanes Saleh, An introduction to probability and statistics, vol. 1.  John Wiley & Sons, 2015.
  * [30] P. Drazin and R. Johnson, Solitons: an Introduction, 2nd ed. Cambridge University Press, 1989.
  * [31] I. S. Gradshteyn, I. M. Ryzhik and R. H. Romer, Tables of integrals, series, and products, 1988.
  * [32] M. Girotti, T. Grava, R. Jenkins, K. McLaughlin, A. Minakov, Soliton v. the gas: Fredholm determinants, analysis, and the rapid oscillations behind the kinetic equation, arXiv:2205.0260, 2022.
  * [33] P. F. Byrd and M. D. Friedman, Handbook of elliptic integrals for engineers and physicists,  Springer-Verlag, 1954.
  * [34] A. M. Kamchatnov, Nonlinear periodic waves and their modulations: an introductory course.  World Scientific, 2000.
  * [35] G. El and M. Hoefer, “Dispersive shock waves and modulation theory,” Physica D: Nonlinear Phenomena, vol. 333, pp. 11–65, Oct 2016.
  * [36] S. Dyachenko, D. Zakharov, and V. Zakharov, “Primitive potentials and bounded solutions of the KdV equation,” Physica D: Nonlinear Phenomena, vol. 333, pp. 148–156, 2016.
  * [37] B. Dubrovin, Functionals of the Peierls-Frölich type and variational principle for Whitham equations, Amer. Math. Soc. Transl. 179 (1997) 35-44
  * [38] G. A. El, A. L. Krylov, S. Venakides, Unified approach to KdV modulations, Comm. Pure Appl. Math. 54 (2001) 1243–1270.
  * [39] T. Grava, F.-R. Tian, The generation, propagation, and extinction of multiphases in the KdV zero-dispersion limit, Comm. Pur. Appl. Math. 55 (12) (2002) 1569–1639.
  * [40] G. B. Whitham, “Non-linear dispersive waves,” Proc. Roy. Soc. Ser. A, vol. 283, pp. 238–261, 1965.
  * [41] P. D. Lax, Hyperbolic systems of conservation laws and the mathematical theory of shock waves.  SIAM, 1973.
  * [42] M. J. Ablowitz, Nonlinear dispersive waves: asymptotic analysis and solitons.  Cambridge University Press, 2011.
  * [43] P. Sprenger, and M. A. Hoefer, “Discontinuous shock solutions of the Whitham modulation equations as zero dispersion limits of traveling waves.” Nonlinearity 33 no. 10 3268 (2020).
  * [44] S. Gavrilyuk, B. Nkonga, K. M. Shyue and L. Truskinovsky, “Stationary shock-like transition fronts in dispersive systems”, Nonlinearity, 33 no. 10 5477 (2020).
  * [45] Nian-Ning Huang.  “Darboux transformations for the Korteweg-de-Vries equation”.  Journal of Physics A: Mathematical and General, 25(2):469, 1992\.
  * [46] A.A. Gelash and D.S. Agafontsev.  Strongly interacting soliton gas and formation of rogue waves.  Phys. Rev. E, 98(4):042210, 2018.
  * [47] A. V. Gurevich, N. G. Mazur, K. P. Zybin, “Statistical limit in a completely integrable system with deterministic initial conditions,” Journal of Experimental and Theoretical Physics, 90, pp. 695–713, 2000.
  * [48] P. J. Prins and S. Wahls, “An accurate $\mathcal{O}(N^{2})$ floating point algorithm for the Crum transform of the KdV equation”, Commun Nonlinear Sci Numer Simul, 102, 105782 (2021)
  * [49] F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds., NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/, Release 1.1.6 of 2022-06-30.
  * [50] A. Gelash, D. Agafontsev, P. Suret, and S. Randoux, “Solitonic model of the condensate”, Phys. Rev. E, 104(4):044213, 2021.
  * [51] P. J. Prins and S. Wahls, “An accurate $\mathcal{O}(N^{2})$ floating point algorithm for the Crum transform of the KdV equation”, Commun Nonlinear Sci Numer Simul, 102, 105782 (2021)
  * [52] B. Doyon, “Lecture notes on Generalised Hydrodynamics”, SciPost Physics Lecture Notes, 018 (2020)
  * [53] T. Bonnemain, G. El, and B. Doyon, “Generalized hydrodynamics of the KdV soliton gas”, J. Phys. A: Math. Theor., 55, 374004, 2022.
  * [54] E. Bettelheim, ‘The Whitham approach to the $c\to 0$ limit of the Lieb–Liniger model and generalized hydrodynamics”, J. Phys. A: Math. Theor., 53, 205204, 2020.
  * [55] A. Gelash, D. Agafontsev, V. Zakharov, G. El, S. Randoux, and P. Suret, “Bound State Soliton Gas Dynamics Underlying the Spontaneous Modulational Instability”, Phys. Rev. Lett., 123, 234102, 2019.
  * [56] D. Agafontsev and V. Zakharov, “Integrable turbulence and formation of rogue waves”, Nonlinearity, 28, 2791–2821, 2015.
  * [57] D. Agafontsev, S. Randoux, and P. Suret, “Extreme rogue wave generation from narrowband partially coherent waves”, Phys. Rev. E, 103(3):032209, 2021.
  * [58] E. D. Belokolos, A. I. Bobenko, V. Z. Enolski, A. R. Its, and V. B. Matveev, Algebro-geometric approach to nonlinear integrable equations.  New York: Springer, 1994.
  * [59] V. B. Matveev, “30 years of finite-gap integration theory,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 366, no. 1867, pp. 837–875, 2008.
  * [60] G. Roberti, G. El, A. Tovbis, F. Copie, P. Suret, and S. Randoux, “Numerical spectral synthesis of breather gas for the focusing nonlinear Schrodinger equation”, Phys. Rev. E, 103(4):042205, 2021.