Soliton Gas: Integrable Kinetic Theory

• Soliton Gas is a macroscopic ensemble of localized solitons whose elastic, factorized scattering preserves amplitudes while imparting velocity-dependent position shifts.
• Kinetic theory and generalized hydrodynamics yield nonlocal transport equations that describe the evolution of the soliton density of states in spectral space.
• Experimental and numerical studies in water tanks and quantum systems validate the theory, revealing extreme events and novel integrable turbulence phenomena.

Searching arXiv for recent and foundational papers on soliton gas to ground the article in primary sources. Soliton gas denotes a macroscopic ensemble of solitons or closely related localized integrable excitations whose individual parameters are random but whose mutual scattering remains elastic and factorized. In classical integrable equations such as KdV and the focusing or defocusing nonlinear Schrödinger equation, and in quantum integrable many-body systems through Bethe quasi-particles, the defining microscopic feature is the velocity-dependent position shift acquired in two-body collisions while amplitudes and asymptotic velocities remain unchanged. At mesoscopic and macroscopic scales, this many-body state is described not by individual trajectories but by a density of states in spectral space and by kinetic or hydrodynamic equations for its transport (Doyon et al., 2017, Suret et al., 2023).

1. Concept and historical placement

The modern notion of a soliton gas goes back to Zakharov’s 1971 proposal of an infinite collection of weakly interacting KdV solitons, and later developments extended the concept from dilute gases of almost non-overlapping solitons to dense gases in which interactions are strong and continuous. In this usage, the “particles” are not molecules but soliton modes of an integrable PDE, each labeled by an inverse-scattering spectral parameter and a phase or norming constant. The central statistical object is therefore a density of states, usually written as $f(\eta;x,t)$ or $u(\lambda;x,t)$, specifying the expected number of solitons in a mesoscopic cell of spectral and physical space (Suret et al., 2023, El, 2021).

Two distinctions are fundamental. First, a dilute gas is collision-sparse and is well approximated by a random collection of nearly isolated solitons, whereas a dense gas is not pointwise decomposable into visibly separate pulses even though its inverse-scattering data still consist of discrete solitonic modes. Second, zero-background gases of solitons must be distinguished from breather gases on nonzero background in focusing NLS, where the presence of a Stokes band changes both the spectral construction and the macroscopic transport theory (El et al., 2019).

The concept is not restricted to classical wave equations. In generalized hydrodynamics of one-dimensional quantum integrable models, Bethe quasi-particles carry the same scattering data encoded by the two-body differential phase $\varphi$, and at the Euler scale their hydrodynamics coincides with that of a classical gas of solitons whose collisions implement the corresponding position shifts (Doyon et al., 2017). This quantum-classical equivalence places soliton gas theory at the interface of integrable dispersive hydrodynamics, inverse scattering, thermodynamic Bethe ansatz, and nonequilibrium many-body transport.

2. Kinetic and generalized-hydrodynamic formulations

In the kinetic description, soliton number is conserved in each spectral layer. For KdV and, in suitable form, for broader classes of integrable models, one writes

$$\partial_t f(\eta,x,t)+\partial_x\bigl[s(\eta;x,t)\,f(\eta,x,t)\bigr]=0,$$

where $f(\eta,x,t)\,d\eta\,dx$ is the expected number of solitons with parameter in $[\eta,\eta+d\eta]$ inside $[x,x+dx]$, and $s(\eta;x,t)$ is the effective velocity of a tracer soliton in the gas. The closure is nonlocal because $s$ is renormalized by the cumulative phase shifts generated by pairwise collisions. In a general collision-rate form,

$$s(\eta)=s_0(\eta)+\int G(\eta,\mu)\,f(\mu)\,[s(\eta)-s(\mu)]\,d\mu,$$

with $u(\lambda;x,t)$0 the isolated-soliton velocity and $u(\lambda;x,t)$1 the two-body shift kernel (Suret et al., 2023, El, 2021).

