Statistical properties of extreme soliton collisions

The soliton turbulence is a challenging topic of recent research which is highly interesting from both, the mathematical and the physical points of view. In physics it describes strongly coherent wave states when the assumptions of a wave (weak) turbulence are violated, and thus the quasi-Gaussian probabilistic description completely fails. In mathematics, it represents deterministic chaos-like behaving systems which are formally solvable by means of the Inverse Scattering Transform, but its machinery is often too much complicated to write down explicit solutions or relations.

The kinetic theory for a soliton gas [1, 2] describes the transport of the soliton spectral density. Due to the failure of a linear superposition property, the kinetic theory cannot help to compute wave fields as such. In particular, the occurrence and probability of extreme events cannot be evaluated. It was shown in Refs. 3, 4 within the Korteweg–de Vries-type equations, that interactions of solitons of the same sign lead to some reduction of the wave field amplitude, whereas collisions of solitons of different signs (which can co-exist in the modified Korteweg – de Vries or the Gardner equations of the focusing types) can produce very high waves (so-called rogue waves). Meanwhile, the question of a meaningful description of the amplitude probability distribution in essentially coherent wave ensembles remains unclear.

It was pointed out recently [5, 6], that the standard definition of the variance

for a sign-defined wave field $u(x)\geq 0$ imposes a formal limit on the quantity which has the meaning of the soliton density, $\rho=N/L>0$, $\rho\leq\rho_{cr}$, where $N$ is the number of solitons per the spatial interval $L$. Here the overline means averaging over the interval $L$, which may be converted to the averaging in the spectral plane of the associated scattering problem to the given integrable equation when $L$ tends to infinity, see Ref. 5. The distribution of eigenvalues may provide information on a chance of large-amplitude wave generation [7].

Let us consider an example of the Korteweg – de Vries (KdV) equation

Its exact $N$-soliton solution $u_{N}(x,t)$ can be obtained via consecutive Darboux transformations [8] which allow a compact representation

Here $W(\cdot)$ denotes the Wronskian for $N$ “seed” functions $\psi_{2s-1}=\cosh{\theta_{2s-1}}$, $\psi_{2s}=\sinh{\theta_{2s}}$ for integer $s\geq 1$, where the phases are $\theta_{j}=k_{j}(x-V_{j}t-x_{j})$, $j=1,2,...,N$. The eigenvalues $k_{j}^{2}$ of the associated scattering problem specify the soliton amplitudes $A_{j}=2k_{j}^{2}$ and velocities $V_{j}=4k_{j}^{2}$, while the constants $x_{j}$ are responsible for the respective positions of solitons at a given time. The solution $u_{N}(x,t)$ is always positive [9].

The averaged quantities in the relation (1) can be calculated explicitly [10, 11] for the $N$-soliton solution (3) in the asymptotic limit,

where $x_{j}^{\pm}$ are reference locations of the solitons in the limits $t\to\pm\infty$. They read: $\overline{u}=4\rho\langle k\rangle$ and $\overline{u^{2}}=\frac{16}{3}\rho\langle k^{3}\rangle$, where the angle brackets denote averaging in the spectral domain. According to these relations, the critical soliton density reads [6]

Strictly speaking, solitons are not compact solutions, and hence the scale $L$ implied by the quantity $\rho$ is not well-defined in this illustrative example. Remarkably, the same expression (5) was obtained for interacting solitons in the thermodynamic limit state in Ref. 5, where the soliton density characteristic was introduced in a rigorous manner.

In the existing literature, most of researchers consider the limit of a small density which corresponds to the “rarefied” gas, when solitons interact in some sense weakly, and when approximate approaches can be efficient. In the present work we are not confined by the assumption of a small density of solitons. Quite the opposite, we consider the extreme soliton superposition, which occurs in synchronous soliton collisions, when the reference locations of all the solitons at $t=0$ coincide with the coordinate origin. This property can be formalized through the following symmetry condition, $u_{N}(-x,-t)=u_{N}(x,t)$; it is fulfilled when all the parameters $x_{j}$ in (3) are put equal to zero: $x_{j}=0$, $j=1,...,N$. An example of such an interaction is shown in Fig. 1.

In what follows, the soliton amplitudes are set decaying exponentially, so that they form a geometric series with the ratio $d>1$:

Within the KdV framework this choice is in fact rather generic as it corresponds to the distribution of eigenvalues of the scattering problem (represented by the stationary Schrödinger equation) for a parabolic potential, which may serve as a first approximation to any hump-like disturbance.

