Early stage of integrable turbulence in 1D NLS equation: the semi-classical approach to statistics

TL;DR

This paper investigates the early statistical evolution of integrable turbulence in the 1D nonlinear Schrödinger equation using a semi-classical approach, deriving explicit formulas and validating them with numerical simulations.

Contribution

It introduces a semi-classical method to analytically describe early-stage statistical properties of 1D-NLSE turbulence, focusing on the evolution of the fourth moment.

Findings

Excellent agreement between analytical formulas and numerical simulations.

Identification of heavy tails in focusing regime and low tails in defocusing regime.

Explicit formula for early evolution of wave amplitude statistics.

Abstract

We examine statistical properties of integrable turbulence in the defocusing and focusing regimes of one-dimensional small-dispersion nonlinear Schrodinger equation (1D-NLSE). Specifically, we study the 1D-NLSE evolution of partially coherent waves having Gaussian statistics at time t=0. Using short time asymptotic expansions and taking advantage of the scale separation in the semi-classical regime we obtain a simple explicit formula describing an early stage of the evolution of the fourth moment of the random wave field amplitude, a quantitative measure of the "tailedness" of the probability density function. Our results show excellent agreement with numerical simulations of the full 1D-NLSE random field dynamics and provide insight into the emergence of the well-known phenomenon of heavy (resp. low) tails of the statistical distribution emerging in the focusing (resp. defocusing) regime of 1D-NLSE.

Equations

$$i\epsilon \frac{\partial \psi }{\partial t}+\frac{{\epsilon }^2}{2}\frac{{\partial }^2\psi }{\partial x^2}+\sigma |\psi {|}^2\psi =0,$$

$$\psi (x,t)=\sum _k{\psi }_k(t)e^{\frac{2i\pi }{L}kx}withk\in \mathbb{Z},$$

$${\psi }_k(t)=\frac{1}{L}{\int }_0^L\psi (x,t)e^{-2i\pi kx}dx.$$

$$N=\frac{1}{L}{\int }_0^L|\psi {|}^{2}dx=\sum _k|{\psi }_k{|}^2,$$

$$P=\frac{1}{L}\int_{0}^{L}\psi_{x}\psi^{*}dx=\sum_{k}\Big{(}\frac{2\pi ik}{L}\Big{)}|\psi_{k}|^{2}.$$

$$H={\epsilon }^2{H}_{\mathrm{L}}+H_{\mathrm{NL}}$$

$$\varepsilon^{2}H_{\rm L}(t)=\frac{\varepsilon^{2}}{2L}\int_{0}^{L}|\psi_{x}|^{2}dx=\frac{\varepsilon^{2}}{2}\sum_{k}\Big{(}\frac{2\pi k}{L}\Big{)}^{2}|\psi_{k}|^{2}$$

$$H_{\mathrm{NL}}(t)=\frac{\sigma }{2L}{\int }_0^L|\psi {|}^{4}dx.$$

$${\kappa }_4(t)=\frac{⟨\frac{1}{L}{\int }_0^L|\psi (x,t){|}^{4}dx⟩}{⟨\frac{1}{L}{\int }_0^L|\psi (x,t){|}^{2}dx{⟩}^2},$$

$$\frac{dH_{\rm NL}}{dt}=\frac{\sigma}{L}\int_{0}^{L}|\psi|^{2}\,\Big{[}\psi\frac{\partial\psi^{*}}{\partial t}+\psi^{*}\frac{\partial\psi}{\partial t}\Big{]}\;dx.\\$$

$$\frac{d{H}_{\mathrm{NL}}}{dt}=\frac{\sigma \epsilon }{L}{\int }_0^{L}Im[({\psi }_x{\psi }^{\ast }{)}^2]dx.$$

$$\psi =\sqrt{\rho }{e}^{i\frac{\varphi }{\epsilon }},u=\frac{\partial \varphi }{\partial x},$$

$$\frac{d{H}_{\mathrm{NL}}}{dt}=\frac{\sigma }{L}{\int }_0^L\rho u{\rho }_{x}dx.$$

$${\kappa }_4(t)=\frac{2⟨H_{\mathrm{NL}}⟩}{\sigma ⟨N{⟩}^2}$$

$$\frac{d{\kappa }_4}{dt}=\frac{2}{\sigma ⟨N{⟩}^2}\frac{d⟨H_{\mathrm{NL}}⟩}{dt}=\frac{2}{⟨N{⟩}^{2}L}{\int }_0^L⟨\rho u{\rho }_x⟩dx.$$