For KdV, where a one-soliton state has amplitude $u(\lambda;x,t)$2 and free speed $u(\lambda;x,t)$3, the logarithmic phase shift determines the specific kernel. In the generalized-hydrodynamic formulation of the KdV soliton gas, the differential scattering phase is

$u(\lambda;x,t)$4

with bare derivatives $u(\lambda;x,t)$5 and $u(\lambda;x,t)$6, and the effective velocity is the dressed ratio

$u(\lambda;x,t)$7

The slowly varying density of states $u(\lambda;x,t)$8 then obeys

$u(\lambda;x,t)$9

exactly as in the Zakharov–El kinetic equation (Bonnemain et al., 2022).

Generalized hydrodynamics expresses the same structure in the rapidity language of integrable many-body systems. The fundamental variable is the quasi-particle density $\varphi$0, satisfying

$\varphi$1

with

$\varphi$2

Here $\varphi$3 is the occupation function and $\varphi$4 the density of states. Conserved densities and currents are then

$\varphi$5

so once $\varphi$6 is known, local observables follow directly (Doyon et al., 2017).

A particularly concrete realization of the quantum-classical equivalence is the “flea gas” or generalized hard-rod model. In this classical gas, point particles move ballistically and, upon collision of rapidities $\varphi$7, jump by

$\varphi$8

The resulting Euler-scale hydrodynamics coincides with GHD, and the associated molecular-dynamics algorithm, using on the order of $\varphi$9 quasi-particles with ballistic evolution plus collision jumps, provides a numerically efficient simulator for inhomogeneous dynamics in integrable chains and the Lieb–Liniger model (Doyon et al., 2017).

3. Spectral theory, thermodynamic limits, and position variables

Inverse scattering gives the natural microscopic state space for a soliton gas. For localized wavefields in focusing NLS, the scattering data consist of a discrete set of eigenvalues $$\partial_t f(\eta,x,t)+\partial_x\bigl[s(\eta;x,t)\,f(\eta,x,t)\bigr]=0,$$0 in the upper half-plane, each corresponding to a soliton, plus a continuous spectrum represented by a reflection coefficient. A soliton gas is the limit of large $$\partial_t f(\eta,x,t)+\partial_x\bigl[s(\eta;x,t)\,f(\eta,x,t)\bigr]=0,$$1 in which the eigenvalues densely populate a region or contour in spectral space with a prescribed density of states $$\partial_t f(\eta,x,t)+\partial_x\bigl[s(\eta;x,t)\,f(\eta,x,t)\bigr]=0,$$2, or equivalently a normalized probability density $$\partial_t f(\eta,x,t)+\partial_x\bigl[s(\eta;x,t)\,f(\eta,x,t)\bigr]=0,$$3 (Suret et al., 2020).

A complementary and more structural construction starts from finite-gap solutions. In focusing NLS, a thermodynamic limit of genus-$$\partial_t f(\eta,x,t)+\partial_x\bigl[s(\eta;x,t)\,f(\eta,x,t)\bigr]=0,$$4 finite-gap fields is taken by collapsing bands while their centers fill a contour $$\partial_t f(\eta,x,t)+\partial_x\bigl[s(\eta;x,t)\,f(\eta,x,t)\bigr]=0,$$5 with macroscopic density. The resulting nonlinear dispersion relations become singular integral equations for the spectral density $$\partial_t f(\eta,x,t)+\partial_x\bigl[s(\eta;x,t)\,f(\eta,x,t)\bigr]=0,$$6 and its temporal counterpart $$\partial_t f(\eta,x,t)+\partial_x\bigl[s(\eta;x,t)\,f(\eta,x,t)\bigr]=0,$$7. For a zero-background soliton gas, the ratio $$\partial_t f(\eta,x,t)+\partial_x\bigl[s(\eta;x,t)\,f(\eta,x,t)\bigr]=0,$$8 satisfies