In the following Secs. III–V the statistical moments of synchronously interacting solitons are examined within the KdV framework. The applicability of the obtained results to a much broader class of equations is discussed in Sec. VI. We relate the fields of focused solitons to the states of the critical density in Sec. VII. Some final conclusions and perspectives are given in Sec. VIII.

The following integrals will be used to measure the statistical moments, similar to Ref. 10:

The choice of the soliton amplitude series (6) ensures that the asymptotic values $\mu_{n}(t\to\pm\infty)$ remain finite for any $n$ when $N\to\infty$.

The use of an ultra-high-precision procedure made it possible to compute the exact $N$-soliton solutions (3) and to calculate the statistical moments (7) for them, when $N$ is large, see details in Ref. 12. We found that the moments $\mu_{3}(t)$ and $\mu_{4}(t)$, which characterize the wave field asymmetry and kurtosis respectively, decrease in all situations when KdV solitons interact. This observation is in line with the limit of a two-soliton interaction [10] and with the direct numerical simulations of rarefied soliton gases (see e.g. Refs. 6, 13). Representative examples of the evolution of $\mu_{3}(t)$ and $\mu_{4}(t)$ are given in Fig. 2 for two choices of the parameter $d$ and different numbers of interacting solitons $N$. It follows from the figures that the curves converge to some limiting ones when $N\gg 1$. For the values of $d$ closer to $1$ the convergence occurs at larger numbers $N$.

When $d$ is sufficiently close to $1$, and $N$ is large, the statistical moments remain approximately constant within long time spans, when the solitons are located most densely. This effect is readily seen in the upper panels in Fig. 2 ($d=1.1$) for $N\gtrsim 10$, but vanishes for smaller numbers $N$ or for significantly larger values of $d$ (see the bottom panel in Fig. 2 for $d=1.6$). When $d$ is large, the soliton amplitudes in the sequence (6) decay quickly, so that apparently only a few of the largest solitons play a role. Obviously, the observation of “plateaux” in the dependencies $\mu_{3}(t)$, $\mu_{4}(t)$ should mean some degeneracy of the soliton state, which we explore below.

It is well-known that the integrability property of the KdV equation is deeply related to the existence of an infinite series of conserved quantities $I_{m}$,

Hereafter we will use the densities $T_{m}$ in the form given in Ref. 14. The two first statistical moments $\mu_{1}$ and $\mu_{2}$ are proportional to the integrals $I_{1}=6\int_{-\infty}^{\infty}{udx}$ and $I_{2}=18\int_{-\infty}^{\infty}{u^{2}dx}$, and thus do not change in time, while the higher-order statistical moments may vary.

The third conservation integral may be presented in the form

With the use of the asymptotic solution (4), the terms in (9) can be explicitly computed for $t\to\pm\infty$,

and then

It was noted in Ref. 10 that since $\mu_{3}$ and $\chi_{3}$ are strictly positive quantities, the decrease of $\mu_{3}$ in the course of the collision leads to the decrease of the second integral, $\chi_{3}$. This should correspond to “smoothing” of the solution in some sense. According to (11), the minimal possible value of $\mu_{3}$, $\mu_{3,min}/\mu_{3}(\infty)=3/4$ is achieved when $\chi_{3}$ totally vanishes. Note that this limit well agrees with the minimum values of $\mu_{3}$ observed in Fig. 2 (left column).

In order to elucidate the nature of the “smoothing” of the wave field, it is instructive to introduce the similarity parameter $\epsilon(t)$ (the so-called Ursell number), which estimates the ratio of the nonlinear term versus the term of the wave dispersion in the evolution equation (2) for a given solution $u(x,t)$ at a given instant of time. Then the conservation integrals (8) may be presented in the form

see Ref. 15. In these notations, the “smoothing” should correspond to the occurrence of very large values of $\epsilon$ (in other words, the small-dispersion limit) at $t\approx 0$. From this, we obtain the estimate for $\mu_{m}(0)$, and then the general formula for any integer $n\geq 1$ assuming $\epsilon\gg 1$:

To obtain the final expression in (13), the integrals $I_{n}$ and $\mu_{n}(\infty)$ should be calculated for the asymptotic solution (4):

The integrals $I_{n}$ for a single soliton may be found in Ref. 15.