$$\{\begin{array}{l}{\rho }_t+(\rho u{)}_x=0,\\ u_t+u{u}_x-\sigma {\rho }_x=0.\end{array}$$

$$\rho (x,t)={\rho }_0(x)+{\rho }_1(x)t+{\rho }_2(x)t^2+{\rho }_3(x)t^3+\mathcal{O}(t^4)$$

$$u_0(x,t)=u_1(x)t+u_2(x)t^2+u_3(x)t^3+\mathcal{O}(t^4).$$

$$\rho (x,t)={\rho }_0-\frac{1}{4}\sigma [{\rho }_0^2{]}_{xx}{t}^2+\mathcal{O}(t^4),$$

$$u(x,t)=\sigma {\rho }_{0x}t-(\frac{1}{12}[{\rho }_0^2{]}_{xxx}+\frac{1}{3}{\rho }_{0x}{\rho }_{0xx})t^3+\mathcal{O}(t^4).$$

$${\kappa }_4(t)-{\kappa }_4(0)=$$

$$-\frac{t^4}{2⟨N{⟩}^{2}L}{\int }_0^L⟨\frac{2}{3}{\rho }_0^2{\rho }_{0x}{\rho }_{0xxx}$$

$$+\frac{17}{6}{\rho }_0{\rho }_{0x}^2{\rho }_{0xx}+\frac{1}{2}{\rho }_{0x}^4⟩dx+\mathcal{O}(t^6).$$

$$\big{<}\psi_{k}\psi_{k^{\prime}}^{*}\big{>}=n_{k}\,\delta_{k}^{k^{\prime}},$$

$$⟨{\psi }_{k_1}{\psi }_{k_2}{\psi }_{k_3}{\psi }_{k_4}^{\ast }{\psi }_{k_5}^{\ast }{\psi }_{k_6}^{\ast }⟩=n_{k_1}{n}_{k_2}{n}_{k_3}$$

$$[{\delta }_{k_4}^{k_1}{\delta }_{k_5}^{k_2}{\delta }_{k_6}^{k_3}+{\delta }_{k_4}^{k_1}{\delta }_{k_5}^{k_3}{\delta }_{k_6}^{k_2}+{\delta }_{k_4}^{k_2}{\delta }_{k_5}^{k_1}{\delta }_{k_6}^{k_3}+$$

$${\delta }_{k_4}^{k_2}{\delta }_{k_5}^{k_3}{\delta }_{k_6}^{k_1}+{\delta }_{k_4}^{k_3}{\delta }_{k_5}^{k_1}{\delta }_{k_6}^{k_2}+{\delta }_{k_4}^{k_3}{\delta }_{k_5}^{k_2}{\delta }_{k_6}^{k_1}]$$

$$\frac{\sigma }{⟨N{⟩}^{2}L}$$

$$\times ⟨{\psi }_{0k_1}{\psi }_{0k_2}{\psi }_{0k_3}{\psi }_{0k_4}^{\ast }{\psi }_{0k_5}^{\ast }{\psi }_{0k_6}^{\ast }⟩(k_2-k_5)(k_3-k_6),$$

$${\kappa }_4(t)-{\kappa }_4(0)$$

$$\times [-2(k_2-k_3{)}^2]t^2+\mathcal{O}(t^4),$$

$${\kappa }_4(t)-{\kappa }_4(0)=8\sigma ⟨H_{\mathrm{L}}(0)⟩t^2+\mathcal{O}(t^4).$$

$${\kappa }_4(t)-$$

$$\displaystyle+4\langle N\rangle\langle\sum_{k}\Big{(}\frac{2\pi k}{L}\Big{)}^{4}|\psi_{0k}|^{2}\rangle\Bigg{]}t^{4}+\mathcal{O}(t^{6}).$$

$$|{\psi }_{0k}{|}^2=n_0{e}^{-\frac{k^2}{(\Delta k{)}^2}}.$$

Full text

I Introduction

Predicting the statistical properties of nonlinear random waves is at the core of research on turbulence Nazarenko (2011). Generally speaking, nonlinear wave propagation in media with self-focusing nonlinearity tends to produce heavy-tailed deviations from the initial Gaussian statistics, observed in the probability density function (PDF) of the random wave field amplitude. In recent years, the question of the emergence of heavy-tailed statistical distributions has been extensively studied in relation to the occurrence of extreme events such as Rogue Waves (RW), mainly in the physical contexts of fluid dynamics Kharif et al. (2009); Onorato et al. (2004, 2005a, 2013) and optics Solli et al. (2007); Mussot et al. (2009); Walczak et al. (2015); Pierangeli et al. (2015); Suret et al. (2016); Närhi et al. (2016); Safari et al. (2017); Tikan et al. (2018). Even though statistical properties of nonlinear defocusing media have been less extensively examined, several experiments have shown that defocusing nonlinearities tend to produce low-tailed deviations from the initial Gaussian statistical distribution Randoux et al. (2014); Bromberg et al. (2010).