$$\partial_t f(\eta,x,t)+\partial_x\bigl[s(\eta;x,t)\,f(\eta,x,t)\bigr]=0,$$9

with

$f(\eta,x,t)\,d\eta\,dx$0

and $f(\eta,x,t)\,d\eta\,dx$1 is precisely the two-soliton position shift. In a slowly varying inhomogeneous gas, this produces the transport equation

$f(\eta,x,t)\,d\eta\,dx$2

for the density of states (El et al., 2019).

This framework distinguishes several regimes. In the ideal or dilute limit, the secular term dominates and the velocity approaches the free value $f(\eta,x,t)\,d\eta\,dx$3. In the opposite sub-exponential scaling, the secular term drops out and one obtains a condensate density from a criticality condition. For example, the bound-state condensate on $f(\eta,x,t)\,d\eta\,dx$4 has

$f(\eta,x,t)\,d\eta\,dx$5

while for a circular condensate $f(\eta,x,t)\,d\eta\,dx$6,

$f(\eta,x,t)\,d\eta\,dx$7

These formulas formalize the distinction between an ideal gas of weakly interacting objects and a limiting state whose properties are set entirely by pairwise interactions (El et al., 2019).

Hydrodynamic reductions provide another exact entry point. Under the $f(\eta,x,t)\,d\eta\,dx$8-functional ansatz

$f(\eta,x,t)\,d\eta\,dx$9

the kinetic equation reduces to a $[\eta,\eta+d\eta]$0-component quasilinear system whose coefficient matrix consists of $[\eta,\eta+d\eta]$1 Jordan blocks. This reduced system possesses a hierarchy of commuting hydrodynamic flows and can be solved by a generalized hodograph method extending Tsarev’s theory beyond the diagonalizable case (Ferapontov et al., 2021).

A persistent issue in dense gases is that the field may no longer reveal where the constituent solitons are. For KdV, this has been addressed by introducing “magnifying-glass positions,” effective positions $[\eta,\eta+d\eta]$2, and a fluid-cell projection that removes solitons outside a mesoscopic interval while leaving the field unchanged inside that interval. On large scales these effective positions satisfy semi-classical Bethe equations with the two-body shifts $[\eta,\eta+d\eta]$3, and a non-rigorous derivation from these equations reproduces El’s kinetic equation (Doyon, 18 May 2026).

4. Statistical regimes, dense interactions, and extreme events

The statistical behavior of a soliton gas depends strongly on density and velocity spread. For the focusing one-dimensional NLS, numerically generated $[\eta,\eta+d\eta]$4-soliton solutions with $[\eta,\eta+d\eta]$5 show that in the rarefied limit $[\eta,\eta+d\eta]$6, for equal amplitudes $[\eta,\eta+d\eta]$7 and velocity variance $[\eta,\eta+d\eta]$8,

$[\eta,\eta+d\eta]$9

and the kurtosis obeys

$[x,x+dx]$0

As density $[x,x+dx]$1 or velocity spread $[x,x+dx]$2 increases, collisions become frequent, both $[x,x+dx]$3 and $[x,x+dx]$4 rise above the dilute-gas formulas, the wave-action spectrum broadens, and the intensity PDF moves closer to the exponential Rayleigh law (Gelash et al., 2018).

The same simulations identify rogue-wave events with multi-soliton collisions rather than with a distinct non-solitonic mechanism. Some of these events have spatial profiles close to Peregrine solutions of various orders, and a carefully tuned three-soliton collision can reproduce the temporal evolution of the second-order Peregrine maximum with high accuracy. The point is not that the gas ceases to be solitonic, but that rational-breather-like waveforms can emerge as collision geometries within a purely solitonic ensemble (Gelash et al., 2018).