In Fig. 3 the evolution of the Ursell number for the case $d=1.1$ is shown for different $N$. The Ursell number is calculated using the integral estimator $\epsilon=\int_{-\infty}^{\infty}{u^{3}dx}/\int_{-\infty}^{\infty}{(u_{x})^{2}dx}=\mu_{3}/{(2\chi_{3})}$. This quantity has the value $\epsilon=2$ for free solitons of any amplitude. One can see from the figure that indeed, the values of $\epsilon$ become very large when the solitons interact. They further grow if $N$ increases. A large number of solitons with close velocities leads to a long time of interaction, what explains the appearance of plateaux in the dependencies $\mu_{3}(t)$ and $\mu_{4}(t)$ in Fig. 2 (top). For $d=1.6$ the parameter $\epsilon$ increases considerably as well, but due to a shorter time of interaction a quasi-stationary behavior of the statistical moments is not observed (Fig. 2, bottom).

One can easily check that the relation (13) aligns with constant moments $\mu_{1}$ and $\mu_{2}$, and provides the limiting values $3/4$ and $1/2$ for the drops of the third and fourth moments, which agree with the numerical solutions shown in Fig. 2. It follows from the figure, that the estimates (13) are applicable to long time spans when the solitons focus, if $d$ is close to $1$ and $N\gg 1$. According to (13), statistical moments of any order decrease when KdV solitons interact, and this effect becomes stronger for higher moments. We have checked numerically up to the order $n=7$ that the estimate (13) is remarkably accurate, see Table 1.

Examples of the numerical solutions $u_{N}$ for $d$ close to $1$ and different finite number of solitons $N$ are shown in Fig. 4 for the instant $t=0$. When $N$ grows, the wave smoothing becomes apparent through a decrease in local maxima and an increase in the amplitude of the solution within the valleys between the humps. Accounting for a greater number of solitons in the series (6) leads to further smoothing of the solution and modification of its peripheral part. Presumably, in the limit of a large number $N$ of solitons with amplitudes distributed according to the geometric progression they focus into a visually one-hump profile for any $d>1$. A similar effect in reconstructing a box potential by solitons of the nonlinear Schrödinger equation was found in Ref. 16.

There is a well-known particular shape of the initial perturbation which gives purely discrete spectrum of the associated scattering problem for the KdV equation [17]:

Its evolution leads to the emergence of $N$ solitons with zero dispersive tail. The spectral parameters are given by the series

which yields the following relation between soliton amplitudes,

Thus, it is close to the sequence (6) with the ratio $d=1+2/N$, when $N$ is large and the soliton sequence numbers $j$ are small (i.e., for the tallest solitons).

For the initial-value problem specified by (16), the ratio for an arbitrary $n$-th statistical moment can be evaluated in a few steps using the asymptotic solution (4) and assuming that the number of solitons is large:

Thus, we obtain exactly the same result (13) provided by the qualitative description in the previous section.

Since the formulae (14) and (15) are general, the obtained result (13) is valid for any distribution of the soliton amplitudes, as soon as the solution does not change the sign and the similarity parameter $\epsilon$ is large; the considered example of a $\text{sech}^{2}$ perturbation (16) supports this conclusion. This dynamic scenario may be considered as a particular realization of a dense soliton gas. The discovered statistical property of ensembles of a large number of unipolar solitons actually has a much broader application beyond the KdV framework. Below we give a few examples of this fact.

First, the relation (13) is applicable to all members of the hierarchy of integrable KdV equations with the scattering problem represented by the stationary Schrödinger equation. Indeed, in this case the $N$-soliton solutions may be constructed using the same Darboux transformation [8]. The difference will be in the time-dependence of solitons only (i.e., other expressions for the soliton velocity $V_{j}(k_{j})$), but this difference vanishes at the focusing moment $t=0$ and has no effect on the statistical moments when the solitons separate at $t\to\pm\infty$. Consequently, the moments $\mu_{n}(0)$ and $\mu_{n}(\infty)$ will have exactly the same values as before. The obtained results obviously can be extended to higher dimensions, such as the Kadomtsev – Petviashvili equation.

The focusing modified Korteweg – de Vries (mKdV) equation

has a different from the classic KdV associated scattering problem, but as other integrable systems, possesses an infinite number of conserved quantities $I_{m}^{(mKdV)}$, $m=1,2,...$. Soliton solutions of the mKdV equation have two branches of solitons on a zero background which differ in a sign (polarity).

