# Modulation instability and convergence of the random phase approximation for
stochastic sea states

Agissilaos G. Athanassoulis Corresponding author email:
<EMAIL_ADDRESS>Department of Mathematics, University of Dundee
Irene Kyza  Department of Mathematics, University of Dundee

###### Abstract

The nonlinear Schrödinger equation is widely used as an approximate model for
the evolution in time of the water wave envelope. In the context of simulating
ocean waves, initial conditions are typically generated from a measured power
spectrum using the random phase approximation, and periodized on an interval
of length $L$. It is known that most realistic ocean waves power spectra do
not exhibit modulation instability, but the most severe ones do; it is thus a
natural question to ask whether the periodized random phase approximation has
the correct stability properties. In this work we specify a random phase
approximation scaling so that, in the limit of $L\to\infty,$ the stability
properties of the periodized problem are identical to those of the continuous
power spectrum on the infinite line. Moreover, it is seen through concrete
examples that using a too short computational domain can completely suppress
the modulation instability.

Keywords:  Stochastic sea state, random phase approximation, Nonlinear
Schrödinger equation, Modulation instability, Alber equation

##  1 Introduction

Ocean wave models can be classified as phase-resolved or phase-averaged. A
phase-resolved model aims to describe each individual wave crest and trough –
for example the first principles, fully nonlinear, free-boundary problem for
potential flow, or asymptotic approximations thereof (Boussinesq, KdV,
Zakharov, nonlinear Schrödinger etc) are all phase-resolved models. On the
other hand, phase-averaged models (such the CSY, Alber and Hasselmann
equations) aim to describe a stochastic sea state on a macroscopic and/or
coarse-grained level. In that sense any phase-averaged model is a statistical
model of waves, and typically uses some additional assumptions in order to
come up with closed equations for a phase-averaged representation. Phase-
resolved models can be very accurate in terms of hydrodynamics. However, they
still have to be initialized, and any prediction with them ultimately is as
good as the initial conditions used. Since the vast majority of data available
for real-world ocean waves is phase-averaged (typically power spectra /
autocorrelations), the crisp exactness of ideal hydrodynamics is out of reach
anyway when modelling the ocean. On the other hand, phase-averaged models have
the advantage of being designed with real-world data in mind. They are well
suited for the large space and time scales required in ocean engineering, and
are known to provide valuable insights for the qualitative behaviour of
stochastic wavefields. Phase-resolved simulations can also be used in
conjunction with data to study stochastic sea states in a Monte Carlo sense,
i.e. when used with a large population of initial conditions (realizations of
the sea state) generated from phase-averaged measurements. This is more
accurate, but also much more computationally expensive than phase-averaged
models – and thus reserved for validation or cases of special interest.
Typically, a measured power spectrum (itself a phase-averaged representation
of the sea surface elevation) is used to generate a population of initial
conditions on a finite computational domain of size $L$ through the random-
phase approximation.

In this work we study the convergence for the random phase approximation in
terms of qualitative behaviour, i.e. whether it correctly captures the
presence or not of modulation instability in the original problem. We work in
a regime where it is known that both modulation instability and its absence
are possible [1, 2, 3, 24, 31, 51].

###  1.1 Phase-averaged models and the modulation instability

In 1962, Nobel-prize winner K. Hasselmann published what is now known as the
Hasselmann equation, a phase-averaged equation describing the energy transfer
between various wavenumbers [35]. The starting point of [35] is a stochastic
sea state, where both the background and perturbations are gaussian and
homogeneous random processes. Both gaussianity and homogeneity are widely
considered to be fundamental features of deep-water gravity waves [46]; this
is largely consistent with empirical observation, and is furthermore supported
intuitively by an informal application of the Central Limit Theorem. Today
many state of the art approaches still use the Hasselmann equation, or similar
moment equations (often collectively called “wave kinetic equations”, WKE) to
study different aspects of real-world ocean waves, including long-distance
propagation [36, 52, 37], extreme events [26, 28, 27] etc. Moreover, it has
inspired a number of rigorous mathematical studies not necessarily closely
tied to real-world ocean waves. These include the well-posedness theory of
generalized wave kinetic equations [30], wave turbulence [15, 53, 16], as well
as the rigorous investigation of the closure conditions involved in the
derivation of the WKE [23, 22, 21].

A few years after the introduction of Hasselmann’s equation, what is now
called modulation instability (MI) entered the picture, in the study of the
nonlinear Schrödinger equation (NLS) as an approximate model for ocean waves
[56, 7] (groups working on different problems were also encountering aspects
of the MI around the same time, see [57] for a historical overview). The
presence of MI in focusing dispersive waves means that plane waves are
inherently unstable. It is now known that the MI is a pretty universal feature
of focusing, dispersive waves [14, 12]. Moreover, [14, 12] point out that the
MI can be seen as the tendency of inhomogeneities to initially grow, and
subsequently settle in a particular pattern. It is natural to ask whether
stochastic wavefields with sufficiently narrow power spectra behave like plane
waves in that respect. The fact that ocean waves are in many cases narrowband
makes this a very natural question.

A first attempt to investigate how the MI might show up in a phase-averaged
moment equation (analogous to Hasselmann’s equation) was made by Longuet-
Higgins [42]. In that work Davey-Stewartson dynamics were used and it was
found that the MI for a plane wave would indeed shows up in the moment
equations. The investigation for a wavefield with a narrow but realistic power
spectrum (not a delta-function) was first performed in [1]. A sufficient
condition for exponential growth of small inhomogeneities was derived, in the
form of a nonlinear system of equations depending on the power spectrum (what
I. E. Alber called “the integral stability eigenvalue relation” [1]). We will
call this kind of instability of a stochastic sea state generalized modulation
instability (gMI), as opposed to the classical MI (MI) which is the growth of
inhomogeneities around a plane wave only. It turns out that the MI can be seen
as a special case of the gMI for a delta-function power spectrum (which
practically amounts to a plane wave) [4]. The notable feature of Alber’s
analysis is that sufficiently narrow and intense spectra can be seen to be
unstable (meaning that inhomogeneities would spontaneously grow), while
sufficiently broad and low-intensity spectra are expected to be stable
(meaning that inhomogeneities do not grow; in fact later it was proved that
inhomogeneities disperse in the stable case [3]). The quantitative resolution
of the stability condition proved to require the development of novel
numerical and analytical techniques, and took several decades to mature [2, 3,
51, 31] – more on that below.

The possibility of realistic sea states for which gMI is present (i.e.
homogeneity is not stable) has profound consequences in the modelling of ocean
waves. This was first articulated in the seminal paper [17], where a second
moment equation is derived for water waves with an assumption of gaussianity
only – not homogeneity. This is called the CSY equation, and it is a
broadband, 2-dimensional two-space moment equation, derived from the Zakharov
equation using only a gaussian closure. It is the most accurate moment
equation for water waves in the literature. Assuming homogeneity is a robust
feature of the sea state, then one can go to sixth order effects for the
homogeneous wavefield, and essentially reproduce Hasselmann’s equation [17].
When a narrowband approximation is used in [17], it is found that the question
of whether homogeneity is preserved can be investigated quantitatively, and
Alber’s condition is essentially recovered. Crucially, if homogeneity is not
stable, then the effects of the gMI are second order. In [2] this point was
further strengthened: assuming that the initial sea state in the CSY equation
consists of an $O(1)$ homogeneous part and an $o(1)$ inhomogeneous part, a CSY
stability condition was derived, controlling whether the inhomogeneous part
grows exponentially or not. This was done without any narrowband
approximation, and it yields the broadband, 2-dimensional counterpart of
Alber’s “stability relation”.

So [1, 17, 2] collectively present a fundamental dichotomy: starting with the
a priori assumptions of gaussianity and quasi-homogeneity of a given sea state
(meaning the initial presence of a small inhomogeneity), an Alber-type
stability condition for that particular sea state can be produced. If it is
found to be stable, then homogeneity is robust, the inhomogeneity is going to
disperse, and one can safely go to higher-order effects in order to derive
homogeneous kinetic wave equations - as in, e.g., [35, 38]. In that sense the
Alber stability of the sea state is a compatibility condition, allowing the
self-consistent use of homogeneity.

If the spectrum is found to be unstable however, then any inhomogeneity (which
in the field could be due to wave breaking, gusts of wind, localized extreme
events, presence of ships etc) will be spontaneously amplified, and this is a
leading order effect. Thus homogeneity cannot be consistently used as a
permanent structural feature of an unstable sea state: the homogeneous sea
state is basically an unstable equilibrium. These sea states where gMI is
present have to be studied with different methods, and are expected to behave
very differently from the stable ones. This pumping of energy in
inhomogeneities for unstable spectra has been reproduced in many works,
including [24, 25, 33, 49, 44, 6, 48, 51]. In the presence of gMI, large
localized extreme events emerge which clearly break homogeneity. It is worth
noting that working out the fundamental scalings of the gMI can give a
qualitative description of these localized extreme events; see Figure 1 for a
comparison of asymptotic analysis of the Alber equation [4] with state of the
art Monte Carlo results [20].