Among nonlinear random waves systems, those that are described by the one-dimensional nonlinear Schrödinger Equation (1D-NLSE) play a unique and fundamental role. First of all, 1D-NLSE describes various physical systems at leading order in the amplitude and secondly, 1D-NLSE is integrable and exhibit many remarkable exact solutions, including solitons (both for the focusing and the defocusing cases) Yang (2010). The salient features of focusing regime of 1D-NLSE are modulational instability of plane waves and the existence of breather solutions localized both in space and in time and thus the focusing 1D-NLSE plays a central role in the study of RWs  Peregrine (1983); Henderson et al. (1999); Akhmediev et al. (2009a, b); Kibler et al. (2010); Chabchoub et al. (2011); Kibler et al. (2011). The study of the statistics of nonlinear wave systems described by 1D-NLSE enters within the framework of integrable turbulence proposed by Zakharov Zakharov (2009); Zakharov and Gelash (2013); Agafontsev and Zakharov (2015). Generally, integrable turbulence arises in nonlinear random wave systems described by integrable equations such as the 1D-NLSE, the Korteweg de Vries (KdV) equation and other physically important equations. Given the absence of resonances in integrable systems, the mechanisms underlying integrable turbulence Walczak et al. (2015); Agafontsev and Zakharov (2015); Suret et al. (2016); Randoux et al. (2016); Soto-Crespo et al. (2016); Akhmediev et al. (2016); Randoux and Suret (2016); Suret et al. (2017); Tikan et al. (2018) are of profoundly different nature than those found in the standard wave turbulence Zakharov et al. (1992); Nazarenko (2011); Picozzi et al. (2014). Since many nonlinear wave systems can be described by partial differential equations having an integrable core part, integrable turbulence has become an active field of research both from the theoretical and experimental perspective  Zakharov (2009); Zakharov and Gelash (2013); Agafontsev and Zakharov (2015); Walczak et al. (2015); Suret et al. (2016); Randoux et al. (2016); Soto-Crespo et al. (2016); Akhmediev et al. (2016); Randoux and Suret (2016); Suret et al. (2017); Tikan et al. (2018).

Random wave fields composed of a linear superposition of a large number of independent Fourier modes having random phases and/or amplitudes play a fundamental role in wave turbulence Zakharov et al. (1992); Nazarenko (2011) and in integrable turbulence  Randoux et al. (2016); Randoux and Suret (2016). As an immediate consequence of the central limit theorem, these nonlinear random waves are characterized by Gaussian statistics. Numerical simulations and optical fiber experiments show that the long-term evolution of such Gaussian random wave fields in the focusing regime of 1D-NLSE is accompanied by the emergence of a heavy-tailed PDF  Walczak et al. (2015); Suret et al. (2016); Tikan et al. (2018). Contrastingly, in the defocusing regime of 1D-NLSE, the long-term statistics of random waves with the same initial distributions, is characterized by a low-tailed PDF  Randoux et al. (2014, 2016); Randoux and Suret (2016). The theoretical description of statistical changes that occur in the course of nonlinear propagation of such wave fields in 1D-NLSE systems is a very active current field of research  Randoux et al. (2016); Randoux and Suret (2016).

The first statistical theory of wave systems described by the integrable 1D-NLSE has been established within the framework of wave turbulence theory Janssen (2003). The wave turbulence theory fundamentally relies on a scale separation where the linear dispersion time $T_{L}$ is supposed to be much shorter than the nonlinear time $T_{NL}$. Wave turbulence theory is principally applicable to sets of random dispersive waves that are engaged in a *weakly * nonlinear resonant interaction Nazarenko (2011); Zakharov et al. (1992) (we note, howevere, that resonant energy exchanges are not possible in wave systems ruled by the integrable 1D-NLSE Janssen (2003); Suret et al. (2011); Zakharov (2009)). A modified treatment of the wave turbulence theory based on the derivation of quasi-kinetic equations for the lowest order moments of the wave field has been proposed in ref. Janssen (2003); Soh et al. (2010); Suret et al. (2011); Picozzi et al. (2014). By using this theoretical approach, the evolution of the kurtosis (the fourth standardized moment of the PDF) of partially coherent waves has been derived both in the focusing and defocusing regime of the 1D-NLSE Janssen (2003). The kurtosis quantifies the “heaviness” of the tail in the PDF of the field amplitude distribution and this theory confirmed that the kurtosis increases in the focusing case and decreases in the defocusing case while the initial field is characterized by Gaussian statistics. However, we stress that this kind of analytical treatment is inherently limited to the weakly nonlinear propagation regime.

The question of integrable turbulence (i.e. statistical theory of random nonlinear wave systems governed by integrable partial differential equations) has been recently approached from a completely different perspective. For a certain class of wave systems, one can take advantage of the mathematical framework of dispersive hydrodynamics – the semi-classical theory of nonlinear dispersive waves Biondini et al. (2016) – in order to analyse the wave evolution asymptotically, by introducing two distinct spatio-temporal scales: the long scale specified by initial conditions and the short scale by the internal coherence length (i.e. the typical size of the coherent soliton-like structures).

The semi-classical, dispersive hydrodynamic approach describes the propagation regimes of a completely opposite nature compared to the regimes considered in the framework of wave turbulence theory. This approach can be applied to the 1D-NLSE propagation if the initial scale of the fluctuations of the power of the complex field $|\psi|^{2}$ are much larger than the one corresponding to the balance between nonlinearity and dispersion. In most of the standard cases, this separation of scales correspond to situations where the nonlinear part of the energy is much greater that the linear (kinetic) part of the energy at the initial time. As shown in ref. Randoux et al. (2017), this scale separation permits one to split the development of integrable turbulence into two distinct stages characterized by qualitatively different dynamical and statistical features. At the initial (we shall call it “pre-breaking”) stage of the evolution nonlinear effects dominate linear dispersion and the wave fronts of the random initial field experience gradual steepening leading to the formation of gradient catastrophes that are subsequently regularized through the generation of dispersive shock waves in the defocusing regime El and Hoefer (2016) and of Peregrine-like breather sequences in the focusing regime Tikan et al. (2017). As shown in ref. Randoux et al. (2017), the dynamical and statistical features that occur at the pre-breaking stage of the defocusing 1D-NLSE can be interpreted in terms of the evolution of random Riemann waves.

