Nonstationary distributions of wave intensities in wave turbulence

Wave turbulence (WT) theory describes random weakly nonlinear wave systems with broadband spectra (see e.g. [1]). The main object in this theory is the wave action spectrum which is the second-order moment of the wave amplitude evolving according to the so-called wave-kinetic equation. Specific attention in the past was given to the stationary power-law solutions of this equation which are analogous to the Kolmogorov spectrum of hydrodynamic turbulence—the Kolmogorov–Zakharov spectra. However, as it was shown in [1, 3–6], the WT approach can also be extended to describing the evolution of the higher-order moments and even to the entire probability density function (PDF) of the wave amplitude. A formal justification of such an extension, based on a rigorous statistical formulation, was presented in [7]. An introduction to the WT, as well as a summary of recent developments, can be found in [1] and in an older text [2].

It has been widely believed that the statistics of random weakly nonlinear wave systems is close to being Gaussian. Derivation of the evolution equation for the PDF of the wave intensities presented in [5] has made it possible to examine this belief. It was shown that such an equation does have a stationary solution corresponding to the Gaussian state. It was also noted that the typical evolution time of the PDF is the same as the one for the wave spectrum. Thus, for non-stationary wave systems one can expect significant deviations from Gaussianity if the initial wave distribution is non-Gaussian. Note that non-Gaussian (usually deterministic) initial conditions for the wave intensity are typical in numerical simulations of WT. Also, there is no reason to believe that initial waves excited by natural phenomena, e.g. sea waves excited by wind, should be Gaussian. Therefore, the study of the evolution of the wave statistics is important for both understanding the fundamental nonlinear processes, as well as for practical predictions such as, for instance, the wave weather forecast.

In the present paper, we will present the full general solution for the PDF equation derived in [5]. Based on this solution we will formulate a condition under which the wave statistics relaxes to the Gaussian state.

Consider a weakly nonlinear wave system dominated by the four-wave interactions bounded by an L-periodic cube in a d-dimensional physical space. (The four-wave system is considered here as an illustrative example only. All results of this paper hold for general N-wave systems for any N. The only difference will be in the expressions for $\gamma_{\bf k}$ and $\eta_{\bf k}$ below; see [1].) The system evolves according to the following Hamiltonian equations for the Fourier coefficients,

$$\begin{align} \displaystyle {\rm i}\dot{a_{\bf k}} = \frac{\partial \mathcal{H}}{\partial {a}_{\bf k}^*}, \quad \mathcal{H}=\sum_{\bf k} \omega_{\bf k} \vert {a}_{\bf k}\vert ^2 + \frac{1}{2}\sum_{{\bf k}_1, {\bf k}_2,{\bf k}_3,{\bf k}_4} W^{{\bf k}_1,{\bf k}_2}_{{\bf k}_3,{\bf k}_4}{a}_{{\bf k}_1}^* {a}_{{\bf k}_2}^* a_{{\bf k}_3} a_{{\bf k}_4} \label{4wH}, \nonumber \end{align} \tag{ 1 }$$

where ${\bf k}, {\bf k}_1, {\bf k}_2, {\bf k}_3, {\bf k}_4 \in \frac L {2\pi} \mathbb{Z}^d$ are wave vectors, $a_{\bf k}(t) \in \mathbb{C}$ is the wave action variable, $W^{{\bf k}_1, {\bf k}_2}_{{\bf k}_3, {\bf k}_4} \in \mathbb{R}$ is an interaction coefficient which is a model-specific function of ${{\bf k}_1, {\bf k}_2, {\bf k}_3, {\bf k}_4}$ (e.g. $W^{{\bf k}_1, {\bf k}_2}_{{\bf k}_3, {\bf k}_4} \equiv 1$ for the Gross–Pitaevskii equation) and $\omega_{\bf k} \in \mathbb{R}$ corresponds to the frequency of mode $\bf k$.

Let us consider the PDF ${{\mathcal P}}(t, s_{\bf k})$ of the wave intensity $J_{\bf k}= (L/2\pi){\hspace{0pt}}^d \vert a_{\bf k}\vert ^2$ defined in a standard way so that the probability for $J_{\bf k}$ to be in the range from $s_{\bf k}$ to $s_{\bf k} +ds_{\bf k}$ is ${{\mathcal P}}(t, s_{\bf k}) {\rm d} s_{\bf k}$. In symbolic form this is represented as