The remaining question is the quantitative resolution of the stability
condition: when is a given power spectrum expected to be stable or unstable?
The papers [2, 3, 51, 31] all investigate the onset of gMI for JONSWAP
parametric spectra. The relevant parametrization here is in $\alpha,$ a
parameter directly affecting the severity of the sea state, and $\gamma,$ a
parameter primarily affecting how narrow the spectrum is. In practice, more
severe sea states (storms) also tend to have narrower spectra – i.e. both
$\alpha$ and $\gamma$ tend to increase or decrease together. While different
techniques (and in some cases different equations) are used, the findings are
broadly in agreement: the vast majority of realistic ocean wave spectra are
stable, but the most extreme sea states on record cross over to the unstable
region. On one hand, this is consistent with the widespread success of
homogeneous approaches and the absence of MI-like effects in mild and moderate
seas [27, 41]. On the other hand, it is consistent with large, dangerous rare
events requiring a separate investigation [47, 46, 39]. It is worth noting
that investigation of non-parametric spectra [5] is also possible, and yields
similar results.

(a)

(b)

Figure 1: Left: Inhomogeneity produced by the unstable modes of the Alber
equation for a JONSWAP spectrum [4]. Right: Average profile of extreme events
for a JONSWAP spectrum; Monte Carlo simulation using modified NLS dynamics
(image provided by T. Grafke) [20].

###  1.2 Onset of instability and Monte Carlo setup

Understanding these unstable sea states is crucial for several reasons: on one
hand, it is likely that they play a key role in real-world oceanic rogue
waves; on the other hand they represent a mathematical regime that is very
poorly understood and in fact most of the time we simply assume them away.
Because of the instability of homogeneity, any phase-averaged model that
assumes permanent homogeneity cannot be used as a first-principles model. For
that we should use Monte Carlo simulations of the underlying phase-resolved
problem.

That is, we have to use phase-averaged data (power-spectra) to create a large
number of realizations of the sea surface elevation. These initial
realizations will be evolved in time over a finite computational (or
experimental) domain of length $L$. Very often periodic boundary conditions
are used (although different kinds of boundary conditions also exist in the
literature, e.g. non-reflective absorbing etc). It is standard practice to do
this with the random-phase approximation [46, 40], basically a superposition
of plane waves with random phases and amplitudes consistent with the power
spectrum.

Let us assume nonlinear Schrödinger (NLS) dynamics governing the evolution in
time of a complex wave envelope $u(x,t).$ The original physical problem is
typically framed in infinite space [43],

$i\partial_{t}U+\frac{p}{2}\Delta U+\frac{q}{2}|U|^{2}U=0,\quad
x\in\mathbb{R},\quad\mathop{\sup}\limits_{x\in\mathbb{R}}|U(x,t)|<\infty,$ (1)

with initial data $U_{0}(x)$ being a realization of a homogeneous random
process with power spectrum $S(k).$ On the other hand, the truncated problem
that we can readily simulate is the periodized NLS

$\begin{array}[]{c}i\partial_{t}u+\frac{p}{2}\Delta
u+\frac{q}{2}|u|^{2}u=0,\qquad x\in[-\frac{L}{2},\frac{L}{2}],\\\
u(-\frac{L}{2},t)=u(\frac{L}{2},t),\qquad\partial_{x}u(-\frac{L}{2},t)=\partial_{x}u(\frac{L}{2},t),\end{array}$
(2)

with initial data of the form $u_{0}(x)=\sum A_{j}e^{2\pi i(k_{j}x+\phi_{j})}$
where the $\phi_{j}$ are i.i.d. random variables uniformly distributed in
$[0,1).$

The main question that we address in this paper is the following:

  * •

Under what conditions can we guarantee that the presence or not of gMI in the
original problem (1) is accurately captured in the periodized problem (2)?

We answer this question by deriving a periodized Alber equation for the
truncated problem, along with its stability relation which contains explicitly
the lengthscale $L.$ Then we can see that the instability condition for the
truncated problem converges to that of the original problem if appropriate
scalings are used (roughly speaking as long as $L$ is large enough and $O(L)$
plane waves are used – see Theorem 2.1). More specifically, we find a way to
create realizations of the sea state that leads to the correct stability
condition. To highlight the role of $L$ in the qualitative behaviour of the
solution, we investigate its role in the fully nonlinear MI in Section 5.

Recent theoretical breakthroughs suggest that different preparation of the
initial phase-resolved sea state may lead to completely different regimes [23,
16]. Thus, the preparation of the initial data deserves more detailed
attention. This is already underway on the engineering side, where
initialization and calibration is now starting to be seen as an iterative
workflow in the context of high quality data [34] rather than a standard
formula.

##  2 Main results

######  Theorem 2.1 (Periodized Alber equation).

Let $S(k)$ be a smooth power spectrum with
$\mathop{supp}S(k)\subseteq[0,k_{max}],$ and consider equation (2) with
stochastic initial data $u(x,0)=u_{0}(x)$ where

$u_{0}(x)=\sum\limits_{j=1}^{M}A_{j}e^{2\pi i(\frac{jx}{L}+\phi_{j})},\qquad
A_{j}=\sqrt{\frac{1}{L}S(\frac{j}{L})},\qquad\phi_{j}\sim\mbox{ i.i.d.
}U(0,1),\qquad M=k_{max}\cdot L.$ (3)

The Alber equation for the evolution in time of inhomogeneity $\rho(x,y,t)$ of
this problem is

$\begin{array}[]{c}i\partial_{t}\rho+\frac{p}{2}\left(\Delta_{x}-\Delta_{y}\right)\rho+q\left[\Gamma(x-y)+\rho(x,y,t)\right]\left[\rho(x,x,t)-\rho(y,y,t)\right]=0,\quad
x,y\in[\frac{L}{2},\frac{L}{2}],\end{array}$

equipped with periodic boundary conditions, where

$\Gamma(x)=\sum\limits_{n}\frac{1}{L}S(\frac{n}{L})e^{2\pi i\frac{xn}{L}},$

and the linear instability condition for this problem is

$\exists\xi\in\mathbb{Z}\qquad\exists\omega\in\\{z\in\mathbb{C}:\operatorname{Re}(z)>0\\}\qquad\qquad\widetilde{h}_{L}(\xi,\omega)=1$
(4)

where

$\widetilde{h}_{L}(\xi,\omega)=\frac{iq}{L}\sum\limits_{k\in\mathbb{Z}}\frac{S(\frac{k-\xi}{L})-S(\frac{k}{L})}{\omega+ip\frac{2\pi^{2}\xi(2k-\xi)}{L^{2}}}.$
(5)

Moreover, in the limit $L\to\infty,$ $\xi\to\infty,$ $\xi/L\to X$

$\lim\limits_{L\to\infty}h_{L}(LX,\omega)=iq\int\limits_{k=-\infty}^{+\infty}\frac{S(k-\frac{X}{2})-S(k+\frac{X}{2})}{\omega+4\pi^{2}ipkX}dk.$
(6)

######  Corollary 2.2 (Convergence of the instability condition as
$L\to\infty$).

Under the assumptions of Theorem 2.1, the linear instability condition (4) of
the periodized problem converges to the linear instability condition for the
infinite line problem, [1, 3],

$\exists
X\in\mathbb{R}\qquad\exists\omega\in\\{z\in\mathbb{C}:\operatorname{Re}(z)>0\\}\qquad\qquad
iq\int\limits_{k=-\infty}^{+\infty}\frac{S(k-\frac{X}{2})-S(k+\frac{X}{2})}{\omega+4\pi^{2}ipkX}dk=1,$
(7)

as $L\to\infty.$

######  Remark 2.3.

We don’t claim that equation (3) is the only way to recover the correct
stability condition. Other random-phase-type strategies to create realizations
of the sea state are also used in the literature, and it seems likely that
they can also be calibrated lead to the same stability condition.

The rest of the paper is organized as follows: Sections 3.1 and 3.2 go over
the derivation of the periodized Alber equation and its stability relation.
They are written in an accessible way that helps motivate and frame the proof,
which can be found in Section 3.3. Methods to check the stability condition
are recalled in Section 4. We explore the role of $L$ in the MI and gMI in
Section 5 \- in particular, it is seen that a short computational domain can
completely stabilize the MI as well as the gMI. Ramifications are discussed in
Section 6.

