2.1 Kinetic equation

In the field of ocean-wave forecasting one tries to predict a number of key properties of the sea state in a domain which is much larger then the typical wavelength of the sea state in a point of interest. The domain may be as large as the whole globe. Ideally one would like to predict the surface elevation $\unicode[STIX]{x1D702}(\boldsymbol{x},t)$ but the prediction of the phases of individual waves is a hopeless venture, because (i) there are no observations of the initial phases of the waves while (ii) long-time integrations of the Zakharov equation exhibit features of chaotic behaviour (Annenkov & Shrira Reference Annenkov and Shrira2001), i.e. the phases have a sensitive dependence on the initial conditions and are therefore in practice not predictable. This is the main reason that in ocean-wave forecasting one tries to predict the two-dimensional wave variance spectrum $F(\boldsymbol{k},\boldsymbol{x},t)$ which is a function of the wavenumber vector $\boldsymbol{k}$ and a slowly varying function of the spatial coordinates $\boldsymbol{x}$ and time $t$ . The wave spectrum represents the ensemble average of the wave variance as function of wavenumber and the main challenge now is whether it is possible to incorporate aspects of the deviations from the ensemble mean sea state so that statements on the probability of extreme events can be made.

It turns out that this is possible by following the statistical approach to nonlinear interactions which was studied extensively in the 1960s and 1970s and the key developments have been described by Hasselmann (Reference Hasselmann1962) and Davidson (Reference Davidson1972). The application of this approach to water waves for slowly varying spectra is extensively discussed in Janssen (Reference Janssen2003), and in chapter 4 of Janssen (Reference Janssen2004). For a similar approach in nonlinear optics see Suret, Picozzi & Randoux (Reference Suret, Picozzi and Randoux2011). Annenkov & Shrira (Reference Annenkov and Shrira2006) have extended this approach by including effects of rapid variations in the wave spectrum, but we will not need to use their result because we concentrate on the statistical properties of a stationary spectrum.

Because the phase of the waves is not predictable, one can hope at best to predict quantities that remove the phase dependence as much as possible. An example is the second moment $\langle a_{1}a_{2}^{\ast }\rangle$ , where the angle brackets denote an ensemble average. Assuming a zero mean value of the amplitude $a_{1}$ , i.e. $\langle a_{1}\rangle =0$ one can obtain the rate of change in time of the second moment from the Zakharov equation. However, because of nonlinearity, the evolution of the second moment is determined by the fourth moment and so on, resulting in an infinite hierarchy of equations (see e.g. Davidson Reference Davidson1972). A meaningful closure of this hierarchy of equations is obtained by making two assumptions, namely (i) that the ensemble of waves is spatially homogeneous and stationary, hence the second moment satisfies

(2.5) $$\begin{eqnarray}\displaystyle \langle a_{1}a_{2}^{\ast }\rangle =N_{1}\unicode[STIX]{x1D6FF}(\boldsymbol{k}_{1}-\boldsymbol{k}_{2}), & & \displaystyle\end{eqnarray}$$

where $N_{1}$ is the spectral action density, and (ii) the probability distribution for the complex amplitude of the waves $a_{1}$ is close to a Gaussian (random phase approximation), which implicitly implies that the ocean waves are weakly nonlinear. It also has implications for the wave spectrum, which should be sufficiently broad. Using these assumptions one may determine then the rate of change of the action density owing to the four-wave interactions and also the rate of change of the fourth moment, so that one finds an expression for a measure for the deviations from the Gaussian distribution, namely excess kurtosis (cf. Janssen Reference Janssen2003, Reference Janssen2004). Assuming that the action density varies slowly in time the evolution equation for the action density becomes

(2.6) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}N_{1}=4\int \text{d}\boldsymbol{k}_{2,3,4}|T_{1,2,3,4}|^{2}\unicode[STIX]{x1D6FF}_{1+2-3-4}G_{i}(\unicode[STIX]{x0394}\unicode[STIX]{x1D714},t)[N_{3}N_{4}(N_{1}+N_{2})-N_{1}N_{2}(N_{3}+N_{4})],\quad & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D714}_{2}-\unicode[STIX]{x1D714}_{3}-\unicode[STIX]{x1D714}_{4}$ measures the distance to the resonant surface, which is given by the resonance condition $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}=0$ and $\boldsymbol{k}_{1}+\boldsymbol{k}_{2}=\boldsymbol{k}_{3}+\boldsymbol{k}_{4}$ . A key role in this problem is played by the complex resonance function $G$ defined as

(2.7) $$\begin{eqnarray}\displaystyle G(\unicode[STIX]{x0394}\unicode[STIX]{x1D714},t)=\frac{1-\exp (-\text{i}\unicode[STIX]{x0394}\unicode[STIX]{x1D714}t)}{\unicode[STIX]{x0394}\unicode[STIX]{x1D714}} & & \displaystyle\end{eqnarray}$$

because the real part will determine the deviations from Gaussian statistics while the imaginary part determines the rate of change of the action density. Thus, $G_{i}=\text{Im}(G)$ is given by

(2.8) $$\begin{eqnarray}\displaystyle G_{i}(\unicode[STIX]{x0394}\unicode[STIX]{x1D714},t)=\frac{\sin (\unicode[STIX]{x0394}\unicode[STIX]{x1D714}t)}{\unicode[STIX]{x0394}\unicode[STIX]{x1D714}} & & \displaystyle\end{eqnarray}$$