In this paper, we extend the analysis of the previous works based on the semi-classical treatment of the 1D-NLSE with random initial data by calculating the short-time evolution of the the normalized fourth moment $\kappa_{4}$ of the amplitude of the field. Similar to the standard kurtosis, the quantity $\kappa_{4}$ describes the degree of the deviation from the initial statistical distribution which is often assumed to be Gaussian Onorato et al. (2016). Using the semi-classical Madelung transform and performing the zero dispersion limit, we derive a general analytical expression for the short-time evolution of the fourth moment of the random 1D-NLSE wave field in terms of hydrodynamic variables, and show that this expression can be further simplified for the wave field having Gaussian statistics at initial time. Our analytical asymptotic results are shown to be in excellent agreement with numerical simulations of the evolution of partially coherent initial data in 1D-NLSE.

This paper is organised as follows :

• Sec. II : using the semi-classical approximation, we identify the initial stage of the 1D-NLSE development of partially coherent waves with the nonlinearity dominated, dispersionless regime, and derive the general expression for the short-time evolution of the fourth order moment $\kappa_{4}$ as a power series expansion in time $t$;

• Sec III : we apply the derived formula for $\kappa_{4}$ to the fundamental case of random waves characterized by Gaussian statistics at time $t=0$;

• Sec IV : we provide a comparison between our semi-classical analytical results and numerical simulations of 1D-NLSE.

II The dispersionless limit of the 1D-NLSE and the time evolution of the fourth-order moment of a random wave field

We consider the 1D-NLSE in the normalized form

[TABLE]

where $\psi$ is a complex field, $\varepsilon$ is the dispersion parameter, $\sigma=-1$ in the defocusing regime and $\sigma=+1$ in the focusing regime.

The 1D-NLSE (1) is considered in a periodic box of size $L$, $\psi(x+L,t)=\psi(x,t)\,\forall t$. The field $\psi$ then can be represented as a Fourier series:

[TABLE]

where the Fourier coefficients are given by

[TABLE]

The “density of particles” $N$ and the momentum $P$ represent integrals of motion and are expressed in terms of Fourier coefficients:

[TABLE]

[TABLE]

The Hamiltonian, that we represent in the form

[TABLE]

is also integral of motion, which is naturally split into two parts: the linear (kinetic energy) part

[TABLE]

and the nonlinear part

[TABLE]

We now assume that the Fourier modes at initial time $\psi_{0k}=\psi_{k}(t=0)=|{\psi_{0k}}|e^{i\phi_{0k}}$ are complex random variables. The complex field (2) is then a random periodic solution of the 1D-NLSE. No particular hypothesis about statistical properties of $\psi_{0k}$ needs to be introduced at this step but we will show in Sec. III that the main result of our analysis can be simplified if the initial statistics of the random wave field is assumed to be Gaussian. We consider random initial conditions for which $N$, $H_{\rm L}(t=0)$ and $H_{\rm NL}(t=0)$ are all $\mathcal{O}(1)$. This is typically achieved by taking the initial power spectrum $n_{0k}=|\psi_{0k}|^{2}$ with the characteristic width $\Delta k\simeq 1$, which implies that the typical spatial size of the initial random fluctuations is also of the order of unity and much larger than the internal coherence length (that is $\epsilon$). Such random waves are often called partially coherent, particularly in the statistical optics context Mandel and Wolf (1995).

Given the 1D-NLSE evolution of individual realizations of the random field $\psi(x,t)$ the challenge is to determine the associated evolution of its statistical characteristics such as the probability density function (PDF) of the amplitude $|\psi|$, the power spectrum $|\psi_{k}|^{2}$ etc. The particular objective of this paper is to determine the short-time evolution of the normalized fourth moment $\kappa_{4}(t)$ defined as

[TABLE]

where the brackets $\langle\dots\rangle$ denote ensemble average that is made over a large number of realizations of the random process $\psi(x,t)$. The independence of Fourier modes composing the signal implies statistical homogeneity, hence the average over space $\lim\limits_{L\to\infty}\frac{1}{L}\int_{0}^{L}f(\psi)\,dx$ of an arbitrary function $f(\psi)$ of the random complex field $\psi$ is equivalent to the ensemble average over a large number of realizations. Nevertheless, it is convenient for our purposes in this paper to use the double (space and ensemble) averaging in the definition (9) of $\kappa_{4}$. We also note that the limit of a large box ($L\to\infty$) corresponds to the situation when a large number of random Fourier components are included inside the core part of the Fourier spectrum, thus fulfilling the pre-requisites of the central limit theorem for the initial condition.

The fourth moment (9) is an important characteristic of the PDF of a random process that quantifies the “heaviness“ of its tail. In particular, it can be used to characterize the deviation from Gaussianity in the course of evolution, when the initial statistics is Gaussian, in which case $\kappa_{4}$ is known to be equal $2$ Onorato et al. (2016); El Koussaifi et al. (2018). The determination of $\kappa_{4}$ is particularly relevant to the rogue wave studies as the formation of a “heavy tail” of the PDF is associated with the frequent appearance of large-amplitude events in the random process’ realizations Onorato et al. (2013); Walczak et al. (2015); Onorato et al. (2016); Randoux et al. (2016).