##  3 Stability analysis for periodized stochastic sea states

###  3.1 Derivation of the stability condition

Let us briefly recall the derivation of the Alber equation in the infinite
space case. Assuming NLS dynamics for the wave envelope on infinite space (1)
and following [1, 3], we can derive the second moment equation

$i\partial_{t}R+\frac{p}{2}\left(\Delta_{x}-\Delta_{y}\right)R+qR(x,y,t)\left[R(x,x,t)-R(y,y,t)\right]=0$
(8)

for the two-space autocorrelation,

$R=R(x,y,t)=E[u(x,t)\overline{u}(y,t)].$ (9)

The standard gaussian closure

$E[u(x,t)\bar{u}(x,t)u(x,t)\bar{u}(y,t)]=2E[u(x,t)\bar{u}(y,t)]E[u(x,t)\bar{u}(x,t)]$

is used to arrive to equation (8), cf. Appendix B of [3] for more details.
Moreover, by separating the wavefield into a homogeneous and non-homogeneous
part

$R(x,y,t)=\Gamma(x-y)+\rho(x,y,t),$ (10)

we derive the standard (infinite space) Alber equation for the evolution of
the inhomogeneous part,

$i\partial_{t}\rho+\frac{p}{2}\left(\Delta_{x}-\Delta_{y}\right)\rho+q\left[\Gamma(x-y)+\rho(x,y,t)\right]\left[\rho(x,x,t)-\rho(y,y,t)\right]=0.$
(11)

Observe that if $\rho=0,$ equation (11) is automatically satisfied. In other
words, any homogeneous autocorrelation $R(x,y,t)=\Gamma(x-y)$ is automatically
a solution of equation (8); the point is to investigate if a particular
homogeneous solution is stable under inhomogeneous perturbations (in which
case an initially small inhomogeneity $\rho$ is guaranteed to stay small and
even disperses) or unstable (in which case any inhomogeneity is expected to
grow rapidly).

Now let us consider the periodization of the NLS on a torus of length $L.$
Starting from equation (2) and performing the same steps as before, the same
Alber equation (11) can still be derived – now equipped with periodic boundary
conditions instead of non-growth at infinity. Moreover, the periodicity of $u$
is inherited by its second moments, and thus we can express moments in terms
of Fourier series on the interval $[-\frac{L}{2},\frac{L}{2}],$

$\Gamma(y)=\sum\limits_{n\in\mathbb{Z}}P_{n}e^{2\pi
i\frac{ny}{L}},\qquad\rho(x,y,t)=\sum\limits_{k,l\in\mathbb{Z}}r_{k,l}(t)e^{2\pi
i\frac{kx+ly}{L}}.$ (12)

For now we treat the Fourier coefficients $P_{n}$ as known. We will revisit
their relationship to the continuous power spectrum $S(k)$ in the next
subsection. Now assuming quasi-homogeneity, i.e.

$\rho(x,y,0)=o(1),$ (13)

we linearise around homogeneity (i.e. drop $O(\rho^{2})$ terms) and express
everything in terms of the Fourier coefficients introduced in (12):

$\begin{array}[]{c}i\sum\limits_{k,l\in\mathbb{Z}}\partial_{t}r_{k,l}(t)e^{2\pi
i\frac{kx+ly}{L}}+\frac{p}{2}\sum\limits_{k,l\in\mathbb{Z}}r_{k,l}(t)\left((2\pi
i\frac{k}{L})^{2}-(2\pi i\frac{l}{L})^{2}\right)e^{2\pi i\frac{kx+ly}{L}}+\\\
+q\left(\sum\limits_{n\in\mathbb{Z}}P_{n}e^{2\pi
i\frac{n(x-y)}{L}}\right)\left(\sum\limits_{k,l\in\mathbb{Z}}r_{k,l}(t)\left[e^{2\pi
i\frac{k+l}{L}x}-e^{2\pi i\frac{k+l}{L}y}\right]\right)=0.\end{array}$ (14)

By taking the inner product with $e^{2\pi i\frac{k^{\prime}x+l^{\prime}y}{L}}$
(and suppressing the primes) we obtain

$\partial_{t}r_{k,l}(t)+i\frac{2\pi^{2}p}{L^{2}}r_{k,l}(t)(k-l)(k+l)-iq\left[P_{-l}-P_{k}\right]\sum\limits_{K\in\mathbb{Z}}r_{K,k+l-K}(t)=0.$
(15)

Denote for brevity

$\vec{r}(t):=\\{r_{k,l}(t)\\}_{k,l\in\mathbb{Z}},$ (16)

now equation (15) is linear and can be summarised as

$\frac{d}{dt}\vec{r}(t)=A\vec{r}(t),$ (17)

in other words the growth of $|r(t)|$ in time can be inferred by the
eigenvalues of the matrix $A.$ This can provide a powerful black-box method
(see [2] for a similar approach in the CSY equation), but would not allow us,
e.g., to show convergence as $L\to\infty.$ Thus, we will proceed to work with
the Laplace transform as in [3].

To help understand better what affects the eigenvalues, first of all we bring
equation (15) to mild form

$\begin{array}[]{c}r_{k,l}(t)=e^{-ip\frac{2\pi^{2}(k+l)(k-l)}{L^{2}}t}r_{k,l}(0)+\int\limits_{\tau=0}^{t}e^{-ip\frac{2\pi^{2}(k+l)(k-l)}{L^{2}}(t-\tau)}iq\left[P_{-l}-P_{k}\right]\sum\limits_{K\in\mathbb{Z}}r_{K,k+l-K}(\tau)d\tau,\end{array}$
(18)

and proceed to define

$f(\xi,t)=\sum\limits_{K\in\mathbb{Z}}r_{K,\xi-K}(t),\qquad\xi\in\mathbb{Z};$
(19)

this corresponds to the $\xi$’th Fourier coefficients of the inhomogeneous
part of the position density $|u(x,t)|^{2}$. (The point is that $f(\xi,t)$
really drives the problem, in the sense that if we now $f(\xi,t)$ then
$r_{k,l}(t)$ can be computed by substitution in equation (18).)

Now set

$k+l=\xi\implies l=\xi-k,\qquad k-l=2k-\xi;$

thus, equation (18) becomes

$r_{k,\xi-k}(t)=e^{-ip\frac{2\pi^{2}\xi(2k-\xi)}{L^{2}}t}r_{k,\xi-k}(0)+\int\limits_{\tau=0}^{t}e^{-ip\frac{2\pi^{2}\xi(2k-\xi)}{L^{2}}(t-\tau)}iq\left[P_{k-\xi}-P_{k}\right]\sum\limits_{K\in\mathbb{Z}}r_{K,\xi-K}(\tau)d\tau.$
(20)

By summing equation (20) in $k$ it follows that

$f(\xi,t)=\phi_{L}(\xi,t)+\int\limits_{\tau=0}^{t}h_{L}(\xi,t-\tau)f(\xi,\tau)d\tau,$
(21)

where

$\phi_{L}(\xi,t)=\sum\limits_{k}e^{-ip\frac{2\pi^{2}\xi(2k-\xi)}{L^{2}}t}r_{k,\xi-k}(0),\qquad
h_{L}(\xi,t-\tau)=\sum\limits_{k}e^{-ip\frac{2\pi^{2}\xi(2k-\xi)}{L^{2}}(t-\tau)}iq\left(P_{k-\xi}-P_{k}\right).$
(22)

This is a simplified form of (15) in the sense that we have one discrete
variable $\xi$ in the new unknown function $f$ instead of two, $k,l$ in the
old variable $\vec{r}$, while we still can recover the full $\vec{r}(t)$
through

$r_{k,l}(t)=e^{-ip\frac{2\pi^{2}(k+l)(k-l)}{L^{2}}t}r_{k,l}(0)+\int\limits_{\tau=0}^{t}e^{-ip\frac{2\pi^{2}(k+l)(k-l)}{L^{2}}(t-\tau)}iq\left[P_{-l}-P_{k}\right]f(k+l,\tau)d\tau$
(23)

(in particular, if $f(t)$ does not grow exponentially, neither does
$\vec{r}(t)$).

By Laplace transform equation (21) becomes

$\widetilde{f}(\xi,\omega)=\widetilde{\phi}_{L}(\xi,\omega)+\widetilde{h}_{L}(\xi,\omega)\widetilde{f}(\xi,\omega)\implies\widetilde{f}(\xi,\omega)=\frac{\widetilde{\phi}_{L}(\xi,\omega)}{1-\widetilde{h}_{L}(\xi,\omega)}.$
(24)