$$\begin{align} \displaystyle {{\mathcal P}}(t,s_{\bf k}) = \langle \delta (s_{\bf k}-J_{\bf k}) \rangle. \label{pasdelta} \nonumber \end{align} \tag{ 2 }$$

Suppose that the waves are weakly nonlinear, so that the quartic part of the Hamiltonian is much less than its quadratic part. Further suppose that the initial complex wave amplitudes $a_{\bf k}$ are independent random variables for each $\bf k$ and that the initial phases of $a_{\bf k}$ are random and equally probable in the range from 0 to $2\pi$. Fields with such statistical properties are called random phase and amplitude (RPA) fields—they are the main objects of WT theory (see [1, 5]). An important fact pointed out in [5] is that the RPA property does not imply Gaussianity and, therefore, it makes sense to study the evolution of the PDF. Assuming the initial conditions of a RPA type, taking the infinite-box limit $L \to \infty$ and the weak amplitude limit, the WT derivation leads to the following evolution equation for ${{\mathcal P}}(t, s_{\bf k})$ (see [1, 5]):

$$\begin{align} \displaystyle \frac{\partial {{\mathcal P}}(t, s_{\bf k})}{\partial t} +\frac{\partial F}{\partial s_{\bf k}}=0,\label{main} \nonumber \end{align} \tag{ 3 }$$

where

$$\begin{align} \displaystyle \label{F} F = -s_{\bf k} \Big(\gamma_{\bf k} {{\mathcal P}} +\eta_{\bf k} \frac{\partial {{\mathcal P}} }{\partial s_{\bf k}}\Big), \nonumber \end{align} \tag{ 4 }$$

with, for the four-wave systems,

$$\begin{align} \displaystyle \eta_{\bf k}(t) = 4\pi \int\vert W^{{\bf k},{\bf k}_1}_{{\bf k}_2,{\bf k}_3}\vert ^2\delta({{\bf k} + {\bf k}_1 -{\bf k}_2 -{\bf k}_2})\delta(\omega_{{\bf k}} + \omega_{{\bf k}_1} - \omega_{{\bf k}_2} -\omega_{{\bf k}_3}) n_{{\bf k}_1} n_{{\bf k}_2} n_{{\bf k}_3} {\rm d}{\bf k}_1{\rm d}{\bf k}_2{\rm d}{\bf k}_3, \label{eta} \nonumber \end{align} \tag{ 5 }$$

$$\begin{align} \displaystyle \gamma_{\bf k}(t) = \,\,&8\pi \int\vert W^{{\bf k},{\bf k}_1}_{{\bf k}_2,{\bf k}_3}\vert ^2\delta({{\bf k} + {\bf k}_1 -{\bf k}_2 -{\bf k}_2})\delta(\omega_{{\bf k}} + \omega_{{\bf k}_1} - \omega_{{\bf k}_2} -\omega_{{\bf k}_3}) \nonumber \\ &\times \Big[n_{{\bf k}_1}(n_{{\bf k}_2}+n_{{\bf k}_3})-n_{{\bf k}_2}n_{{\bf k}_3}\Big]{\rm d}{\bf k}_1{\rm d}{\bf k}_2{\rm d}{\bf k}_3, \label{gam} \nonumber \end{align} \tag{ 6 }$$

where $n_{\bf k} = \langle J_{\bf k} \rangle$ is the wave action spectrum. The infinite-box limit results in passing to the continuous wave number description; each wave number integration in the above equations is over $\mathbb{R}^d$.

Note that the PDF equation was formally derived for the moment when $t=0$. However, it is assumed to be valid over a typical evolution time during which ${{\mathcal P}}$ may experience significant changes. This assumption is rather common for the kinetic theory—it is called the 'propagation of chaos'. Discussion of the validity of such an assumption in WT can be found in [1].

In this paper, we will find the general time-dependent solution of the PDF equation (3). The problem can be split into two independent parts. In part I, we find the spectrum $n_{\bf k} \equiv \int_0^\infty s {{\mathcal P}}(t, s_{\bf k}) \, {\rm d} s$ by solving the wave-kinetic equation

$$\begin{align} \displaystyle \dot n_{\bf k} = \eta_{\bf k} - \gamma_{\bf k} n_{\bf k} \label{eq:ndot} \nonumber \end{align} \tag{ 7 }$$

