A Graphs of the evolution of wave-action spectrum and spatial correlation function.

# Integrable turbulence and formation of rogue waves

## Abstract

In the framework of the focusing Nonlinear Schrödinger (NLS) equation we study numerically the nonlinear stage of the modulation instability (MI) of the condensate. The development of the MI leads to formation of “integrable turbulence” [Zakharov V. E., Stud. Appl. Math. 122, 219-234 (2009)]. We study the time evolution of it’s major characteristics averaged across realizations of initial data – the condensate solution seeded by small random noise with fixed statistical properties.

We observe that the system asymptotically approaches to the stationary integrable turbulence, however this is a long process. During this process momenta, as well as kinetic and potential energies, oscillate around their asymptotic values. The amplitudes of these oscillations decay with time as , the phases contain the nonlinear phase shift that decays as , and the frequency of the oscillations is equal to the double maximum growth rate of the MI. The evolution of wave-action spectrum is also oscillatory, and characterized by formation of power-law region in the small vicinity of the zeroth harmonic with exponent close to 2/3. The corresponding modes form “quasi-condensate”, that acquire very significant wave action and macroscopic potential energy.

The probability density function (PDF) of wave amplitudes asymptotically approaches to Rayleigh distribution in oscillatory way. Nevertheless, in the beginning of the nonlinear stage the MI slightly increases the occurrence of rogue waves. This takes place at the moments of potential energy modulus minima, where the PDF acquires “fat tales” and the probability of rogue waves occurrence is by about two times larger than in the asymptotic stationary state.

Presented facts need theoretical explanation.

## I Introduction.

Today the total amount of experimental evidence of rogue waves emergence on the surface of fluid and in optical fibers is huge [[1], [2], [3], [4], [5]]. Thus, the development of a consistent theory of these events is urgently needed. The simplest nonlinear mathematical model for the description of rogue waves phenomenon is the modulation instability (MI) developing from the condensate solution in the framework of the focusing one-dimensional Nonlinear Schrodinger (NLS) equation [[1], [2], [3]]. Without loss of generality we will use the NLS equation in following form:

 iΨt−Ψ+Ψxx+|Ψ|2Ψ=0. (1)

The simplest ”condensate” solution of this equation is unstable. If we consider modulations to the condensate as

 Ψ=1+κexp(ikx+iΩt),|κ|≪1, (2)

and linearize Eq. (1) against the condensate, we obtain

 Ω2=k4−2k2. (3)

The modulations with turn out to be unstable, and the maximum growth rate of the instability,

 γ0=maxkImΩ=1, (4)

is realized at . Thus, the characteristic length of the instability is , and the characteristic time is .

To study the nonlinear stage of the MI, one has to solve Eq. (1) with the initial data in the form

 Ψ|t=0=1+ϵ(x),|ϵ(x)|≪1. (5)

It should be noted that the problem of the MI development on the background of any condensate solution , , and for the focusing NLS equation

 iΨt+BΨxx+G|Ψ|2Ψ=0

with arbitrary dispersion and nonlinearity coefficients renormalizes to Eqs. (1), (5), as can be seen after the scaling and gauge transformations , , and .

If at , then the problem can be solved analytically with the help of the inverse scattering transformation [[6], [7], [8], [9]]. The MI in this case leads to formation of different types of solitonic solutions. The scenario of the MI essentially depends on the fine details of the initial perturbation . The pure real perturbation leads to formation of homoclinic solutions of Peregrine type [[10]], while the pure imaginary one generates ”superregular” solitonic solutions described in the publications of V.E. Zakharov and A.A. Gelash [[8], [9]]. The general case when both imaginary and real parts of the perturbation are present is not properly studied yet.

In spite of the apparent significance of these results, they do not answer to the main question – what happens if the perturbation is not localized? To study it, we solve the NLS equation (1) numerically in the box with periodic boundary. Theoretically speaking, this problem can also be solved analytically. Any periodic solution of the NLS equation can be expressed explicitly in terms of Jacobi theta-functions over a certain hyperbolic curve [[7]]. However, this beautiful mathematical result can hardly be used for practical purposes. In our numerical experiments is a small random noise, and we use the number of harmonics of order . Then, to model it’s evolution in terms of Jacobi functions, we have to make the genus of the curve of order . It is unrealistic so far to follow this evolution by the use of the exact analytical methods.

Therefore, we rely completely on numerical experiments. We use integrability of the NLS equation only in the weakest sense. Integrability implies conservation of infinite number of integrals of motion. The first three of these invariants are wave action,

 N=1L∫L/2−L/2|Ψ(x,t)|2dx, (6)