Thus, the existence of poles in the right half plane,

$\exists\omega^{*}\,\,:\operatorname{Re}(\omega^{*})>0\,\,\mbox{ and
}\widetilde{h}_{L}(\xi,\omega^{*})=1$

means that $f(\xi,t)$ exhibits exponential growth. In terms of working out
explicitly the function $\widetilde{h}_{L},$ observe that for any $\omega$
with $\operatorname{Re}\omega>0$ we have

$\begin{array}[]{c}\widetilde{h}_{L}(\xi,\omega)=\int\limits_{t=0}^{\infty}e^{-\omega
t}\left[\sum\limits_{k\in\mathbb{Z}}e^{-ip\frac{2\pi^{2}\xi(2k-\xi)}{L^{2}}t}iq\left(P_{k-\xi}-P_{k}\right)\right]dt=iq\sum\limits_{k\in\mathbb{Z}}\frac{P_{k-\xi}-P_{k}}{\omega+ip\frac{2\pi^{2}\xi(2k-\xi)}{L^{2}}}\end{array}$
(25)

and similarly

$\widetilde{\phi}_{L}(\xi,\omega)=\sum\limits_{k\in\mathbb{Z}}\frac{r_{k,\xi-k}(0)}{\omega+ip\frac{2\pi^{2}\xi(2k-\xi)}{L^{2}}},$
(26)

so that finally

$\begin{array}[]{c}\widetilde{f}(\xi,\omega)=\frac{\sum\limits_{k\in\mathbb{Z}}\frac{r_{k,\xi-k}(0)}{\omega+ip\frac{2\pi^{2}\xi(2k-\xi)}{L^{2}}}}{1-iq\sum\limits_{k\in\mathbb{Z}}\frac{P_{k-\xi}-P_{k}}{\omega+ip\frac{2\pi^{2}\xi(2k-\xi)}{L^{2}}}}.\end{array}$
(27)

###  3.2 Discretizing continuous power spectra

In the previous subsection we took the Fourier coefficients $P_{n}$ for
granted. Here we consider how they are related to the original continuous
power spectrum $S(k).$ Given a spectrum $S(k)$ a standard (periodic)
realization of the envelope corresponding to it can be generated through

$u_{0}(x)=\sum\limits_{j=1}^{M}A_{j}e^{2\pi i({k_{j}x}+\phi_{j})},\qquad
A_{j}=\sqrt{\,\delta k\,S(k_{j})},\qquad\phi_{j}\sim U(0,1),\qquad
k_{j}=j\cdot\delta k,$ (28)

so that $A_{j}^{2}\approx\int\limits_{k_{j-1}}^{k_{j}}S(k)dk$ [46]. Observe
that in general there is some flexibility in the wavenumber discretization,

$k_{j}=j\cdot\delta k,\quad j=1,2,\dots,M,\qquad\delta k=\frac{m}{L},\qquad
m,M\in\mathbb{N}.$ (29)

Now the autocorrelation of the random-phase process (28) is

$\begin{array}[]{c}\Gamma(x-y)=E[u(x)\overline{u}(y)]=\sum\limits_{j,j^{\prime}=1}^{M}E[A_{j}\bar{A}_{j^{\prime}}e^{2\pi
i[k_{j}x-k_{j^{\prime}}y+\phi_{j}-\phi_{j^{\prime}}]}]=\sum\limits_{j,j^{\prime}=1}^{M}A_{j}\bar{A}_{j^{\prime}}e^{2\pi
i[k_{j}x-k_{j^{\prime}}y]}E[e^{2\pi i(\phi_{j}-\phi_{j^{\prime}})}]=\\\
=\sum\limits_{j,j^{\prime}=1}^{M}A_{j}\bar{A}_{j^{\prime}}e^{2\pi
i[k_{j}x-k_{j^{\prime}}y]}\delta_{j,j^{\prime}}=\sum\limits_{j=1}^{M}|A_{j}|^{2}e^{2\pi
ik_{j}(x-y)}.\end{array}$

Therefore