(the latter follows from taking the first moment of equation (3)), and in part II we use the solution for $n_{\bf k}(t)$ in the PDF equation (3) and solve the latter. Below, we will deal with part II only thereby assuming that the solution of the wave-kinetic equation exists and that the equation itself remains valid for a non-zero period of time $[0, t_{\rm max})$ (with $t_{\rm max} \in \mathbb{R}_+$ or $t_{\rm max} =+\infty$) and, therefore, the integrals defining $\eta_{\bf k}$ and $\gamma_{\bf k}$ are well-defined and convergent.

The initial function ${{\mathcal P}}(0, s_{\bf k})$ must satisfy the following properties:

• 1.  
  $\int_0^\infty {{\mathcal P}}(0, s_{\bf k}) \, {\rm d} s =1.$ This is a probability normalization property which is preserved by equation (3).
• 2.  
  $n_{\bf k}(0) = \int_0^\infty s {{\mathcal P}}(0, s_{\bf k}) \, {\rm d} s <\infty.$ This condition is implied by the existence of the solution of the wave-kinetic equation. The higher moments of ${{\mathcal P}}(0, s_{\bf k})$ are allowed to be, in principle, infinite.

Let us introduce the generating function

$$\begin{align} \displaystyle \mathcal{Z}(t,\lambda_{\bf k}) = \langle {\rm e}^{-\lambda_{\bf k} J_{\bf k}(t)} \rangle = \int^{\infty}_{0}{{\mathcal P}}(\lambda_{\bf k},t) {\rm e}^{-\lambda_{\bf k} s_{\bf k}} {\rm d}s_{\bf k} \label{eq:Z} \nonumber \end{align} \tag{ 8 }$$

where $\lambda_{\bf k}$ is a real parameter. Existence of $\mathcal{Z}(t, \lambda_{\bf k})$ is guaranteed by the integrability of $\mathcal{P}$. Note that this definition is different from the one used in [5] by the sign of the exponent. Here, we have changed the sign in order to comply with the standard relation between $\mathcal{P}$ and $\mathcal{Z}$ via the Laplace transform, as expressed in equation (8).

In what follows we will drop subscript $\bf k$ for brevity whenever it does not cause ambiguity.

The inverse Laplace transformation of $\mathcal{Z}(t, \lambda)$ gives:

$$\begin{align} \displaystyle {{\mathcal P}}(t,s) = \frac{1}{2\pi {\rm i}}\lim_{T\to \infty} \int^{T+i\infty}_{T-i\infty}{{\mathcal Z}}(\lambda,t) {\rm e}^{s\lambda} {\rm d}\lambda. \nonumber \end{align} \tag{ 9 }$$

Given ${{\mathcal Z}}$, one can easily calculate the moments of the wave intensity,

$$\begin{align} \displaystyle M^{(\,p)}_{\bf k} \equiv \langle J_{\bf k}^{\,p}\rangle =(-1)^{\,p} {{\mathcal Z}}_{\lambda\cdots\lambda}\vert _{\lambda=0} = (-1)^{\,p} \langle J^{\,p} {\rm e}^{\lambda J^{\,p}}\rangle \vert _{\lambda=0}, \nonumber \end{align} \tag{ 10 }$$

where $p \in \mathbb{N}$ is the order of the moment and subscript λ means differentiation with respect to λ. In particular, for the waveaction spectrum we have

$$\begin{align*} \displaystyle n_{\bf k} = - {{\mathcal Z}}_{\lambda}\vert _{\lambda=0}. \nonumber \end{align*}$$

The evolution equation for $\mathcal{Z}$ can be obtained by taking the Laplace transform of the equation (8), which gives

$$\begin{align} \displaystyle \dot{\mathcal{Z}} =- \lambda\eta \mathcal{Z} -(\lambda^2\eta+\lambda\gamma)\mathcal{Z}_{\lambda}.\label{g\,f} \nonumber \end{align} \tag{ 11 }$$

Note that the sign differences in this equation with respect to the corresponding equation in [5] is due to the sign difference in our definition of $\mathcal{Z}$.

As we mentioned before, the moments $M^{(\,p)}_{\bf k}$ with $p\geqslant 2$ are allowed to be infinite. Thus, we require the function $\mathcal{Z}$ to be differentiable with respect to λ only once. If this property holds at $t=0$ then it will hold over the relevant period of time $[0, t_{\rm max})$ over which the solution of the wave-kinetic equation exists.

Previously in [5], the general steady state solution of equation (11) was presented:

$$\begin{align} \displaystyle {\mathcal{Z}} = \frac 1 {1+ \lambda_{\bf k} n_{\bf k}}. \label{g\,fst} \nonumber \end{align} \tag{ 12 }$$

This solution corresponds to Gaussian statistics of the wave field (Rayleigh distribution for the wave intensity, respectively).

Below, we will study the evolving system, in which the parameters $\eta, \gamma$ are time-dependent. The goal of this paper is first to find the solution of equation (11) and then obtain the respective time-dependent PDF.

We can find the general solution of equation (11) by using the method of characteristics. Rewriting this equation in the characteristic form we have:

$$\begin{align} \displaystyle \frac{{\rm d}\lambda}{{\rm d}t}= \Big({\gamma}+\lambda \eta \Big)\lambda, \quad \frac{{\rm d}\mathcal{Z}}{{\rm d}t}=- \lambda\eta\mathcal{Z}.\label{cm} \nonumber \end{align} \tag{ 13 }$$

Changing variable to $\mu = \lambda {\rm e}^{-\int_0^t {\gamma(t^{\prime})}{\rm d}t^{\prime}}$ in the first of these equations, we transform it into

$$\begin{align} \displaystyle \frac{{\rm d}\mu(t)}{{\rm d}t} = \eta \mu \lambda =\eta \mu^2 {\rm e}^{\int_0^t {\gamma(t^{\prime})}{\rm d}t^{\prime}}. \label{cm1} \nonumber \end{align} \tag{ 14 }$$

Solving this equation we have

$$\begin{align} \displaystyle -\frac 1 {\mu(t)}+\frac 1 {\mu_0} = \int^t_0 \eta(t^{\prime}) {\rm e}^{\int^{t^{\prime}}_0 \gamma(t^{\prime\prime}){\rm d}t^{\prime\prime}}{\rm d}t^{\prime}, \nonumber \end{align} \tag{ 15 }$$

where $\mu_0 = \mu(0) = \lambda_0$. This gives for $\lambda(t)$:

$$\begin{align} \displaystyle \lambda(t) = \frac{\lambda_0 {\rm e}^{\int^t_0 \gamma(t^{\prime}){\rm d}t^{\prime}}}{1 - \lambda_0 \int^t_0 \eta(t^{\prime}) {\rm e}^{\int^{t^{\prime}}_0 \gamma(t^{\prime\prime}){\rm d}t^{\prime\prime}}{\rm d}t^{\prime}}.\label{slt} \nonumber \end{align} \tag{ 16 }$$

This relation has a simpler form in terms of n rather than η. Indeed, n satisfies equation (7), integrating which we have

$$\begin{align} \displaystyle n(t) =n(0) \, {\rm e}^{-\int^t_0 \gamma(t^{\prime}){\rm d}t^{\prime}} + {\rm e}^{-\int^t_0 \gamma(t^{\prime}){\rm d}t^{\prime} } \int^t_0 \eta(t^{\prime}) {\rm e}^{\int^{t^{\prime}}_0 \gamma(t^{\prime\prime}){\rm d}t^{\prime\prime}}{\rm d}t^{\prime}. \label{eq:n} \nonumber \end{align} \tag{ 17 }$$

Using this identity, we find

$$\begin{align} \displaystyle \lambda(t) = \frac{\lambda_0 }{{\rm e}^{-\int^t_0 \gamma(t^{\prime}){\rm d}t^{\prime}} - \lambda_0 \left(n(t)-n(0) {\rm e}^{-\int^t_0 \gamma(t^{\prime}){\rm d}t^{\prime}}\right)}.\label{slt1} \nonumber \end{align} \tag{ 18 }$$

Conversely, we have

$$\begin{align} \displaystyle \lambda_0 = \frac{\lambda {\rm e}^{-\int^t_0 \gamma(t^{\prime}){\rm d}t^{\prime}} }{1+ \lambda \left(n(t)-n(0) {\rm e}^{-\int^t_0 \gamma(t^{\prime}){\rm d}t^{\prime}}\right)}.\label{slt1i} \nonumber \end{align} \tag{ 19 }$$

Now, from the second of equations (13) and from the first equality in equation (14) we see that the log derivative of $\mathcal{Z}$ is equal to the negative log derivative of μ. Thus,