If the similarity parameter of the mKdV equation $\delta$, which may be estimated as $\delta^{2}=\int_{-\infty}^{\infty}{w^{4}dx}/\int_{-\infty}^{\infty}{w_{x}^{2}dx}$ becomes large (i.e., the small-dispersion regime takes place), then to the leading order the conservation integrals are proportional to the statistical moments, $I_{m}^{(mKdV)}\propto\mu_{2(m-1)}^{(mKdV)}\left(1+O(\delta^{-2})\right)$, $\mu_{n}^{(mKdV)}=\int_{-\infty}^{\infty}{w^{n}dx}$. The evolution of the parameter $\delta$ in two examples of interaction of mKdV solitons of the same sign is shown in Fig. 5 with red curves. A similar approach based on the Darboux transform is used to construct $N$-soliton solutions of the mKdV equation, see Ref. 4. The soliton amplitudes are specified according to the same law (6) as in the case of the KdV equation. During interaction, the parameter greatly exceeds the value $\delta^{2}=2$ which corresponds to free soliton solutions, if $d$ is close to $1$ and $N$ is sufficiently large.

Therefore, when the nonlinear term in the equation dominates over the term of dispersion, one can straightforwardly estimate the relative decrease of the statistical moments with even orders $n$. Remarkably, these estimates exactly coincide with the ones for the KdV equation (13). Even more, we have checked numerically that odd orders of the statistical moments satisfy the relation (13) too,

Note the perfect agreement between numerical results and the analytical estimates in Table 1.

In fact, equations (2) and (20) are related by the complex Miura transformation

which maps solutions of the mKdV equation $w(x,t)$ to complex-valued solutions $q(x,t)$ of the complex Korteweg – de Vries (cKdV) equation on $q(x,t)$ which has the form identical to (2) (see e.g. Ref. 18).

The balance of terms in the Miura transformation is controlled by the similarity parameter $\delta$ of the mKdV equation, which quantifies the contribution of the nonlinear term in (22) with respect to the imaginary term of dispersion, so that the real part effectively dominates at the focal point $t\approx 0$. That provides an estimate $\int_{-\infty}^{\infty}{q^{n}(x,t\approx 0)dx}\approx\int_{-\infty}^{\infty}{w^{2n}(x,t\approx 0)dx}$. At asymptotically large times $|t|\gg 1$ solitons of the mKdV equation $w_{j}(x,t)=2\lambda_{j}\text{sech}{\left[2\lambda_{j}(x-4\lambda_{j}^{2}t)\right]}$ are mapped to $q_{j}(x,t)=2\lambda_{j}^{2}\text{sech}^{2}{\left[\lambda_{j}(x-4\lambda_{j}^{2}t)-i\pi/4\right]}$, where real $\lambda_{j}$ is now the spectral parameter of the scattering problem for the mKdV equation. The expression for $q_{j}(x,t)$ coincides with the single soliton solution for the KdV equation (see (4) putting $k_{j}=\lambda_{j}$) except for the imaginary phase shift. The imaginary phase shift does not influence the integrals $\int_{-\infty}^{\infty}{q^{n}dx}$, they remain real-valued. Therefore, one can straightforwardly obtain that $\int_{-\infty}^{\infty}{q^{n}(x,t\to\pm\infty)dx}=2^{1-n}\int_{-\infty}^{\infty}{w^{2n}(x,t\to\pm\infty)dx}$. As a result, the statistical moments for the solution $q(x,t)$ of the complex KdV equation, $\mu_{n}^{(cKdV)}=\int_{-\infty}^{\infty}{q^{n}dx}$, may be also expressed through the ratio (21):

During the intervals of almost constant reduced statistical moments the interacting KdV solitons are located most closely and are characterized by large Ursell numbers. It seems reasonable to associate these extreme quasi-stationary states with the situation of the critical soliton density.

Importantly, the focused solitons represent essentially non-uniform in space wave fields, thus the adopted formulas, which relate the averaging in the physical space to the averaging in the spectral domain and eventually lead to the formula (5), become questionable. One can straightforwardly check that the definition (5) for the soliton amplitudes taken in the form of a geometric series (6) yield the values proportional to $N$, hence the quantity $\rho_{cr}$ grows infinitely if $N\to\infty$. Bearing this in mind, let us consider the quantity $l_{cr}=N/\rho_{cr}$ which has the meaning of some characteristic size, assuming that $\rho_{cr}$ is still determined by (5). For the exponential distribution of the soliton amplitudes (6) the values of $l_{cr}$ are finite:

To estimate the actual size of the spatial domain $l_{foc}$ occupied by the solution at $t\approx 0$, we consider the coordinate shift which solitons experience when pass through all the other solitons in the train (see Fig. 1). The shift for the $p$-th soliton, $\Delta x_{p}=|x_{p}^{+}-x_{p}^{-}|$, is given by the classic formula