which is time dependent and selects the quartets that contribute to the rate of change of the action spectrum. It is stressed that the action density changes due to non-resonant and resonant transfer. This follows from two limits of $G_{i}$ . For small times we have

(2.9) $$\begin{eqnarray}\displaystyle G_{i}(\unicode[STIX]{x0394}\unicode[STIX]{x1D714},t)\approx t & & \displaystyle\end{eqnarray}$$

while for large times we have

(2.10) $$\begin{eqnarray}\displaystyle \lim _{t\rightarrow \infty }G_{i}(\unicode[STIX]{x0394}\unicode[STIX]{x1D714},t)=\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}(\unicode[STIX]{x0394}\unicode[STIX]{x1D714}). & & \displaystyle\end{eqnarray}$$

Hence, according to (2.6), for short times the evolution of the action density $N$ is caused by both resonant and non-resonant four-wave interactions, while for large times, when the resonance function evolves towards a $\unicode[STIX]{x1D6FF}$ -function, only resonant interactions contribute to spectral change. Note, that in the standard treatment of resonant wave–wave interactions (cf. for example Hasselmann Reference Hasselmann1962 and Davidson Reference Davidson1972) it is argued that the resonance function $G_{i}(\unicode[STIX]{x0394}\unicode[STIX]{x1D714},t)$ may be replaced by its time-asymptotic value (2.10), because the action density spectrum is assumed to be a slowly varying function of time. However, the time required for the resonance function to evolve towards a delta function may be so large that, in the meantime, considerable changes in the statistics of ocean waves have occurred. For this reason the full expression for the resonance function is kept.

At the same time it is emphasized that non-resonant interactions typically occur on a much shorter time scale than the resonant interactions. With $\unicode[STIX]{x1D716}$ a typical wave steepness and $\unicode[STIX]{x1D714}_{0}$ a typical angular frequency of the wave field, one finds from the kinetic equation (2.6) that for non-resonant interactions, corresponding to the short-time limit (2.9), $\unicode[STIX]{x1D70F}_{NL}=O(1/\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D714}_{0})$ , while for resonant interactions $\unicode[STIX]{x1D70F}_{NL}=O(1/\unicode[STIX]{x1D716}^{4}\unicode[STIX]{x1D714}_{0})$ . Therefore, since the non-resonant interactions occur on such a short time scale they determine the transient behaviour of action density and wave statistics.

Connected to this, it is remarked that the random phase approximation has implications for spectral shape. In order for it to be valid the autocorrelation time $\unicode[STIX]{x1D70F}_{ac}=1/\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D714}}$ (here $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D714}}$ is the width of the frequency spectrum) should be much shorter than the nonlinear evolution time scale of the non-resonant interactions $\unicode[STIX]{x1D70F}_{NL}=1/\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D714}_{0}$ . As a consequence, one finds that $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D714}}/\unicode[STIX]{x1D714}_{0}\gg \unicode[STIX]{x1D716}^{2}$ hence wave spectra should be sufficiently broad.

It is also important to note the following difference between the one-dimensional (1-D) and the 2-D kinetic equations. In one dimension only non-resonant interactions give rise to nonlinear transfer (Dyachenko & Zakharov Reference Dyachenko and Zakharov1994). Now, the resonance function $G_{i}$ becomes progressively narrower in the course of time to such an extent that, for large times only, resonant wave–wave interactions are selected. In that event in one dimension there is no change of the wave spectrum possible anymore so that for large times a steady state is achieved. And, according to the numerical simulations and statistical theory also excess kurtosis evolves towards a finite value. In sharp contrast, for the case of two-dimensional propagation, resonant wave–wave interactions do contribute to spectral change, and no steady state solution is achieved. It is therefore of interest to study the time evolution of excess kurtosis in the two-dimensional case.

Finally, it is remarked that the spectral action density $N$ is equivalent to a number density because $\unicode[STIX]{x1D714}N/g$ is the spectral energy density, while $\boldsymbol{k}N/g$ is the spectral momentum density. The kinetic equation (2.6) admits, just as the deterministic Zakharov equation, conservation of total action, wave momentum and energy. In wave forecasting the key parameter of interest is, however, not the spectral action density, but it is the wave variance spectrum $F(\boldsymbol{k})$ (for a thorough discussion of this see Janssen (Reference Janssen2004), which also quotes a number of historical references). It is defined as

(2.11) $$\begin{eqnarray}\displaystyle F(\boldsymbol{k})=\frac{\unicode[STIX]{x1D714}N(\boldsymbol{k})}{g}, & & \displaystyle\end{eqnarray}$$

and it is normalized in such a way that the integral of $F$ over wavenumber space equals the variance $\langle \unicode[STIX]{x1D702}^{2}\rangle$ of the sea surface. Alternatively, one uses the two-dimensional frequency-directional spectrum. Using the wavenumber spectrum $F$ it is defined as

(2.12) $$\begin{eqnarray}\displaystyle E(\unicode[STIX]{x1D714},\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D714}\,\text{d}\unicode[STIX]{x1D703}=F(\boldsymbol{k})\,\text{d}\boldsymbol{k}, & & \displaystyle\end{eqnarray}$$

where the direction $\unicode[STIX]{x1D703}$ is implicitly defined through the wave vector $\boldsymbol{k}:k_{x}=k\cos \unicode[STIX]{x1D703},k_{y}=k\sin \unicode[STIX]{x1D703}$ .

2.2 Deviations from normality