$$\begin{align} \displaystyle \mathcal{Z}(t, \lambda)=\mathcal{Z}_0 \frac {\mu_0} \mu= \frac {\mathcal{Z}_0 \lambda_0} \lambda {\rm e}^{\int_0^t {\gamma(t^{\prime})}{\rm d}t^{\prime}} = \frac {\mathcal{Z}_0 }{1+ \lambda \left(n(t)-n(0) {\rm e}^{-\int^t_0 \gamma(t^{\prime}){\rm d}t^{\prime}}\right)}.\label{solnZ} \nonumber \end{align} \tag{ 20 }$$

where $\mathcal{Z}_0 = \mathcal{Z}(0, \lambda_0)$ and $\lambda_0$ must be substituted in terms of λ solving for it from equation (19).

Equations (20) and (19) allow us to prove the following theorem.

Theorem. Let a WT system be described by the characteristic function ${\mathcal{Z}}(t, \lambda)$ satisfying equation (11) over a non-zero period of time $[0, t_{\rm max})$ (with $t_{\rm max} \in \mathbb{R}_+$ or $t_{\rm max} =+\infty$) with functions $\eta(t) \in C^1[0, t_{\rm max})$ and $\gamma(t)\in C^1[0, t_{\rm max})$ related to the spectrum $n(t)\in C^1[0, t_{\rm max})$ via equations (5) and (6) and $n(t)$ satisfying the wave-kinetic equation (7). We seek a solution for ${\mathcal{Z}}(t, \lambda)$ with an initial condition $\mathcal{Z}(0, \lambda) \in C^1(\mathbb{R})$. Then:

• 1.  
  Solution ${\mathcal{Z}}(t, \lambda) \in C^1\left([0, t_{\rm max}), \mathbb{R}\right)$ exists and is given by equation (20).
• 2.  
  Expression (12) is a stationary solution for $t\in [0, t_{\rm max})$. In other words, wave fields which are Gaussian initially will remain Gaussian for all time (for which a solution to the wave-kinetic equation exists).
• 3.  
  If $t_{\rm max} = \infty$, the wave statistics will asymptotically become Gaussian if

  $$\begin{align} \displaystyle \lim_{t \to \infty} \frac {n(0) {\rm e}^{-\int_0^t {\gamma(t^{\prime})}{\rm d}t^{\prime} }}{n(t)} = 0. \label{cond} \nonumber \end{align} \tag{ 21 }$$

To prove the first part, it suffices to note that the characteristics (18) are not self-intersecting and that expression (20) is differentiable in t and in λ.

To prove the second part we simply substitute $\mathcal{Z}_0 = 1/(1+ \lambda_0 n_0)$ into equation (20) and, after using (19), obtain $\mathcal{Z} = 1/(1+ \lambda n)$, which corresponds to Gaussian statistics.

To prove the third part we notice that if condition (21) is satisfied from equations (19), (20) and (8) we have:

$$\begin{align*} \displaystyle \lim_{t \to \infty} \lambda_0 (t, \lambda) &= \frac{\lambda }{{\rm e}^{\int^t_0 \gamma(t^{\prime}){\rm d}t^{\prime}} + \lambda n(0) \left[\frac{n(t)}{n(0)} {\rm e}^{\int^t_0 \gamma(t^{\prime}){\rm d}t^{\prime}} -1 \right] } = 0, \Longrightarrow \lim_{t \to \infty} \mathcal{Z}_0 \nonumber \\ &= \mathcal{Z}(0, 0) =1, \Longrightarrow \lim_{t \to \infty} \mathcal{Z}(t, \lambda) = \frac {1 }{1+ \lambda n}. \nonumber \end{align*}$$

Remarks 

• 1.  
  Condition (21) is satisfied for the inertial range modes in forced-dissipated systems converging towards the steady state. Indeed, in the case $\gamma \to \eta/n$, which is a positive constant (mode $\bf k$ being fixed), the time integral of this quantity diverges as $t \to \infty$.
• 2.  
  In the absence of forcing and dissipation, wave spectrum $n_{\bf k}$ decays to zero for any mode $\bf k$ as $t \to \infty$, and so does $\gamma_{\bf k}$. Thus the integral of $\gamma_{\bf k}(t)$ may converge as $t \to \infty$, which means that non-Gaussianity of some (or all) wave modes may persist as $t \to \infty$.
• 3.  
  In general, function $\gamma_{\bf k}(t)$ is not sign-definite, and there may be transient time periods where $\gamma_{\bf k}(t)<0$. The deviation from Gaussianity of some (or all) wave modes may increase during these periods.
