{'fraction_non_alphanumeric': 0.035779905231602356, 'fraction_numerical': 0.013731747413209554, 'mean_word_length': 4.689958734525447, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical simulation of a solitonic gas in KdV and KdV-BBM equations

The collective behaviour of soliton ensembles (i.e. the solitonic gas) is studied using the methods of the direct numerical simulation. Traditionally this problem was addressed in the context of integrable models such as the celebrated KdV equation. We extend this analysis to non-integrable KdV–BBM type models. Some high resolution numerical results are presented in both integrable and nonintegrable cases. Moreover, the free surface elevation probability distribution is shown to be quasi-stationary. Finally, we employ the asymptotic methods along with the Monte–Carlo simulations in order to study quantitatively the dependence of some important statistical characteristics (such as the kurtosis and skewness) on the Stokes–Ursell number (which measures the relative importance of nonlinear effects compared to the dispersion) and also on the magnitude of the BBM term.


Introduction
Solitary wave solutions play the central rôle in various nonlinear sciences ranging from hydrodynamics to solid and plasma physics [44,34,30]. These solutions can propagate without changing its shape. However, the most intriguing part consists in how these solutions interact with each other. The binary interactions of solitary waves have been studied in the context of various nonlinear wave equations [44,23,39,25,9,6]. It is well known that in integrable models the collision of two solitons is elastic, i.e. they interact without emitting any radiation. In non-integrable models usually the interactions are nearly elastic [4].
The collective behaviour of soliton ensembles is much less understood nowadays. When a large number of solitary waves are considered simultaneously the researchers usually speak about the so-called solitonic turbulence or a solitonic gas. The literature on this topic is abundant. Some recent studies on solitonic gas turbulence in the KdV framework include [31,35,32]. The solitonic turbulence in nonintegrable NLS-type equations was studied in [47,12] and the authors showed that in conservative nonintegrable systems the solitonic gas is a statistical attractor whose dimension decreases with time. Recently it was shown both numerically and experimentally that solitonic ensembles appear in the laminar-turbulent transition in a fibre laser [41], modeled by a non-integrable nonlinear Schrödinger-type equation. However, the dominant number of studies is based on integrable models. This apparent contradiction motivated mainly our investigation to quantify the non-integrable effects onto the collective behaviour of solitons.
An approximate theoretical description of solitonic gases was proposed by V. Zakharov (1971) [45] using the kinetic theory. Later this research direction has been successfully pursued by G. El & A. Kamchatnov [14,15] who used the Inverse Scattering Technique (IST) [1] limited only to the integrable models. In this study the problem of solitonic gases will be investigated using the methods of direct numerical simulation. The evolution of random wave fields including solitonic gases was simulated numerically in [31,35,10] using symplectic, multi-symplectic and pseudo-spectral methods. However, previous investigators considered only a limited number of solitons (a few dozens) to simulate a solitonic gas. In this study we will adopt the pseudo-spectral method since it provides the high accuracy and computational efficiency necessary to handle large computational domains. Our goal will consist in: • investigate the influence of soliton interactions on statistical characteristics of the wave field, • construct the Probability Density Function (PDF) and compute the first four statistical moments of the solitonic turbulence, • study the rôle of non-integrable terms on the characteristics of soliton ensembles.
The present manuscript is organized as follows. In Section 2 we derive the governing equation used in this study and in Section 3 numerical results on a solitonic gas dynamics are presented. Finally, the main conclusions of this study are outlined in Section 4.

Mathematical model
As the starting point we choose the celebrated Korteweg-de Vries equation [22,24,20,31] (in dimensional variables) which models the undirectional propagation (here in the rightwards direction) of weakly nonlinear and weakly dispersive waves: where η(x, t) is the vertical excursion of the free surface above the still water level, h is the uniform undisturbed water depth and c = √ gh is the speed of linear gravity waves (g being the gravity acceleration).
The KdV equation (2.2) is known to be integrable [18,27]. However, the full water wave problem is known to be a non-integrable system, since the interaction of solitary waves is inelastic [7,16,8,9]. Moreover, we will modify the original equation (2.2) in order to include an additional dispersive term of the BBM-type [2]. Thus, the resulting KdV-BBM equation will not be integrable [17,11,21].
Consider the following scaled dependent and independent variables: where a 0 is the characteristic wave amplitude and l is the characteristic wavelength. In dimensionless variables KdV equation (2.2) reads: where parameter ε := a 0 /h measures the nonlinearity and µ 2 := h /l 2 is the dispersion parameter. The relative importance of these two effects is measured by the so-called Stokes-Ursell number [42] (sometime denoted as Ur): The last equation can be further simplified if we perform an additional change of variables: which yields the following simple equation including explicitly the Stokes-Ursell number S: 2) The last scaled KdV equation can be further generalized by using the low-order asymptotic relations in order to alternate higher-order terms as it was proposed by Bona & Smith (1976) [5] and Nwogu (1993) [28]. This step is rather standard and we do not provide here the details of the derivation:  [21]. Below we will study the solitonic gas behaviour 1 under the dynamics of the KdV (2.2) (δ = 0) and KdV-BBM (2.3) (δ = 0) equations. In the absence of the KdV term, we recover the celebrated BBM equation derived in [33,2].
Remark 1. We note that for a particular value of the Stokes-Ursell number S ≡ 1 another simpler scaling is possible when all the lengthes (x and η) are scaled by the mean water depth h. However, we do not adopt it in this study since below in Section 3.1 the dependence of some statistical characteristics on the Stokes-Ursell number S is investigated.