momentum,

 P=i2L∫L/2−L/2(Ψ∗xΨ−ΨxΨ∗)dx, (7)

and total energy,

 E=H2+H4,Hd=1L∫L/2−L/2|Ψx|2dx,H4=−12L∫L/2−L/2|Ψ|4dx. (8)

Here is kinetic and is potential energy. We define these integrals with the prefactor for further convenience. We use method of numerical simulations that conserves very well the first 12 invariants.

Our study has two main goals. First, we expect that after a very long evolution the result of the MI of the condensate should be the stationary “integrable turbulence” [[11]] – thermodynamically equilibrium state defined by infinite number of invariants. In our experiments we indeed observe that the system asymptotically approaches towards it’s stationary turbulent state. The investigation of this state has fundamental importance. Note that the similar research has recently been made for the focusing NLS equation, but with incoherent wave field initial conditions [[12]] (see also [[13], [14], [15], [16], [17]] for the dependence of the final state on the Benjamin-Feir index and [[18]] for the defocusing NLS equation). Since we study integrable system which “remembers” it’s initial state through infinite number of integrals of motion, it is not surprising that our asymptotic turbulent state differs from that of [[12]]. Second, we examine the beginning of the nonlinear stage of the MI and the subsequent evolution towards the asymptotic turbulent state in order to understand the characteristic features of rogue waves emergence in the framework of the focusing NLS equation. In this sense our study is in line with the intensive modern research on the statistics of waves in different nonlinear systems, and many of these systems fall under the category of generalized one-dimensional NLS equation [[4], [5], [19], [20], [21], [22], [23], [24], [25], [26]].

Since we perform our simulations in the finite box , after a very long time we will encounter with Fermi - Pasta - Ulam (FPU) recurrence phenomenon [[27], [28]]. This means that at some point of time the evolution towards the stationary integrable turbulence will stop, and the system will move back towards the condensate state. Thus, the system of finite size does not have the asymptotic stationary turbulent state in it’s true sense. However, it has quasi-asymptotic state in which it spends most of the time between the development of the MI and the FPU recurrence. Since the time of the FPU recurrence tends to infinity with the system size , the quasi-asymptotic state approaches to the stationary integrable turbulence as .

We stop our simulations before we observe tendency towards the FPU recurrence. Technically, we perform convergence study comparing our results obtained in the computational boxes and . As soon as we observe deviations between these results, we stop simulations in the smaller box and switch to simulations in the box , comparing the results with that from the box . We repeat this procedure until our computational resources allow. In this sense our results can be considered as the subsequent approximations of the stationary integrable turbulence. Taking this into account we will continue to use term ”asymptotic state” in it’s initial meaning.

One of the important characteristics of the turbulence is wave-action spectrum

 Ik(t)=⟨|Ψk(t)|2⟩. (9)

Here and below stands for arithmetic average across ensemble of initial data and is Fourier transform of . We define forward and backward Fourier transformations as

 Ψk(t) = F[Ψ(x,t)]=1L∫L/2−L/2Ψ(x,t)e−ikxdx, (10) Ψ(x,t) = F−1[Ψk(t)]=∑kΨk(t)eikx, (11)

where is wavenumber and is integer. In our simulations we use where is integer, so that our spectral band contains exact wavenumbers where the maximum growth rate of the MI is achieved. Wave-action spectrum is the spectral density of wave action, since

 ⟨N⟩=∑kIk(t). (12)