A different route to dense states appears in the adiabatic-growth scenario for focusing NLS integrable turbulence. Starting from homogeneous Gaussian noise and evolving the system with weak linear pumping until a target intensity is reached, one observes a sequence of statistically stationary NLS states: first a dilute gas, then a soliton-dominated regime, and eventually a dense bound-state soliton gas in which the real parts of the eigenvalues collapse toward a common value. In a representative bound-state regime with $[x,x+dx]$5, the reported statistics are $[x,x+dx]$6, $[x,x+dx]$7, $[x,x+dx]$8, and rogue-wave probability

$[x,x+dx]$9

about two orders of magnitude above the linear value $s(\eta;x,t)$0 (Agafontsev et al., 2022).

By contrast, high soliton content does not automatically imply isolated extreme crests. Deep-ocean measurements analyzed by nonlinear Fourier transform define the soliton energy ratio

$s(\eta;x,t)$1

with

$s(\eta;x,t)$2

The identified high-$s(\eta;x,t)$3 open-ocean events have short periods, relatively small wave heights, steepness $s(\eta;x,t)$4, and $s(\eta;x,t)$5, yet skewness is approximately zero, kurtosis is unremarkable, and the abnormality index is low, so these states are soliton-dominated without exhibiting a single conspicuous “rogue” extreme (Lee et al., 6 Oct 2025).

5. Experimental realizations and field observations

Several experimental platforms now realize or detect soliton gases directly, spanning shallow-water flumes, deep-water wave tanks, and open-ocean measurements.

System and paper	Setup	Central result
Shallow-water bidirectional gas (Redor et al., 2019)	34 m flume, depth $s(\eta;x,t)$6 cm, sinusoidal forcing near $s(\eta;x,t)$7 Hz	Stationary dense bidirectional gas; collisions agree with Kaup–Boussinesq two-soliton theory
Deep-water spectral synthesis (Suret et al., 2020)	148 m $s(\eta;x,t)$8 5 m $s(\eta;x,t)$9 3 m tank, 20 gauges, $s$0-soliton signals with $s$1 up to 128	First controlled synthesis of a dense hydrodynamic soliton gas with prescribed DOS
Counter-propagating gas jets (Fache et al., 2023)	140 m tank, two $s$2 soliton beams, $s$3	Measured densities and velocities agree with spectral kinetic theory
Deep open-ocean detection (Lee et al., 6 Oct 2025)	5 977 Eluanbi records, NFT-based $s$4, directional filtering	First field evidence of deep-ocean soliton gas sea states
2D wave-turbulence transition (Leduque et al., 30 Sep 2025)	27 m $s$5 30 m basin, $s$6 m, stereoscopic profilometry	Transition from dispersive turbulence to soliton gas with threshold near $s$7

In the shallow-water laboratory realization, a sinusoidal wavetrain undergoes fission into thousands of solitons that propagate back and forth in a 34 m flume. Despite viscous damping with decay time $s$8 s, the gas reaches a statistically stationary state in which head-on and overtaking collisions can be isolated, fitted, and compared with exact Kaup–Boussinesq two-soliton solutions. The space-time Fourier spectrum resolves a pair of nearly straight soliton branches, a weakly excited linear-radiation branch, and standing-mode signatures. Reported one-point statistics for the right-going component are skewness $s$9 and kurtosis $$s(\eta)=s_0(\eta)+\int G(\eta,\mu)\,f(\mu)\,[s(\eta)-s(\mu)]\,d\mu,$$0, in agreement with integrable-system simulations at comparable Ursell number $$s(\eta)=s_0(\eta)+\int G(\eta,\mu)\,f(\mu)\,[s(\eta)-s(\mu)]\,d\mu,$$1 (Redor et al., 2019).