Linear well-posedness
In order to ensure the linear well-posedness of equation (2.3), the free parameter δ has to satisfy the following constraint δ ≥ 0. In the following we will consider only nonnegative values of this parameter. Recall that for δ = 0 we recover an scaled version of the classical KdV equation (2.2).

Ivariants
Assuming that the solution η(x, t) has either the compact support or decays sufficiently fast at the infinity along with its first derivative (η → 0, η x → 0 as x → ±∞), one can easily show that the following quantities are conserved [11]: In other words, I 1 (t) ≡ I 1 (0) and I 2 (t) ≡ I 2 (0), ∀t > 0. The invariant I 1 (t) is related to the mass conservation property, while the integral I 2 (t) can be assimilated to the generalized kinetic energy. The conservation of these quantities has not only the theoretical importance, but also the practical one. For example, it will allow us to check the global accuracy of the employed numerical scheme. We note also that the same invariants hold also in finite, but periodic domains (below we use periodic boundary conditions). The conservation laws to the BBM equation can be found in [29].

Solitary wave solutions
Equation (2.3) admits an exact localized (solitary) travelling wave solution which can be found analytically: The dependence of the solitary wave shape on parameters a, S and δ is shown on Figure 1. For instance, one can see that solitary waves become thiner when the amplitude a (and hence the speed) and/or the Stokes-Ursell number S are increased (see Figure 1(a,c)). On the other hand, the increase of the BBM coefficient δ leads to the growth of the tail (see Figure 1(b)). If both parameters S and δ are increased simultaneously, the 'thinning' effect of the Stokes-Ursell number dominates (see Figure 1(d)). The solitary waves interact elastically only in the integrable KdV case (δ ≡ 0) [26], while some dispersive tails are generated after the interaction in the general KdV-BBM model. We note that some authors came to the wrong conclusion about the elasticity of solitary wave interactions in the BBM equation [13] on the basis of low accuracy numerical simulations.