In the past most attention has been devoted to understanding the evolution of the wave spectrum caused by four-wave interactions. However, if there is interest in extreme events it may also be relevant to obtain information regarding the deviations from the Gaussian distribution as caused by the weakly nonlinear four-wave interactions. Because of the symmetries of the Zakharov equation, the first moment of interest is then the fourth moment and the related excess kurtosis $C_{4}$ , defined as

(2.13) $$\begin{eqnarray}\displaystyle C_{4}=\frac{\langle \unicode[STIX]{x1D702}^{4}\rangle }{3\langle \unicode[STIX]{x1D702}^{2}\rangle ^{2}}-1, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D702}$ the surface elevation. It turns out that excess kurtosis is an important measure for the probability of extreme sea states. Since the fourth moment is known it is then straightforward to determine excess kurtosis with the result (Janssen Reference Janssen2003)

(2.14) $$\begin{eqnarray}\displaystyle C_{4}=\frac{4}{g^{2}\langle \unicode[STIX]{x1D702}^{2}\rangle }\int \text{d}\boldsymbol{k}_{1,2,3,4}T_{1,2,3,4}\unicode[STIX]{x1D6FF}_{1+2-3-4}(\unicode[STIX]{x1D714}_{1}\unicode[STIX]{x1D714}_{2}\unicode[STIX]{x1D714}_{3}\unicode[STIX]{x1D714}_{4})^{1/2}\times G_{r}(\unicode[STIX]{x0394}\unicode[STIX]{x1D714},t)N_{1}N_{2}N_{3}, & & \displaystyle\end{eqnarray}$$

with $G_{r}$ the real part of the resonance function $G$ , or

(2.15) $$\begin{eqnarray}\displaystyle G_{r}(\unicode[STIX]{x0394}\unicode[STIX]{x1D714},t)=\frac{1-\cos (\unicode[STIX]{x0394}\unicode[STIX]{x1D714}t)}{\unicode[STIX]{x0394}\unicode[STIX]{x1D714}}. & & \displaystyle\end{eqnarray}$$

Note that for large times $G_{r}\rightarrow P/\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ (with $P$ denoting the principle value) and therefore, unlike the evolution of the action density, the kurtosis does not involve a Dirac $\unicode[STIX]{x1D6FF}$ -function. Therefore the kurtosis is determined by both the resonant and non-resonant interactions.

It is instructive to apply (2.14) to the case of a narrow-band wave spectrum in one dimension, because this allows us to make a connection with the necessary condition for instability, discussed in the previous section (see (2.6)). Hence, performing the usual Taylor expansion around the carrier wavenumber $k_{0}$ to lowest significant order, one finds for large times (Janssen Reference Janssen2003)

(2.16) $$\begin{eqnarray}\displaystyle C_{4}=\frac{8\unicode[STIX]{x1D714}_{0}^{2}}{g^{2}\langle \unicode[STIX]{x1D702}^{2}\rangle }\frac{T_{0}}{\unicode[STIX]{x1D714}_{0}^{\prime \prime }}\int \text{d}p_{1,2,3,4}\frac{\unicode[STIX]{x1D6FF}_{1+2-3-4}}{p_{1}^{2}+p_{2}^{2}-p_{3}^{2}-p_{4}^{2}}N_{1}N_{2}N_{3}, & & \displaystyle\end{eqnarray}$$

where $p=k-k_{0}$ is the wavenumber with respect to the carrier and $\unicode[STIX]{x1D714}_{0}^{\prime \prime }$ is the second derivative of angular frequency with respect to wavenumber. It is seen that the sign of the kurtosis is determined by the ratio $T_{0}/\unicode[STIX]{x1D714}_{0}^{\prime \prime }$ , which is the same parameter that determines whether a wave train is stable or not to sideband perturbations. It should be remarked that, numerically, the integral is found to be negative, at least for bell-shaped spectra. Hence, from (2.16) it is immediately evident that for an unstable wave system which has negative $T_{0}/\unicode[STIX]{x1D714}_{0}^{\prime \prime }$ the kurtosis will be positive and thus will result in an increased probability of extreme events, because of the presence of heavy tails. On the other hand, for a stable wave system a reduction in the probability of extreme events is expected, because excess kurtosis is negative.

These findings have been confirmed by means of numerical simulations of an ensemble of waves of the one-dimensional NLS equation for the focusing case (positive $T_{0}$ ) and the defocussing case (negative $T_{0}$ ), where the initial probability distribution function is a Gaussian (Janssen Reference Janssen2003). Observational and further numerical evidence of the generation of a heavy tail in the probability distribution was found for ocean waves by Onorato et al. (Reference Onorato, Osborne, Serio, Cavaleri, Brandini and Stansberg2004) and in nonlinear optics by Walczak, Randoux & Suret (Reference Walczak, Randoux and Suret2015) and by Suret et al. (Reference Suret, El Koussaifi, Tikan, Evain, Randoux, Szwaj and Bielawski2016). In fact, a quantitative experimental comparison between hydrodynamics and optics which focusses on the emergence of heavy tails has been recently presented by El Koussaifi et al. (Reference El Koussaifi, Tikan, Toffoli, Randoux, Suret and Onorato2018). Evidence for a reduced tail in the probability distribution for the defocusing NLS equation has been found by Randoux et al. (Reference Randoux, Walczak, Onorato and Suret2014).

