{'fraction_non_alphanumeric': 0.03802770327368924, 'fraction_numerical': 0.018370871146709777, 'mean_word_length': 4.608489594065526, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Macroscopic dynamics of incoherent soliton ensembles: soliton-gas kinetics and direct numerical modeling

We undertake a detailed comparison of the results of direct numerical simulations of the integrable soliton gas dynamics with the analytical predictions inferred from the exact solutions of the relevant kinetic equation for solitons. We use the KdV soliton gas as a simplest analytically accessible model yielding major insight into the general properties of soliton gases in integrable systems. Two model problems are considered: (i) the propagation of a `trial' soliton through a one-component `cold' soliton gas consisting of randomly distributed solitons of approximately the same amplitude; and (ii) collision of two cold soliton gases of different amplitudes (soliton gas shock tube problem) leading to the formation of an incoherend dispersive shock wave. In both cases excellent agreement is observed between the analytical predictions of the soliton gas kinetics and the direct numerical simulations. Our results confirm relevance of the kinetic equation for solitons as a quantitatively accurate model for macroscopic non-equilibrium dynamics of incoherent soliton ensembles.


Introduction
Dynamics of incoherent nonlinear dispersive waves have been the subject of very active research in nonlinear physics for several decades, most notably in the contexts of ocean wave dynamics and nonlinear optics (see e.g. [16,20,23]). Two major areas where statistical properties of random ensembles of nonlinear waves play essential role are wave turbulence and rogue wave studies (see [15,19] and references therein).
A very recent direction in the statistical theory of nonlinear dispersive waves introduced by V.E. Zakharov (2009) [34] is turbulence in integrable systems. It was suggested in [34] that many questions pertinent to a turbulent motion can be formulated in the framework of completely integrable systems. Such an 'integrable turbulence' theory has two natural premises: (i) nonlinear wave phenomena are often so complex that they must be described in statistical terms; (ii) integrable systems capture essential properties of many wave processes occurring in the real-world systems. Physical relevance of integrable turbulence theory was recently demonstrated in the fibre optics experiments [24,31].
Solitons play the key role in the characterisation of nonlinear wave fields in dispersive media, therefore the theory of soliton gases in integrable comprises an important part of the general theory of integrable turbulence [34]. The very recent observations of dense statistical ensembles of solitons in shallow water wind waves in the ocean [2] well modelled by the KdV equation and in the laminar/turbulent transition in fibre lasers [29] described by the defocusing Nonlinear Schrödinger (NLS) equation provide further physical motivation for the development of the theory of soliton gas/soliton turbulence in integrable systems. In this paper we shall be using the KdV soliton gas as a simplest analytically accessible model yielding a major insight into the general properties of soliton gases in integrable systems.
Macroscopic dynamics of a Korteweg-de Vries (KdV) soliton gas are determined by the fundamental 'microscopic' properties of two-soliton interactions [3]: (i) soliton collisions are elastic, i.e. the interaction does not change the soliton amplitudes (or, more precisely, the discrete spectrum levels in the associated linear spectral problem for the quantum-mechanical Schrödinger operator); (ii) after the interaction, each soliton acquires an additional phase shift; (iii) the total phase shift of a 'trial ' soliton acquired during a certain time interval can be calculated as a sum of the 'elementary' phase shifts in pairwise collisions of this soliton with other solitons during this time interval. Thus, the macroscopic dynamics of a soliton gas are essentially determined by two-soliton interactions. This fact enabled Zakharov in 1971 to introduce the kinetic equation for a 'diluted ', small-density, gas of solitons for the KdV equation. The generalisation of Zakharov's equation to finite densities derived in [6] has required the consideration of the thermodynamic-type limits for finite-gap potentials and the associated Whitham modulation equations [11,32]. A straightforward, physical derivation of the kinetic equation was made in [8]. We stress that kinetic description of soliton gas makes an emphasis on the particle-like nature of solitons. At the same time, solitons represent nonlinear coherent wave structures so the total random nonlinear wave field associated with a soliton gas can be naturally interpreted as soliton turbulence [10]. The effect of two-soliton collisions on the properties of the statistical moments of soliton turbulence was studied in [22]. An effective method for the numerical computation of soliton gas was developed in [4] and was applied to the numerical modelling of soliton gas in the KdV and the KdV-BBM (Benjamin-Bona-Mahoni) equations in [5].
The kinetic equation for solitons derived in [6,8] was shown in [9,21] to possess some remarkable mathematical properties. In particular, it was shown that it has an infinite number of integrable hydrodynamic reductions which is a strong indication of integrability of the full kinetic equation. Integrability of hydrodynamic reductions opens a broad perspective for obtaining various exact solutions to the kinetic equation. However, the quantitative confirmation of the relevance of the predictions of the soliton kinetic theory to the actual macroscopic dynamics of soliton gases still remains an open problem. Indeed, as was already mentioned, the formal derivation of the kinetic equation involves certain singular limiting transition of the thermodynamic type. The mathematical conditions required for this transition are not necessarily applicable to physically (or even numerically) accessible soliton systems. This is why it is vitally important to have a direct numerical confirmation of the validity of the kinetic equation. The main goal of the pesent paper is thus to test the relevance of the soliton gas kinetics to the actual 'particle dynamics' of incoherent soliton ensembles. With this aim in view we compare some model exact solutions of the kinetic equation with the results of the high accuracy direct numerical simulations of the KdV soliton gas.
The paper is organised as follows. Section 2 provides a brief account on the kinetic equation for solitons. In Section 3 hydrodynamic reductions of the kinetic equation are considered for one-and two-component soliton gases and two model problems are considered: the propagation of a trial soliton through a one-complonent cold soliton gas and the collision of two cold gases (the soliton gas shock tube problem). Section 4 is devoted to the direct numerical modelling of the problems considered in Section 2 and detailed comparisons of the analytical and numerical solutions. The main conclusions of this study are outlined in Section 5.

Kinetic equation for a soliton gas
We consider the KdV equation in the canonical form We introduce soliton gas as an infinite collection of KdV solitons randomly distributed on the line with non-zero density. This intuitive definition lacks precision but it is sufficient for the purposes of this paper. A mathematically consistent definition of a soliton gas as the thermodynamic limit of finite-gap potentials will be briefly described below (see [10] for details). Let each soliton in the gas be 'labeled ' by the spectral parameter η i ⩾ 0 so that λ i = −η 2 i is the corresponding discrete eigenvalue in the spectral problem for the linear Schrödinger operator associated with (2.1) in the inverse scattering transform (IST) formalism. We assume that the discrete values η i are distributed with certain density on some finite interval, say [0, 1] and replace η i by the continuous variable η ∈ [0, 1] (see the details in [6,9,10]). As is well known (see e.g. [3]) the amplitude of an isolated KdV soliton with the spectral parameter η is a = 2η 2 and its speed is S = 4η 2 . It is clear that in the soliton gas, the mean speed of the soliton with the same parameter η will differ from 4η 2 . Indeed, due to the pairwise collisions with other solitons (each leading to a 'phase-shift') the distance covered by this 'trial ' soliton over some time interval ∆t ≫ 1 will be different from 4η 2 ∆t. This intuitive reasoning was used in the original Zakharov paper [33] for the derivation of the approximate kinetic equation describing macroscopic dynamics of a 'rarefied soliton gas'. The full, non-perturbative equation for a 'dense' gas of KdV solitons was derived in [6] by considering a singular, thermodynamic type limiting transition for the modulation Whitham equations describing slow evolution of the finite-gap solutions of the KdV equation. The key property of the thermodynamic limit is the special spectral band-gap distribution (scaling) that preserves finiteness of the integrated density of states [14,17] in the infinite-band limit. It was shown in [6,10] that in the thermodynamic limit the integrated density of states yields the spectral measure f (η) dη of the soliton gas. In a spatially inhomogeneous soliton gas f (η) ≡ f (η; x, t) so that f (η 0 ; x, t) dη dx is the number of solitons with the spectral parameter η ∈ (η 0 , η 0 + dη) and located in the spatial interval (x, x + dx) at the moment t. The integral is the total physical (as opposed to spectral) density of the soliton gas, i.e. the number of solitons per unit length. For the KdV equation (2.1) the evolution of the spectral density f (η, x, t) is described by the integro-differential kinetic equation [6] where we use the shorthand notations f (η) ≡ f (η, x, t) and s(η) ≡ s(η, x, t), the latter having the meaning of the mean soliton gas velocity (or you can view it as the velocity of a 'trial ' soliton with the spectral parameter η placed in the soliton gas characterised by the distribution function f (µ)). Zakharov's approximate kinetic equation for a rarified soliton gas [33] is obtained from (2.3), (2.4) by assuming κ ≪ 1 and retaining only the first order correction in (2.4). It is important to stress that the typical scale of variations of x and t in the kinetic equation (2.3) is much larger than in the KdV equation (2.1), governing the primitive, 'microscopic' evolution. It was noted in [8] that the formal procedure of the thermodynamic limit of the Whitham equations can be viewed as the justification of a simple straightforward derivation of the kinetic equations for solitons using the original Zakharov's reasoning [33] based on the phase shift expressions for two-soliton collisions. The key difference is that one should use the full, unapproximated mean velocity function s(η) (rather than unperturbed free soliton velocity 4η 2 ) in the calculation of the correction to the velocity of free soliton which leads to a self-consistent definition of the local mean velocity of solitons with the spectral parameter η.
In the kinetic description (2.3), (2.4) of a soliton gas the solitons are viewed as particles moving with certain speeds and interacting with each other according to the nonocal closure equation (2.4). On the other hand, one is also interested in the nonlinear wave field u(x, t) associated with the soliton gas dynamics (an integrable soliton turbulence). As is well known (see e.g. [18,19]) a turbulent wave field is usually characterised by the moments ⟨u n ⟩ over the statistical ensemble which, due to ergodicity of soliton turbulence [10], can be computed as spatial averages It was shown in [6,7,10] that the two first moments in the KdV soliton turbulence are calculated in terms of the spectral distribution function f (η, x, t) as A fundamental restriction imposed on the distribution function f (η) follows from nonnegativity of the variance The consequences of this restriction have been explored in [7]. Here it will inform the choice of the soliton gas parameters for the numerical modelling.

Hydrodynamic reductions and exact solutions
To get a better insight into the properties of the soliton gas dynamics we consider the hydrodynamic reductions of the kinetic equation (2.3), (2.4). Such hydrodynamic reductions enable one to derive some simple, physically interesting exact solutions which could then be compared with the results of direct numerical modelling of the KdV equation.
The family of the simplest N-component (N-beam) 'cold gas' reductions is selected by the delta-function multiflow ansatz, well known in plasma physics [26]: Physically, the i-th complonent of the soliton gas described by the distribution (3.1) consists of an infinite sequence of nearly identical solitons having the spectral parameter η distributed in a narrow ε-vicinitty of η = η i such that ε η i ≪ 1, and distributed by Poisson on x ∈ (−∞, ∞) with the density f i (x, t) which can slowly vary in space in time [7,21]. The hydrodynamic reductions obtained by (3.1) for arbitrary N have been thoroughly analysed in [9,21] where their integrability in the sense of the generalised hodograph transform was proven. Here we shall be mostly looking at the case of a two-component gas yielding the simplest nontrivial results that can be verified numerically. We consider two model problems: (i) the propagation of a 'trial ' soliton through a one-component 'cold ' soliton gas and (ii) collision of two one-component soliton gases -the shock tube problem.

Propagation of a trial soliton through one-component soliton gas
We consider a 'trial ' soliton with the spectral parameter η = η 1 moving through a onecomponent soliton gas with the distribution function where the density f 0 (x, t) is found by the substitution of the distribution (3.2) into (2.3), (2.4), yielding the linear transport equation Substituting the distribution function (3.2) into the expressions (2.5) for the moments we obtain from (2.6) the restriction for the soliton gas density [7] The mean velocity of the trial soliton s(η 1 ; x, t) can then be found from formula (2.4), where s 1 (x, t) ≡ s(η 1 ; x, t)). Expressing s 1 from (3.5) we obtain is the classical phase-shift expression for the two-soliton collision [3].