Thus, the right-hand side of Eq. (12) is conserved by the motion. According to (10), all wave-action of the condensate is concentrated in the zeroth harmonic ,

 Ik={1,k=0,0,k≠0. (13)

During the development of the MI we observe that wave action disperse across other harmonics. This happens in the form of the oscillatory exchange of wave action between the zeroth harmonic from one hand, and the rest of the spectrum from the other. In the result of this process wave-action spectrum approaches towards the asymptotic spectrum. In the beginning of the MI the spectrum has discontinuity at in the form of a high peak occupying the zeroth harmonic only. This peak appears from the initial data (5). The remarkable result of our experiments is that the peak does not disappear with the arrival to the nonlinear stage of the MI, but instead decays in oscillatory way and remains detectable for a long time after the beginning of the nonlinear stage. After the peak finally disappears, the singularity in the spectrum at transforms to power-law behavior at with exponent close to 2/3. The corresponding modes have very large scales in the physical space and can be called ”quasi-condensate”, that in the asymptotic turbulent state has about 40% of wave action, less than 1% of kinetic energy and about 10% of potential energy. The asymptotic spectrum decays monotonically as ; this decay is slower at , very fast near and , and close to exponential from .

Another important characteristic of the turbulence is the (simultaneous) spatial correlation function,

 g(x,t)=⟨1L∫L/2−L/2Ψ(y,t)Ψ∗(y−x,t)dy⟩. (14)

It is connected with wave-action spectrum by the relation

 g(x,t)=F−1[Ik(t)],

that follows from Eqs. (10)-(11). Spatial correlation function is fixed to unity at , since

 g(0,t)=⟨N⟩, (15)

and for the ensemble of initial data (5) wave action almost coincides with unity (for our experiments ).

In the nonlinear stage of the MI we observe that spatial correlation function also evolves in oscillatory way approaching towards the asymptotic correlation function. While the peak at in wave-action spectrum is present, decays with to some nonzero level that is determined by the magnitude of the peak. When the peak disappears, the correlation function decays to zero as . At lengths the asymptotic spatial correlation is close to Gaussian. Here is it’s full width at half maximum.

We also measure the probability density function (PDF) of wave amplitudes . Let us suppose that the current state of a system consists of multitude of uncorrelated linear waves,

 Ψ(x)=∑k|Ψk|ei(kx+ϕk). (16)

If phases are random and uncorrelated, the number of waves is large enough, and amplitudes fall under the conditions of central limit theorem, then real and imaginary parts of field are Gaussian-distributed, and the PDF of wave amplitudes coincides with Rayleigh distribution [[29]] (see example on FIG. 1),

 PR(|Ψ|)=2|Ψ|σ2e−|Ψ|2/σ2. (17)

For a system that conserves wave action and has Rayleigh PDF, the parameter can be readily calculated as

 ⟨N⟩=⟨|Ψ|2⟩=∫+∞0|Ψ|2PR(|Ψ|)d|Ψ|=σ2.

Here is ensemble and space average of squared amplitude. For the ensemble of initial data (5) this leads to conclusion .

Since , the PDF of squared amplitudes is exponential if the corresponding amplitude PDF is Rayleigh one, and vice versa; for this leads to

 PR(|Ψ|2)=e−|Ψ|2. (18)

It is more convenient to examine exponential dependencies than Rayleigh ones, and thereby we measure the PDF of squared amplitudes and compare the results with exponential dependency (18) that we call Rayleigh one for simplicity. Throughout the publication we use term ”PDF” only in relation to PDF of wave amplitudes. We measure the PDF for entire field , in contrast to PDFs for local maximums or absolute maximums, and use normalization,

 ∫+∞0P(|Ψ|2)d|Ψ|2=1.

The knowledge of the PDF gives us the probability of occurrence of waves exceeding certain threshold ,

 W(Y,t)=∫+∞YP(|Ψ|2,t)d|Ψ|2. (19)

In case of Rayleigh PDF (18) this probability takes the simple form

 WR(Y,t)=e−Y. (20)

In addition to the PDF we measure ensemble average kinetic and potential energies, and also the moments

 M(n)(t)=⟨1L∫+L/2−L/2|Ψ(x,t)|ndx⟩. (21)

The moments are connected to the PDF as

 M(n)(t)=∫+∞0|Ψ|nP(|Ψ|,t)d|Ψ|. (22)

Thus, for a system with Rayleigh PDF (18) the moments can be easily calculated,

 M(n)R=Γ(n2+1), (23)

where is gamma-function. The moment does not change with time since .

In the nonlinear stage of the MI we observe that kinetic and potential energies, as well as the moments , , oscillate with time around their asymptotic values. The amplitudes of these oscillations decay with time as , the phases contain the nonlinear phase shift that decays as , and the period of the oscillations is equal to . Thus, the frequency of the oscillations is equal to the double maximum growth rate of the MI. The asymptotic values of kinetic and potential energies are 0.5 and -1 respectively, while the asymptotic moments coincide with Rayleigh predictions (23).

The PDF in the asymptotic turbulent state coincides with Rayleigh one (18). Thus, the probability of occurrence of waves turns out to be the same as in a wave field described by linear equations (20). The level of nonlinearity of the turbulence can be estimated by the parameter

 Q=|⟨H4⟩||⟨Hd⟩|.

For weak turbulence the Rayleigh PDF would be a natural result. However, for the NLS equation we observe ”moderately strong” turbulence with in the asymptotic state.

Nevertheless, we confirm that in the beginning of the nonlinear stage the MI moderately increases the occurrence of rogue waves. According to the standard definition [[1], [2], [3]], a rogue wave is a wave that exceeds at least two times the significant wave height . The significant wave height is calculated as the average wave height of the largest 1/3 waves. It is easy to calculate that for Rayleigh PDF (18) the significant wave amplitude is , and rogue waves must exceed in amplitude or in squared amplitude.

In our experiments we observe that in the beginning of the nonlinear stage of the MI the PDF evolves significantly with the oscillations of kinetic and potential energies. At the points of time corresponding to local maximums and minimums of potential energy modulus , the PDF acquires ”fat tails” and significantly exceeds Rayleigh PDF (18) in the two regions of squared amplitudes and respectively. It is interesting that the evolution of the PDF goes in such a way that in the beginning of the nonlinear stage the ”standard” rogue waves appear even less frequently than predicted by Rayleigh PDF (20).

The waves from the first region are ”imperfect” rogue waves, since they do not match the criterion for the ”standard” rogue waves . The ”imperfect” rogue waves are the typical outcome of the MI, and can be seen at the first several local maximums of . In space these waves form a modulated lattice of large waves with distance between them close to the characteristic length of the MI. In the beginning of the nonlinear stage the probability of occurrence of such waves with is by about three times larger than Rayleigh one (20). The crests of the ”imperfect” rogue waves are mostly composed of the imaginary part of wave field , . At the first, third, and so on, local maximums of it is positive , and at the second, fourth, and so on, local maximums – negative .

The similar scenario is realized in case of the Akhmediev breather [[30], [31], [32]] that corresponds to the maximum growth rate of the MI. At the time of it’s maximal elevation this solution is purely imaginary, and at it’s maximums the imaginary part is positive . After the decay this solution changes the phase of the condensate by . Thus, the following Akhmediev breather – if it appears – should have negative imaginary part at it’s crests, the third Akhmediev breather – positive, and so on.

However, it is unclear how these solutions may appear from random statistically homogeneous in space noise one after another with the short interval between them. This interval is equal to the period of the oscillations of the moments, as well as kinetic and potential energies, is close to 4 in the beginning of the nonlinear stage of the MI and approaches to with time. Also, spatial correlation function of the ”imperfect” rogue waves significantly decreases after a few characteristic lengths of the MI, and takes (locally in time) minimal values. For the Akhmediev breather it remains periodic. Wave-action spectrum of the ”imperfect” rogue waves also is not very similar to that of the Akhmediev breather.

In the beginning of the nonlinear stage of the MI and at the local minimums of potential energy modulus the waves from the second region appear by about two times more frequently for than predicted by Rayleigh PDF (20). These waves are very rare events, and represent in space a singular high peak with full width at half maximum of about and duration in time of about . These ”large” rogue waves appear on the background of perturbed wave field that is usually less than in amplitude. Statistically at this time wave field is significantly correlated, spatial correlation function takes (locally in time) maximal values, and wave-action spectrum has (locally in time) maximal zeroth harmonic with the rest of the spectrum minimally excited.

The crests of the ”large” rogue waves are composed mostly of real part of wave field , . At the first, third, and so on, local minimums of it is negative , and at the second, fourth, and so on, local minimums – positive . It is interesting, that the Peregrine solution [[10]] has similar property: at the time of it’s maximal elevation this solution is purely real and negative at it’s maximum amplitude. However, the Peregrine solution has slightly smaller maximal squared amplitude .

We also observe extremely large waves with amplitudes of up to . However, the accuracy of our simulations is insufficient to study the time evolution of the PDF and the probability of occurrence of such waves.

The paper is organized as follows. In the next Section we describe the numerical methods that we used in the framework of the current study. Section 3 is devoted to the investigation of the asymptotic state of the integrable turbulence, while in Section 4 we describe how the temporal evolution towards this state is arranged. Section 5 contains conclusions and acknowledgements. In the Appendix A we provide detailed graphs for wave-action spectrum and spatial correlation function in the beginning of the nonlinear stage of the MI. Detailed graphs for the evolution of the probability of large waves occurrence are given in the Appendix B.

## Ii Numerical methods.

We integrate Eq. (1) numerically on the time interval in the box , , with periodic boundary. We use , where is integer, in order to have in our spectral band exact wavenumbers where the maximum growth rate of the MI is achieved. We continue integration for so long time in order to approach to the asymptotic turbulent state as close as our computational resources allow. We do not integrate beyond since starting from we observe tendency towards the FPU recurrence: kinetic and potential energies, as well as the moments , significantly deviate from their asymptotics, as well as from the results obtained on larger computational box (for that large we cannot allow ourselves comparison with the box ). We also performed simulations on smaller boxes and found the same phenomenon starting from for , for , and so on.

We would like to stress, that we have a very good quantitative agreement of our results obtained on different computational boxes before the tendency towards the FPU phenomenon appears. Thus, results for coincide with that obtained on larger boxes up to , for – up to and for – up to . Note that in order to compare wave-action spectrum obtained on different boxes , we need to use it in the form , since for finite wavenumber actually models area in the spectrum. Here is the distance between the subsequent wavenumbers. Below we will continue to use the spectrum in the definition (9), where the spectrum of the condensate takes the convenient form (13).

We use Runge-Kutta 4th-order method, and calculate spatial derivatives and wave-action spectrum with the help of Fast Fourier transformations (FFT) routines. We perform our simulations on uniform grid with the spatial grid size , where is the number of nodes on the grid. Thus, all summation over wavenumbers in (11), (12), and so on, where is integer, is performed in the spectral band .

We change the spatial grid size adaptively after the analysis of Fourier components of the solution : we reduce when at large wave numbers exceed and increase when this criterion allows. The distance between the subsequent wavenumbers is fixed by the length of the computational box , and the range of wavenumbers is determined by . Thereby, the modification of only adds or removes harmonics with large wavenumbers. This allows us to perform interpolation from one uniform grid to another by simply transferring the shared part of the spectrum to the new grid. We checked that the error of such interpolation is comparable with the round-off errors. In order to prevent appearance of numerical instabilities, time step also changes with as , .

We start our simulations on the grid with nodes. In order to calculate the ensemble average characteristics, we interpolate the solution from the current grid (determined by the spectrum at the current time ) to fixed grid with nodes. We checked that such grids are sufficient for our computational box, comparing our results with that obtained on larger grids. We start from the initial data (5) and use statistically homogeneous in space initial noise that can be written symbolically as

 ϵ(x)=A0(√8πθL)1/2F−1[e−k2/θ2+iξk]. (24)

Here is noise amplitude, is noise width in k-space and are arbitrary phases for each and each noise realization within the ensemble of initial data. The average squared amplitude of noise in x-space can be calculated as

 ¯¯¯¯¯¯¯|ϵ|2=1L∫L/2−L/2|ϵ(x)|2dx=√8πθLA20L∑k1,k2e−(k21+k22)/θ2+i(ξk1−ξk2)∫+L/2−L/2ei(k1−k2)xdx= =√8πθLA20∑ke−2k2/θ2≈√8πθLA20(2πL)−1∫+∞−∞e−2k2/θ2dk=A20. (25)

We performed several experiments for several ensembles of initial data that differed from each other by noise parameters and . We didn’t find significant dependence of our results on noise amplitude , except that the system arrives to the nonlinear stage of the MI faster for larger . We also changed noise width in k-space in the broad range from to , and obtained the same results for these experiments. In the next two Sections we present the results obtained for the most broad in k-space noise distribution , , that we studied. Inside the range of the instability such noise can be treated as a white noise.

We also tested the following initial noise distribution,

 ϵ2(x)=A0(√8πθL)1/2F−1[10−vk×e−k2/θ2+iξk], (26)

where is uniformly distributed over random value for each , and . The multiplier introduces the detuning between the amplitudes of noise in -space by up to 10 orders of magnitude. However, we came to very similar results, though the oscillatory evolution of the system in the nonlinear stage of the MI, that we will demonstrate below, became slightly less regular. In our opinion this means that our results should be visible for a very wide variety of statistical distributions of noise.

The NLS equation (1) has an infinite number of integrals of motion [[7]]. The first three of these integrals are wave action (6), momentum (7) and total energy (8), the fourth one is

 c4[Ψ(x)]=1L∫+∞−∞[ΨΨ∗xxx+12Ψddx(|Ψ|2Ψ∗)+|Ψ|2ΨΨ∗x]dx,

and so on. Our scheme provides very good conservation of the first 12 invariants with accuracy better than . We measure absolute errors for integrals with even numbers , and relative errors for integrals with odd numbers , since for our initial data integrals with even numbers are very close to zero. The first three invariants are conserved by our scheme with accuracy better than .

In our experiments we use ensembles of 1000 initial distributions each. We checked our statistical results against the size of the ensembles, the parameters of our numerical scheme and the implementation of other numerical methods (Runge-Kutta 5th-order, Split-Step 2nd- and 4th-order methods [[33], [34]]), and found no difference.

We also compared our results with that for the MI of the condensate in the framework of the Ablowitz–Ladik (AL) equation [[35], [36]],

 idΨndt+Ψn+1−2Ψn+Ψn−1h2−Ψn+|Ψn|2Ψn+1+Ψn−12=0. (27)

The AL system is defined on the grid with periodic boundary, where is (integer) node number and is the total number of nodes, and is analogous to the NLS equation. As shown in [[36]], the problem of the MI of the condensate for the AL system has one free parameter that has the meaning of the constant of coupling between the nodes. The NLS equation can be obtained after the substitution and in the limit . Therefore, for the AL system (27) may be considered as the scheme of numerical integration of the NLS equation with fixed discretization along spatial dimension. Using this scheme together with Runge-Kutta 4th order method, we arrived to exactly the same results as in case of the described above numerical scheme for the integration of the NLS equation.

Comparison of our results with that for the AL system is also important from the point of view of integrability. Indeed, the AL system is also completely integrable with the help of the inverse scattering transformation. The schemes of numerical integration of the NLS equation break the integrability due to discretization along both spatial and temporal dimensions. However, during numerical simulations of the AL system the integrability is broken due to discretization along temporal dimension only. Therefore, since our results for the NLS equation and the AL system with coincide, we may conclude that violation of integrability due to spatial discretization does not affect our results.

## Iii Integrable turbulence: the asymptotic state.

In the linear stage of the MI, when perturbations to the condensate are small, wave-action spectrum almost coincides with the spectrum of the condensate (13), spatial correlation function is almost indistinguishable from unity , and the PDF represents a very high and thin peak around . All of the moments (21) are unity , kinetic energy is zeroth , and potential energy is equal to .

With the development of the MI this situation changes; for our initial data the system arrives to the nonlinear stage approximately at (the characteristic time of the instability is ). In the nonlinear stage we observe that the moments with exponents , and also kinetic and potential energies oscillate with time around their asymptotic values, and the amplitudes of these oscillations decay as time increases. The evolution of wave-action spectrum, spatial correlation function and the PDF is more complex, as in addition to oscillations they simultaneously change their forms. These functions also evolve towards some asymptotic forms that they take at late times. Thus, it is possible to say that during the nonlinear stage of the MI the system evolves from the condensate state towards some asymptotic turbulent state. The asymptotic state is characterized by independent on time wave-action spectrum, spatial correlation function, the PDF, the moments, and also kinetic and potential energies, that we below call asymptotic ones. We study these asymptotic characteristics in the current Section, while the temporal evolution towards the asymptotic state - in the next Section.

We find asymptotic characteristics by averaging the corresponding functions over time , where the deviations of these functions with time are sufficiently small. Nevertheless, their temporal evolution is still visible even at . However, this evolution is oscillatory-like, and averaging over sufficiently long time interval should provide us very good approximation to the asymptotic characteristics. We would like to stress that such time-averaging procedure gives very similar results already starting from time interval . Therefore, we believe that integration beyond on larger computational boxes should not provide us very different behavior.

FIG. 2a,b shows asymptotic wave-action spectrum . The spectrum decays monotonically as . This decay is slower at , very fast near and , and close to exponential from . The remarkable property of the spectrum is that it has singularity at zeroth harmonic . As will be shown in the next Section, in the beginning of the nonlinear stage of the MI this singularity represents a high peak occupying the zeroth harmonic only. The peak decays in oscillatory way, remaining detectable up to . At later times we observe the formation of power-law dependence for with exponent close to 2/3 (see FIG. 2b). The asymptotic wave-action spectrum for these wavenumbers is very well approximated by the function

 Ik≈b|k|−α, (28)

where and . At the asymptotic spectrum has finite value .

The power-law behavior of the asymptotic spectrum means that the corresponding modes near are relatively large. Wave action concentrated in modes can be calculated as

 ⟨N(k0)⟩=∑|k|≤k0Ik. (29)

It coincides with the ensemble average wave action (6) in the limit . It turns out that approximately 41% of wave action is concentrated in modes from the power-law region , , since , and about 7% of wave action is concentrated in just 3 modes and (see also FIG. 6a). Note that integration of the asymptotic (28) over modes , , gives very similar result:

 N(k0)=∑|k|≤k0Ik≈2Δk∫k00Ikdk≈2b(1−α)Δkk1−α0≈0.47,

where is the distance between the subsequent wavenumbers. The modes have scales in the physical space comparable with the length of the integration box , and thus can be called ”quasi-condensate”.

The macroscopic wave action concentration into quasi-condensate from one hand, and the FPU phenomenon from the other are the reasons why the numerical simulations that we perform must be implemented on computational boxes with very large lengths. For smaller computational boxes the distance between wavenumbers is too large, and the region of wavenumbers containing about 40% of wave action cannot be carefully resolved.

We also detect another region of power-law dependence of the asymptotic wave-action spectrum at wavenumbers , where the spectrum decays close to (see FIG. 3b). The maximum growth rate (4) of the MI is realized at ; the corresponding modulations have characteristic scale in the physical space. Thus, the second power-law region corresponds to modes with scales , or 1–2.5 characteristic scales of the MI. These modes acquire about 25% of wave-action. We think that after a very long evolution computed on a very large computational box the two power-law regions in the spectrum might merge in one , , with some shared exponent . However, the computational resources that we have at our disposal are insufficient to check this hypothesis.

Modes with wavenumbers , corresponding to scales smaller than in the physical space, decay in the asymptotic spectrum close to exponential law , . These modes have about 5% of wave-action.

The asymptotic spatial correlation function is shown on FIG. 4a,b. It’s characteristic scale, defined as full width at half maximum, is . At lengths it is close to Gaussian (see Eq. (15): ),

 g(x)≈exp[−4ln2(xxcorr)2]. (30)

At the asymptotic correlation function decays very slowly to about as . In the region this decay is very well approximated as

 g(x)≈b1|x|+b2, (31)

with the coefficients and (see FIG. 4b). In the region the decay is even slower, but we believe that this is the effect of the finiteness of the computational box. Indeed, by construction spatial correlation function is periodic , and should be even , therefore it’s derivatives at the borders of the computational box should be zeroth . Thus, near the borders the behavior of the correlation function should deviate from (31).

The squared amplitude asymptotic PDF coincides with Rayleigh PDF (18), as shown on FIG. 5a. Thus, we come to surprising conclusion, that despite the nonlinearity of the NLS equation it’s asymptotic PDF is the same that would be for a wave field described by linear equations. We additionally checked this conclusion by calculating the asymptotic values of the moments (21) , and found that these values coincide with their Rayleigh predictions (23), at least for exponents (see FIG. 5b). The asymptotic moment allows us to calculate the ensemble average potential energy in the asymptotic state,

 ⟨H4⟩=−12M(4)A≈−1. (32)

Combined with the conservation of total energy , this yields the ensemble average kinetic energy ,

 ⟨Hd⟩≈0.5. (33)

We also calculated independently with the same result. Thus, in the asymptotic state we have ”moderately strong” turbulence with . This makes the conclusion of Rayleigh statistics for amplitudes even more surprising.

In fact, the relation (32) is itself truly remarkable. According to (11),

 −⟨H4⟩=12∑k1,k2,k3,k4⟨Ψk1Ψk2Ψ∗k3Ψ∗k4⟩δ(k1+k2−k3−k4), (34)

where is Kronecker delta

 δ(k)={1,k=0,0,k≠0.

The four-wave momentum in (34) can be represented as

 ⟨Ψk1Ψk2Ψ∗k3Ψ∗k4⟩=Ik1Ik2(δ(k1−k3)δ(k2−k4)+δ(k1−k4)δ(k2−k3))+Jk1,k2,k3,k4, (35)

where is the cumulant. Since , we obtain

 −⟨H4⟩≈1+12∑k1,k2,k3,k4Jk1,k2,k3,k4δ(k1+k2−k3−k4). (36)

Together with the relation (32) this yields

 ∣∣∣∑k1,k2,k3,k4Jk1,k2,k3,k4δ(k1+k2−k3−k4)∣∣∣≪1. (37)

The latter result might mean that the cumulant in the asymptotic turbulent state is zeroth. Thus, the relation can be considered as an indication that the stationary integrable turbulence is purely Gaussian. The calculation of the cumulant is a cumbersome problem, and we will study it in a separate publication.

It is interesting to examine the spectral distribution of the ensemble average kinetic and potential energies. The spectral density of kinetic energy can be obtained from Eq. (8) and Eq. (10)-(11). For convenience we divide it by the distance between the subsequent wavenumbers , and use it in the following form:

 T(k)=k2(Ik+I−k)Δk. (38)

Kinetic energy concentrated within modes can be calculated as

 ⟨Hd(k0)⟩=∑|k|≤k0k2Ik=Δk∑0≤k≤k0T(k). (39)

Calculation of the spectral distribution of potential energy is more complex. For this purpose we introduce the new function ,

 ~Ψ(x,t)=F−1[~Ψk(t)], (40)

that contains only modes with from the original solution of the NLS equation ,

 ~Ψk(t)={Ψk(t),|k|≤k0,0,|k|>k0, (41)

and is Fourier components of original solution. Then we find potential energy concentrated in these modes and average it across the realizations of initial data,

 ⟨H4(k0)⟩=−⟨12L∫L/2−L/2|~Ψ(x,t)|4dx⟩, (42)

Hence, the spectral density of potential energy can be calculated as

 U(k)=⟨H4(k+Δk)⟩−⟨H4(k)⟩Δk. (43)

In the limit the quantities (39) and (42) coincide with the ensemble average kinetic and potential energies (8) respectively.

The distribution across wavenumbers of kinetic and potential energies, as well as their spectral densities and , are shown on FIG. 6a,b. The spectral density of kinetic energy monotonically increases from zero at zeroth harmonic to at , where it achieves maximum, sharply decreases at , and then decays to zero starting from . The spectral density of potential energy modulus has sharp maximum at , monotonically decreases to at , where it achieves local minimum, and then behaves similar to , acquiring another maximum at .

It is interesting that is always larger than , so that in the asymptotic turbulent state all modes are essentially nonlinear. The ”most nonlinear modes” are those of the quasi-condensate , that contain less than 1% of kinetic energy, about 10% of potential energy, and about 40% of wave action. Modes with contain about 60% of kinetic and potential energies and about 55% of wave action, while the exponentially decaying modes have about 5% of wave action, about 40% of kinetic energy and 30% of potential energy.

## Iv The nonlinear stage of the modulation instability: evolution towards the asymptotic state.

It will be convenient for us to study the evolution towards the asymptotic turbulent state on the example of ensemble average kinetic and potential energies, and also the moments (see FIG. 7a,b). Up to the perturbations to the condensate are small, so that the moments and the energies do not change substantially from their initial values , and respectively. At the MI arrives to it’s nonlinear stage; the moments start to oscillate around their asymptotic Rayleigh values (23), kinetic energy - around 0.5, and potential energy – around -1. The moment oscillates in-phase with potential energy , and antiphase with the moments , , and kinetic energy , so that the positions in time of local maximums and minimums of and coincide with the positions of local minimums and maximums of , , and respectively.

We study the time dependence of the oscillations on the example of moment . FIG. 8a shows that the amplitude of the oscillations of , that we measure as the modulus of deviations of local maximums and minimums of from it’s asymptotic value , is very well approximated by the function with the prefactor . The period of the oscillations changes from at to at . We think that this is the effect analogous to the nonlinear phase shift. Indeed, one can search for the approximation of in the form

 M(1)(t)≈M(1)A+pt3/2sin(Φ(t)),Φ(t)=st+ϕnl(t)+Φ0, (44)

where is phase, is constant frequency, is constant phase, and the nonlinear phase shift should be proportional to the amplitude of the oscillations multiplied by time , or with constant . Then, the phases at the local maximums of should be equal to

 Φ(tmax)=stmax+q√tmax+Φ0=π2+2πm,

and at the local minimums – to

 Φ(tmin)=stmin+q√tmin+Φ0=3π2+2πm,

where is integer number. We find all the subsequent extremums and of from one hand, and their phases from the other hand by setting for the first maximum, for the second maximum, and so on. Then, with the help of the least squares method we determine the coefficients , and . After that we check that the nonlinear phase shift

 Φ(t)−st−Φ0, (45)

calculated at the extremums of , indeed is very well approximated by the function , as shown on FIG. 8b.

We observe that anzats (44) fits very well to the experimental time dependence of the moments, and also kinetic and potential energies as well. The example of such fit for is shown on FIG. 9. The phases for the moment and potential energy coincide, and differ by from the phases for the moments , , and kinetic energy . We checked that anzats (44) without the nonlinear phase shift, or with the exponent of the nonlinear phase shift significantly different from –0.5, fits significantly worse to the experimental data. It is interesting to note that the period of the oscillations is almost equal to , and their frequency is almost equal to the double maximum growth rate of the MI (4), . We think that the frequency should coincide with 2, and the period should accordingly coincide with , and we measure the frequency almost the same for all of our experiments irrespective of the statistics of initial noise. However, the nature of such correspondence is unclear for us yet.

Wave-action spectrum, spatial correlation function and the PDF of squared amplitudes also evolve in oscillatory way with time, approaching to their asymptotic forms at late times. The ”turning points” for the evolution of these functions – the points in time where the motion of , and at fixed , and respectively changes to roughly the opposite – approximately coincide with the local maximums and minimums of the moments, and also kinetic and potential energies. For definiteness, below we will refer to such points in time on the example of extremums of potential energy modulus . This choice has also straightforward physical sense: at the local maximums of the effect of nonlinearity is the largest, and at the local minimums – the smallest. Note, that the spectrum, the correlation function and the PDF do not evolve exactly as