The initial conditions in these experiments assume a finite bandwidth spectrum, and give, under stationary conditions, rise to significant deviations from normality. Agafontsev & Zakharov (Reference Agafontsev and Zakharov2015) have reached a different conclusion when, as initial condition, a plane wave is chosen, which is a stationary solution of the NLS equation. For this reason they perturbed the plane wave with a small amount of random noise, so that the initial wave spectrum is almost a Dirac-delta function, hence very narrow. It turns out that the stationary statistics for large time become Gaussian. This is quite a surprise because the numerical results show that the system is fairly nonlinear as the nonlinear part of the Hamiltonian is of the same order of magnitude as the linear part. It remains an open question at the moment why there is such an apparent sensitive dependence on the initial condition. Note that in the context of the present statistical approach this question cannot be answered because the autocorrelation scales of the initial condition used by Agafontsev & Zakharov (Reference Agafontsev and Zakharov2015) are very large and most likely violate the condition that the spectrum should be sufficiently broad.

It is also possible to obtain the narrow-band approximation for the case of two-dimensional propagation. This derivation is, however, quite involved and details may be found in Janssen & Bidlot (Reference Janssen and Bidlot2009) which are summarized in Mori, Onorato & Janssen (Reference Mori, Onorato and Janssen2011) and in Fedele (Reference Fedele2015). As already mentioned, wave forecasting is about predicting the angular frequency-directional spectrum. For this reason, the two-dimensional, narrow-band version of the excess kurtosis is expressed in terms of $E$ , defined in (2.12). Excess kurtosis is then given as function of dimensionless short crestedness $R$ and dimensionless time $\unicode[STIX]{x1D70F}$ in terms of a six-dimensional integral

(2.17) $$\begin{eqnarray}\displaystyle C_{4}=J(R,\unicode[STIX]{x1D70F})BFI^{2}, & & \displaystyle\end{eqnarray}$$

where the coefficient $J(R,\unicode[STIX]{x1D70F})$ has a cubic dependence on the wave spectrum $E$ , i.e.

(2.18) $$\begin{eqnarray}\displaystyle J(R,\unicode[STIX]{x1D70F})=2\int \text{d}\unicode[STIX]{x1D708}_{1,2,3}\,\text{d}\unicode[STIX]{x1D719}_{1,2,3}\hat{E}_{1}\hat{E}_{2}\hat{E}_{3}G_{r}(\unicode[STIX]{x0394}\unicode[STIX]{x1D714},\unicode[STIX]{x1D70F}). & & \displaystyle\end{eqnarray}$$

Here, the wave spectrum has been normalized with the wave variance, i.e. $\hat{E}_{1}(\unicode[STIX]{x1D708}_{1},\unicode[STIX]{x1D719}_{1})=E/\langle \unicode[STIX]{x1D702}^{2}\rangle$ .

Note that the narrow-band approximation is basically an expansion of the interaction coefficients and the dispersion relation $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}(k)$ around the peak angular frequency $\unicode[STIX]{x1D714}_{0}=\unicode[STIX]{x1D714}(k_{0})$ and dominant wave direction $\unicode[STIX]{x1D703}_{0}$ of the two-dimensional frequency spectrum, having a frequency width $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D714}}$ and a directional width $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D703}}$ . Here, the dimensionless frequency $\unicode[STIX]{x1D708}$ and the scaled direction $\unicode[STIX]{x1D719}$ are defined through

(2.19a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{0}(1+\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D714}}\unicode[STIX]{x1D708}),\quad \unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{0}+\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D703}}\unicode[STIX]{x1D719}. & & \displaystyle\end{eqnarray}$$

Furthermore,

(2.20) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}=\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D714}}^{2}\unicode[STIX]{x1D714}_{0}t & & \displaystyle\end{eqnarray}$$

is dimensionless time, and

(2.21) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D714}=(\unicode[STIX]{x1D708}_{3}-\unicode[STIX]{x1D708}_{1})(\unicode[STIX]{x1D708}_{3}-\unicode[STIX]{x1D708}_{2})-R(\unicode[STIX]{x1D719}_{3}-\unicode[STIX]{x1D719}_{1})(\unicode[STIX]{x1D719}_{3}-\unicode[STIX]{x1D719}_{2}) & & \displaystyle\end{eqnarray}$$

is the frequency mismatch, where the short crestedness $R$ , defined as

(2.22) $$\begin{eqnarray}\displaystyle R={\textstyle \frac{1}{2}}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D703}}^{2}/\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D714}}^{2}, & & \displaystyle\end{eqnarray}$$

measures the importance of directional aspects of the wave spectrum with respect to the frequency aspects. In particular, for vanishing $R$ , corresponding to the case of infinitely long crests, the envelope wave dynamics is determined by the balance between nonlinearity and the linear dispersion of the waves, while for finite $R$ such a balance does not exist. Finally, $BFI$ is the Benjamin–Feir index given by $BFI=\unicode[STIX]{x1D716}\sqrt{2}/\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D714}}$ with $\unicode[STIX]{x1D716}=k_{0}m_{0}^{1/2}$ the wave steepness, which is the product of the peak wavenumber $k_{0}$ and the variance of the surface elevation $m_{0}$ .

Not much progress has been made with this general expression for excess kurtosis for narrow-band waves. The only general result found thus far assumes that the spectrum has the same form in frequency and direction, i.e. $E(\unicode[STIX]{x1D708},\unicode[STIX]{x1D719})=E(\unicode[STIX]{x1D719},\unicode[STIX]{x1D708})$ . Under this condition it can be shown, by means of interchanging integration variables, that the following relation holds for $C_{4}$ :