Soliton gas shock tube problem
We now consider a two-component soliton gas by introducing the distribution function in the form (3.9) Substitution of (3.8) into (2.3), (2.4) leads to the system of two conservation laws 10) where s 1,2 (x, t) ≡ s(η 1,2 , x, t) and It follows from (3.9), (3.11), (3.12) that if s 1 > 4η 2 1 then s 1 > s 2 and s 2 < 4η 2 2 . Similarly, if s 1 < 4η 2 1 then s 1 < s 2 and s 2 > 4η 2 2 . We note that, in the degenerate case η 1 = η 2 ≡ η 0 the ansatz (3.8) reduces to the onecomponent distribution (3.2) with f 0 = f 1 + f 2 as expected. Expressing s 1,2 in terms of f 1,2 from (3.11) we obtain Note that by setting f 1 ≡ 0 we recover the expression (3.6) for the speed of the trial soliton s 1 (to establish the correspondence with (3.6) the index '2' in the first expression (3.13) should be replaced with 0, also α 2 becomes α). The density of the two-component soliton gas (2.2) is It is not difficult to show that equations (3.10), (3.11) assume Riemann invariant form in variables s 1 , s 2 [8], System (3.14) is linearly degenerate [25] which implies: (i) the absence of the nonlinear wave-breaking effects in a two-component soliton gas (see [9] for the relevant account of the properties of linearly degenerate hydrodynamic type systems) and (ii) unavailability of simple-wave solutions (indeed one can easily see that s 2 (s 1 ) implies that s 1,2 are constants).
The component densities f 1,2 (x, t) are expressed in terms of the velocities s 1,2 by the expressions derived from (3.11) Using the distribution function (3.8) we obtain for the two first moments (2.5) of the wave field in the two-component gas: : We now consider the Riemann problem for the two-component soliton gas characterized by the spectral distribution function (3.8) corresponding to the shock tube problem: an initial contact dscontinuity separating gases of different density where f 10 , f 20 > 0 are some constants satisfying f i0 ⩽ η i 3 (see (3.4)). We also assume that η 1 > η 2 . Note that, unlike in the classical gas-dynamics shock tube problem, the initial velocity of the soliton gases is not zero but is fully determined by the density distribution (3.18) via relations (3.13), We note that the values of s 2 for x < 0 and s 1 for x > 0 corresponds to the respective zero-density components (cf. (3.18)) and thus, are fictitious initial parameters providing consistency with the equation to be solved. Since the governing equations (3.14) are quasilinear the solution of the Riemann problem must depend on x t alone. Since due to linear degeneracy system (3.14) does not have nonconstant simple wave solutions one has to look for weak similarity solutions of the original conservation laws (3.10). The required solution represents a combination of three constant states separated by two contact discontinuities (see [8]). For the total density κ = f 1 + f 2 we have (3.20) The schematic representation of this solution is shown in Figure 1. The values f 1c and f 2c as well as the velocities c ± of the discontinuities are found from the Rankine-Hugoniot jump conditions Here c − and c + are the velocities of the left and right discontinuities respectively, and f 1c , f 2c and s 1c , s 2c are the densities and velocities of the soliton gas components in the interaction region x ∈ [c − t, c + t]. The velocities s 1c and s 2c are expressed in terms of f 1c , f 2c by relations (3.13). Solving (3.21) and (3.22) we obtain: (3.23) The speeds c ± of the boundaries of the interaction region are given by The expanding interaction region in the soliton gas shock tube problem can be viewed as an incoherent dispersive shock wave, a stochastic counterpart of the traditional, coherent dispersive shock wave (DSW) forming due to a dispersive regularisation of the Riemann initial data in the KdV equation [13] (we need to make a clear disctinction between the studied here incoherent DSWs, which are generated in the Riemann problems for soliton gases and the incoherent DSWs recently observed in Fourier spectra evolution of random waves [12]). In contrast with the coherent DSWs, the incoherent DSW generated in the collision of two soliton gases does not have a disctinct structure of a slowly modulated wavetrain but is characterised by the increased itensity of fluctuations A 2 (2.6) of the random nonlinear wave field, compared to the values of A 2 in the colliding soliton gases at t = 0. It is not difficult to show that the density of the two-component soliton gas in the incoherent DSW region κ c = f 1c + f 2c > f 10 , f 20 but κ c < f 10 + f 20 (see Figure 1).
In conclusion of this section we note that the above results for a two-component gas cannot be in principle derived from the approximate kinetic equation obtained in [33] for a rarefied gas. Indeed, it was shown in [33] that the kinetic equation for a rarified soliton gas prescribes linear instability for a two-component gas so that any density perturbation would grow without bound making the underlying equation inapplicable. The full kinetic equation used here yields hyperbolic hydrodynamic reduction (3.14) for a two-component gas implying stability for the full range of admissible densities. This conclusion will be confirmed by direct numerical simulations in the next section.

