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

The nonlinear Schr\"odinger 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.


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, freeboundary 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].

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 deltafunction) 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 Op1q homogeneous part and an op1q 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 α, a parameter directly affecting the severity of the sea state, and γ, 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 α and γ 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.

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 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].
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 upx, tq.The original physical problem is typically framed in infinite space [43], with initial data U 0 pxq being a realization of a homogeneous random process with power spectrum Spkq.
On the other hand, the truncated problem that we can readily simulate is the periodized NLS with initial data of the form u 0 pxq " ř A j e 2πipk j x`ϕ j q where the ϕ j are i.i.d.random variables uniformly distributed in r0, 1q.
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 OpLq 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.

Main results
Theorem 2.1 (Periodized Alber equation).Let Spkq be a smooth power spectrum with supp Spkq Ď r0, k max s, and consider equation (2) with stochastic initial data upx, 0q " u 0 pxq where The Alber equation for the evolution in time of inhomogeneity ρpx, y, tq of this problem is iB t ρ `p 2 p∆ x ´∆y q ρ `q rΓpx ´yq `ρpx, y, tqs rρpx, x, tq ´ρpy, y, tqs " 0, x, y P r L 2 , L 2 s, equipped with periodic boundary conditions, where Γpxq " and the linear instability condition for this problem is where (5) Corollary 2.2 (Convergence of the instability condition as L Ñ 8).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], as L Ñ 8.
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

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 iB t R `p 2 p∆ x ´∆y q R `qRpx, y, tq rRpx, x, tq ´Rpy, y, tqs " 0 (8) for the two-space autocorrelation, R " Rpx, y, tq " Erupx, tqupy, tqs.
Observe that if ρ " 0, equation ( 11) is automatically satisfied.In other words, any homogeneous autocorrelation Rpx, y, tq " Γpx ´yq 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 ρ 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 r´L 2 , L 2 s, Γpyq " For now we treat the Fourier coefficients P n as known.We will revisit their relationship to the continuous power spectrum Spkq in the next subsection.Now assuming quasi-homogeneity, i.e. ρpx, y, 0q " op1q, we linearise around homogeneity (i.e. drop Opρ 2 q terms) and express everything in terms of the Fourier coefficients introduced in (12): By taking the inner product with e 2πi k 1 x`l 1 y L (and suppressing the primes) we obtain Denote for brevity ⃗ rptq :" tr k,l ptqu k,lPZ , now equation ( 15) is linear and can be summarised as in other words the growth of |rptq| 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 Ñ 8. 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 r k,l ptq " e ´ip 2π 2 pk`lqpk´lq and proceed to define f pξ, tq " this corresponds to the ξ'th Fourier coefficients of the inhomogeneous part of the position density |upx, tq| 2 .(The point is that f pξ, tq really drives the problem, in the sense that if we now f pξ, tq then r k,l ptq can be computed by substitution in equation (18).)Now set k `l " ξ ùñ l " ξ ´k, k ´l " 2k ´ξ; thus, equation (18) becomes By summing equation (20) in k it follows that where This is a simplified form of (15) in the sense that we have one discrete variable ξ in the new unknown function f instead of two, k, l in the old variable ⃗ r, while we still can recover the full ⃗ rptq through r k,l ptq " e ´ip 2π 2 pk`lqpk´lq pt´τ q iq rP ´l ´Pk s f pk `l, τ qdτ (in particular, if f ptq does not grow exponentially, neither does ⃗ rptq).By Laplace transform equation ( 21) becomes Thus, the existence of poles in the right half plane, Dω ˚: Repω ˚q ą 0 and r h L pξ, ω ˚q " 1 means that f pξ, tq exhibits exponential growth.In terms of working out explicitly the function r h L , observe that for any ω with Re ω ą 0 we have and similarly so that finally r f pξ, ωq "

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 Spkq.Given a spectrum Spkq a standard (periodic) realization of the envelope corresponding to it can be generated through A j e 2πipk j x`ϕ j q , A j " b δk Spk j q, ϕ j " U p0, 1q, k j " j ¨δk, Spkqdk [46].Observe that in general there is some flexibility in the wavenumber discretization, Now the autocorrelation of the random-phase process ( 28) is Γpx ´yq " Erupxqupyqs " M ř j,j 1 "1 ErA j Āj 1 e 2πirk j x´k j 1 y`ϕ j ´ϕj 1 s s " M ř j,j 1 "1 A j Āj 1 e 2πirk j x´k j 1 ys Ere 2πipϕ j ´ϕj 1 q s " " M ř j,j 1 "1 A j Āj 1 e 2πirk j x´k j 1 ys δ j,j 1 " M ř j"1 |A j | 2 e 2πik j px´yq . Therefore " m L P p n L q, n " mj for some j 0, n mod m ‰ 0 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 ¨δk " M m L .
Since we are interested in gravity waves it is reasonable to select a fixed maximum wavenumber, e.g.k max « 10.Thus if k max " Op1q and L " 1, this leads to which in principle could allow for different wavenumber discretization strategies.Now if we return to the choice m " 1, that was made in the statement of Theorem 2.1, equation (30) implies that