(2.23) $$\begin{eqnarray}\displaystyle C_{4}(R,\unicode[STIX]{x1D70F})=-\frac{1}{R}C_{4}\left(\frac{1}{R},R\unicode[STIX]{x1D70F}\right). & & \displaystyle\end{eqnarray}$$

This is a straightforward extension of a result found in Janssen & Bidlot (Reference Janssen and Bidlot2009), who obtained a relation for infinite time $\unicode[STIX]{x1D70F}$ . This is a powerful relation because once one knows $C_{4}$ for $R<1$ , equation (2.23) immediately gives the kurtosis for $R>1$ . Thus, in the remainder of this paper we will restrict our attention to the range $0<R<1$ . From this relation it is also clear that $R=1$ plays a special role. Substituting $R=1$ into (2.23) one finds at once that $C_{4}$ vanishes for all time,

(2.24) $$\begin{eqnarray}\displaystyle C_{4}(1,\unicode[STIX]{x1D70F})=0. & & \displaystyle\end{eqnarray}$$

It turns out that $R=1$ is the demarcation line between nonlinear focusing (for $R<1$ ) and nonlinear defocusing for $R>1$ .

Note that (2.17) and (2.18) give the evolution in time of the excess kurtosis subject to the initial condition that kurtosis vanishes, i.e. that initially the probability density function of surface elevation is normally distributed. Therefore, although initially the ocean waves follow the normal law, the nonlinear interactions force the system out of the normal state, and the question is of interest how far the system is forced out of the normal state and how long this condition lasts. In addition, it will turn out that the answer to this question depends in a sensitive manner on the role of the directional distribution, i.e. on the value of the parameter $R$ . For finite $R$ it will be found that for large times the probability distribution evolves towards a Gaussian. The exception is the limit of vanishing $R$ , corresponding to the case of the one-dimensional NLS equation. This equation is integrable, because it can be solved by means of the inverse scattering transform, its kinetic equation only allows non-resonant interactions and for large times excess kurtosis is found to evolve towards a non-zero value.

The first part of this question, regarding how far the system is forced out of the normal state, has been extensively studied in the laboratory by Onorato et al. (Reference Onorato, Waseda, Toffoli, Cavaleri, Gramstad, Janssen, Kinoshita, Monbaliu, Mori, Osborne, Serio, Stansberg, Tamura and Trulsen2009) who recorded the maximum value of excess kurtosis as a function of the directional width of the spectrum. Additional information on the dependence of excess kurtosis on the directional width follows from Mori et al. (Reference Mori, Onorato and Janssen2011) who performed numerical simulations for an ensemble of waves using the two-dimensional NLS equation. Time evolution and the maximum of excess kurtosis were studied in some detail.

In order to further analyse this complicated problem, consider the special case of a Gaussian spectrum, which is symmetrical with respect to the origin,

(2.25) $$\begin{eqnarray}\displaystyle \hat{E}_{1}=\frac{1}{2\unicode[STIX]{x03C0}}\text{e}^{-(1/2)(\unicode[STIX]{x1D708}_{1}^{2}+\unicode[STIX]{x1D719}_{1}^{2})}, & & \displaystyle\end{eqnarray}$$

and $J(R,\unicode[STIX]{x1D70F})$ becomes

(2.26) $$\begin{eqnarray}\displaystyle J(R,\unicode[STIX]{x1D70F})=\frac{2}{(2\unicode[STIX]{x03C0})^{3}}\int _{-\infty }^{\infty }\text{d}\unicode[STIX]{x1D708}_{1,2,3}\,\text{d}\unicode[STIX]{x1D719}_{1,2,3}\text{e}^{-(1/2)(\unicode[STIX]{x1D708}_{1}^{2}+\unicode[STIX]{x1D719}_{1}^{2}+\unicode[STIX]{x1D708}_{2}^{2}+\unicode[STIX]{x1D719}_{2}^{2}+\unicode[STIX]{x1D708}_{3}^{2}+\unicode[STIX]{x1D719}_{3}^{2})}G_{r}(\unicode[STIX]{x0394}\unicode[STIX]{x1D714},\unicode[STIX]{x1D70F}). & & \displaystyle\end{eqnarray}$$

The time behaviour of excess kurtosis is now determined by (2.17) and (2.26). Numerical evaluation of the integral in (2.26) is possible for relatively short times, but convergence problems, caused by the rapidly varying resonance function $G_{r}$ , are evident for large times. Therefore, an analytical approach was clearly desirable. The first author realized that, based on the work of Fedele et al. (Reference Fedele, Cherneva, Tayfun and Guedes Soares2010), this type of integral may perhaps be evaluated for arbitrary time, and he requested Francesco Fedele to investigate this possibility. This was indeed feasible (see Fedele Reference Fedele2015), and he was able to express the evolution in time of the excess kurtosis in terms of a complex one-dimensional integral. The main result found was that for $0<R<1$ the function $J$ will reach a maximum at $\unicode[STIX]{x1D70F}_{max}=1/\sqrt{3R}$ after which it decays for large times to zero. For $R>1$ the function will reach a minimum before it decays to zero.

In this paper, we show by analytical means some further interesting properties of the analytical result. However, first, a concise derivation of Fedele’s elegant integral expression for the evolution in time of the excess kurtosis is given.

4.1 Proof that for $R>0$ the integral $J$ vanishes for $\unicode[STIX]{x1D70F}\rightarrow \infty$