Numerical experiments
In this section we perform direct numerical simulations of the KdV soliton gas and compare the numerical solutions with the corresponding solutions of the kinetic equation obtained in Section 3.

Numerical method
In order to solve numerically the KdV equation we employ the standard pseudo-spectral Fourier collocation technique [1,28]. This method is briefly explained below.
Denote byû(k, t) = F {u} the Fourier transform of u(x, t) in x, where k is the wavenumber. Then, by Fourier-transforming the KdV equation (2.1) yieldŝ u t − ik 3û = −3ik (u 2 ). (4.1) The most computationally efficient way consists in computing the spatial derivatives in spectral space while the nonlinear product is computed in real space and de-aliased using the classical 3 2th rule. The overall implementation is very efficient thanks to the FFT algorithm. In order to improve the time-stepping we will use the so-called integrating factor technique. This consists of the exact integration of the linear terms of (4.1), viz.
where the linear and nonlinear operators L and N are defined through their symbols as This allows to increase substantially the accuracy and the stability region of the time marching scheme (see, for example [28]). Finally, the resulting system of ODEs is discretized in time by the Verner's embedded adaptive 9(8) Runge-Kutta scheme [30]. The time step is adapted automatically according to the H211b digital filter approach [27]. In order to estimate the accuracy of the numerical solution, one can follow the values of quantities which are conserved during the evolution. For the KdV equation we check the values of the first three invariants during all numerical computations. A typical evolution of these invariants is represented in Figure 2.