Fig. 1 shows a typical initial evolution of a random wave in the regime where the cubic (Kerr) nonlinearity dominates linear dispersive effects, which corresponds to the semiclassical regime described by Eq. (1) with $\varepsilon\ll 1$ (in the numerics we took $\varepsilon=0.1$). As shown in Fig. 1, the self-focusing dynamics tends to produce bright peaks while the self-defocusing dynamics leads to a decrease of the peak amplitudes but is accompanied by steepening of slopes in the random amplitude profile. While only the short-time evolution of the wave system is shown in Fig. 1, a longer development leads to the formation of gradient catastrophes – the explosion of the first derivatives of the wave’s profile. These gradient catastrophes have qualitatively different geometrical nature in the defocusing regime (the wave-breaking singularity Pomeau et al. (2008)) and the focusing regime (the elliptic umbilic singularity Dubrovin et al. (2008)). In both cases the gradient catastrophes are regularized by dispersive effects via the generation of nonlinear short wavelength oscillations: breather structures in the focusing regime Bertola and Tovbis (2013) and dispersive shock waves in the defocusing regime (see El and Hoefer (2016) and references therein). For convenience, we shall call the initial nonlinear evolution preceding the formation of gradient catastrophes, the “pre-breaking stage” in both defocusing and focusing regimes. The advantage of the semi-classical, dispersive-hydrodynamic approach employed in this paper is that it enables one to asymptotically separate the pre-breaking and post-breaking stages of the evolution, which exhibit qualitatively different behaviors and require very different analytical methods for their descriptions.

The starting point of our analysis is the evolution of the nonlinear part $H_{\rm NL}$ of the Hamiltonian. Differentiating (8) we obtain

[TABLE]

Using Eq. (1) and integrating by parts, one readily finds

[TABLE]

Now, using the Madelung transformation

[TABLE]

Eq. (11) can be rewritten as

[TABLE]

Noticing from (9), (4), (8) that

[TABLE]

one obtains

[TABLE]

We now derive an analytical expression for $\kappa_{4}(t)$ for short evolution times, $t\ll 1$. If the dispersion parameter is small, $\varepsilon\ll 1$, the initial dynamics are dominated by nonlinearity. To describe these dynamics analytically, we consider the semi-classical limit of the 1D-NLSE (1) which is found by applying the Madelung transform (12) and letting $\varepsilon\to 0$. Assuming smooth evolution of $\rho(x,t)$ and $u(x,t)$ in the pre-breaking regime, we obtain in the limit $\varepsilon\to 0$ the following well-known set of nonlinear geometric optics equations Forest et al. (2009); Wabnitz et al. (2013); Fatome et al. (2014); Kodama and Wabnitz (1995); Randoux et al. (2017)

[TABLE]

If $\sigma=-1$, equations (16) are identical to the shallow-water equations for an incompressible fluid with $\rho>0$ and $u$ interpreted as the fluid depth and the depth-averaged horizontal fluid velocity, respectively. In the nonlinear fiber optics context, $\rho$ represents the instantaneous optical power and $u$ represents the instantaneous frequency (or chirp) Wetzel et al. (2016).

Rigorous proofs of the point-wise convergence, as $\varepsilon\to 0$, of solutions of the 1D-NLSE (1) to the solutions of the dispersionless system (16) with the same initial data, prior to the formation of gradient catastrophe, can be found for certain classes of initial data in Jin et al. (1999) (defocussing) and in Kamvissis et al. (2003); Tovbis et al. (2004) (focussing). Some important particular exact solutions of system (16) for the focusing case have been found as early as in 1960-70s (see Talanov (1965); Akhmanov et al. (1966); Gurevich and Shvartsburg (1970)). A detailed mathematical analysis of the pre-breaking dynamics in the defocusing case can be found in Forest et al. (2009) (see also Moro and Trillo (2014) for the special case of the wave breaking into vacuum).

It follows from the above consideration that, to study the pre-breaking dynamics of partially coherent waves in 1-D NLSE (1) we need to be able to describe random solutions of system (16) obtained by evolving initial data $\rho(x,0)$, $u(x,0)$ with given statistics (e.g. corresponding to the Gaussian statistics of $\psi$). The study of such solutions has been recently initiated for the defocusing case in the context of the interaction of random Riemann waves in fiber optics Randoux et al. (2017) (see also Suret et al. (2017)). In this connection one must stress that the term “pre-breaking dynamics” is understood here in the probabilistic sense, as for random initial data there is always a non-zero probability of having gradient catastrophe at any, arbitrarily small, moment of time. However, due to initial data having typical size $\Delta x=\mathcal{O}(1)$, we assume that for small $\varepsilon$ the contribution of such early gradient catastrophes to the statistics is negligibly small.

To this end, with the short-time, pre-breaking evolution in mind, we look for the solutions of Eqs. (16) in the form of the time power series expansions for the realizations $\rho(x,t)$ and $u(x,t)$:

[TABLE]

In (17) we assumed that initially, $u_{0}=0$, which agrees with typical physical condition $u_{0}\ll\rho_{0}$ satisfied in standard realistic experimental conditions. Indeed timescales of amplitude and phase in partially coherent waves are generally similar ($\mathcal{O}(1)$ here). Considering the normalizations given by the Eq. (12), this means that the derivative of the phase $\frac{\partial(\phi_{0}/\epsilon)}{\partial x}=\mathcal{O}(1)$ and thus $u_{0}=\frac{\partial\phi_{0}}{\partial x}=\mathcal{O}(\epsilon)$ whereas $\rho_{0}=\mathcal{O}(1)$. This assumption is for example satisfied in the experiments on the propagation of partially coherent light through optical fibers, see Randoux et al. (2017).

Substituting Eqs. (17) into Eqs. (16) we obtain

[TABLE]

Next, substituting Eqs. (18) into Eq. (15) and integrating in time, we obtain the following expression for the time evolution of the normalized fourth moment of the field amplitude:

[TABLE]

Eq. (19) is our main general result. One can make now two important observations. The first one is that Eq. (19) shows that the normalized fourth-order moment $\kappa_{4}(t)$ of the field evolves quadratically with time at leading order for $t\ll 1$. The second observation is that Eq. (19) shows that the increasing or decreasing nature of the time evolution of $\kappa_{4}(t)$ is determined by the value taken by $\sigma$. In the focusing regime ($\sigma=+1$), $\kappa_{4}(t)$ is an increasing function of time which means that the nonlinear evolution of the wave field is characterized by PDFs that exhibit heavy tailed deviations from the initial statistical distribution. On the other hand, in the defocusing regime ($\sigma=-1$) $\kappa_{4}(t)$ becomes a decreasing function of time which implies low-tailed deviations from the initial statistics occurring in this regime. The statistical features described by Eq. (19) are in full qualitative agreement with the results that have been recently obtained in numerical and experimental investigations of integrable turbulence Randoux et al. (2016, 2014); Walczak et al. (2015); Randoux and Suret (2016); Suret et al. (2016, 2017).

Let us emphasize that the decreasing or increasing nature of the time evolution of $\kappa_{4}$ has also been shown to be determined by the defocusing or focusing nature of the propagation regime for weakly nonlinear dispersive random waves that are described by the 1D-NLSE Janssen (2003). Theoretical approaches that have been used in the weakly nonlinear regime are based on the wave turbulence theory and they consist in deriving quasi-kinetic equations for the lowest order moments of the wave field Suret et al. (2011); Janssen (2003). Dispersion plays crucial role in that consideration. Our work is based on a completely different, dispersive-hydrodynamic approach, where dispersive effects are initially not of dominant but of perturbative nature.

III Initial conditions with Gaussian statistics

Eq. (19) represents a general result that is derived with the only assumption that $\varepsilon\ll 1$. As we already stressed, it is valid before the typical time of the gradient catastrophe occurrence, i.e. for $t\ll 1$ (for random initial conditions with typical scales for $\rho$- and $x$-variations at $\mathcal{O}(1)$). Importantly, Eq. (19) is derived without any assumption on the nature of the initial statistics of the random wave field. In this section we show that Eq. (19) can be further simplified if the random wave field taken as initial condition has Gaussian statistics. To this end, we assume that the random initial field $\psi(x,0)$ is composed of a linear superposition of a large number of independent random Fourier modes $\psi_{k}(t=0)=\psi_{0k}=|{\psi_{0k}}|e^{i\phi_{0k}}$, so that by the central limit theorem $\psi(x,0)$ is a Gaussian random field  Nazarenko (2011).

In the random phase and amplitude (RPA) model, $|\psi_{0k}|$ and $\phi_{0k}$ are both taken as randomly-distributed variables Nazarenko (2011). Here, we will mainly use the so-called random phase (RP) model in which only the phases $\phi_{0k}$ of the Fourier modes are considered as being random Nazarenko (2011). In this model, the phase of each Fourier mode is randomly and uniformly distributed between $-\pi$ and $\pi$. Moreover, the phases of separate Fourier modes are assumed to be uncorrelated so that $\langle e^{i\phi_{0k}}e^{i\phi_{0k^{\prime}}}\rangle=\delta_{k}^{k^{\prime}}$. In the above expression, the brackets, as usual, represent the averaging over an ensemble of many realizations of the random process; $\delta_{k}^{k^{\prime}}$ is the Kronecker symbol defined by $\delta_{k}^{k^{\prime}}=1$ if $k=k^{\prime}$ and $\delta_{k}^{k^{\prime}}=0$ if $k\neq k^{\prime}$. With the assumptions of the RP model described above, the statistics of the initial field is homogeneous, which means that all statistical moments of the initial complex field $\psi(x,t=0)=\psi_{0}(x)$ do not depend on $x$ Picozzi (2007); Picozzi et al. (2014). This RP description of the initial random field has been shown to describe in a satisfatory way many experiments performed in the field of integrable turbulence Randoux et al. (2016, 2014); Walczak et al. (2015); Randoux and Suret (2016); Suret et al. (2016, 2017); Onorato et al. (2004, 2005a, 2001, 2005b).