It is important to prove that for large $\unicode[STIX]{x1D70F}$ the integral $J$ vanishes, $J(R,\unicode[STIX]{x1D70F}=\infty )=0$ . At least two proofs of this are known to us. The argument used in Fedele (Reference Fedele2015), p. 28, to prove this gives that $\lim _{\unicode[STIX]{x1D70F}\infty }J(R,\unicode[STIX]{x1D70F})$ indeed exists, but it is not obvious that this limit indeed equals 0; also see the discussion on this matter below.

The first proof starts from (3.12) which is regarded for now as an integral along the real axis. We have seen that the integral $J$ has a maximum at $\unicode[STIX]{x1D70F}=1/\sqrt{3R}$ and we regard this as a special point along the integration path. Therefore, we write

(4.1) $$\begin{eqnarray}\displaystyle J(R,\unicode[STIX]{x1D70F}=\infty )={\textstyle \frac{2}{3}}\text{Im}[I_{1}+I_{2}], & & \displaystyle\end{eqnarray}$$

where

(4.2a,b ) $$\begin{eqnarray}\displaystyle I_{1}=\int _{0}^{1/\sqrt{3R}}\frac{\text{d}z}{N(z)},\quad I_{2}=\int _{1/\sqrt{3R}}^{\infty }\frac{\text{d}z}{N(z)}, & & \displaystyle\end{eqnarray}$$

with $N(z)=(z-\text{i})^{1/2}(z+(1/3)\text{i})^{1/2}(Rz+\text{i})^{1/2}(Rz-(1/3)\text{i})^{1/2}$ . The next step is to introduce a transformation of variables in the second integral in such a way that the integration range becomes identical to the integration range of the first integral. The transformation

(4.3) $$\begin{eqnarray}\displaystyle z=\frac{1}{3Rw}, & & \displaystyle\end{eqnarray}$$

with Jacobian $=-1/3Rw^{2}$ , hence $\text{d}z=-\text{d}w/3Rw^{2}$ , indeed changes the integration range from $[1/\sqrt{3R},\infty ]$ to $[0,1/\sqrt{3R}]$ . Now, the important point to note is that under this transformation the function $N(z)$ is invariant in the sense that, apart from the Jacobian of the transformation, it equals the complex conjugate of $N(w)$ , i.e.

(4.4) $$\begin{eqnarray}\displaystyle N(z)=\frac{1}{3Rw^{2}}N^{\ast }(w). & & \displaystyle\end{eqnarray}$$

Therefore, using the transformation (4.3) in the integral $I_{2}$ , we find that $I_{2}=I_{1}^{\ast }$ so that

(4.5) $$\begin{eqnarray}\displaystyle J(R,\unicode[STIX]{x1D70F}=\infty )={\textstyle \frac{2}{3}}\text{Im}[I_{1}+I_{1}^{\ast }]=0, & & \displaystyle\end{eqnarray}$$

i.e. the integral $J$ vanishes for large times. The restriction to this result is that it only holds for finite $R$ since for vanishing $R$ there is no maximum of $J$ for finite time.

It is emphasized that the vanishing of the integral $J$ , and hence kurtosis for large time, is not a property for arbitrary spectral shape. We think that this property is closely connected to the ‘invariance’ property (4.4) of the denominator $N(z)$ which in the end is determined by our choice of a Gaussian spectrum (2.25). Numerical computations for spectra that are not symmetrical with respect the peak frequency do give a non-zero value of kurtosis for large time.

The second proof is rather involved, but in the end it utilizes the same invariance property for $N(z)$ to show the vanishing of $J$ . This proof regards the integral in (3.12) as an integral in the complex domain, and one applies the theorem that if there are no singularities inside the domain of a closed contour then the resulting integral vanishes. The contour was chosen to be (a) $0<x<\infty$ , (b) the circle with large radius $\unicode[STIX]{x1D70C}$ and phase $\unicode[STIX]{x1D703}$ with $0<\unicode[STIX]{x1D703}<\unicode[STIX]{x03C0}/2$ , and then back to the origin from $y=+\text{i}\infty$ to $y=0$ . Now, because of the branch points in (3.12), branch cuts need to be introduced, and basically this is done along the imaginary axis from the largest positive branch point to the smallest negative branch point. The phases of all relevant functions are then well defined. In addition, for $R>0$ one can show that the integral along the circle with large radius vanishes as $\unicode[STIX]{x1D70C}\rightarrow \infty$ , and therefore the integral along the real axis can be expressed in terms of the integral along the imaginary axis. The advantage is, namely, that the integral along the imaginary axis is much simpler. In fact, one finds for $R<1/3$

(4.6) $$\begin{eqnarray}\displaystyle J(R,\infty )=\frac{1}{3}\left\{\int _{0}^{1}\frac{\text{d}y}{N}-\int _{1/3R}^{\infty }\frac{\text{d}y}{N}\right\}, & & \displaystyle\end{eqnarray}$$

where $N=\{|y-1||y+1/3||Ry+1||Ry-1/3|\}^{1/2}$ . For $R>1/3$ a similar result is found. In the final step it can be shown that the two integrals cancel each other, therefore $J(\unicode[STIX]{x1D70F}\rightarrow \infty )$ vanishes. The easiest way to show this is by means of a transformation of the integration variable in the second integral, i.e. $y=1/3Rx$ , which is identical to the transformation (4.3)! For completeness details of this calculation are given in appendix A.

The contour integration approach, as employed in the second proof, can be used to establish further analytic results. This is done in all detail in appendix C, where the behaviour of $\max _{\unicode[STIX]{x1D70F}}J(R,\unicode[STIX]{x1D70F})$ near $R=0$ is obtained. The proof technique in appendix C can, furthermore, be used to obtain an expression of $J(R,1/\sqrt{3R})$ in terms of the incomplete elliptic integral of the first kind.