Due to the factor $k_{p}^{-1}$, the shift of small-amplitude solitons is not bounded from above when $k_{p}\to 0$. On the other hand, the shift of the largest soliton $\Delta x_{1}$ is the least affected by an addition of small-amplitude solitons. Therefore we represent the characteristic size of the focused solution through the shift of the largest soliton, $l_{foc}=\Delta x_{1}$. The blue dashed lines in Fig. 1 correspond to the plots $(-\Delta x_{1}/2+V_{1}t,t)$ for $t<0$ and $(\Delta x_{1}/2+V_{1}t,t)$ for $t>0$. They follow the paths of the largest solitons when $|t|>>1$ and confine the interval $-l_{foc}/2\leq x\leq l_{foc}/2$ when cross the axis $t=0$. For the soliton amplitudes distributed according to (6), this quantity can be calculated assuming $N$ is large and $d$ is close to $1$,

The estimated size of the focused wave has the same dependence on the governing parameter $d$ as the critical size (24), with the proportionality coefficient $l_{cr}/l_{foc}\approx 3.65$.

For the shape (16) the characteristic size $l_{foc}$ of the focused solution can be estimated using (25) as $l_{foc}=\Delta x_{1}\approx 2\ln{2}$. One can also calculate the critical width $l_{cr}=N/\rho_{cr}$ for this example directly using the formula for the critical density (5) and the distribution of soliton amplitudes according to (17), $l_{cr}=3$; so for this case $l_{cr}/l_{foc}\approx 2.16$.

For the examples of multisoliton solutions $u_{50}$ and $u_{90}$ shown in Fig. 4 the characteristic sizes are $l_{cr}\approx 263$, $\Delta x_{1}\approx 104$ and $l_{cr}\approx 408$, $\Delta x_{1}\approx 126$ correspondingly. It may be seen from the figure that the scales $l_{foc}$ and $l_{cr}$ are in good compliance with the visible size of the focused wave train. One may speculate that the extreme compression of solitons which occurs in the considered scenarios of soliton collisions corresponds to the states with the minimum allowed (critical) size of the soliton ensemble. The scale $l_{foc}=\Delta x_{1}$ underestimates the actual visible size of the focused solution, and can be used as the estimate from below.

In this work we present a general idea that dense ensembles of solitons of the same sign (say, KdV-type solitons) can be considered as the strongly-nonlinear / small-dispersion wave states, what allows to express the statistical moments in terms of the spectral parameters of the associated scattering problem. Synchronous collisions of many solitons with slowly decaying amplitudes have been considered as a particular case when the dense soliton state can occur. We made a qualitative assumption (without a rigorous justification) on the relation between the critical soliton density and the minimal size of the focused solitons, which shows a reasonable agreement.

A broad applicability of the employed concept to other integrable systems is anticipated, we have given several confirming examples. Since weakly non-integrable systems typically inherit the soliton dynamics with little alteration (see e.g. Ref. 11), we may expect that the obtained results can be applied to non-integrable generalizations of the equations too.

At the same time, the approach apparently does not allow one to estimate the statistical moments in the case of simultaneous collisions of solitons with different signs, when they can cause wave amplification, which is not limited from above with an increase in the number of solitons [3, 4]. It follows from the direct numerical calculations that in the course of collision of mKdV solitons with alternating signs, the similarity parameter associated with the integral ratio, $\delta^{2}=\int_{-\infty}^{\infty}{w^{4}dx}/\int_{-\infty}^{\infty}{w_{x}^{2}dx}$, decreases from the value $\delta^{2}=2$ which corresponds to free soliton solutions, but not much. It never falls below the value $\delta=1$, see the blue curves in Fig. 5. Hence, a discrimination of terms which constitute the conservation integrals does not occur. Since phases of soliton solutions are not affected by their signs, the synchronous collision of solitons with alternating polarities formally leads to the creation of an equally dense soliton state, as in the case when solitons have the same signs. However, the statistical moments of the focusing solitons of different signs change rapidly, in contrast to the case of unipolar solitons.

In that way, a counterintuitive situation takes place, when interacting solitons of the same sign produce a strongly nonlinear state in terms of the similarity parameter, but which does not contain high waves. At the same time, the collision of solitons of different signs may cause extremely high wave amplitudes but is characterized by a relatively small ratio of nonlinearity versus dispersion.

All authors contributed equally to this work.

Data available on request from the authors.