$\begin{array}[]{c}P_{n}=\frac{1}{L}\int\limits_{y=-\frac{L}{2}}^{\frac{L}{2}}e^{-2\pi
i\frac{ny}{L}}\sum\limits_{j=1}^{M}|A_{j}|^{2}e^{2\pi
ik_{j}y}dy=\frac{1}{L}\sum\limits_{j}|A_{j}|^{2}\int\limits_{y=-\frac{L}{2}}^{\frac{L}{2}}e^{2\pi
i(k_{j}-\frac{n}{L})y}dy=\sum\limits_{j}|A_{j}|^{2}\delta_{k_{j},\frac{n}{L}}=\\\
=\sum\limits_{j}\frac{m}{L}P(j\delta k)\delta_{j\delta
k,\frac{n}{L}}=\left\\{\begin{array}[]{cl}\frac{m}{L}P(\frac{n}{L}),&n=mj\mbox{
for some }j\\\ 0,&n\mod m\neq 0\end{array}\right.\end{array}$ (30)

Now let us go back to the parameters $m,M:$ the parameter $M$ controls how
many distinct frequencies/wavenumbers are used, while $m$ controls how close
these wavenumbers are taken. Thus the largest wavenumber included is

$k_{max}=M\cdot\delta k=\frac{Mm}{L}.$

Since we are interested in gravity waves it is reasonable to select a fixed
maximum wavenumber, e.g. $k_{max}\approx 10.$ Thus if $k_{max}=O(1)$ and $L\gg
1,$ this leads to

$M\cdot m=O(L),$ (31)

which in principle could allow for different wavenumber discretization
strategies. Now if we return to the choice

$m=1,\qquad M=k_{max}\cdot L,$ (32)

that was made in the statement of Theorem 2.1, equation (30) implies that

$\begin{array}[]{c}P_{n}=\frac{1}{L}S(\frac{n}{L}).\end{array}$ (33)

###  3.3 Proof of theorem 2.1

In section 3.1 the expression for $\widetilde{h}_{L}(\xi,\omega)$ was derived
in terms of the Fourier coefficients $P_{n},$ cf. equation (25). In section
3.2 it was shown how the scaling of equation (3) for the initial data leads to
the Fourier coefficients $P_{n},$ cf. equation (33). Combining equations (25)
and (33) leads to equation (5).

Now, observe that with the change of variables $k^{\prime}=2k-\xi$ (which
means that $k^{\prime}\in 2\mathbb{Z}-\xi,$ $k=\frac{k^{\prime}+\xi}{2},$
$k-\xi=\frac{k^{\prime}-\xi}{2},$) we have

$\begin{array}[]{c}\displaystyle\widetilde{h}_{L}(\xi,\omega)=\frac{iq}{L}\sum\limits_{k\in\mathbb{Z}}\frac{S(\frac{k-\xi}{L})-S(\frac{k}{L})}{\omega+ip\frac{k2\pi^{2}\xi(2k-\xi)}{L^{2}}}=\par\frac{iq}{2}\sum\limits_{k^{\prime}\in
2\mathbb{Z}-\xi}\frac{S(\frac{k^{\prime}-\xi}{2L})-S(\frac{k^{\prime}+\xi}{2L})}{\omega+4\pi^{2}ip\frac{k^{\prime}}{2L}\frac{\xi}{L}}\frac{2}{L}.\end{array}$
(34)

Now denote $X\in\mathbb{R}$ a fixed number independent from $L$ and take the
joint limit

$L\to\infty,\quad\xi\to\infty\qquad\mbox{ so that }\quad\frac{\xi}{L}\to X.$

In that regime $\widetilde{h}_{L}(\xi,\omega)$ can be seen as a Riemann sum
converging to

$\begin{array}[]{c}\lim\widetilde{h}_{L}=\frac{iq}{2}\int\limits_{s=-\infty}^{+\infty}\frac{S(\frac{s}{2}-\frac{X}{2})-S(\frac{s}{2}+\frac{X}{2})}{\omega+4\pi^{2}ip\frac{s}{2}X}ds.\end{array}$
(35)

Where $k^{\prime}/L\to s$ and $2/L\to ds.$ This double-spacing is because
$k^{\prime}$ in equation (34) spans the odd or even numbers only, depending on
$\xi;$ that’s why $2/L\to ds.$ Equation (6) follows with the change of
variables $k=s/2$ in equation (35).

With regard to the convergence of the Riemann integral in equation (35), note
that as along as $\operatorname{Re}(\omega)\neq 0,$ the denominator does not
become zero, and the integrand inherits the smoothness and decay properties of
the power spectrum $S(k).$ ∎

##  4 Resolving the stability condition

The instability conditions both for the periodized problem and the infinite
line problem (equation (4) and (7) respectively) are of the form

$\exists\xi\qquad\exists\omega\in\\{z\in\mathbb{C}:\operatorname{Re}(z)>0\\}\qquad\qquad\widetilde{h}(\xi,\omega)=1$
(36)

for some known function $\widetilde{h}(\xi,\omega)$ which is by construction
analytic in $\omega$ as long as $\operatorname{Re}(\omega)>0.$ This is a
system of two equations (the real and imaginary part of $\widetilde{h}$) in
three unknowns ($\xi,$ $\operatorname{Re}(\omega),$
$\operatorname{Im}(\omega)$) and as such determining the existence or not of
solutions is not trivial.

A technique that can be used to determine the existence or not of solutions is
based on the argument principle, and was analyzed in detail in [3]. In a
nutshell, it is based on the idea that, given a closed subset
$D\subseteq\\{z\in\mathbb{C}:\operatorname{Re}(z)>0\\}$ with piecewise smooth
boundary $\Gamma=\partial
D\subseteq\\{z\in\mathbb{C}:\operatorname{Re}(z)>0\\},$ the curve
$\widetilde{h}(\xi,\Gamma)=\\{\widetilde{h}(\xi,s)\in\mathbb{C},\,\,s\in\Gamma\\}$
circumscribes
$\widetilde{h}(\xi,D)=\\{\widetilde{h}(\xi,s)\in\mathbb{C},\,\,s\in D\\}.$ A
natural choice for $D$ is a right semicircle of radius $R\frac{1}{\epsilon}$
with the strip $\\{\operatorname{Re}(Z)<\epsilon\\}$ excluded. As $\epsilon\to
0,$ this set exhausts the open right half plane. So now we have $\Gamma$ being
a D-shaped contour consisting of the straight line
$\\{\epsilon+is\\}_{s=-\frac{1}{\epsilon}}^{\frac{1}{\epsilon}}$ and the arc
$\\{\frac{1}{\varepsilon}(cos\theta+i\sin\theta)\\}_{\theta=\mathop{arccos}\varepsilon^{2}}^{-\mathop{arccos}\epsilon^{2}}.$
Computing $\widetilde{h}(\xi,s)$ for $s\in\Gamma$ is straightforward;
moreover, for small enough epsilon the values of $\widetilde{h}$ on the arc
are going to be very close to zero, hence the nontrivial part is the images of
the straight line – this is analogous to Nyquist plots in control theory.

A more elementary technique is just plotting all the values attained on a fine
enough grid in the domain $D.$ Finally, when the existence of an unstable mode
is detected, plotting the distance $|\widetilde{h}-1|$ can help locate its
approximate location and work locally if necessary. Indeed, a number of
numerical choices are required (numerical tolerances etc), it is often helpful
to use several different methods in conjunction.

##  5 MI and gMI can be completely suppressed on a short interval

In this section we go over in some more detail the role of $L,$ the length of
the computational domain, in the MI and gMI. We start by recalling in detail
the MI for the real line in Section 5.1, and for the interval of length $L$ in
section 5.2. Then we investigate numerically the nonlinear phase of the MI in
Section 5.3. We plot the inhomogeneity (in this case the difference from the
plane wave background) and discuss its qualitative behaviour for various
values of $L.$ In Section 5.4 we investigate numerically the nonlinear
evolution of a small inhomogeneity added to a homogeneous background solution
which serves as a simple model realization of a sea state with a very narrow
spectrum. By plotting the inhomogeneity we see a picture completely analogous
to the MI: a short computational domain $L$ can completely stabilize the
problem, i.e. make the inhomogeneity stay small – while on larger
computational domains the inhomogeneity grows by two orders of magnitude. In
all cases the inhomogeneity is numerically zero at the endpoints of the
interval, while the background solution satisfies periodic boundary
conditions.

Figure 2: Space-time plot of the modulus of the inhomogeneity, $|\delta_{2}(x,t)|,$ when propagated on intervals of different length $L$ according to equations (51) – (53). Growth of the inhomogeneity to $O(1)$ is evidence of MI. The bifurcation length $L_{c}$ is reported in equation (49). Top left: $L=0.98L_{c}.$ Top right: $L=1.3L_{c}.$ Bottom left: $L=3L_{c}.$ Bottom right: $L=10L_{c}.$ | $N=0.98$ | $N=1.3$ | $N=3$ | $N=10$
---|---|---|---|---
$j=1$ | $0.0359$ | $2.59$ | $3.08$ | $3.03$
$j=2$ | $0.0393$ | $2.59$ | $3.09$ | $3.03$
$j=3$ | $0.149$ | $2.65$ | $3.17$ | $3.12$
$j=4$ | $0.25$ | $2.69$ | $3.22$ | $3.18$
$j=5$ | $0.03$ | $2.57$ | $3.05$ | $2.96$
Table 1: Maximum value of the inhomogeneity
$\mathop{\mathrm{max}}\limits_{x,t}|\delta(x,t)|$ for $x\in[-L/2,L/2]$ and
$t\in[0,10],$ with $L=NL_{c}$ (cf. equation (49)) and initial inhomogeneity
$\delta_{j}(x,0)$ (cf. equation (52)).

###  5.1 MI on the real line

We briefly go through the standard MI calculation on the real line, for
equation (1). This will prepare the ground for the subsequent analysis on a
bounded interval. Equation (1) is known to be well-posed in Zhidkov spaces
[29, 58]; this is the natural framework for this discussion. Also, since we
are dealing with the focusing NLS, we will assume without loss of generality
that $p,q>0.$ Moreover, equation (1) admits the exact solution

$w(x,t):=Ae^{iqA^{2}t},$ (37)

which is the simplest plane wave solution with amplitude $A>0$. To study the
stability of this plane wave solution one can consider whether small
perturbations grow. This leads to the perturbed initial value problem

$\begin{array}[]{c}i\partial_{t}u+p\Delta u+q|u|^{2}u=0,\qquad
u(x,0)=A(1+\delta_{0}(x))\end{array}$ (38)

for some initial perturbation $\delta_{0}(x)$ which is small in an appropriate
sense, $\delta_{0}=o(1)$. Using the substitution

$u(x,t)=Ae^{iqA^{2}t}(1+\delta(x,t))$ (39)

one readily computes that problem (38) is equivalent to

$\begin{array}[]{c}i\delta_{t}+p\Delta\delta+qA^{2}(\delta+\bar{\delta})+qA^{2}(\delta+\bar{\delta})\delta+qA^{2}|\delta|^{2}(1+\delta),\qquad\delta(x,0)=\delta_{0}(x).\end{array}$
(40)

Dropping higher order terms, we obtain the linearized problem for the
perturbation, namely

$\begin{array}[]{c}i\delta_{t}+p\Delta\delta+qA^{2}(\delta+\bar{\delta}),\qquad\delta(x,0)=\delta_{0}(x).\end{array}$
(41)

By expanding equation (41) into its real and imaginary parts, and denoting

$\delta(x,t)=\alpha(x,t)+i\beta(x,t),$ (42)

we eventually obtain the system

$\begin{array}[]{c}\partial_{tt}\beta+\left(p^{2}\Delta\Delta+2pqA^{2}\Delta\right)\beta=0,\qquad\partial_{t}\alpha+p\Delta\beta=0.\end{array}$
(43)

This now can be solved explicitly with separation of variables, leading to the
construction of the modes

$\begin{array}[]{c}\zeta\in\mathbb{R},\qquad\beta_{\zeta}(x,t)=e^{i[\zeta
x+\omega(\zeta)t]}+c.c.,\qquad\omega^{2}(\zeta)=\zeta^{2}[p^{2}\zeta^{2}-2pqA^{2}].\end{array}$
(44)

(Here c.c. stands for complex conjugate.) More general solutions can be formed
by superpositions of these modes, $\beta(x,t)=\int\limits
M(\zeta)\beta_{\zeta}(x,t)d\zeta.$ The instability is due to the fact that,

$\begin{array}[]{l}|\zeta|<A\sqrt{2\frac{q}{p}}\implies\omega(\zeta)=\pm
i|\omega(\zeta)|\end{array}$ (45)

and thus the corresponding modes $\beta_{\zeta}(x,t)$ of equation (44) contain
an exponentially growing component. Therefore the solution $\delta$ of the
linearized equation (41) generally grows exponentially in time, i.e. the plane
wave solution of the NLS is linearly unstable. Moreover the unstable
wavenumbers and their rate of growth follow from this analysis.

In this problem there are three parameters, $p,q,A.$ However by rescaling the
problem according to $\tau=qA^{2}\,t,$ $\chi=A\sqrt{\frac{q}{p}}\,x,$
$U(\chi,\tau)=\frac{1}{A}\,u(x,t),$ the equation is mapped to the “canonical”
NLS $i\partial_{\tau}U+\partial_{\chi\chi}U+|U|^{2}U=0,$ which has the plane
wave solution $W(\chi,\tau)=e^{i\tau}.$ That is, the exact values of $p,q,A$
don’t play an important role and the nature of the modulation instability is
the same for any $p,q,A>0.$ This widely understood universality of the MI may
have contributed to an implicit expectation that the periodized problem
automatically enjoys similar properties.

###  5.2 MI on an interval of length $L$

Now let us consider the NLS equation on an interval of length $L$ equipped
with periodic boundary conditions, namely (2). The steps described above in
equations (37)-(43) were carried out with boundary conditions of boundedness
at infinity; however, one readily checks that each step is equally valid with
the periodic boundary conditions on $[-L/2,L/2].$ Indeed, the plane wave is
still a solution, the algebra leading to (40) is the same, and finally the
linearization and separation of real and imaginary parts are the same. So now
we have to solve the problem

$\begin{array}[]{c}\partial_{tt}\beta+\left(p^{2}\Delta\Delta+2pqA^{2}\Delta\right)\beta=0,\quad\beta(-\frac{L}{2},t)=\beta(\frac{L}{2},t),\quad\partial_{x}\beta(-\frac{L}{2},t)=\partial_{x}\beta(\frac{L}{2},t).\end{array}$
(46)

This is the first time that the parameter $L$ really comes into play, and
separation of variables now leads to the discrete modes

$\begin{array}[]{c}\displaystyle\beta_{n}(x,t)=e^{2\pi
i(\frac{nx}{L}+\omega_{n}t)}+c.c.,\qquad(2\pi\omega_{n})^{2}=(\frac{2\pi
n}{L})^{2}[p^{2}(\frac{2\pi n}{L})^{2}-2pqA^{2}].\end{array}$ (47)

While this is analogous to equation (44), there is a very important
difference: depending on the values of $p,q,A,L,$ there may not be even one
unstable mode. Indeed, $\omega_{0}=0,$ and for any $0\neq n\in\mathbb{Z}$ we
have

$\omega_{n}^{2}<0\iff L>\frac{2\pi|n|}{A}\sqrt{\frac{p}{2q}}.$ (48)

That is, the computational domain $L$ has to be larger than

$L_{c}:=\frac{2\pi\sqrt{p}}{A\sqrt{{2q}}},$ (49)

otherwise the MI is completely suppressed. Observe moreover that the
lengthscale $L_{c}$ does not depend on the initial inhomogeneity
$\delta_{0}(x),$ and becomes larger when the nonlinearity becomes weaker (i.e.
when $q,A>0$ decrease). For water waves with typical wavenumber $k_{0}$
(equivalently typical wavelength $\lambda_{0}=2\pi/k_{0}$), the coefficients
are $p=\sqrt{g}/(8k_{0}^{3/2}),$ $q=\sqrt{g}k_{0}^{5/2}/2$ [43] leading to

$L_{c}=\frac{1}{4\pi\sqrt{2}}\frac{\lambda_{0}^{2}}{A}.$ (50)

Note that the waves that would be dangerous for ships in the ocean and carry
most of the surface wave energy have wavelengths in the hundreds of meters.

Figure 3: Space-time plot of the modulus of the inhomogeneity,
$|\delta(x,t)|,$ when propagated on intervals of different length $L$
according to equations (54) – (58). The reference length $L_{0}$ is reported
in equation (54). Top left: $L=L_{0}.$ Top right: $L=2L_{0}.$ Bottom left:
$L=3L_{0}.$ Bottom right: $L=10L_{0}.$

###  5.3 Numerical investigation of the fully nonlinear MI

In what follows we present numerical solutions of the problem

$\begin{array}[]{c}i\partial_{t}u+p\Delta u+q|u|^{2}u=0,\qquad
u(x,0)=A(1+\delta_{j}(x))\\\
u(-\frac{L}{2},t)=u(\frac{L}{2},t),\quad\partial_{x}u(-\frac{L}{2},t)=\partial_{x}u(\frac{L}{2},t),\end{array}$
(51)

for the initial inhomogeneities

$\begin{array}[]{c}\delta_{1}(x)=0.03\cdot
A\cdot\mathrm{sech}(15x)\cos(5x),\qquad\delta_{2}(x)=0.03\cdot
A\cdot\mathrm{sech}(15x),\\\ \delta_{3}(x)=0.03\cdot A\cdot
e^{-3x^{2}},\quad\delta_{4}(x)=0.03\cdot A\cdot
e^{-x^{4}},\quad\delta_{5}(x)=0.06\cdot A\cdot x\cdot e^{-x^{4}},\end{array}$
(52)

and $p=q=A=1.$ The bifurcation length in this case is $L_{c}\approx 4.45;$ the
solutions are computed for $L=0.98L_{c},$ $L=1.3L_{c},$ $L=3L_{c}$ and
$L=10L_{c}.$ The computation time is $t\in[0,10].$ The numerical scheme used
is a relaxation in time with second-order finite differences in space, and it
satisfies mass and energy conservation on the discrete level [9, 10]. The
inhomogeneity $\delta(x,t)$ for $t\in[0,10]$ is defined as

$\delta(x,t):=(u(x,t)-Ae^{iqA^{2}t})\frac{1}{A}e^{-iqA^{2}t}$ (53)

consistently with equation (39). Numerical results can be found in Figure 2 as
well as in Table 1.

The initial inhomogeneities $\delta_{1}$ and $\delta_{2}$ are localized
$\mathrm{sech}$ bumps, and are quite similar to each other. $\delta_{3}$ and
$\delta_{4}$ are also bumps with different profiles, and $\delta_{5}$ is a
localised wave. In Table 1 the maximum moduli of the inhomogeneities are
recorded when propagated on intervals of different lengths. There is a clear
confirmation of the abrupt bifurcation, namely an abrupt change in behaviour
from $L=0.98L_{c}$ to $L=1.3L_{c}.$

On the larger intervals $L=3L_{c},$ $L=10L_{c},$ coherent structures emerge.
These are a sign that the MI has been resolved in a manner qualitatively
similar to what would happen on the real line (that is, until the structures
reach close to the boundary of the computational domain). These structures all
involve localised maxima supported on neighbourhoods of size roughly $L_{c}$
each. In the larger domain, the well-known space-time cone [14, 13, 12] is
clearly visible.

The presence of a space-time cone in the nonlinear phase of the MI, reported
for a wide array of dispersive problems [14, 13, 12], means that there is
generically a finite speed of propagation of the boundary between the plane
wave region (for $|x|\gg 1$) and the “inhomogeneous” region. Moreover, the
“inhomogeneous” region seems to settle in a pattern determined by the equation
and not by the specific initial condition (see also Figures 1 and A1 of [6]).
These patterns that appear within the cone seem to be related to so-called
dispersive shock solutions [11].

Crucially, when the length of the domain becomes $L<L_{c},$ we don’t merely
see a periodised or slightly off version of the cone, but something completely
different: the amplification is small or non-existent and there is no spatial
pattern at all, see left graph in Figure 2. Indeed, for all initial data, when
$L=0.98L_{c},$ the inhomogeneity initially disperses rapidly and then slowly
grows in homogeneous way over the whole computational domain.

###  5.4 An instance of gMI on an interval of length $L$

We set

$L_{0}=4.4518,\qquad p=q=1$ (54)

and consider an initial condition

$u_{0}(x)=0.018e^{-2\pi i\frac{x}{L_{0}}}+0.899+0.1252e^{2\pi
i\frac{x}{L_{0}}}.$ (55)

This $u_{0}$ can be thought of as a perturbation of the plane wave solution
(37) at $t=0,$ and the lengthscale $L_{0}$ is numerically close to the $L_{c}$
of problem (51) – (53). The goal of this computation is to see whether a
phenomenology analogous to what was seen in Section 5.3 arises beyond plane
wave backgrounds. Thus the initial condition $u_{0}(x)$ plays the role of a
realization of a very narrow spectrum.

In this setting the “homogeneous background solution” $u(x,t)$ is the solution
of

$\begin{array}[]{c}i\partial_{t}u+p\Delta u+q|u|^{2}u=0,\qquad
u(x,0)=u_{0}(x)\\\
u(-\frac{L}{2},t)=u(\frac{L}{2},t),\quad\partial_{x}u(-\frac{L}{2},t)=\partial_{x}u(\frac{L}{2},t),\end{array}$
(56)

and the inhomogeneously perturbed solution $V(x,t)$ is the solution of

$\begin{array}[]{c}i\partial_{t}v+p\Delta v+q|v|^{2}v=0,\qquad
v(x,0)=u_{0}(x)+0.07\mathrm{sech}(15x)\cos(5x)\\\
v(-\frac{L}{2},t)=v(\frac{L}{2},t),\quad\partial_{x}v(-\frac{L}{2},t)=\partial_{x}v(\frac{L}{2},t).\end{array}$
(57)

The inhomogeneity $\delta(x,t)$ is thus defined as

$\delta(x,t)=v(x,t)-u(x,t),$ (58)

and the main question is whether $\delta(x,t)$ grows significantly or not for
different values of $L.$

Immediately when we consider background solutions other than plane waves, we
observe that the length $L$ and the initial conditions $u_{0},$ $\delta$ have
to be compatible: the homogeneous $u_{0}$ has to satisfy periodic BCs on $L,$
hence in this case we will consider only integer multiples of $L_{0}.$ On the
other hand, the localized $\delta(x,0)$ has to be numerically zero at the ends
of the domain, up to satisfactory accuracy. By taking $L=NL_{0},$ for
$N\in\\{1,2,3,10\\}$ both of these requirements are satisfied. The final time
is taken to be $T=10.3.$

The numerical scheme used is the same as in Section 5.3 (a relaxation in time
with second-order finite differences in space, and it satisfies mass and
energy conservation on the discrete level [9, 10]). The time step used is
$dt=4\cdot 10^{-3}$ and the space mesh size is $dx=4\cdot 10^{-3}.$ The
modulus $|\delta(x,t)|$ can be found in Figure 3, and it is clear that for the
short enough computational domain $L=L_{0}$ the growth of the inhomogeneity is
completely suppressed (in particular the inhomogeneity in the $L^{\infty}$
norm is not even doubled). On the other hand, the space-time cone, inside
which the inhomogeneity grows by more than an order of magnitude is very
clearly visible when $L=10L_{0}.$

To provide some more context, the background solution and the inhomogeneity
are plotted in the intitial and final time in Figures 4 and 5. Observe how the
background solution qualitatively behaves similarly to the plane wave
solution.

Figure 4: Plots of the real part (blue) and imaginary part (red) of the
background solution $u,$ as in equation (56), on the interval with length
$L=10L_{0}.$ Left: for $t=0.$ Right: for $t=T.$

Figure 5: Plots of the real part (blue) and imaginary part (red) of the
perturbed solution $v,$ as in equation (57), on the interval with length
$L=10L_{0}.$ Left: for $t=0.$ Right: for $t=T.$

##  6 Conclusions

In this paper we presented for the first time a construction of the random
phase approximation in the NLS equation that comes with a convergence result:
using the scaling of equation (3), for large enough $L,$ the instability
condition of the truncated problem (2) converges to that of the original
problem (1).

In the proof it is clear that equation (3) is not the only way to prepare
initial data that would have this property. There are other ways to prepare
initial data in the literature, e.g. with randomized amplitudes as well as
phases of the modes [5], and it is possible that these could also lead to a
consistent stability condition. On the other hand, recent results in different
regimes [16, 15, 23] highlight that different ways to prepare the initial
conditions may lead to genuinely different outcomes.

The importance of using large enough computational domains is well understood
for the classical MI. This was further elaborated in Section 5, where a
lengthscale of $L_{c}=O(\lambda_{0}^{2})$ was found to be necessary for the MI
to be present even for plane waves. If $\lambda_{0}$ is in the hundreds of
meters, as in ocean waves, this could mean $L_{c}$ hundreds of wavelengths
long.

The role of $L$ in the gMI is less well understood. In works where the onset
of gMI was actively investigated, it has been reported that long computational
domains are needed; indeed in [51], it was reported that a computational
domain of $15L_{c}$ was needed for a JONSWAP sea state (where $L_{c}$ is the
corresponding critical length for the plane wave with the same wavelength and
wave height as the peak values for the sea state). Around $130$ wavelengths
were found to be required in [6], while in [48] the computational domain is
about $100$ wavelengths. On the other hand, in many papers the question of the
length of computational domain seems to not be addressed in any detail. While
more work is needed to conclusively resolve this question for realistic ocean
wave spectra, we have demonstrated here (in Section 5.3 and Figure 3) that the
capacity of short computational domains to suppress the growth of
inhomogeneities extends beyond plane wave backgrounds.

This work was carried out assuming NLS dynamics. This is of course an
approximation for ocean waves, but it should be noted that recent rigorous
works have shown that the MI for the fully nonlinear water waves problem is
accurately captured by the nonlinear Schrödinger equation [8, 45].
Furthermore, the analysis carried out here for the truncated NLS problem can
be carried out for the truncated Zakharov problem in the spirit of [2], thus
producing a stability relation for the periodized CSY equation (which is
2-dimensional and broadband, since it assumes Zakharov dyamics). This would
also allow the investigation of crossing seas scenarios, either in fully 2D
through the aforementioned CSY equation, or using an Alber system for the
interaction of two unidirectional wavetrains at an angle [32, 5]. Crossing
seas are very important with regard to real-world marine safety, and much more
can be done to investigate them [19, 50, 54]. More broadly, incorporating more
physics into the model, such as wind input [4, 55] and vorticity [18] seems to
be a very promising direction with regard to real-world ocean waves.

Acknowledgment. AA would like to thank Prof. G. A. Athanassoulis, Prof. T.
Sapsis, Dr. O. Gramstad and Dr. T. Tang for helpful discussions.

## References

  * [1] I. E. Alber, The Effects of Randomness on the Stability of Two-Dimensional Surface Wavetrains, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 363 (1978), pp. 525–546.
  * [2] D. Andrade and M. Stiassnie, New solutions of the C.S.Y. equation reveal increases in freak wave occurrence, Wave Motion, 97 (2020), p. 102581.
  * [3] A. Athanassoulis, G. Athanassoulis, M. Ptashnyk, and T. Sapsis, Strong solutions for the Alber equation and stability of unidirectional wave spectra, Kinetic and Related Models, 13 (2020), pp. 703–737.
  * [4] A. Athanassoulis, G. Athanassoulis, and T. Sapsis, Localized instabilities of the Wigner equation as a model for the emergence of Rogue Waves, Journal of Ocean Engineering and Marine Energy, 3 (2017), pp. 353–372.
  * [5] A. Athanassoulis and O. Gramstad, Modelling of Ocean Waves with the Alber Equation: Application to Non-Parametric Spectra and Generalisation to Crossing Seas, Fluids, 6 (2021), p. 291.
  * [6] A. G. Athanassoulis, Phase Resolved Simulation of the Landau-Alber Stability Bifurcation, Fluids, 8 (2023).
  * [7] T. B. Benjamin and J. E. Feir, The disintegration of wave trains on deep water, Journal of Fluid Mechanics, 27 (1967), pp. 417–430.
  * [8] M. Berti, A. Maspero, and P. Ventura, Full description of Benjamin-Feir instability of Stokes waves in deep water, Inventiones Mathematicae, 230 (2022), pp. 651–711.
  * [9] C. Besse, A relaxation scheme for the nonlinear Schrodinger equation, SIAM Journal on Numerical Analysis, 42 (2004), pp. 934–952.
  * [10] C. Besse, S. Descombes, G. Dujardin, and I. Lacroix-Violet, Energy-preserving methods for nonlinear Schrödinger equations, IMA Journal of Numerical Analysis, 41 (2021), pp. 618–653.
  * [11] G. Biondini, Riemann problems and dispersive shocks in self-focusing media, Physical Review E, 98 (2018).
  * [12] G. Biondini, S. Li, D. Mantzavinos, and S. Trillo, Universal Behavior of Modulationally Unstable Media, SIAM Review, 60 (2018), pp. 888–908.
  * [13] G. Biondini and D. Mantzavinos, Long-Time Asymptotics for the Focusing Nonlinear Schrödinger Equation with Nonzero Boundary Conditions at Infinity and Asymptotic Stage of Modulational Instability, Communications on Pure and Applied Mathematics, 70 (2017), pp. 2300–2365.
  * [14] G. Biondini and D. D. Mantzavinos, Universal Nature of the Nonlinear Stage of Modulational Instability, Physical Review Letters, 116 (2016), pp. 1–5.
  * [15] T. Buckmaster, P. Germain, Z. Hani, and J. Shatah, Onset of the wave turbulence description of the longtime behavior of the nonlinear Schrödinger equation, Inventiones Mathematicae, 225 (2021), pp. 787–855.
  * [16] T. Congy, G. A. El, G. Roberti, A. Tovbis, S. Randoux, and P. Suret, Statistics of extreme events in integrable turbulence.  jul 2023.
  * [17] D. R. Crawford, P. G. Saffman, and H. C. Yuen, Evolution of a random inhomogeneous field of nonlinear deep-water gravity waves, Wave Motion, 2 (1980), pp. 1–16.
  * [18] C. W. Curtis and M. Murphy, Evolution of Spectral Distributions in Deep-Water Constant Vorticity Flows, Water Waves, 2 (2020), pp. 361–380.
  * [19] S. Davison, A. Benetazzo, F. Barbariol, G. Ducrozet, J. Yoo, and M. Marani, Space-time statistics of extreme ocean waves in crossing sea states, Frontiers in Marine Science, 9 (2022), p. 1002806.
  * [20] G. Dematteis, T. Grafke, and E. Vanden-Eijnden, Rogue Waves and Large Deviations in Deep Sea, Proceedings of the National Academy of Sciences, (2018), p. 201710670.
  * [21] Y. Deng and Z. Hani, On the derivation of the wave kinetic equation for NLS, Forum of Mathematics, Pi, 9 (2021).
  * [22] , Propagation of chaos and the higher order statistics in the wave kinetic theory, (2021), pp. 1–47.
  * [23] , Full derivation of the wave kinetic equation, Inventiones Mathematicae, 233 (2023), pp. 543–724.
  * [24] K. Dysthe, H. E. Krogstad, and P. Müller, Oceanic Rogue Waves, Annual Review of Fluid Mechanics, 40 (2008), pp. 287–310.
  * [25] K. B. Dysthe, K. Trulsen, H. E. Krogstad, and H. Socquet-Juglard, Evolution of a narrow-band spectrum of random surface gravity waves, Journal of Fluid Mechanics, 478 (2003), pp. 1–10.
  * [26] F. Fedele, On the kurtosis of deep-water gravity waves, Journal of Fluid Mechanics, 782 (2015), pp. 25–36.
  * [27] F. Fedele, J. Brennan, S. Ponce De León, J. Dudley, and F. Dias, Real world ocean rogue waves explained without the modulational instability, Scientific Reports, 6 (2016), pp. 1–11.
  * [28] F. Fedele, Z. Cherneva, M. A. Tayfun, and C. Guedes Soares, Nonlinear Schrödinger invariants wave statistics, Physics of Fluids, 22 (2010), pp. 1–9.
  * [29] C. Gallo, Schrödinger group on Zhidkov spaces, Advances in Differential Equations, 9 (2004), pp. 509–538.
  * [30] P. Germain, A. D. Ionescu, and M.-B. Tran, Optimal local well-posedness theory for the kinetic wave equation, Journal of Functional Analysis, 279 (2020), p. 108570.
  * [31] O. Gramstad, Modulational Instability in JONSWAP Sea States Using the Alber Equation, in ASME 2017 36th International Conference on Ocean, Offshore and Arctic Engineering, 2017.
  * [32] O. Gramstad, E. Bitner-Gregersen, K. Trulsen, and J. C. Nieto Borge, Modulational instability and rogue waves in crossing sea states, Journal of Physical Oceanography, 48 (2018), pp. 1317–1331.
  * [33] O. Gramstad and K. Trulsen, Influence of crest and group length on the occurrence of freak waves, Journal of Fluid Mechanics, 582 (2007), pp. 463–472.
  * [34] E. Gubesch, N. Abdussamie, I. Penesis, and C. Chin, Physical and numerical modelling of extreme wave conditions, Ocean Engineering, 283 (2023), p. 115055.
  * [35] K. Hasselmann, On the non-linear energy transfer in a gravity-wave spectrum: Part 1. General theory, Journal of Fluid Mechanics, 12 (1962), pp. 481–500.
  * [36] K. Hasselmann, S. Hasselmann, E. Bauer, P. A. E. M. Janssen, L. Bertotti, P. Lionello, V. C. Cardone, J. A. Greenwood, M. Reistad, L. Zambresky, and J. A. Ewing, The WAM Model - A Third Generation Ocean Wave Prediction Model, Journal of Physical Oceanography, 18 (1988), pp. 1775–1810.
  * [37] P. A. Janssen, Progress in ocean wave forecasting, Journal of Computational Physics, 227 (2008), pp. 3572–3594.
  * [38] P. A. E. M. Janssen, Nonlinear Four-Wave Interactions and Freak Waves, Journal of Physical Oceanography, 33 (2003), pp. 863–884.
  * [39] C. Kharif and E. Pelinovsky, Physical mechanisms of the rogue wave phenomenon, European Journal of Mechanics, B/Fluids, 22 (2003), pp. 603–634.
  * [40] G. J. Komen, L. Cavaleri, M. Donelan, K. Hasselmann, S. Hasselmann, and P. A. E. M. Janssen, Dynamics and Modelling of Ocean Waves, Cambridge University Press, 1994.
  * [41] Y.-C. Lee, M. Brühl, D.-J. Doong, and S. Wahls, Nonlinear Fourier Classiication of 663 Rogue Waves Measured in the Philippine Sea Nonlinear Fourier Classification of 663 Rogue Waves Measured in the Philippine Sea, (2022).
  * [42] M. S. Longuet-Higgins, On the nonlinear transfer of energy in the peak of a gravity-wave spectrum: a simplified model, Proceedings of the Royal Society A, 347 (1976), pp. 311–328.
  * [43] C. C. Mei, M. Stiassnie, and D. K.-P. Yue, Theory and Applications of Ocean Surface Waves, vol. 23 of Advanced Series on Ocean Engineering, World Scientific, jul 2005.
  * [44] N. Mori and T. Yasuda, Effects of high-order nonlinear wave-wave interactions on gravity waves, in Rogue Waves 2000, G. Athanassoulis and M. Olagnon, eds., Brest, 2000.
  * [45] H. Q. Nguyen and W. A. Strauss, Proof of Modulational Instability of Stokes Waves in Deep Water, Communications on Pure and Applied Mathematics, 76 (2023), pp. 1035–1084.
  * [46] M. K. Ochi, Ocean Waves : the Stochastic Approach, Cambridge University Press, 1998.
  * [47] M. Olagnon and G. A. Athanassoulis, eds., Rogue Waves 2000, Ifremer, 2000.
  * [48] M. Onorato, A. Osborne, R. Fedele, and M. Serio, Landau damping and coherent structures in narrow-banded 1 + 1 deep water gravity waves, Physical Review E, 67 (2003), p. 046305.
  * [49] M. Onorato, A. R. Osborne, M. Serio, L. Cavaleri, C. Brandini, and C. T. Stansberg, Extreme waves, modulational instability and second order theory: wave flume experiments on irregular waves, European Journal of Mechanics, B/Fluids, 25 (2006), pp. 586–601.
  * [50] M. Onorato, D. Proment, and A. Toffoli, Freak waves in crossing seas, European Physical Journal: Special Topics, 185 (2010), pp. 45–55.
  * [51] A. Ribal, A. V. Babanin, I. Young, A. Toffoli, and M. Stiassnie, Recurrent solutions of the Alber equation initialized by Joint North Sea Wave Project spectra, Journal of Fluid Mechanics, 719 (2013), pp. 314–344.
  * [52] W. E. Rogers, Phase-Averaged Wave Models, in Ocean Wave Dynamics, WORLD SCIENTIFIC, apr 2020, pp. 163–204.
  * [53] J. Skipp, J. Laurie, and S. V. Nazarenko, An effective semilocal model for wave turbulence in two-dimensional nonlinear optics, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 479 (2023).
  * [54] J. N. Steer, M. L. Mcallister, A. G. L. Borthwick, and T. S. V. D. Bremer, Experimental Observation of Modulational Instability in Crossing Surface Gravity Wavetrains, Fluids, 4 (2019), pp. 1–15.
  * [55] A. Toffoli, A. Alberello, H. Clarke, F. Nelli, A. Benetazzo, F. Bergamasco, B. N. Ntamba, M. Vichi, and M. Onorato, Observations of rogue seas in the Southern Ocean, (2023), pp. 1–7.
  * [56] V. E. Zakharov, Stability of periodic waves of finite amplitude on the surface of a deep fluid, Journal of Applied Mechanics and Technical Physics, 9 (1968), pp. 190–194.
  * [57] V. E. Zakharov and L. A. Ostrovsky, Modulation instability: The beginning, Physica D: Nonlinear Phenomena, 238 (2009), pp. 540–548.
  * [58] P. Zhidkov, The Cauchy problem for the nonlinear Schrödinger equation, tech. rep., Joint Inst. for Nuclear Research No. JINR-R–5-87-373, 1987.