Figure 3.

Comparison of exact normalized integral $J/N_{J}$ with the asymptotic expansion (4.8) for $R=0.5$ , the red line showing the first term and the green line showing the sum of the two terms in (4.8).

4.2 Decay of $J(R,\unicode[STIX]{x1D70F})$ as $\unicode[STIX]{x1D70F}\rightarrow \infty$

Having established that the integral $J$ of (3.12) vanishes for large dimensionless time $\unicode[STIX]{x1D70F}$ , we now study how $J$ approaches this limit. Starting from (3.12) we have

(4.7) $$\begin{eqnarray}\displaystyle J(R,\unicode[STIX]{x1D70F})=\frac{2}{3}\text{Im}\int _{0}^{\unicode[STIX]{x1D70F}}\frac{\text{d}z}{N(z)}=-\frac{2}{3}\text{Im}\int _{\unicode[STIX]{x1D70F}}^{\infty }\frac{\text{d}z}{N(z)}, & & \displaystyle\end{eqnarray}$$

where again $N(z)=(z-\text{i})^{1/2}(z+(1/3)\text{i})^{1/2}(Rz+\text{i})^{1/2}(Rz-(1/3)\text{i})^{1/2}$ . The large time behaviour of $J$ is now found by expanding the integrand of the second integral of (4.7) into powers of $1/z$ and integrating the expansion. The details of this calculation are given in appendix B, and the final result is

(4.8) $$\begin{eqnarray}\displaystyle J(R,\unicode[STIX]{x1D70F})=\frac{1}{R}\left(\frac{1}{R}-1\right)\left[\frac{1}{9\unicode[STIX]{x1D70F}^{2}}-\left(\frac{7}{R^{2}}+\frac{4}{R}+7\right)\frac{1}{162\unicode[STIX]{x1D70F}^{4}}+\cdots \!\right], & & \displaystyle\end{eqnarray}$$

an expansion which is valid for $\unicode[STIX]{x1D70F}>1/R$ , as is verified in figure 3, which compares the exact result for $R=0.5$ with the above asymptotic expansion. The first term in (4.8) agrees with what can be obtained from Fedele (Reference Fedele2015), (3.5). The second-order term has a surprisingly large impact for moderate values of dimensionless time $\unicode[STIX]{x1D70F}$ .

Figure 4.

Comparison of numerical maximum value of normalized integral $J/N_{J}$ with the first- and second-order approximation (4.10). Also, the expansion (4.14) is shown.

4.3 Behaviour of $\max _{\unicode[STIX]{x1D70F}}J(R,\unicode[STIX]{x1D70F})$ as $R\downarrow 0$

It is important to have an estimate of the severity of the sea state during extreme events. Such an indicator may be provided by the maximum of the excess kurtosis as a function of time. It is therefore of interest to get an estimate of the maximum of the integral $J$ , which occurs at time $\unicode[STIX]{x1D70F}=1/\sqrt{3R}$ . For $0<R<1$ we have

(4.9) $$\begin{eqnarray}\displaystyle \max _{\unicode[STIX]{x1D70F}\geqslant 0}J(R,\unicode[STIX]{x1D70F})=J\left(R,\frac{1}{\sqrt{3R}}\right), & & \displaystyle\end{eqnarray}$$

where $J(R,\unicode[STIX]{x1D70F})$ is given by the last expression in (4.7). The detailed analysis is given in appendix C. Taylor expansion around $R=0$ then gives the eventual result

(4.10) $$\begin{eqnarray}\displaystyle J\left(R,\frac{1}{\sqrt{3R}}\right)=\frac{\unicode[STIX]{x03C0}}{3\sqrt{3}}\left(1-\frac{4}{\unicode[STIX]{x03C0}}\sqrt{3R}+\left(\frac{2\sqrt{3}}{\unicode[STIX]{x03C0}}+\frac{1}{3}\right)R+O(R^{3/2})\right). & & \displaystyle\end{eqnarray}$$

The first term of this Taylor expansion agrees with what can be found from Fedele (Reference Fedele2015), (4.3), (4.5). In order to assess the accuracy of the approximation (4.10), we compare in figure 4 the values of the first-order approximation and the second-order approximation with the numerically computed value of $J(R,1/\sqrt{3R})$ , normalized with $\unicode[STIX]{x03C0}/3\sqrt{3}$ . Although the first-order approximation deviates quite quickly from the ‘exact’ result, it is rewarding to see how much the second-order term improves the agreement with ‘exact’ up to values of $R=0.4$ . This is surprising as the expansion in (4.10) has a radius of convergence of only $1/3$ .

Here it is mentioned that once more, we use complex analysis and Cauchy’s theorem to simplify the integral. The contour integration approach can be used to classify the integral $J(R,1/\sqrt{3R})$ in terms of the incomplete elliptic integral of the first kind that is available in many mathematical packages. This result is also instrumental for checking numerical work. It can then be shown that