• 4.  
  Recall that in this paper we assume that the wave-kinetic equation remains valid over a certain period of time. Thus, remarks 1 and 2 only apply to the systems for which all assumptions of the WT theory survive for an infinite time. We will discuss this point further in the conclusions section.

Now let us analyse the PDF corresponding to evolving systems. Let us think of a simple case with a deterministic initial wave intensity, $\mathcal{P}(0, s) = \delta(s - J)$. We will call such a solution $\mathcal{P}_J(t, s)$. Then $\mathcal{Z}(0, \lambda) = {\rm e}^{-\lambda J}$. In fact, since the inverse Laplace transform is a linear operation, the considered solution is nothing but the Green's function for the general problem with an arbitrary initial condition $P(0, s)$:

$$\begin{align} \displaystyle \mathcal{P}(t,s) =\int_0^\infty \mathcal{P}(0,J) \mathcal{P}_J(t,s) {\rm d}J. \label{solnPgen} \nonumber \end{align} \tag{ 22 }$$

Let us take the inverse Laplace transform of $\mathcal{Z}(t, \lambda)$ given by equation (20) to obtain $\mathcal{P}_\delta(t, s)$ at $t>0$:

$$\begin{align} \displaystyle \mathcal{P}_J(t,s) = \frac{1}{2\pi {\rm i}}\lim_{T\to +\infty} \int^{T+i\infty}_{T-i\infty} {\rm e}^{\lambda s}\mathcal{Z}(\lambda) {\rm d}\lambda = \frac{1}{2\pi {\rm i}} \lim_{T\to +\infty} \int^{T+i\infty}_{T-i\infty} \frac {{\rm e}^{\lambda s-\lambda_0 J} }{1+ \lambda \tilde n } {\rm d}\lambda, \label{solnPst} \nonumber \end{align} \tag{ 23 }$$

where

$$\begin{align} \displaystyle \tilde n = n(t)-J {\rm e}^{-\int^t_0 \gamma(t^{\prime}){\rm d}t^{\prime}} \nonumber \end{align} \tag{ 24 }$$

(note that $n(0) =J$). Substituting $\lambda_0$ from (18) and changing the integration variable to $\rho = \lambda +1/\tilde n$, we have:

$$\begin{align} \displaystyle \mathcal{P}_J(t,s) = \frac{{\rm e}^{-\frac{s}{\tilde n} - a \tilde n}}{2\pi {\rm i} {\tilde n}} \lim_{T\to +\infty} \int^{T+i\infty}_{T-i\infty} \frac{{\rm e}^{s \rho +\frac a \rho} }\rho {\rm d} \rho =\frac 1 {{\tilde n}} {{\rm e}^{-\frac{s}{\tilde n} - a \tilde n}} I_0(2\sqrt{as}). \label{solnPst1} \nonumber \end{align} \tag{ 25 }$$

where $a = \frac{J}{{\tilde n}^2} {\rm e}^{-\int^t_0 \gamma(t^{\prime}){\rm d}t^{\prime}}$ and $I_0(x)$ is the zeroth modified Bessel function of the first kind. Note that $I_0(0)=1$, so we recover $\mathcal{P}_\delta \to \mathcal{P}_G = \frac 1 n {\rm e}^{-s/n}$ as $t \to \infty$ if condition (21) is satisfied provided that s is not too large, $as \ll 1$.

Now let us suppose that condition (21) is satisfied and let us consider the asymptotic behaviour of the probability distribution at large s and large t with $as \gg 1$ (i.e. s is much larger than $1/a$, which is itself large). Taking into account that $I_0(x) \mathop{\longrightarrow}\limits^{x\rightarrow \infty} \frac{{\rm e}^x}{\sqrt{2\pi x}}$, we have:

$$\begin{align} \displaystyle \mathcal{P}_J(t,s) \to \frac{\mathcal{P}_G}{(2\pi)^{1/2}(as)^{1/4}}{\rm e}^{2\sqrt{as} - as} \ll \mathcal{P}_G \quad {\rm for} \quad as\gg 1, \; \int^t_0 \gamma(t^{\prime}){\rm d}t^{\prime} \gg 1. \label{front} \nonumber \end{align} \tag{ 26 }$$

Thus, we see a front at $s \sim s^{*}(t) = 1/a$ moving toward large s as $t \to \infty$. The PDF ahead of this front is depleted with respect to the Gaussian distribution, whereas behind the front it asymptotes to the Gaussian. Obviously, the same kind of behaviour will be realised for any solution (22) arising from initial data having finite support in s.