In deep water, inverse-scattering methods made controlled synthesis possible. A 148 m wave flume was driven by a signal constructed with the dressing algorithm of Gelash and Agafontsev, producing $$s(\eta)=s_0(\eta)+\int G(\eta,\mu)\,f(\mu)\,[s(\eta)-s(\mu)]\,d\mu,$$2-soliton solutions with eigenvalues chosen uniformly in a spectral rectangle and random norming phases. Nonlinear spectral analysis at successive gauges reconstructs the empirical DOS $$s(\eta)=s_0(\eta)+\int G(\eta,\mu)\,f(\mu)\,[s(\eta)-s(\mu)]\,d\mu,$$3. In the ideal 1D-NLSE this DOS would be invariant, but the experiment observed a slow deformation attributed to weak non-integrable effects such as third-order dispersion, induced mean flow, and Hilbert-transform terms described by the extended Dysthe equation (Suret et al., 2020).

The interaction of two soliton-gas “jets” provides a more directly hydrodynamic test. In a 140 m tank, two monochromatic spectral clouds centered at $$s(\eta)=s_0(\eta)+\int G(\eta,\mu)\,f(\mu)\,[s(\eta)-s(\mu)]\,d\mu,$$4 and $$s(\eta)=s_0(\eta)+\int G(\eta,\mu)\,f(\mu)\,[s(\eta)-s(\mu)]\,d\mu,$$5, with $$s(\eta)=s_0(\eta)+\int G(\eta,\mu)\,f(\mu)\,[s(\eta)-s(\mu)]\,d\mu,$$6, were realized as two $$s(\eta)=s_0(\eta)+\int G(\eta,\mu)\,f(\mu)\,[s(\eta)-s(\mu)]\,d\mu,$$7 soliton ensembles of opposite velocities. The reduced two-beam kinetic theory predicts explicit densities $$s(\eta)=s_0(\eta)+\int G(\eta,\mu)\,f(\mu)\,[s(\eta)-s(\mu)]\,d\mu,$$8 and velocities $$s(\eta)=s_0(\eta)+\int G(\eta,\mu)\,f(\mu)\,[s(\eta)-s(\mu)]\,d\mu,$$9 in the interaction region, controlled by

$u(\lambda;x,t)$00

For nine values of $u(\lambda;x,t)$01 between 0.2 and 0.9, the measured density and velocity changes lie on the theoretical curves within experimental error bars (Fache et al., 2023).

Field evidence in the deep open ocean comes from Taiwan waters. Among 20 523 records that satisfied quality control, unimodality, and deep-water criteria, eleven had $u(\lambda;x,t)$02, corresponding to an overall occurrence rate of $u(\lambda;x,t)$03. Because directional interference can artificially inflate the apparent soliton content, a probabilistic directional filtering method was applied. The mean $u(\lambda;x,t)$04 decreases from approximately $u(\lambda;x,t)$05 for the full spectrum to approximately $u(\lambda;x,t)$06 for $u(\lambda;x,t)$07 and approximately $u(\lambda;x,t)$08 for $u(\lambda;x,t)$09, yet three Eluanbi records remain above $u(\lambda;x,t)$10 after truncation, confirming genuine one-directional soliton-gas states (Lee et al., 6 Oct 2025).

A broader transition from weak turbulence to soliton gas has also been observed in two-dimensional random water waves. In a 27 m $u(\lambda;x,t)$11 30 m basin at depth $u(\lambda;x,t)$12 m, reducing the dispersion parameter by lowering the forcing frequency transforms a dispersive regime with energy concentrated on the linear dispersion surface into a shallow-water regime with a straight soliton ridge $u(\lambda;x,t)$13, heavier-than-Tayfun elevation statistics, and a nearly flat low-frequency spectrum followed by exponential decay. The reported onset of the soliton-dominated state occurs near $u(\lambda;x,t)$14, corresponding to $u(\lambda;x,t)$15, $u(\lambda;x,t)$16, and $u(\lambda;x,t)$17 (Leduque et al., 30 Sep 2025).

6. Extensions, non-integrable analogues, and current directions