Test 1: Propagation of a trial soliton through a one-component gas
We now present the results of the numerical simulation of the propagation of a trial soliton with given spectral parameter η = η 1 through the one-component soliton gas with η = η 0 to compare with the theoretical results of Section 3.1. The value for the comparison is the mean (i.e. averaged over a large interval) velocity of the trial soliton is given by formula (3.6) in which f 0 = const.
In the simulations, the initial condition for the one-component KdV soliton gas is composed of finite but sufficiently large number M of solitons (in our experiments M = 200) with the random amplitude a = 2η 2 chosen from the normal distribution of the spectrum η with mean η 0 and fixed standard deviation σ = 2 × 10 −2 , separated by a space lag ∆ 0 whose  value is directly related to the gas density: The parameters ℓ and ǫ i are respectively the starting point of the random lattice and a random (uniform) perturbation to the i-th soliton position, taken in the interval ǫ i ∈ [−1, 1]. So that, in order to increase (or decrease) the density κ 0 , it is only required to change the value of the space lag ∆ 0 .
With an added trial η 1 -soliton the initial-boundary conditions for the KdV equation (2.1) assume the form The snapshots of the trial soliton evolution are shown in Figure 3. One can see that the trial soliton undergoes a noticeable acceleration as predicted by the theory. The quantitative comparisons of the numerically found values for the averaged speed of the trial soliton with the formula (3.5) are shown in Figure 4 for three different sets of parameters of the soliton gas.
The comparisons show an excellent agreement between the results of direct numerical simulations and the predictions of the kinetic theory. In all simulations the condition (3.4) for the soliton gas density is satisfied.