In this paper we have obtained the general solutions for the generating function and for the PDF of wave intensities in WT, equations (20), (19) and (22), (25), respectively. This allowed us to prove a theorem stating that wave fields which are Gaussian initially will remain Gaussian for all time and that the wave statistics will asymptotically become Gaussian if condition (21) is satisfied. We have also found, subject to condition (21), an asymptotic solution for the PDF (26) where the Gaussian distribution forms behind a front propagating toward large wave intensities.

Condition (21) is satisfied for the inertial range modes in forced-dissipated systems approaching a steady state. Thus, the Gaussian statistics will form at large times for such modes in these systems. An interesting subclass of solutions in forced-dissipated systems is when the spectrum is in a steady state regime from the initial moment of time (i.e. it is a stationary solution of the wave-kinetic equation), while the PDF is not initially Rayleigh (i.e. the initial wave field is not Gaussian). For example, the initial wave intensities can be deterministic, i.e. their PDFs are delta-functions, as is often taken in numerical simulations of WT. In this case, equation (25) simplifies with $\int^t_0 \gamma(t^{\prime}){\rm d}t^{\prime} = \gamma t$.

Since the characteristic evolution times for the spectrum $n_{\bf k}$ and for the PDF are the same, the latter will remain non-Gaussian over a substantial time if the initial field is non-Gaussian. Such situations should be considered typical rather than being an exception in natural conditions (where initial waves arise, e.g. from an instability which does not necessarily produce Gaussian statistics) and in numerical simulations (where typically the wave intensities are taken to be deterministic).

Moreover, in the absence of forcing and dissipation, spectrum $n_{\bf k}$ decays to zero for any mode $\bf k$ as $t \to \infty$, and so does $\gamma_{\bf k}$. Thus the integral $\int^t_0 \gamma_{\bf k}(t^{\prime}){\rm d}t^{\prime}$ may converge as $t \to \infty$, which means that non-Gaussianity of some (or all) wave modes may persist as $t \to \infty$. Furthermore, since $\gamma_{\bf k}(t)$ is not sign definite, there may be transient time periods where $\gamma_{\bf k}(t)<0$. The deviation from Gaussianity of some (or all) wave modes may increase during these periods.

The present paper considers a four-wave system as an illustrative example, but it is clear that the obtained solutions are more general and apply to wave systems with resonances of any order (one would simply have to use different expressions for the integrals $\gamma_{\bf k}(t)$ and $\eta_{\bf k}(t)$ corresponding to the resonance of the considered order; see e.g. [1]). Note that our solution for the PDF (22) and (25) is expressed in terms of the initial PDF and the spectrum $n_{\bf k}(t)$ (recall that $\gamma_{\bf k}(t)$ depends on $n_{\bf k}(t)$ via equation (6)). On the other hand, $n_{\bf k}(t)$ obeys the wave-kinetic equation (7) which is not easy to solve for non-stationary systems. However, it is quite straightforward to solve the wave-kinetic equation numerically, after which the resulting $n_{\bf k}(t)$ can be used in the analytical formula for the PDF's Green function (25). Since the latter formula is very simple, we believe that it can be very effective in practical calculations especially in the situations where non-Gaussianity is important, e.g. in wave weather forecasts including prediction of emerging strong waves—so-called freak waves.

Finally, we would like to mention an additional source of non-Gaussianity which was discussed in [5]. It is related to the situations in which the weak WT description based on the wave-kinetic equation breaks down in finite time. (Recall that in the present paper we are only concerned with the periods of time over which the WT description remains valid.) Physically, this corresponds to formation of strongly nonlinear structures and wave breaking. Often it leads to the formation of a wave spectrum such that the linear and the nonlinear time scales are of the same order over a certain range of wave numbers—the so-called critical balance (CB) state (see e.g. [1]). However, there may be situations where the CB state at large wave numbers may coexist with weak WT at lower wave numbers. In this case the effect of wavebreaking occurring in the CB range onto the weak turbulence range is to produce a finite flux in probability space; see equation (4). This leads to power-law PDF tails; see [5].

Yeontaek Choi's work is supported by the National Institute for Mathematical Sciences (NIMS) and the Korean Union of Public sector Research and Professional workers (KUPRP). The work of Young-Sam Kwon is supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2013R1A1A2057662).