Proof of theorem 2.1
In section 3.1 the expression for r h L pξ, ωq 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 1 " 2k ´ξ (which means that Now denote X P R a fixed number independent from L and take the joint limit In that regime r h L pξ, ωq can be seen as a Riemann sum converging to Where k 1 {L Ñ s and 2{L Ñ ds.This double-spacing is because k 1 in equation ( 34) spans the odd or even numbers only, depending on ξ; that's why 2{L Ñ 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 Repωq ‰ 0, the denominator does not become zero, and the integrand inherits the smoothness and decay properties of the power spectrum Spkq.
for some known function r hpξ, ωq which is by construction analytic in ω as long as Repωq ą 0. This is a system of two equations (the real and imaginary part of r h) in three unknowns (ξ, Repωq, Impωq) 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 Ď tz P C : Repzq ą 0u with piecewise smooth boundary Γ " BD Ď tz P C : Repzq ą 0u, the curve r hpξ, Γq " t r hpξ, sq P C, s P Γu circumscribes r hpξ, Dq " t r hpξ, sq P C, s P Du.A natural choice for D is a right semicircle of radius R 1 ϵ with the strip tRepZq ă ϵu excluded.As ϵ Ñ 0, this set exhausts the open right half plane.So now we have Γ being a D-shaped contour consisting of the straight line tϵ `isu s"´1 ϵ and the arc t 1 ε pcosθ `i sin θqu ´arccos ϵ 2 θ"arccos ε 2 .Computing r hpξ, sq for s P Γ is straightforward; moreover, for small enough epsilon the values of r 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 | r 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.

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.

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 wpx, tq :" Ae iqA 2 t , ( 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 iB t u `p∆u `q|u| 2 u " 0, upx, 0q " Ap1 `δ0 pxqq (38) for some initial perturbation δ 0 pxq which is small in an appropriate sense, δ 0 " op1q.Using the substitution upx, tq " Ae iqA 2 t p1 `δpx, tqq one readily computes that problem ( 38) is equivalent to iδ t `p∆δ `qA 2 pδ `δq `qA 2 pδ `δqδ `qA 2 |δ| 2 p1 `δq, δpx, 0q " δ 0 pxq.
By expanding equation ( 41) into its real and imaginary parts, and denoting δpx, tq " αpx, tq `iβpx, tq, we eventually obtain the system This now can be solved explicitly with separation of variables, leading to the construction of the modes (Here c.c. stands for complex conjugate.)More general solutions can be formed by superpositions of these modes, βpx, tq " ş M pζqβ ζ px, tqdζ.The instability is due to the fact that, and thus the corresponding modes β ζ px, tq of equation ( 44) contain an exponentially growing component.Therefore the solution δ 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 τ " qA 2 t, χ " A b q p x, U pχ, τ q " 1 A upx, tq, the equation is mapped to the "canonical" NLS iB τ U Bχχ U `|U | 2 U " 0, which has the plane wave solution W pχ, τ q " e iτ .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.

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 r´L{2, L{2s.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 This is the first time that the parameter L really comes into play, and separation of variables now leads to the discrete modes 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, ω 0 " 0, and for any 0 ‰ n P Z we have That is, the computational domain L has to be larger than otherwise the MI is completely suppressed.Observe moreover that the lengthscale L c does not depend on the initial inhomogeneity δ 0 pxq, 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 λ 0 " 2π{k 0 ), the coefficients are p " ?g{p8k 3{2 0 q, q " ?gk 5{2 0 {2 [43] leading to 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.

Numerical investigation of the fully nonlinear MI
In what follows we present numerical solutions of the problem iB t u `p∆u `q|u| 2 u " 0, upx, 0q " Ap1 `δj pxqq up´L 2 , tq " up L 2 , tq, B x up´L 2 , tq " B x up L 2 , tq, for the initial inhomogeneities δ 1 pxq " 0.03 ¨A ¨sechp15xq cosp5xq, δ 2 pxq " 0.03 ¨A ¨sechp15xq, δ 3 pxq " 0.03 ¨A ¨e´3x 2 , δ 4 pxq " 0.03 ¨A ¨e´x 4 , δ 5 pxq " 0.06 ¨A ¨x ¨e´x 4 , and p " q " A " 1.The bifurcation length in this case is L c « 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 P r0, 10s.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 δpx, tq for t P r0, 10s is defined as δpx, tq :" pupx, tq ´Ae iqA 2 t q 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 δ 1 and δ 2 are localized sech bumps, and are quite similar to each other.δ 3 and δ 4 are also bumps with different profiles, and δ 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| " 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.

An instance of gMI on an interval of length L
We set L 0 " 4.4518, p " q " 1 and consider an initial condition u 0 pxq " 0.018e 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 pxq plays the role of a realization of a very narrow spectrum.In this setting the "homogeneous background solution" upx, tq is the solution of and the inhomogeneously perturbed solution V px, tq is the solution of The inhomogeneity δpx, tq is thus defined as δpx, tq " vpx, tq ´upx, tq, and the main question is whether δpx, tq 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 , δ 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 δpx, 0q has to be numerically zero at the ends of the domain, up to satisfactory accuracy.56), on the interval with length L " 10L 0 .Left: for t " 0. Right: for t " T.  57), on the interval with length L " 10L 0 .Left: for t " 0. Right: for t " T.
By taking L " N L 0 , for N P t1, 2, 3, 10u 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 ¨10 ´3 and the space mesh size is dx " 4 ¨10 ´3.The modulus |δpx, tq| 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 8 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.

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 " Opλ 2 0 q was found to be necessary for the MI to be present even for plane waves.If λ 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 2dimensional 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.

Figure 2 :
Figure 2: Space-time plot of the modulus of the inhomogeneity, |δ 2 px, tq|, when propagated on intervals of different length L according to equations (51) -(53).Growth of the inhomogeneity to Op1q 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 .

Figure 4 :
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 :
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.