Numerical results
In order to solve numerically equation (2.3) we use a Fourier-type pseudo-spectral method with 3/2-antialiasing rule [40]. embedded adaptive 9(8) Runge-Kutta scheme [43]. The time step is chosen adaptively using the so-called H211b digital filter [37,38] to meet some prescribed error tolerance (generally of the order of machine precision ∼ 10 −15 ). The number of Fourier modes, the length of the computational domain and other numerical parameters are specified in Table 1.
In long time simulations presented below the invariants I 1,2 were conserved within the numerical accuracy 10 −11 and 10 −9 correspondingly. For the sake of illustration the error of the invariant I 1 (t) conservation in KdV and KdV-BBM simulations is shown on Figure 2. This accuracy is satisfactory to draw some robust conclusions on a solitonic gas statistical behaviour.
The initial condition for the KdV equation is composed of a finite number N s of solitons with random ampmlitudes a i ∼ N(a, σ), i = 1, . . . , N s separated by quasi-uniform distance ∆x s which is randomized to improve the ergodocity of the initial state: where ∆x s denotes the mean value reported in Table 1 and N(µ, σ) is the normal distribution with the mean µ and the variance σ. The solitonic gas initial state generated in this way is depicted on Figure 3(a). We use the same parameters for the initial solitonic gas state for the simulation with the KdV-BBM equation except for the domain size and the spacing between two solitons (see Table 1). They are larger in the KdV-BBM case since the soliton width increases with the parameter δ ≥ 0 (see Figure 1). Consequently, we can say that two initial conditions are approximatively isomorphic up to the horizontal coordinate stretching transformation. The simulation times T 1 (KdV) and T 2 (Kdv-BBM) are chosen so that an average soliton has enough time to go around the whole computational domain. The final states of the simulations are shown on Figures 3(b,c). One can notice how the initially quasi-uniform distribution of solitons is mixed by forming instantaneous soliton clusters as long as some void spaces due to the mass conservation property. The complete space-time dynamics simulated using the KdV and KdV-BBM equations is depicted on Figures 4(a,b). Individual lines correspond to solitons trajectories. The convergence of these lines corresponds to solitons collisions. It might appear on Figure 4 that collisions involve multiple solitons, however it is not the case. A zoom on a portion of the space-time domain is shown on Figure 5. One can see that interactions are only binary in agreement with [45]. It is interesting to estimate the total number of collisions in our simulation. The models under consideration are unidirectional. Thus, we can experience only overtaking collisions. Since the simulation time is chosen so that almost every soliton has enough time to travel across the whole computational domain, the number of collisions scales with O(m + m − ), where m ± is the number of solitons which travel with the speed above (below) the average. By construction of the initial condition we have m ± = O( 1 2 N s ). Consequently, the total number of collisions scales with O( 1 4 N 2 s ). Since the initial conditions are approximatively self-similar, after an appropriate rescaling of the spatial and time variables, the space-time dynamics is similar in both simulations (see Figures 4(a,b)). The difference between two simulations can be noticed if one makes a zoom on solitons background in order to see small radiating oscillations due to the inelasticity of collisions in the KdV-BBM case. This zoom on a portion of the computational domain is shown on Figure 6. In contrast, Figure 6(a) shows the absence of phonon modes in the KdV simulation.
It is custom to use the statistical methods to describe random wave fields [3,19]. The probability distribution of the normalized free surface elevation η 0 (x, t) := η(x, t) − η / η 2 1/2 at times t = 0 and t = T 1 is shown on Figure 7 (under the KdV dynamics). One can see that this distribution is quasi-stationary which is a direct consequence of the fact that solitons preserve perfectly their shape during the interactions. We note that the KdV-BBM numerical result shows the same invariance property since the inelasticity   Table 1 for the values of T 1,2 ). Parameter τ 0 denotes the time needed for an average soltion to go over the whole computational domain.   is too weak to modify significantly the probability distribution. Moreover, for comparison we plot also on the same Figure the standard normal (Gaussian) distribution. One can see that numerical results show much heavier tails than the standard distribution depicted on the same plot. The probability distribution can be characterized by several parameters. Perhaps two most important characteristics are listed hereinbelow: , which measures the heaviness of the spectrum tail.
These quantities are defined in terms of the normalized free surface elevation moments µ n := η n 0 . We note that κ = 3 and ς = 0 for the normal (Gaussian) distribution N(0, 1). The evolution of these quantities is shown on Figure 8. The kurtosis κ is shown on top panels (a,b) and the skewness ς on the bottom (c,d). The KdV simulation results are represented on the left images (a,d) and the KdV-BBM on the right (b,d). One can see that the qualitative behaviour of these quantities is similar in integrable and nonintegrable cases. After a rapid initial transient period both quantities κ and ς enter in a quasistationary regime which consists of fast small amplitude (±1.6%) oscillations around the mean value. Initial values of the moments are higher, but then it drops down quickly due to soliton interactions. We note also that corresponding values of κ and ς are slightly lower in the KdV-BBM case due to the differences in solitary wave shapes.