Given the delta-correlation of the random phases, the second moment of a field composed from the linear superposition of a large number of independent Fourier components having Gaussian statistics is readily evaluated as

[TABLE]

and the sixth moment can be factored into products of the second moments by using Wick’s decomposition Nazarenko (2011)

[TABLE]

Now, using Eq. (21), one can evaluate the coefficient for the $\mathcal{O}(t^{2})$ term in the expansion (19):

[TABLE]

where we have used the notation $\psi_{0k_{i}}=\psi_{k_{i}}(0)$ for the Fourier component at $t=0$. Using Eq. (21), we obtain the following expression for the short time evolution of the fourth-order moment of a random wave field that has Gaussian statistics at initial time, i.e. $\kappa_{4}(0)=2$:

[TABLE]

where $n_{0k_{i}}=n_{k_{i}}(0)$ are the components of the power spectrum at $t=0$. Using Eqs. (5), (4), (7) and taking into account that $P=0$ for our random Gaussian field, we can finally rewrite Eq. (23) as

[TABLE]

(Note that $H_{\rm L}=\mathcal{O}(1)$ and $H_{\rm NL}=\mathcal{O}(1)$ whereas the linear part of the Hamiltonian (6) is $\mathcal{O}(\epsilon^{2})$).

A similar, but somewhat lengthy, calculation permits one to obtain a more accurate expression that includes $\mathcal{O}(t^{4})$ correction (see Appendix VI):

[TABLE]

Eqs. (24) and (25) show that the time evolution of the fourth moment of the initially Gaussian random wave field is determined by the linear part $H_{\rm L}(0)$ of the Hamiltonian computed for the initial condition.

Eq. (25) can be further simplified if we assume that the shape of the Fourier power spectrum of the initial random field is described by a Gaussian

[TABLE]

The amplitude $n_{0}\in\mathbb{R^{+}}$ has to be determined from the normalization condition provided by Eq. (4). The linear energy density determined from Eq. (7) is $\langle H_{\rm L}\rangle=\frac{N(\Delta k)^{2}}{4}$, and we can finally rewrite Eq. (25) only in terms of the density of particles (or optical power) $N$ and of the width $\Delta k$ of the initial Fourier spectrum

[TABLE]

IV Numerical simulations

In this section, we use numerical simulations of Eq. (1) to investigate the range and the degree of validity of the semi-classical approach to the statistics of integrable turbulence presented in Secs. II and III. The initial condition used in our numerical simulations is a random complex field having Gaussian statistics. The amplitudes of the Fourier components are taken to be distributed according to Eq. (26). In our numerical simulations, the spectral phases $\phi_{0k}$ are random, statistically independent real numbers, uniformly distributed between $-\pi$ and $+\pi$. The width $\Delta k$ of the initial spectrum profile (26) is taken to be unity ($\Delta k=1$), and the value of $\varepsilon$ in (1) is taken to be $0.1$. The numerical simulations are performed by using a pseudo-spectral method with the numerical box having size $L=256$ that is discretized by using $2^{16}$ points.

Fig. 2 shows the time evolution of the normalized fourth moment $\kappa_{4}(t)$ of the random wave field in the focusing regime ($\sigma=+1$). The curve plotted with black line represents the result of the numerical simulation of Eq. (1). At time $t=0$, the theoretical value taken by $\kappa_{4}(t)$ is $2$ due to the Gaussian nature of the initial statistics. As can be seen from Fig. 2, in the numerical experiments the value taken by $\kappa_{4}(t)$ at $t=0$ slightly differs from $2$ (see the inset in Fig. (2)) because the conditions of the central limit theorem are not perfectly fulfilled in our numerical simulations. Indeed, because of the finite value of $L$, the number of Fourier modes in the spectrum given by Eq. 26 is finite, in particular, in the Full Width at Half Maximum (FWHM) we count $\frac{\sqrt{2\ln 2}L}{\pi}=95$ modes. Importantly, the deviation of the initial condition from Gaussian statistics affects only the value of $\kappa_{4}(0)$ in Eq. (21) but not the evolution. In the initial (before the formation of a gradient catastrophe) stage of the nonlinear evolution of the random wave, $\kappa_{4}(t)$ is at first an increasing function of time that later reaches a maximum around $t\sim 0.6$. Then $\kappa_{4}(t)$ becomes a decaying function of time that reaches a stationary value around $\sim 4$ at long evolution time. A similar evolution of $\kappa_{4}(t)$ has already been evidenced in numerical simulations presented in ref. Onorato et al. (2016).

The curves plotted with green and red lines in Fig. 2 show monotonic evolutions of $\kappa_{4}(t)$ that are obtained from Eq. (27). In particular the curves plotted in the inset of Fig. 2 clearly reveal a very good quantitative agreement between numerical and theoretical results. In particular, a better agreement between numerics and theory is obtained by including the fourth-order correction term found in the time expansion of the solution, see Eq. (27). A significant quantitative disagreement is found between our theoretical results and the numerical simulation at evolution times greater than $\sim 0.2$. This arises from the fact that our approach is only valid at evolution times that are shorter that the typical wave breaking time (the pre-breaking description). The significant occurrence of wave breakings at evolution times greater than $\sim 0.2$ has strong influence on the wave evolution and subsequently the wave statistics in a way that cannot be accounted for by using our pre-breaking treatment.