Bidirectionality requires a refinement of the scalar theory because overtaking and head-on collisions need not have the same sign of phase shift. In integrable Eulerian systems this leads to two classes: isotropic gases, for which head-on and overtaking shifts have the same sign, and anisotropic gases, for which they have opposite signs. The resulting kinetic equations involve two densities of states $u(\lambda;x,t)$18 and $u(\lambda;x,t)$19 and two dressed velocities $u(\lambda;x,t)$20 and $u(\lambda;x,t)$21. Concrete realizations include the defocusing NLS gas, which is isotropic, and the resonant NLS or Kaup–Boussinesq gas, which is anisotropic (Congy et al., 2020).

In quasi-one-dimensional superfluids, a dark-soliton gas can be described by a phase-space kinetic theory closer to Vlasov than to inverse-scattering hydrodynamics. The microscopic Klimontovich density $u(\lambda;x,t)$22 leads, after ensemble averaging and neglect of correlations, to

$u(\lambda;x,t)$23

with a self-consistent pseudo-potential determined by the pair interaction Hamiltonian. Linearization around a homogeneous state predicts an acoustic branch $u(\lambda;x,t)$24; for a cold Fermi-like equilibrium the first-sound speed is

$u(\lambda;x,t)$25

and direct GP simulations reveal a narrow line below the Bogoliubov phonon band in agreement with this prediction (Pereira et al., 2020).

The term “soliton gas” has also been extended beyond integrable models. In the Schamel equation

$u(\lambda;x,t)$26

random ensembles of exact $u(\lambda;x,t)$27-type solitary waves display a robust gas-like phenomenology. In the unipolar case the statistics remain close to KdV behavior, with only weak radiative corrections, but in the bipolar case opposite-polarity collisions generate high-amplitude events and a steady increase in kurtosis from $u(\lambda;x,t)$28 to $u(\lambda;x,t)$29, unlike the nearly constant-kurtosis behavior of the integrable mKdV gas (Flamarion et al., 2023). A plausible implication is that some macroscopic signatures of soliton turbulence survive modest integrability breaking, while tail statistics can become more sensitive to inelastic focusing mechanisms.

Recent asymptotic work has broadened the model class. For the focusing Ablowitz–Ladik lattice, a dense soliton gas has been constructed as a large-$u(\lambda;x,t)$30 limit of $u(\lambda;x,t)$31-soliton data condensing on two imaginary intervals, with a Fredholm-determinant representation and explicit large-space and large-time asymptotics (Chen et al., 17 Mar 2026). For derivative NLS, the continuum limit of a reflectionless $u(\lambda;x,t)$32-soliton Riemann–Hilbert problem yields a $u(\lambda;x,t)$33-problem, an associated Fredholm determinant $u(\lambda;x,t)$34-function, and a kinetic equation for the effective group velocity of a trial soliton moving through the gas; the long-time asymptotics is stratified into one-phase, two-phase, and three-phase theta-function regions (Wang et al., 23 Jun 2026).

Within GHD, the KdV soliton gas has also acquired a full thermodynamic structure. The free-energy density and flux are written as

$u(\lambda;x,t)$35

with $u(\lambda;x,t)$36 for Maxwell–Boltzmann statistics, and the static covariance matrices of charges and currents are

$u(\lambda;x,t)$37

Numerical tests confirm these static predictions for KdV gases only with classical Maxwell–Boltzmann statistics, not with Fermi–Dirac or Bose–Einstein choices; the same framework also imports Euler-scale dynamical correlations, diffusion, and large-deviation formalisms into soliton-gas theory (Bonnemain et al., 2022).

Soliton gas has thus become a technical umbrella for several closely connected objects: kinetic equations for densities of states, thermodynamic limits of finite-gap spectra, generalized-hydrodynamic transport laws, experimentally realized incoherent soliton ensembles, and, more recently, continuum spectral limits with Fredholm-determinant and Riemann–Hilbert characterizations. Across these settings, the organizing principle remains the same: microscopic factorized scattering with calculable phase shifts determines macroscopic transport, statistics, and asymptotics.