Estimation of statistical moments
In order to estimate efficiently kurtosis and skewness for various values of parameters S and δ without performing a series of direct numerical simulations we will adopt an approximate analytical method of statistical moments estimation employed also in [36].
Let us introduce the average density of a solitonic gas: where N S is the number of solitons and ℓ is the half-length of the computational domain (see also Table 1). We will assume that the solitonic gas is rarefied, i.e. ρ ≪ 1. Under this assumption we can represent approximatively the instantaneous free surface elevation η(x, t) as a linear superposition of distinct solitary waves (the interacting part is neglected): where {a i }, {c si } and {x i } are respectively the amplitudes, speeds and phase shifts of individual solitary waves. By assuming that the supports of solitons do not overlap, we can estimate the statistical moments of any order. For the sake of simplicity let us consider the first order µ 1 , i.e. the mean: where we took the limit ℓ → +∞ in order to compute analytically the integrals. Note also that the solitary wave parameter κ i is a function of the amplitude a i according to the relations given in (2.4).
Higher order moments µ n = η n , n > 1 can be computed in a similar way. In this study we will need the moments up to the fourth order: Using these moments we can estimate the skewness ς and kurtosis κ in the rarefied gas limit ρ → 0: In order to validate these asymptotic expressions we compare their predictions for the solitonic gas described in Section 3. The skewness ς and kurtosis κ are computed from numerical simulations and the time average is then taken. The results of the comparison are provided in Table 2. One can see that the simple analytical model presented in this Section is clearly able to describe rarefied solitonic gases. Moreover, one can infer from (3.1), (3.2) the behaviour of the skewness and kurtosis as the density of the gas decreases ρ → 0. However, as ρ → 0, the jump from the initial value of ς(0) or κ(0) to the stationary one will diminish, since the interactions between solitons become sparser.
The asymptotic formulas (3.1), (3.2) can be used to estimate the skewness ς and kurtosis κ for various values of model parameters δ and S in initial stages of the solitonic gas evolution. However, these formulas contain the statistical averages. The proposed Monte-Carlo approach is briefly summarized here. Namely, we generate M random independent realizations of a solitonic gas, each sample consisting of N s solitons. The numerical values of parameters M and N s used in simulations are given in Table 1. Then we use the assumption of the rarefied gas along with the knowledge of the analytical soliton shape in order to estimate some statistical quantities of the gas (see equations (3.1), (3.2), for example). A big number of Monte-Carlo realisations allows to annihilate local fluctuations and to obtain robust estimations.
The numerical results are presented on Figures 9 -11 where the shadowed areas show the statistical error due to Monte-Carlo simulations (±σ, ±1.96σ, where σ is the estimated variance). From these results one can clearly see that for the fixed density ρ an increase in the Stokes-Ursell number S leads to an increase in ς and κ (see Figure 9). The BBM coefficient δ has the antagonistic effect as it is illustrated on Figure 10. When both parameters S and δ are increased simultaneously 2 , the dependence of statistical quantities is not monotonic anymore (see Figure 11).
2 It means that we introduce an auxiliary homotopy parameter λ ∈ [0, 1] which parametrizes the simultaneous change of parameters S and δ in the following way: In the computations presented below we took the values: S min = 0.1, S max = 6.0, δ min = 0.0 and δ max = 4.0.

Conclusions and perspectives
In this study we presented several numerical experiments on the solitonic gas turbulence in the framework of an integrable KdV and a nonintegrable regularized KdV-BBM equation. The numerical results reported above generalize previous investigations [31,10] where only a limited number (a few dozens) of solitons were used to represent a solitonic gas. Consequently, we reduce the statistical error according to the law of large numbers.
First of all, we showed that the probability distribution for the solitonic gas remains quasi-invariant during the system evolution for both KdV and KdV-BBM cases. The special attention was payed to the statistical characteristics such as kurtosis and skewness which measure the 'heaviness' of tails and the asymmetry of the free surface elevation distribution. In particular, using the asymptotic methods and Monte-Carlo simulations we showed that both skewness and kurtosis increase with the Stokes-Ursell number S and decrease when the BBM term coefficient δ. When both parameters S and δ are increased gradually and simultaneously, these effects are in competition: first we observe the increase of these statistical characteristics, but then, this tendency is inversed and they decrease after reaching their respective maximal values (see Figure 11). We would like to underline that the proposed Monte-Carlo methodology is much less computationally expensive than direct numerical simulations. Despite the small number of Monte-Carlo runs (M = 100) the estimated statistical error is sufficiently small for the purposes of this study. On the other hand, this approach is restricted, strictly speaking, to the situations where the solitons are well separated, since it does not take into account superpositions of solitons.
The present study opens a number of perspectives for future investigations. More general nonlinearities could be included into the model along with some weak dissipative and forcing effects. This could allow us to observe Kolmogorov spectra of a solitonic gas [46]. The nonintegrable effects need some time to be accumulated. Consequently, even longer simulation times are needed. The interaction of a solitonic gas with a random radiation field has to be studied as well.