Fig. 3 shows the comparison between the numerical simulation of Eq. (1) and the theoretical result given by Eq. (27) in the defocusing regime ($\sigma=-1$). In the defocusing regime, $\kappa_{4}(t)$ is a monotonically decreasing function of time, as already evidenced in ref. Onorato et al. (2016). As for the focusing regime, a very good quantitative agreement is obtained between numerics and the theory at short evolution time ($t<0.2$), i.e. before the typical occurrence of gradient catastrophes.

V Conclusion

In this paper we have undertaken an analytical study of the problem of the evolution of a random wave field in the 1D-NLSE for both focusing and defocusing regimes. This has been done from the perspective of dispersive hydrodynamics, a semiclassical theory of nonlinear dispersive waves exhibiting two distinct spatio-temporal scales: the long scale specified by initial conditions and the short scale by the internal coherence length (i.e., the typical size $\varepsilon$ of the coherent structures) Biondini et al. (2016). This scale separation enabled us to split the time evolution of the nonlinear random wave system (integrable turbulence) into the initial, “pre-breaking” stage, preceding the formation of gradient catastrophes, when the evolution of the 1D NLSE wave field is almost everywhere smooth, and the “post-breaking” stage characterised by the generation of short-scale nonlinear oscillations (breathers or dispersive shock waves depending on the focusing vs. defocusing character of the 1D-NLSE).

Our work is concerned with the initial, pre-breaking stage of the semi-classical integrable turbulence, when the dynamical and statistical features can be analytically described in terms of random solutions of the dispersionless (nonlinear geometric optics) system (16). As a result, we have derived a simple asymptotic formula describing the evolution of the normalized fourth moment of the random wave field. This formula, applied to the problem of the 1D-NLSE evolution of random field having initial Gaussian statistics, describes the initial stage of the formation of heavy tails of the PDF of the field amplitude in the focusing case and the formation of low tails in the defocusing case.

Recently, an exact and general identity that relates the changes in the statistical properties of the wave field to the changes of its Fourier spectrum has been derived by using the Hamiltonian structure of 1D-NLSE Onorato et al. (2016). In other words, the knowledge of the fourth order moment also provides the description of spectral properties. The general description of the stationary state of integrable turbulence and the theoretical prediction of the fourth order moment is still an open fundamental question. In the weakly nonlinear regime, the wave turbulence approach provides a statistical description of the nonlinear propagation of random wave fields in 1DNLSE systems Janssen (2003); Mori and Janssen (2006); Soh et al. (2010); Suret et al. (2011). Recently, using an approach based on the so-called large deviation theory, it has been shown that rogue waves obey a large deviation principle, i.e. the heavy tails of the PDF of the random wave field are dominated by single realizations Dematteis et al. (2018a, b). This approach is very promising but does not provide a simple formula for the evolution of the statistics. In this article we have demonstrated that the semi-classical approach is an extremely powerful tool enabling one to describe in a simple way the early stage of integrable turbulence in the strongly nonlinear/ small dispersion regime. The proposed methodology can be applied to the description of partially coherent random nonlinear waves described by other integrable equations, including shallow water waves described by the KdV equation and its extensions. In particular, the pre-breaking statistics of bi-directional random shallow water waves is equivalent to that described by the defocusing 1D-NLSE and studied in this paper.

The semi-classical approach to the statistics of random waves in integrable systems is general and can be used beyond the short-time asymptotic regime. It is known very well that, in the semi-classical limit the evolution of nonlinear dispersive waves after the gradient catastrophe point is described by the so-called Whitham modulation equations Whitham (1974) governing the behaviour of the averaged integrals of motion, and replacing the dispersionless system (16) (see Kamchatnov (2000); El and Hoefer (2016); El et al. (2016); Tovbis and El (2016) and references therein for the application of the Whitham theory to the defocusing and focusing 1D-NLSE). Such an extension of the proposed method to longer times is very promising but also highly nontrivial.

VI Appendix: Computation of $\mathcal{O}(t^{4})$ corrections for the case of Gaussian statistics at $t=0$

Here we provide the simplified expressions for the three terms that are found in the integral giving the coefficient of the $O(t^{4})$ term in Eq. (19). To obtain these expressions, we assume Gaussian statistics at the initial time and use Eq. (20) and Eq. (21) to obtain

[TABLE]

[TABLE]

[TABLE]

Acknowledgements.

This work was supported by EPSRC grant EP/R00515X/1 (GE) and Dstl grant DSTLX-1000116851(GR, GE, SR). It has also been partially supported by the Agence Nationale de la Recherche through the LABEX CEMPI project (ANR-11-LABX-0007) and by the Ministry of Higher Education and Research, Hauts-De-France Regional Council and European Regional Development Fund (ERDF) through the Contrat de Projets Etat-Région (CPER Photonics for Society P4S). GE and GR thank the PhLAM laboratory at the University of Lille for hospitality and partial financial support (GE).