Test 2: Soliton gas shock tube problem
Following the strategy proposed in the previous section, we build the initial condition as a superposition of two distinct populations of solitons separated at t ⋆ = 0 by an empty gap, so that with α represent the single gas component (α = 1, 2). As in the pevious case the amplitudes of the two gas components are Gaussian random values distributed with the means η 1 and η 2 and standard deviations σ 1 = 10 −4 and σ 2 = 2×10 −2 respectively. Again, the respective densities f 10 and f 20 can be easily changed by tuning the parameters ∆ 1 and ∆ 2 in equation (4.3). The numerical solution of the KdV equation with initial condition (4.3) is presented in Figure 5. We now perform the comparison of the parameters of this numerical solution with the weak analytical solution of the soliton gas shock tube problem obtained in Section 3.1. Specifically, we are interested in comparing the total density of solitons κ = f 1c + f 2c in the interaction (incoherent DSW) region and in the speeds c ± of its edges.
The first observation is that the incoherent DSW forming due to the interaction of two cold soliton gases is stable in agreement with the hyperbolic nature of the two-component hydrodynamic reduction (3.14) of the full kinetic equation. The comparisons for the total density as the function of time in the interaction region is presented in Figure 6. One can see three distinct regions in the presented numerical plots. The value of the total density is initially equal to the sum f 10 + f 20 of the component densities and then decreases through the equilibration process to the stationary value κ c (highlighted in all threee plots) which is in excellent agreement with the predictions of the theory based on the weak solution (3.20)  of the two-component hydrodynamic reduction of the kinetic equation. The subsequent decrease of the density seen in the numerical plots is due to an inherent restriction of the numerical experiment involving finite number of solitons so the interaction region is sustained only for a finite interval of time.
The comparisons for the velocities of the edges of the interaction region (incoherent DSW) is presented in Figure 7 and demonstrate excellent agreement between the analytical and numerical results. We shall also compare the values of the two first moments (3.16), which for the interaction region assume the form (u) c = 4 (η 1 f 1c + η 2 f 2c ); (u 2 ) c = 16 3 (η 3 1 f 1c + η 3 2 f 2c ) . where f 1c and f 2c are determined in terms of the initial data by formulae (3.23). The results of the comparison are presented in Fig. 8b. Again, the excellent agreement is observed. One can also see that the condition (3.17), A 2 c = (u 2 ) c − (u) 2 c > 0, is satisfied.

Conclusions
We have undertaken a detailed comparison of the macroscopic dynamics of the KdV soliton gas predicted by the kinetic equation for solitons with the results of direct numerical simulations of the KdV equation. The simulations involved 200 solitons enabling an accurate determination of macroscopic parameters of the soliton gas. Two test problems have been considered: the propagation of a trial soliton through a one-component 'cold ' soliton gas and the shock tube problem involving the interaction of two cold gases with different parameters leading to the formation of an incoherent dispersive shock wave. In both cases the excellent agreement between the asymptotic analytical predictions of the kinetic equation and the direct KdV 'molecular dynamics' numerical simulations has been observed. This confirms validity of the kinetic equation for solitons as a quantitatively accurate model for the description of non-equilibrium dynamics of soliton gases in integrable systems. The challenging problem is now to study the structure and evolution of the definitive statistical characteristics of integrable soliton turbulence (PDF, power spectrum density etc.). This will be the subject of future work.