(4.11) $$\begin{eqnarray}\displaystyle J\left(R,\frac{1}{\sqrt{3R}}\right)=\left\{\begin{array}{@{}l@{}}\displaystyle \frac{k}{\sqrt{R}}(F(x_{1},k)-F(x_{2},k)),\quad 0\leqslant R<\frac{1}{3},\\ \displaystyle \frac{1}{\sqrt{R}}(F(x_{1}^{-1},k^{-1})-F(x_{2}^{-1},k^{-1})),\quad \frac{1}{3}<R\leqslant 1,\end{array}\right. & & \displaystyle\end{eqnarray}$$

where

(4.12a-c ) $$\begin{eqnarray}\displaystyle k=\frac{4\sqrt{3R}}{3(1+R)},\quad x_{1}=\frac{1}{2}\sqrt{3(1+R)},\quad x_{2}=\frac{(3+\sqrt{3R})\sqrt{R+\frac{1}{3}}-(1-\sqrt{3R})\sqrt{R+3}}{\sqrt{32R}}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

and, for $0\leqslant x<1$ , $0\leqslant k\leqslant 1$ ,

(4.13) $$\begin{eqnarray}\displaystyle F(x,k)=\int _{0}^{x}\frac{\text{d}t}{\sqrt{(1-t^{2})(1-k^{2}t^{2})}}, & & \displaystyle\end{eqnarray}$$

the incomplete elliptic integral of the first kind. The proof of this result is quite involved, and we omit it. Using the Keisan Online Calculator we found perfect agreement with the numerical results for $J(R,1/\sqrt{3R})$ in table 1 of appendix C.

4.4 Behaviour of $\max _{\unicode[STIX]{x1D70F}}J(R,\unicode[STIX]{x1D70F})$ as $R\rightarrow 1$

In appendix D the behaviour of $\max _{\unicode[STIX]{x1D70F}}J(R,\unicode[STIX]{x1D70F})=J(R,1/\sqrt{3R})$ is obtained as $R\rightarrow 1$ . The result is

(4.14) $$\begin{eqnarray}\displaystyle J(R,1/\sqrt{3R})=({\textstyle \frac{1}{2}}-{\textstyle \frac{1}{4}}\sqrt{3})\{(1-R)+(1-R)^{2}\}+O((1-R)^{3}). & & \displaystyle\end{eqnarray}$$

This result is derived in a straightforward manner from the expression (3.12) for $J(R,\unicode[STIX]{x1D70F}=1/\sqrt{3R})$ by differentiation. Observe that the coefficients of $(1-R)$ and $(1-R)^{2}$ in (4.14) are identical. The expansion (4.14), without the term $O((1-R)^{3})$ , is shown in figure 4, and shows good agreement with the numerical result in the range $0.6<R<1$ .

4.5 Approximating $\max _{\unicode[STIX]{x1D70F}}J(R,\unicode[STIX]{x1D70F})$ on the range $0<R<1$

We wonder whether it would be possible to use the approximation (4.10) and (4.14) in a compromise fit so as to get a reasonable approximation over the whole range $0<R<1$ . Due to the highly different nature, as a function of $R$ , of $J(R,1/\sqrt{3R})$ near $R=0$ and $R=1$ , a simple convex combination of (4.10) and (4.14) fails to give such an approximation. Instead we propose the approximation

(4.15) $$\begin{eqnarray}\displaystyle J\left(R,\frac{1}{\sqrt{3R}}\right)\approx \frac{\unicode[STIX]{x03C0}}{3\sqrt{3}}(1-\unicode[STIX]{x1D6FC}\sqrt{R}+\unicode[STIX]{x1D6FD}R+\unicode[STIX]{x1D6FF}R^{2}),\quad 0<R<1, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}=4\sqrt{3}/\unicode[STIX]{x03C0}$ and $\unicode[STIX]{x1D6FD}=(1/3+2\sqrt{3}/\unicode[STIX]{x03C0})$ , see (4.10), and $\unicode[STIX]{x1D6FF}$ to be chosen such that both sides of (4.15) agree at $R=1$ . This gives the condition $1-\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FD}+\unicode[STIX]{x1D6FF}=0$ from which $\unicode[STIX]{x1D6FF}=2\sqrt{3}/\unicode[STIX]{x03C0}-4/3=-0.2307$ . The compromise fit (4.15) is compared with the simple fit (3.16) and with the exact result obtained numerically in figure 2. The compromise fit achieves a maximum absolute error in (4.15) of the order $10^{-3}$ , and thus outperforms fit (3.16) by far. Furthermore, the values of the derivatives of both sides in (4.15) at $R=1$ are $-0.0670$ and $-0.0774$ , respectively.

We also tried a compromise fit of the form

(4.16) $$\begin{eqnarray}\displaystyle J\left(R,\frac{1}{\sqrt{3R}}\right)\approx \frac{\unicode[STIX]{x03C0}}{3\sqrt{3}}(1-\unicode[STIX]{x1D6FC}\sqrt{R}+\unicode[STIX]{x1D6FD}R++\unicode[STIX]{x1D6FE}R^{3/2}+\unicode[STIX]{x1D6FF}R^{2}),\quad 0<R<1, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ as in (4.15) and $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D6FF}$ such that (i) both sides of (4.16) and their derivatives agree at $R=1$ , or, (ii) both sides of (4.16) agree at $R=1$ and $R=1/3$ . It was found that these fits do not improve the fit (4.15), with fit (i) systematically too small and fit (ii) systematically too large on a large part of the range $0<R<1$ . The average of the fits (i) and (ii) performs about as good as fit (4.15), with parameter values $\unicode[STIX]{x1D6FE}=-0.0070$ , $\unicode[STIX]{x1D6FF}=-0.2237$ that are remarkably close to the parameter values $\unicode[STIX]{x1D6FE}=0.0000,\unicode[STIX]{x1D6FF}=-0.2307$ pertaining to fit (4.15).