2.1. Equations of motion

We consider nonlinear hydroelastic waves propagating along an elastic sheet (e.g. ocean waves in sea ice) on top of a two-dimensional ideal fluid of infinite depth. In the thin-plate approximation, the ice sheet is assumed to coincide with the fluid surface $\{ y = \eta (x, t)\}$ and to bend in unison with it. The fluid domain is then given by $S(\eta ) := \{ (x,y): ~x\in \mathbb {R}, -\infty < y <\eta (x,t) \}.$ Assuming that the fluid is incompressible, inviscid and irrotational, it is described by a potential flow such that the velocity field $\boldsymbol {u}(x,y,t) = \boldsymbol {\nabla } \varphi$ satisfies

(2.1)\begin{equation} \nabla^2 \varphi = 0 , \end{equation}

in the fluid domain $S(\eta )$, where the symbol $\boldsymbol {\nabla }$ denotes the spatial gradient $(\partial _x, \partial _y)$.

On the surface $y = \eta (x,t)$, we impose the standard kinematic condition

(2.2)\begin{equation} \partial_t \eta = \partial_y \varphi - (\partial_x \eta) (\partial_x \varphi), \end{equation}

as well as the dynamic boundary condition

(2.3)\begin{equation} \partial_t \varphi ={-}g\eta - \frac{1}{2} \left(\partial_x^2 \varphi + \partial_y^2 \varphi \right) - \mathcal{D} \left( \partial_{s}^2 \kappa + \frac{1}{2} \kappa^3 \right) - \mathcal{P} \kappa . \end{equation}

This condition is derived from the Bernoulli balance of forces following Plotnikov & Toland (Reference Plotnikov and Toland2011) for the interaction between an elastic ice sheet and water beneath it, while the compression term is proportional to the curvature of the fluid–ice interface. Here, $g$ is the acceleration due to gravity, $\mathcal {D}$ is the coefficient of flexural rigidity, $\mathcal {P}$ is the coefficient of ice compression, $\kappa$ is the curvature of the fluid–ice interface caused by the plate deflection and $s$ is the arclength along this interface. The curvature $\kappa$, in terms of $\eta$, is given by

(2.4)\begin{equation} \kappa = \frac{\partial_x^2 \eta}{\left( 1 + (\partial_x \eta)^2 \right)^{3/2}} , \end{equation}

so that the part in (2.3) that describes the bending of the ice sheet becomes

(2.5) \begin{align} \partial_s^2 \kappa + \frac{1}{2} \kappa^3 &= \frac{1}{\sqrt{1+(\partial_x \eta)^2}} \partial_x \left( \frac{1}{\sqrt{1+(\partial_x \eta)^2}} \partial_x \left( \frac{\partial_x^2 \eta}{\left( 1 + (\partial_x \eta)^2 \right)^{3/2}} \right) \right)\nonumber\\ &\quad + \frac{1}{2} \left( \frac{\partial_x^2 \eta}{\left( 1 + (\partial_x \eta)^2 \right)^{3/2}} \right)^3 . \end{align}

The boundary condition at the bottom reads

(2.6)\begin{equation} |\boldsymbol{\nabla} \varphi| \to 0, \quad \text{as} \ y \to -\infty . \end{equation}

2.2. Hamiltonian formulation

Following Zakharov (Reference Zakharov1968) and Craig & Sulem (Reference Craig and Sulem1993), the system of (2.2)–(2.6) has a Hamiltonian formulation in terms of the variables $(\eta, \xi )$, where $\xi (x,t) = \varphi (x, \eta (x,t), t)$ denotes the trace of the fluid velocity potential evaluated at the surface (Guyenne & Părău Reference Guyenne and Părău2012).

Indeed, introducing the DNO for the fluid domain, which associates with the Dirichlet data $\xi$ on $y=\eta (x,t)$ the normal derivative of $\varphi$ at the surface with a normalizing factor, namely

(2.7)\begin{equation} G(\eta):\left. \xi \mapsto \sqrt{1+(\partial_x \eta)^2} \partial_{n} \varphi \right|_{y=\eta}, \end{equation}

the equations of motion (2.2)–(2.3) are equivalent to the following canonical Hamiltonian system:

(2.8)\begin{equation} \partial_t \begin{pmatrix} \eta \\ \xi \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} \partial_\eta H \\ \partial_\xi H \end{pmatrix}, \end{equation}

where the Hamiltonian $H(\eta, \xi )$ is the total energy

(2.9) \begin{align} H (\eta, \xi) &= \frac{1}{2} \int_\mathbb{R} \xi G(\eta) \xi \, {{\rm d}\kern0.07em x} + \frac{1}{2} \int_\mathbb{R} \left( g\eta^2 + \mathcal{D} \frac{(\partial_x^2 \eta)^2}{\left( 1+ (\partial_x \eta)^2 \right)^{5/2}}\right.\nonumber\\ &\quad -\left. \vphantom{\frac{(\partial_x^2 \eta)^2}{\left( 1+ (\partial_x \eta)^2 \right)^{5/2}}}2\mathcal{P} \left( \sqrt{1+(\partial_x \eta)^2} -1 \right) \right) \,{{\rm d}\kern0.07em x}. \end{align}

The first integral is the kinetic energy, while the second integral is the potential energy due to gravity, bending (or rigidity) and compression. The Hamiltonian $H$, together with the momentum

(2.10)\begin{equation} I = \int_\mathbb{R} \eta (\partial_x \xi)\, {{\rm d}\kern0.07em x}, \end{equation}

and the volume

(2.11)\begin{equation} V = \int_\mathbb{R} \eta \, {{\rm d}\kern0.07em x}, \end{equation}

are invariants of motion.

In Fourier variables, the system (2.8) preserves its canonical form. Indeed, denoting the Fourier transform of $f(x)$ by $\widehat {f_k} = (1/\sqrt {2{\rm \pi} }) \int _\mathbb {R} {\rm e}^{-{\rm i}\,kx} f(x) \,{{\rm d}\kern0.07em x}$, we have

(2.12)\begin{equation} \partial_t \begin{pmatrix} \hat \eta_{{-}k} \\ \hat \xi_{{-}k} \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} \partial_{\hat \eta_k} H \\ \partial_{\hat \xi_k} H \end{pmatrix}, \end{equation}

where we have used that $(\hat \eta _{-k}, \hat \xi _{-k}) = (\bar {\hat {\eta }}_k, \bar {\hat {\xi }}_k)$ since $\eta (x)$ and $\xi (x)$ are real-valued functions. Due to the conservation of volume (2.11), we can choose $\hat \eta _0 = 0.$ In the following, we drop the hats when denoting Fourier modes if there is no confusion.

2.3. Non-dimensionalization

The ice parameters are defined by

(2.13a–c)\begin{equation} \mathcal{D} = \frac{\sigma}{\rho}, \quad \sigma = \frac{Eh^3}{12(1-\nu^2)}, \quad \mathcal{P} = \frac{Ph}{\rho}, \end{equation}

where $\rho$ is the density of the underlying fluid, $h$ is the thickness of the ice sheet, $\nu$ is Poisson's ratio for ice, $E$ is Young's modulus and $P$ is the compressive stress. Numerical values of these physical parameters are listed in table 1 of Părău & Dias (Reference Părău and Dias2002) for two sets of experimental data. Introducing the characteristic length and time scales

(2.14a,b)\begin{equation} \ell = \left(\frac{\sigma}{\rho g} \right)^{1/4} = \left( \frac{\mathcal{D}}{g} \right)^{1/4}, \quad \tau = \left(\frac{\sigma}{\rho g^5} \right)^{1/8} = \left( \frac{\mathcal{D}}{g^5} \right)^{1/8}, \end{equation}

respectively (Guyenne & Părău Reference Guyenne and Părău2012; Xu & Guyenne Reference Xu and Guyenne2023), the dimensionless equations of motion follow (2.8) with Hamiltonian (2.9) modulo the change

(2.15)\begin{equation} \mathcal{P} \to \frac{Ph}{\sqrt{\sigma \rho g}} = \frac{\mathcal{P}}{\sqrt{g \mathcal{D}}}, \quad \mathcal{D} \to 1, \quad g \to 1, \end{equation}

where the first dimensionless parameter measures the relative importance of compression to gravity and rigidity. In typical physical situations, it is of order 1 (see for example Liu & Mollo-Christensen (Reference Liu and Mollo-Christensen1988) for realistic values of physical parameters and a discussion by Slunyaev & Stepanyants (Reference Slunyaev and Stepanyants2022)).

For convenience, we will use the same notations hereafter but it is understood that all variables and parameters are now dimensionless according to this choice of non-dimensionalization.

2.4. Taylor expansion of the Hamiltonian near equilibrium

It is known that the DNO is analytic in $\eta$ (Coifman & Meyer Reference Coifman and Meyer1985), and admits a convergent Taylor series expansion

(2.16)\begin{equation} G(\eta) = \sum_{m=0}^\infty G^{(m)}(\eta), \end{equation}

about $\eta =0$. For each $m$, the term $G^{(m)}(\eta )$ is homogeneous of degree $m$ in $\eta$, and can be calculated explicitly via recursive relations (Craig & Sulem Reference Craig and Sulem1993). Denoting $D = -{\rm i} \, \partial _x$, the first three terms are

(2.17)\begin{equation} \left. \begin{gathered} G^{(0)}(\eta) = |D|, \\ G^{(1)}(\eta) = D \eta D - G^{(0)}\eta G^{(0)}, \\ G^{(2)}(\eta) ={-}\dfrac{1}{2} \left( |D|^2 \eta^2 G^{(0)}+ G^{(0)} \eta^2 |D|^2 - 2G^{(0)} \eta G^{(0)} \eta G^{(0)} \right). \end{gathered} \right\} \end{equation}

This in turn provides an expansion of the Hamiltonian near the stationary solution $(\eta,\xi ) = (0,0)$,

(2.18)$$\begin{align} H (\eta, \xi) &= \frac{1}{2} \int_\mathbb{R} \left[\vphantom{\left. \xi G^{(2)} (\eta) \xi - \frac{5}{2} \mathcal{D} (\partial_x \eta)^2 (\partial_x^2 \eta)^2 + \frac{1}{4} \mathcal{P} (\partial_x \eta)^4 + \cdots \right]} \xi G^{(0)}(\eta) \xi + g\eta^2 + \mathcal{D} (\partial_x^2 \eta)^2 - \mathcal{P} (\partial_x \eta)^2 + \xi G^{(1)}(\eta) \xi\right.\nonumber\\ &\quad +\left. \xi G^{(2)} (\eta) \xi - \frac{5}{2} \mathcal{D} (\partial_x \eta)^2 (\partial_x^2 \eta)^2 + \frac{1}{4} \mathcal{P} (\partial_x \eta)^4 + \cdots \right]\, {{\rm d}\kern0.07em x}. \end{align}$$

In Fourier variables, the above expansion can be written as

(2.19)\begin{equation} H = H^{(2)} + H^{(3)}+ H^{(4)} + \cdots \end{equation}

where each term $H^{(m)}$ is homogeneous of degree $m$ in $\eta$ and $\xi$. In particular, we have

(2.20)\begin{equation} \left. \begin{aligned} H^{(2)} & = \frac{1}{2} \int \left( |k| |\xi_k|^2 + (g - \mathcal{P} k^2 + \mathcal{D} k^4) |\eta_k|^2 \right) \,{\rm d}k,\\ H^{(3)} & ={-} \frac{1}{2\sqrt{2{\rm \pi}}} \int (k_1 k_3 + |k_1| |k_3|) \xi_1 \eta_2 \xi_3 \delta_{123} \,{\rm d}k_{123} \end{aligned} \right\} \end{equation}

and

(2.21)\begin{align} H^{(4)}& ={-} \frac{1}{8{\rm \pi}} \int |k_1| |k_4| (|k_1| + |k_4| - 2 |k_3+k_4|) \xi_1 \eta_2 \eta_3 \xi_4 \delta_{1234} \,{\rm d}k_{1234}\nonumber\\ &\quad + \frac{1}{4{\rm \pi}} \int \left( \frac{5 \mathcal{D}}{2} k_1^2 k_2^2 k_3 k_4 + \frac{\mathcal{P}}{4} k_1 k_2 k_3 k_4 \right) \eta_1 \eta_2 \eta_3 \eta_4 \delta_{1234} \,{\rm d}k_{1234}, \end{align}

where we have used the compact notations $(\eta _j, \xi _j)= (\eta _{k_j}, \xi _{k_j})$, $\textrm {d}k_{1 \ldots n} = \textrm {d}k_1 \ldots \textrm {d}k_n$ and $\delta _{1 \ldots n} = \delta (k_1+\cdots +k_n)$, where $\delta ({\cdot })$ is the Dirac distribution $\delta (k) = (1/2{\rm \pi} ) \int _\mathbb {R} \textrm {e}^{-\textrm {i}\,kx}\, {\textrm {d} x}$. The domain of integration is now omitted from all integrals and is understood to be $\mathbb {R}$ for each $x_j$ or $k_j$. Notice that the ice parameters $\mathcal {P}$ and $\mathcal {D}$ appear explicitly in the expressions for $H^{(2)}$ and $H^{(4)}$, but not in $H^{(3)}$.

2.5. Dispersion relation

The linearized system for (2.8) around $(\eta, \xi )=(0,0)$ is

(2.22)\begin{equation} \left. \begin{gathered} \partial_t \eta = |D| \xi,\\ \partial_t \xi ={-}g\eta - \mathcal{P} \partial_x^2 \eta - \mathcal{D} \partial_x^4 \eta. \end{gathered} \right\} \end{equation}

For

(2.23)\begin{equation} \omega^2(k) = \omega_k^2 := |k| \left( g - \mathcal{P} k^2 + \mathcal{D} k^4 \right),\end{equation}

equation (2.22) admits periodic plane-wave solutions

(2.24a,b)\begin{equation} \eta(x,t) \propto {\rm e}^{{\rm i}(kx-\omega_k t)} \quad \text{and} \quad \xi(x,t) \propto {\rm e}^{{\rm i}(kx-\omega_k t)}. \end{equation}

Equation (2.23) is known as the linear dispersion relation. By construction, for a particular choice of $\mathcal {P}, \mathcal {D}$ and $k$ the value of $\omega _k^2$ may be either positive or negative. When $\omega _k^2$ is positive, the solutions in (2.24a,b) are travelling waves since $\omega _k \in \mathbb {R}$. On the other hand, when $\omega _k^2<0$, one has $\omega _k \in \textrm {i} \, \mathbb {R}$, and therefore, the solutions in (2.24a,b) represent either evanescent or growing waves due to the factor $\textrm {e}^{\pm |\omega _k| t}$ in the expressions of $\eta$ and $\xi$.

For parameters $\mathcal {P}, \mathcal {D}$ and $k$ with $\omega _k^2>0$, it is convenient to introduce the complex symplectic coordinates which diagonalize the quadratic Hamiltonian $H^{(2)}$ in (2.20). In the Fourier space, these are given by

(2.25)\begin{equation} \begin{pmatrix} z_k \\ \bar z_{{-}k} \end{pmatrix} = \dfrac{1}{\sqrt 2} \begin{pmatrix} a_k & {\rm i} \, a_k^{{-}1} \\ a_k & -{\rm i} \, a_k^{{-}1} \end{pmatrix} \begin{pmatrix} \eta_k \\ \xi_k \end{pmatrix}, \end{equation}

where $a_k^2 := \omega _k/|k|$. The quadratic term $H^{(2)}$ becomes

(2.26)\begin{equation} H^{(2)} = \int \omega_k |z_k|^2 \,{\rm d}k, \end{equation}

while the cubic term $H^{(3)}$ takes the form

(2.27) \begin{equation} H^{(3)} = \frac{1}{8\sqrt{\rm \pi}} \int ( { k}_1 { k}_3 +|k_1| |k_3|) \frac{a_1 a_3}{a_2} (z_1-\bar z_{{-}1})(z_2+\bar z_{{-}2})(z_3-\bar z_{{-}3}) \delta_{123} \,{\rm d}k_{123}, \end{equation}

where $z_{\pm j} := z_{\pm { k}_j}$ and $a_j := a_{{ k}_j}$. Moreover, since the mapping (2.25) is canonical (Guyenne et al. Reference Guyenne, Kairzhan and Sulem2022a), the full Hamiltonian system (2.12) becomes

(2.28)\begin{equation} \partial_t \begin{pmatrix} z_k \\ \bar z_{{-}k} \end{pmatrix} = \begin{pmatrix} 0 & -{\rm i} \\ {\rm i} & 0 \end{pmatrix} \begin{pmatrix} \partial_{z_k} H \\ \partial_{\bar z_{{-}k}} H \end{pmatrix}. \end{equation}

We note that the expressions (2.26) and (2.27) coincide in structure with the corresponding formulae for the quadratic and cubic Hamiltonians in the case of deep-water surface waves (Craig & Sulem Reference Craig and Sulem2016). This is because the parameters of ice rigidity and ice compression, appearing in our problem, are hidden in the definitions of $\omega _k$ and $a_k$.

One finds the values of $\mathcal {P}$ when $\omega _k^2$ is positive for any $k$ by analysing the quartic inequality

(2.29)\begin{equation} \mathcal{D} k^4 - \mathcal{P} k^2 + g > 0. \end{equation}

Recall that $g = 1$ and $\mathcal {D} = 1$ in dimensionless units. For a given $\mathcal {P}$, one has $\omega _k^2 > 0$ for small or large $k$. When $k$ takes moderate values, $\omega _k^2$ could be negative. One can find the zeros of $\omega _k^2$ analytically but it is more convenient to plot the dispersion relation $\omega _k^2$ as a function of $k$ for various values of the parameter $\mathcal {P}$ as shown in figure 1. Note that $\omega _k^2$ remains always positive whenever $\mathcal {P} < 2$. When $\mathcal {P} = 2$, $\omega _k^2 = 0$ for $k = 1$. A detailed analysis of the linear dispersion relation in terms of various physical parameters including ice bending and compression can also be found in Schulkes, Hosking & Sneyd (Reference Schulkes, Hosking and Sneyd1987).

Figure 1. Linear dispersion relation $\omega ^2(k)$ as a function of $k$ for $\mathcal {P} = 0.1$ (blue), $\mathcal {P} = 1$ (red), $\mathcal {P} = 2$ (magenta), $\mathcal {P} = 5$ (green).

We are interested in having $\omega _k^2 > 0$ for all real values of $k$ since we are looking for travelling wave solutions, not evanescent or growing waves. Thus, throughout the paper, we make the assumption

(2.30)\begin{equation} 0 \le \mathcal{P} < 2. \end{equation}

3.1. Existence of resonant triads

We now address the question of existence of non-zero resonant triads $(k_1, k_2,k_3)$ satisfying

(3.1)\begin{equation} k_1 + k_2 + k_3 = 0, \end{equation}

and at least one of the equations

(3.2)\begin{equation} \omega_1 \pm \omega_2 \pm \omega_3 = 0, \end{equation}

where $\omega _j$ stands for $\omega _{k_j}$. To find all the possible triads satisfying (3.1)–(3.2), we consider the following function:

(3.3)\begin{align} d_{123} = d(k_1,k_2,k_3) := (\omega_1 + \omega_2 + \omega_3)(\omega_1 + \omega_2 - \omega_3)(\omega_1 - \omega_2 + \omega_3)(\omega_1 - \omega_2 - \omega_3). \end{align}

From the construction of $d_{123}$, we have that $(k_1,k_2,k_3)$ is a resonant triad if and only if it solves the system

(3.4) \begin{equation} \left. \begin{gathered} k_1+k_2+k_3 = 0, \\ d(k_1,k_2,k_3) = 0, \\ k_j \neq 0, \quad j = \{ 1,2,3 \}. \end{gathered} \right\} \end{equation}

The analogue of the system (3.4) was studied in detail by Craig & Sulem (Reference Craig and Sulem2016) for surface gravity water waves. In that case, due to the absence of physical parameters $\mathcal {P}$ and $\mathcal {D}$, the expression for $d_{123}$ is much simpler and the computations can be performed by hand. However, in the present case, it is difficult to solve (3.4) explicitly due to the complexity of the dispersion relation (2.23) and the high degree of the polynomial function $d_{123}$.

Furthermore, as shown in Craig & Sulem (Reference Craig and Sulem2016), the function $d_{123}$ appears together with the constraint

(3.5)\begin{equation} k_1 k_3+|k_1|\,|k_3| \neq 0. \end{equation}

Therefore, we provide estimates for $d_{123}$ under the constraint (3.5) and the assumption $k_1+k_2+k_3=0$.

Lemma 3.1 Assuming $\mathrm{sgn} (k_1) = \mathrm{sgn} (k_3)$ and $k_1+k_2+k_3=0$, we have

(3.6)\begin{equation} d_{123} = k_1k_3 \tilde{d}(k_1, k_3), \end{equation}

where

(3.7)$$\begin{align} \tilde{d} (k_1, k_3) &:= k_1 k_3 (k_1 + k_3)^2 \left( 3\mathcal{P} -5 \mathcal{D} (k_1^2+k_1k_3+k_3^2) \right)^2\nonumber\\ &\quad - 4 (g-\mathcal{P} k_1^2 + \mathcal{D} k_1^4)(g-\mathcal{P} k_3^2 + \mathcal{D} k_3^4). \end{align}$$

The function $\tilde {d}$ has the following symmetry properties:

(3.8a,b)\begin{equation} \tilde{d} (k_1, k_3) = \tilde{d} (k_3, k_1) \quad \text{and} \quad \tilde{d} (k_1, k_3) = \tilde{d} ({-}k_1, -k_3). \end{equation}

Proof. The results are based on direct computations under the above assumptions. Below, we provide a short overview of steps.

First, rewrite $d_{123}$ as

(3.9)\begin{equation} d_{123} = (\omega_1^2 + \omega_3^2 - \omega_2^2)^2 - 4 \omega_1^2 \omega_3^2. \end{equation}

Observe that the second term on the right-hand side of (3.9) identifies with the second line in (3.7). For the first term on the right-hand side of (3.9), we write

(3.10) \begin{align} \omega_1^2 + \omega_3^2 - \omega_2^2 &= g (|k_1|+ |k_3|-|k_2|) - \mathcal{P} (|k_1|^3+ |k_3|^3-|k_2|^3)\nonumber\\ &\quad + \mathcal{D} (|k_1|^5+ |k_3|^5-|k_2|^5). \end{align}

Given the assumptions, the first term vanishes since

(3.11)\begin{equation} |k_1|+ |k_3|-|k_2| = 0, \end{equation}

and for the cubic terms, we have

(3.12)\begin{equation} |k_1|^3+ |k_3|^3-|k_2|^3 = |k_2| (k_1^2 - k_1k_3 + k_3^2) - |k_2|^3 ={-}3 |k_2| k_1k_3. \end{equation}

Similarly, for fifth-power terms, we have

(3.13)\begin{equation} |k_1|^5+ |k_3|^5-|k_2|^5 ={-}5|k_2| k_1k_3 (k_1^2 + k_1k_3 + k_3^2). \end{equation}

Substituting these expressions back into (3.9), we obtain (3.6). The properties (3.8a,b) are straightforward to verify.

Based on Lemma 3.1, finding resonant triads $(k_1, k_2, k_3)$ in (3.4) under the constraint (3.5) is equivalent to finding the non-zero roots of the equation

(3.14)\begin{equation} \tilde{d} (k_1, k_3) = 0, \quad \text{with} \ {\rm sgn} (k_1) = {\rm sgn} (k_3), \end{equation}

where $k_2 = -k_1-k_3$. Solving (3.14) explicitly seems to be impossible due to the complicated structure of $\tilde {d}(k_1,k_3)$ in (3.7). Therefore, we find solutions to (3.14) numerically. Since $\tilde {d}$ is invariant under sign change of variables, $\tilde {d} (k_1, k_3) = \tilde {d} (-k_1, -k_3)$, we only concentrate on the positive roots of (3.14) and the negative roots are found by symmetry. Setting $\mathcal {P} = 1$ for simplicity, the solution curve $\mathcal {C}^+$ for positive roots of (3.14) is displayed in figure 2. A similar solution curve is present for other values of $\mathcal {P} < 2$. We note that the curve $\mathcal {C}^+$ is unbounded, has the $k_1$- and $k_3$-axes as asymptotes and never intersects these axes. The negative roots of (3.14) are given by the solution curve $\mathcal {C}^-$, which is obtained from $\mathcal {C}^+$ by symmetry with respect to the origin. Combining these curves, all non-zero roots of (3.14) are given by

(3.15)\begin{equation} \mathcal{C} := \mathcal{C}^+\cup \mathcal{C}^-. \end{equation}

Figure 2. Solution curve $\mathcal {C}^+$ (blue curve) for positive roots $(k_1, k_3)$ of (3.14) and its neighbourhood $\mathcal {C}_\mu ^+$ (grey area) for $\mathcal {P} = 1$.

In our analysis later, various integrals involve expressions where the function $\tilde {d}(k_1, k_3)$ appears in the denominator. To avoid this ‘small divisors’ issue, we will restrict the integral to regions away from the solution curve $\mathcal {C}$. To do so, we construct a neighbourhood of $\mathcal {C}$, referred to as $\mathcal {C}_\mu$, for sufficiently small $\mu$, such that

(3.16)\begin{equation} |\tilde{d}(x, y)| \geq \text{const.}_\mu, \quad \text{for} \ \text{all } (x, y) \in \mathbb{R}^2 \backslash \mathcal{C}_\mu, \end{equation}

where $\textrm {const.}_\mu$ is a constant depending on $\mu$ only. More precisely, at every point $(k_1, k_3) \in \mathcal {C}^+$, denote a normal vector to $\mathcal {C}^+$ by $\boldsymbol {n} (k_1, k_3)$. Then, consider the set of points in the $(k_1,k_3)$-plane given by

(3.17)\begin{equation} \left\{ (k_1, k_3) + \mu \, t \, \boldsymbol{n} (k_1,k_3): t \in [{-}1,1] \text{ and } (k_1, k_3) \in \mathcal{C}^+ \right\}, \end{equation}

which geometrically looks like a ‘tube’ around the curve $\mathcal {C}^+$ with a width $\mu$. We then construct a neighbourhood of $\mathcal {C}^+$, denoted by $\mathcal {C}_\mu ^+$, as

(3.18)\begin{equation} \mathcal{C}_\mu^+ := \{ \text{the set } (3.17)\} \cap \{(x,y): x\geq 0, y\geq 0\}. \end{equation}

A sketch of $\mathcal {C}_\mu ^+$ is shown in figure 2 for $\mathcal {P} = 1$. The neighbourhood $\mathcal {C}_\mu ^-$ for the curve $\mathcal {C}^-$ is constructed similarly, and can be seen as the reflection of $\mathcal {C}_\mu ^+$ with respect to the origin. Combining these sets, we get the neighbourhood of the set $\mathcal {C}$ in (3.15) as $\mathcal {C}_\mu := \mathcal {C}_\mu ^+ \cup \mathcal {C}_\mu ^-$.

We now provide formal arguments to show that (3.16) is satisfied for the set $\mathcal {C}_\mu$ constructed above. It suffices to show that (3.16) is satisfied on the boundary of $\mathcal {C}_\mu ^+$ for small values of $k_1$, which corresponds to the far end of the neighbourhood in figure 2 located at the top left-hand corner. Then, the validity of (3.16) for all $(x, y)$ outside $\mathcal {C}_\mu$ follows by continuity arguments and symmetry.

First, we note that, under assumption (2.30), the points $(k_1, k_3) \in \mathcal {C}$ satisfy $k_3 \approx (4g/(25 \mathcal {D} k_1))^{1/3}$ whenever $k_1 \ll 1$. Indeed, from (3.7), it is clear that if $k_1 \ll 1$ then $k_3 \gg 1$. As a result, for $k_1 \ll 1$, we have

(3.19)\begin{equation} 0 = \tilde{d}(k_1, k_3) \approx 25 \mathcal{D}^2 k_1 k_3^7 - 4 g \mathcal{D} k_3^4 \implies k_3 \approx \left(\frac{4g}{25 \mathcal{D} k_1} \right)^{1/3}. \end{equation}

Therefore, the normal vector of $\mathcal {C}^+$ for small values of $k_1$ is approximately given by

(3.20)\begin{equation} \boldsymbol{n} (k_1,k_3) = \frac{\left( k_1, 3k_1^{4/3} \left( \dfrac{25\mathcal{D}}{4g} \right)^{1/3} \right)}{\left\| k_1, 3k_1^{4/3} \left( \dfrac{25\mathcal{D}}{4g} \right)^{1/3} \right\|} \approx (1, 0). \end{equation}

The points on the boundary of $\mathcal {C}_\mu ^+$ with $k_1 \ll 1$ are

(3.21)\begin{equation} (0, k_3) \quad \text{and} \quad (k_1+ \mu, k_3), \end{equation}

where $k_3$ satisfies the estimate (3.19). Applying simple algebraic steps to (3.7), it can be shown that the value of $\tilde {d} (0, k_3)$ under the assumption (2.30) is bounded away from zero by

(3.22)\begin{equation} \left| \tilde{d} (0, k_3) \right| \geq \left| \tilde{d} (0, \sqrt{\mathcal{P}/(2\mathcal{D})}) \right| = \frac{g}{\mathcal{D}} (4g\mathcal{D} - \mathcal{P}^2) >0. \end{equation}

To estimate the value of $\tilde {d}(k_1 + \mu, k_3)$ we use the Taylor expansion

(3.23)\begin{equation} \tilde{d}(k_1 + \mu, k_3) \approx \tilde{d}(k_1, k_3) + \mu \, \partial_{k_1} \tilde{d}(k_1, k_3) = \mu \, \partial_{k_1} \tilde{d}(k_1, k_3). \end{equation}

Taking the derivative of (3.7) with respect to $k_1$, we have

(3.24)\begin{align} \partial_{k_1} \tilde{d} (k_1, k_3) &= k_3 (3k_1+k_3) (k_1 + k_3) \left( 3\mathcal{P} -5 \mathcal{D} (k_1^2+k_1k_3+k_3^2) \right)^2\nonumber\\ &\quad - 10 k_1 k_3 (k_1 + k_3)^2 \left( 3\mathcal{P} -5 \mathcal{D} (k_1^2+k_1k_3+k_3^2) \right) \mathcal{D} \left( 2 k_1 + k_3 \right)\nonumber\\ &\quad + 4 (2 \mathcal{P} k_1 - 4 \mathcal{D} k_1^3)(g-\mathcal{P} k_3^2 + \mathcal{D} k_3^4), \end{align}

which, under (3.19), is approximated by

(3.25)\begin{equation} \partial_{k_1} \tilde{d} (k_1, k_3) \approx 25 \mathcal{D}^2 k_3^7, \end{equation}

whenever $k_1 \ll 1$. As a result, the estimate in (3.23) takes the form

(3.26)\begin{equation} \tilde{d}(k_1 + \mu, k_3) \approx 25 \mu \mathcal{D}^2 k_3^7, \end{equation}

which is bounded away from zero as desired. Later in the analysis, the function $d_{123}$ will appear in the denominator and we will set $\mu = \varepsilon$, where $\varepsilon$ is the small parameter in the modulational regime.

In view of (3.4) and (3.6), the set of resonant triads satisfying (3.1)–(3.2) is given by

(3.27)\begin{equation} \mathcal{R} := \{(k_1, -k_1-k_3, k_3) \in \mathbb{R}^3: (k_1, k_3) \in \mathcal{C} \}. \end{equation}

The neighbourhood of $\mathcal {R}$ is defined by

(3.28)\begin{equation} \mathcal{N}_\mu := \{ (k_1, -k_1-k_3, k_3) \in \mathbb{R}^3: (k_1, k_3) \in \mathcal{C}_\mu\}, \end{equation}

as illustrated in figure 2. We also define a characteristic function $\chi _{\mathcal {N}_\mu } (k_1, k_2, k_3)$ supported in the neighbourhood $\mathcal {N}_\mu$ as

(3.29)\begin{equation} \chi_{\mathcal{N}_\mu} (k_1, k_2, k_3) := \left\{ \begin{array}{@{}ll@{}} 1, & \text{if } (k_1, k_2, k_3) \in \mathcal{N}_\mu, \\ 0, & \text{otherwise}. \end{array} \right. \end{equation}

3.2. Phase and group speeds

Assuming $k > 0$ without loss of generality (since $\omega _k$ is an even function of $k$), the phase speed associated with the linear dispersion relation (2.23) is given by

(3.30)\begin{equation} c(k) = \frac{\omega_k}{k} = \sqrt{\frac{g - {\mathcal{P}} k^2 + {\mathcal{D}} k^4}{k}}, \end{equation}

while the group speed reads

(3.31)\begin{equation} c_g(k) = \partial_k \omega_k = \frac{g - 3 {\mathcal{P}} k^2 + 5 {\mathcal{D}} k^4}{2 \omega_k}. \end{equation}

It can be easily shown that, if the phase speed $c$ has a local minimum at $k = k_{min}$, then the phase and group speeds coincide at this minimum. The equation $c = c_g$ reduces to $g + {\mathcal {P}} k^2 - 3 {\mathcal {D}} k^4 = 0$ which is quadratic in $k^2$. The only possible solution that is real and positive takes the form

(3.32)\begin{equation} k_{min} = \sqrt{\frac{{\mathcal{P}} + \sqrt{{\mathcal{P}}^2 + 12 g {\mathcal{D}}}}{6 {\mathcal{D}}}}. \end{equation}

In the flexural-gravity case (${\mathcal {P}} = 0$), this solution yields the well-known value $k_{min} = (g/(3 {\mathcal {D}}))^{1/4}$ (Guyenne & Părău Reference Guyenne and Părău2012). Figure 3 shows the phase and group speeds for values $\mathcal {P} = \{ 1, 2, 5 \}$.

Figure 3. Phase speed $c(k)$ (blue) and group speed $c_g(k)$ (red) as functions of $k$ for (a) $\mathcal {P} = 1$, (b) $\mathcal {P} = 2$, (c) $\mathcal {P} = 5$.

4.1. Canonical transformation

We first construct a transformation that eliminates non-resonant terms from the cubic Hamiltonian (2.27). More precisely, we are looking for a canonical transformation of the physical variables $(\eta,\xi )$

(4.3)\begin{equation} \tau: w = \begin{pmatrix} \eta \\ \xi \end{pmatrix} \longmapsto w' = \begin{pmatrix} \eta' \\ \xi' \end{pmatrix}, \end{equation}

defined in a neighbourhood of the origin, such that the transformed Hamiltonian satisfies

(4.4a,b)\begin{equation} H'(w') = H(\tau^{{-}1}(w')), \quad \partial_t w' = J \, \boldsymbol{\nabla} H'(w'), \end{equation}

and reduces to

(4.5)\begin{equation} H'(w') = H^{(2)}(w') + Z^{(3)} + Z^{(4)} + \cdots + Z^{(m )} + R^{(m+1)}, \end{equation}

where $Z^{(j )}$ only consists of resonant terms of order $j$ and $R^{(m+1)}$ is the remainder term at order $m$ (Craig & Sulem Reference Craig and Sulem2016; Craig, Guyenne & Sulem Reference Craig, Guyenne and Sulem2021b). The transformation $\tau$ is obtained as a Hamiltonian flow $\psi$ from ‘time’ $s=-1$ to ‘time’ $s=0$ governed by

(4.6)\begin{equation} \partial_s \psi = J \, \boldsymbol{\nabla} K(\psi), \quad \psi(w')|_{s=0} = w' , \quad \psi(w')|_{s={-}1} = w , \end{equation}

and associated with an auxiliary Hamiltonian $K$. Such a transformation is canonical and preserves the Hamiltonian structure of the system. The Hamiltonian $H'$ satisfies $H'(w') = H(\psi (w'))|_{s=-1}$ and its Taylor expansion around $s=0$ is

(4.7)\begin{equation} H'(w') = H(\psi(w'))|_{s=0} - \frac{{\rm d} H}{{\rm d}s}(\psi(w'))|_{s=0} + \frac{1}{2} \frac{{\rm d}^2 H}{{\rm d}s^2}(\psi(w'))|_{s=0} - \cdots. \end{equation}

Abusing notations, we further drop the primes and use $w = (\eta,\xi )^\top$ to denote the new variable $w'$. Terms in this expansion can be expressed using Poisson brackets as

(4.8)\begin{align} H(\psi(w))|_{s=0} & = H(w), \nonumber\\ \frac{{\rm d} H}{{\rm d}s}(\psi(w))|_{s=0} & = \int \left( \partial_\eta H \partial_s \eta + \partial_\xi H \partial_s \xi \right) \,{{\rm d}\kern0.07em x},\nonumber\\ & = \int \left( \partial_\eta H \partial_\xi K - \partial_\xi H \partial_\eta K \right) \,{{\rm d}\kern0.07em x} = \{K, H\}(w), \end{align}

and similar expressions can be obtained for higher-order $s$-derivatives. The Taylor expansion of $H'$ around $s=0$ now has the form

(4.9)\begin{equation} H'(w) = H(w) - \{K, H\}(w) + \frac{1}{2} \{K, \{K, H\}\}(w) - \cdots. \end{equation}

Substituting this transformation into the expansion (2.19) of $H$, we obtain

(4.10)\begin{align} H'(w) & = H^{(2)}(w) + H^{(3)}(w) + \cdots \nonumber\\ & \quad - \{K,H^{(2)} \} (w) -\{K, H^{(3)}\} (w) -\{K, H^{(4)}\} (w) - \cdots \nonumber\\ & \quad + \frac{1}{2}\{K,\{K,H^{(2)} \}\} (w) +\frac{1}{2} \{K,\{K,H^{(3)}\}\} (w) + \cdots . \end{align}

If $K$ is homogeneous of degree $m$ and $H^{(n)}$ is homogeneous of degree $n$, then $\{K,H^{(n)} \}$ is of degree $m+n-2$. Thus, if we construct an auxiliary Hamiltonian $K=K^{(3)}$ that is homogeneous of degree $3$ and satisfies the relation

(4.11)\begin{equation} H^{(3)}- \{K^{(3)},H^{(2)} \} = 0, \end{equation}

we will have eliminated all cubic terms from the transformed Hamiltonian $H'$. We can repeat this process at each order.

4.2. Third-order Birkhoff normal form

We recall that the complex symplectic coordinates $z_{j}$ and $\bar {z}_{-j}$ diagonalize the coadjoint operator $\textrm {coad}_{H^{(2)}} := \{{\cdot }, H^{(2)}\}$, that is, the linear operation of taking Poisson brackets with $H^{(2)}$ (Craig & Sulem Reference Craig and Sulem2016; Guyenne et al. Reference Guyenne, Kairzhan and Sulem2022a). When applying the operator to monomial terms such as $\mathcal {I} := \int z_1 z_2 \bar z_{-3} \delta _{123} \,\textrm {d}{k}_{123},$ we have

(4.12)\begin{equation} \{ \mathcal{I}, H^{(2)}\} = {\rm i}\, \int (\omega_1 + \omega_2 - \omega_3) z_1 z_2 \bar z_{{-}3} \delta_{123} \,{\rm d}{k}_{123}. \end{equation}

We use the diagonal property as in (4.12) to find the auxiliary Hamiltonian $K^{(3)}$ from (4.11). The presence of resonant triads (3.27) does not allow us to eliminate $H^{(3)}$ completely by virtue of (4.11). Instead, we are only able to find $K^{(3)}$ such that

(4.13)\begin{equation} \{K^{(3)},H^{(2)} \} = H_{NoRes}^{(3)}, \end{equation}

where $H_{NoRes}^{(3)}$ stands for the non-resonant part of the third-order Hamiltonian. Explicitly, we define the non-resonant part of $H^{(3)}$ as $H_{NoRes}^{(3)} := H^{(3)} - H_{Res}^{(3)}$, where

(4.14)\begin{equation} H_{Res}^{(3)} = \frac{1}{8\sqrt{\rm \pi}} \int \chi_{\mathcal{N}_\mu} (k_1, k_2, k_3) S_{123} (z_1 \bar{z}_{{-}2} z_3 + \bar{z}_{{-}1} {z}_{2} \bar{z}_{{-}3} ) \delta_{123} \,{\rm d}{k}_{123}, \end{equation}

is a resonant part of (2.27),

(4.15)\begin{equation} S_{123} = S(k_1, k_2, k_3) := (k_1 k_3 + |k_1| |k_3|) \frac{a_1 a_3}{a_2}, \end{equation}

and $\chi _{\mathcal {N}_\mu }$ is the characteristic function defined in (3.29).

Proposition 4.1 The cohomological equation (4.13) has a unique solution $K^{(3)}$ which, in complex symplectic coordinates, is

(4.16)$$\begin{align} K^{(3)} &= \frac{1}{8{\rm i}\,\sqrt{\rm \pi}} \int S_{123} \left[ \frac{z_1 z_2 z_3- \bar z_{{-}1} \bar z_{{-}2} \bar z_{{-}3} }{\omega_1 + \omega_2+ \omega_3} - 2 \frac{z_1 z_2 \bar z_{{-}3}- \bar z_{{-}1} \bar z_{{-}2} z_3}{\omega_1+ \omega_2 -\omega_3} \right.\nonumber\\ &\quad + \left.\left(1-\chi_{\mathcal{N}_\mu} (k_1, k_2, k_3) \right) \frac{ z_1\bar z_{{-}2} z_3- \bar z_{{-}1} z_2 \bar z_{{-}3} }{\omega_1 -\omega_2 +\omega_3} \right] \delta_{123} \,{\rm d}{k}_{123}. \end{align}$$

Alternatively, in the variables $(\eta,\xi )$, $K^{(3)}$ has the form

(4.17)\begin{equation} K^{(3)} = \frac{1}{\sqrt{2{\rm \pi}}} \int \left( P_{123} \eta_1 \eta_2 \xi_3 + Q_{123}~ \eta_1 \xi_2 \eta_3 + R_{123} \xi_1 \xi_2 \xi_3 \right) \delta_{123} \,{\rm d}{k}_{123}, \end{equation}

where the denominator $d_{123}$ is given in (3.3) and the coefficients are

(4.18)\begin{align} \left. \begin{aligned} P_{123} & = \frac{1 + {\rm sgn} (k_1) {\rm sgn} (k_3)}{4 \tilde{d} (k_1, k_3)} a_1^2 \left( 4 \omega_1 (\omega_1^2 - \omega_2^2 - \omega_3^2) - \chi_{\mathcal{N}_\mu} (k_1, k_2, k_3) \varPi_{123} \right),\\ Q_{123} & = \frac{1 + {\rm sgn} (k_1) {\rm sgn} (k_3)}{8 \tilde{d} (k_1, k_3)} \left( \frac{a_1^2 a_3^2}{a_2^2} \right) \left( 8 \omega_1 \omega_2 \omega_3 + \chi_{\mathcal{N}_\mu} (k_1, k_2, k_3) \varPi_{123} \right),\\ R_{123} & = \frac{1 + {\rm sgn} (k_1) {\rm sgn} (k_3)}{8 \tilde{d} (k_1, k_3)} \left( \frac{1}{a_2^2} \right) \left( 4 \omega_2 (\omega_1^2 - \omega_2^2 + \omega_3^2) - \chi_{\mathcal{N}_\mu} (k_1, k_2, k_3) \varPi_{123} \right), \end{aligned} \right\} \end{align}

with $\varPi _{123} := (\omega _1 + \omega _2 + \omega _3)(\omega _1 + \omega _2 - \omega _3) (\omega _1 - \omega _2 - \omega _3)$.

We point out that the coefficients $P_{123}, Q_{123}$ and $R_{123}$ are well defined in the neighbourhood of the resonant triads (3.27). This is due to the construction of $K^{(3)}$ in (4.16) which avoids the resonant triads that make the denominator $\omega _1-\omega _2 + \omega _3$ equal to $0$. This is the principal reason why we solve (4.13) to find $K^{(3)}$ rather than the full relation $\{ K^{(3)}, H^{(2)} \} = H^{(3)}$ in (4.11). The latter equation would lead to an ill-defined $K^{(3)}$ with vanishing denominator $\omega _1-\omega _2 + \omega _3$, which we cannot handle.

The third-order normal form defining the new coordinates is obtained as the solution map at $s=0$ of the Hamiltonian flow

(4.19)\begin{equation} \partial_s \begin{pmatrix} \eta \\ \xi \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} \partial_{\eta}K^{(3)} \\ \partial_{\xi} K^{(3)} \end{pmatrix}, \end{equation}

with initial condition at $s= -1$ being the original variables $(\eta,\xi )$. Equivalently, in Fourier coordinates,

(4.20a,b)\begin{equation} \partial_s \eta_{{-}k} = \partial_{\xi_k} K^{(3)}, \quad \partial_s \xi_{{-}k}={-} \partial_{\eta_k} K^{(3)}, \end{equation}

where

(4.21)\begin{equation} \left. \begin{aligned} \partial_{\xi_k} K^{(3)} & = \frac{1}{\sqrt{2{\rm \pi}}} \int \left( \left( P_{12k} + Q_{1k2} \right) \eta_1 \eta_2 + \left(R_{12k} + R_{2k1} + R_{k12} \right) \xi_1 \xi_2 \right) \delta_{12k} \,{\rm d}{k}_{12},\\ \partial_{\eta_k} K^{(3)} & = \frac{1}{\sqrt{2{\rm \pi}}} \int \left( P_{1k2} + P_{k12} +Q_{12k} + Q_{k21} \right) \eta_1 \xi_2 \delta_{12k} \,{\rm d}{k}_{12}, \end{aligned} \right\} \end{equation}

by virtue of (4.17).

5.1. Explicit computations of $T_0$, $T_{Res}$ and $T_{NoRes}$

Here, we provide precise formulae for the coefficients $T_0$, $T_{Res}$ and $T_{NoRes}$ appearing in (5.5)–(5.6).

Proposition 5.1 We have $T_0 = T_0^{(1)} + T_0^{(2)}$, where

(5.8) \begin{equation} \left. \begin{aligned} T_0^{(1)} & ={-}V_{12({-}3)({-}4)} -V_{({-}4)({-}3)21} -V_{1({-}3)2({-}4)} -V_{({-}4)2({-}3)1}\\ &\quad +V_{1({-}4)({-}3)2} +V_{({-}3)21({-}4)},\\ T_0^{(2)} & = \frac{k_1 k_2 k_3 k_4}{32 {\rm \pi}a_1 a_2 a_3 a_4} \left( 5\mathcal{D} (k_1k_2 -k_1k_3 - k_1k_4 - k_2k_3 - k_2k_4 + k_3k_4) - 3 \mathcal{P} \right), \end{aligned} \right\} \end{equation}

with

(5.9)\begin{equation} V_{1234} = \frac{a_1 a_4}{32 {\rm \pi}a_2 a_3} |k_1|\, |k_4| (|k_1|+|k_4| - 2|k_3+k_4|). \end{equation}

Proposition 5.2 Let $B_{123} := S_{123} (1-\chi _\mathcal {N_\mu } (k_1, k_2, k_3))$ with $S_{123}$ and $\chi _{\mathcal {N}_\mu }$ given in (4.15) and (3.29), respectively. Then, $T_{NoRes} = T_{NoRes}^{(1)} + T_{NoRes}^{(2)} + T_{NoRes}^{(3)}$ with

(5.10)\begin{align} T_{NoRes}^{(1)} &= \frac{1}{64 {\rm \pi}} \left(S_{({-}1-2)12} + S_{2({-}1-2)1} + S_{12({-}1-2)}\right) \left(S_{({-}3-4)34} + S_{4({-}3-4)3} + S_{34({-}3-4)}\right)\nonumber\\ & \quad \times \left( \frac{1}{\omega_{1}+ \omega_{2}+ \omega_{1+2}} + \frac{1}{\omega_{3}+ \omega_{4}+ \omega_{3+4}} \right)\nonumber\\ &\quad + \frac{1}{16{\rm \pi}} \left( S_{(3-1)({-}3)1} + S_{1(3-1)({-}3)} \right) \left( S_{(4-2)2({-}4)} + S_{({-}4)(4-2)2} \right)\nonumber\\ &\quad \times \left( \frac{1}{\omega_{3-1}+ \omega_{3}- \omega_{1}} + \frac{1}{\omega_{4-2}+ \omega_{2}- \omega_{4}} \right)\nonumber\\ &\quad - \frac{1}{16{\rm \pi}} S_{12({-}1-2)} S_{34({-}3-4)} \left( \frac{1}{\omega_{1}+ \omega_{2}- \omega_{1+2}} + \frac{1}{\omega_{3}+ \omega_{4}- \omega_{3+4}} \right), \end{align}

(5.11) \begin{align} T_{NoRes}^{(2)} &= \frac{1}{32 {\rm \pi}} \left( S_{12({-}1-2)} B_{3({-}3-4)4} + S_{34({-}3-4)} B_{1({-}1-2)2} \right)\nonumber\\ &\quad \times \left( \frac{1}{\omega_1 + \omega_2 - \omega_{1+2}} + \frac{1}{\omega_3 + \omega_4 - \omega_{3+4}} \right)\nonumber\\ &\quad - \frac{1}{16{\rm \pi}} B_{(4-2)({-}4)2} \left( S_{(3-1)({-}3)1} + S_{({-}3)(3-1)1} \right)\nonumber\\ &\quad\times \left( \frac{1}{\omega_{3-1} + \omega_3-\omega_1} + \frac{1}{\omega_{4-2} + \omega_2-\omega_4} \right)\nonumber\\ & \quad - \frac{1}{16{\rm \pi}} B_{(3-1)1({-}3)} \left( S_{(4-2)2({-}4)} + S_{2(4-2) ({-}4)} \right)\nonumber\\ &\quad \times \left( \frac{1}{\omega_{3-1} + \omega_3-\omega_1} + \frac{1}{\omega_{4-2} + \omega_2-\omega_4} \right) \end{align}

and

(5.12)\begin{align} T_{NoRes}^{(3)} &={-} \frac{1}{64 {\rm \pi}} B_{1({-}1-2)2} B_{3({-}3-4)4} \left( \frac{1}{\omega_{1} + \omega_2 - \omega_{1+2}} + \frac{1}{\omega_{3} + \omega_4 - \omega_{3+4}} \right)\nonumber\\ &\quad + \frac{1}{16 {\rm \pi}} B_{(3-1)1({-}3)} B_{(4-2)({-}4)2} \left( \frac{1}{\omega_{3-1} + \omega_3 - \omega_1} + \frac{1}{\omega_{4-2} + \omega_2 - \omega_4} \right). \end{align}

Here, based on (4.15), $S_{(3-1)(-3)1}$ reads

(5.13)\begin{equation} S_{(3-1)({-}3)1} = \left( (k_3-k_1)k_1 + |k_3-k_1| |k_1| \right) \frac{a(k_3-k_1) a(k_1)}{a({-}k_3)}, \end{equation}

and $B_{(3-1)(-3)1} = S_{(3-1)(-3)1} ( 1- \chi _{\mathcal {N}_\mu } (k_3-k_1, -k_3, k_1) )$.

Proposition 5.3 Recall that $S_{123}$ is given by (4.15) and $B_{123} := S_{123} (1-\chi _{\mathcal {N}_\mu } (k_1, k_2, k_3))$. Then, $T_{Res} = T_{Res}^{(1)} + T_{Res}^{(2)}$ with

(5.14)\begin{align} T_{Res}^{(1)} &= \frac{1}{32{\rm \pi}} \left( \frac{S_{12({-}1-2)}}{\omega_1 + \omega_2 - \omega_{1+2}} [\chi S]_{3({-}3-4)4} + \frac{S_{34({-}3-4)}}{\omega_3 + \omega_4 - \omega_{3+4}} [\chi S]_{1({-}1-2)2} \right)\nonumber\\ &\quad - \frac{1}{16 {\rm \pi}(\omega_{3-1} + \omega_3 - \omega_{1})} \left( S_{(3-1)({-}3)1} + S_{({-}3)(3-1)1} \right) [\chi S]_{(4-2)({-}4)2}\nonumber\\ &\quad - \frac{1}{16 {\rm \pi}(\omega_{4-2} + \omega_2 - \omega_{4})} \left( S_{(4-2)2({-}4)} + S_{2(4-2)({-}4)} \right) [\chi S]_{(3-1)1({-}3)}, \end{align}

and

(5.15)\begin{align} T_{Res}^{(2)} &= \frac{1}{16{\rm \pi}} \left( \frac{B_{(3-1)1({-}3)}}{\omega_{3-1} + \omega_3 - \omega_{1}} [\chi S]_{(4-2)({-}4)2} + \frac{B_{(4-2)({-}4)2}}{\omega_{4-2} + \omega_2 - \omega_{4}} [\chi S]_{(3-1)1({-}3)} \right)\nonumber\\ &\quad - \frac{1}{64 {\rm \pi}} \left( \frac{B_{1({-}1-2)2}}{ \omega_1 + \omega_2 - \omega_{1+2}} [\chi S]_{3({-}3-4)4} + \frac{B_{3({-}3-4)4}}{ \omega_3 + \omega_4 - \omega_{3+4}} [\chi S]_{1({-}1-2)2} \right), \end{align}

where we define

(5.16)\begin{equation} [\chi S]_{123} := B_{123} - S_{123} = \chi_{\mathcal{N}_\mu} (k_1, k_2, k_3) S_{123}. \end{equation}

Proof. The proof is similar to the proof of proposition 5.2 in Appendix C. In particular, using the decomposition of $K^{(3)}$ given by terms (C1), we write

(5.17)\begin{equation} \{K^{(3)}, H_{Res}^{(3)}\} = \{K^{(3)}_0, H_{Res}^{(3)}\} + \{K^{(3)}_\chi, H_{Res}^{(3)}\}. \end{equation}

Then the first Poisson bracket leads to the coefficient (5.14), while the second bracket implies (5.15).

5.2. Modulational ansatz

We restrict ourselves to solutions in the form of near-monochromatic waves with carrier wavenumber $k_0 > 0$. For this purpose, we introduce the modulational ansatz

(5.18)\begin{equation} k = k_0 + \varepsilon \lambda, \quad \text{where} \ \frac{\lambda}{k_0} = O(1),\quad \varepsilon \ll 1, \end{equation}

and, accordingly, a function $U$ is defined in the Fourier space as

(5.19a,b)\begin{equation} U(\lambda) = z(k_0 + \varepsilon \lambda), \quad \bar{U}(\lambda) = \bar{z}(k_0 + \varepsilon \lambda), \end{equation}

where the time dependence is omitted. In the physical space,

(5.20)\begin{equation} z(x) = \frac{1}{\sqrt{2{\rm \pi}}} \int z(k) {\rm e}^{{\rm i}\,kx} \,{\rm d}k = \frac{\varepsilon}{\sqrt{2{\rm \pi}}} \int U(\lambda) {\rm e}^{{\rm i}\, k_0 x} {\rm e}^{{\rm i}\, \lambda \varepsilon x} \,{\rm d}\lambda = \varepsilon \, u(X) {\rm e}^{{\rm i}\,k_0 x}, \end{equation}

where $u$, as a function of the long scale $X = \varepsilon x$, is the inverse Fourier transform of $U$. After the change of variables (5.18), the following integral identity for any function $C$ holds:

(5.21)\begin{align} & \int C(k_1, k_2, k_3, k_4) z_1 z_2 \bar{z}_3 \bar{z}_4 \delta_{1+2-3-4} \,{\rm d}k_{1234} \nonumber\\ &\qquad = \varepsilon^3 \int \tilde{C}(\lambda_1, \lambda_2, \lambda_3, \lambda_4) U_1 U_2 \bar{U}_3 \bar{U}_4 \delta_{1+2-3-4}^{(\lambda)} \,{\rm d}\lambda_{1234}, \end{align}

where $U_j$ stands for $U(\lambda _j)$, the Dirac distribution on the right-hand side is defined as $\delta _{1+2-3-4}^{(\lambda )} = (1/2{\rm \pi} ) \int \exp ({-\textrm {i} (\lambda _1+ \lambda _2 - \lambda _3 - \lambda _4)X}) \,\textrm {d}X$ and $C(k_1, k_2, k_3, k_4) = \tilde {C}(\lambda _1, \lambda _2, \lambda _3, \lambda _4)$.

5.3. Vanishing of the resonant cubic term $H^{(3)}_{Res}$

We show that, under the modulational regime (5.18), the resonant cubic Hamiltonian $H^{(3)}_{Res}$ given in (4.14) is sufficiently small and can be neglected. The estimates will be done for a cubic term of type

(5.22)\begin{equation} \mathcal{I} := \int T_{k_1, k_2, k_3} z_1 \bar{z}_{2} z_3 \delta(k_1-k_2+k_3) \,{\rm d}k_{123}, \end{equation}

which can be similarly repeated for each term of (4.14). After the change of variables (5.18), we have

(5.23)\begin{equation} \mathcal{I} = \frac{\varepsilon^3}{2{\rm \pi}} \int {\rm e}^{-{\rm i}\,k_0 x} \int \tilde{T}_{\lambda_1, \lambda_2, \lambda_3} U_1 \bar{U}_2 U_3 \exp({-{\rm i}(\lambda_1 - \lambda_2 + \lambda_3) (\varepsilon x)})\, {\rm d}\lambda_{123} \,{{\rm d}\kern0.07em x}, \end{equation}

where $\tilde {T}_{\lambda _1, \lambda _2, \lambda _3} :={T}_{k_0+\varepsilon \lambda _1, k_0+ \varepsilon \lambda _2, k_0+ \varepsilon \lambda _3}$. We identify the inner integral above with the function $f(\varepsilon x)$ as

(5.24)\begin{equation} f(\varepsilon x) := \int \tilde{T}_{\lambda_1, \lambda_2, \lambda_3} U_1 \bar{U}_2 U_3 \exp({-{\rm i}(\lambda_1 - \lambda_2 + \lambda_3) (\varepsilon x)}) \,{\rm d}\lambda_{123}. \end{equation}

As a result, we have

(5.25)\begin{equation} \mathcal{I} = \frac{\varepsilon^3}{2{\rm \pi}} \int {\rm e}^{-{\rm i}\,k_0 x} f(\varepsilon x) \,{{\rm d}\kern0.07em x}. \end{equation}

To show that the integral (5.25) is negligible, we use the scale-separation lemma 4.4 of Guyenne et al. (Reference Guyenne, Kairzhan and Sulem2022a).

5.4. Quartic interactions in the modulational regime

We approximate the coefficients in propositions 5.1 and 5.2 under the modulational regime (5.18).

Lemma 5.4 Under the modulational ansatz (5.18), we have

(5.26)\begin{align} &\int T_0 z_1 z_2 \bar{z}_3 \bar{z}_4 \delta_{1+2-3-4} \,{\rm d}k_{1234}\nonumber\\ &\quad = \varepsilon^3 \int \left(c_0^l + \varepsilon c_0^r (\lambda_2 + \lambda_3) \right) U_1 U_2 \bar{U}_3 \bar{U}_4 \delta_{1+2-3-4}^{(\lambda)} \,{\rm d}\lambda_{1234} + O(\varepsilon^5), \end{align}

where $T_0$ is given in proposition 5.1 and the coefficients are

(5.27)\begin{equation} \left. \begin{aligned} & c_0^l = \frac{k_0^3}{8{\rm \pi}} + \frac{k_0^6}{16 {\rm \pi}\omega_0^2} \left(\frac{3}{2} \mathcal{P} - 5\mathcal{D} k_0^2 \right),\\ & c_0^r = \frac{3 k_0^2}{16{\rm \pi}} + \frac{k_0^6}{16 {\rm \pi}\omega_0^2} \left[ \left({\frac{3}{2} \mathcal{P}} - 5\mathcal{D} k_0^2 \right) \left(\frac{2}{k_0} + \frac{g+\mathcal{P} k_0^2 -3\mathcal{D} k_0^4}{2\omega_0^2} \right) -5\mathcal{D} k_0 \right]. \end{aligned} \right\} \end{equation}

Lemma 5.5 Under the modulational ansatz (5.18), we have

(5.28)$$\begin{align} &\int T_{NoRes}^{(1)} z_1 z_2 \bar{z}_3 \bar{z}_4 \delta_{1+2-3-4} \,{\rm d}k_{1234}\nonumber\\ &\quad = \varepsilon^3 \int \left(c_1^l + \varepsilon c_1^r (\lambda_2 + \lambda_3) \right) U_1 U_2 \bar{U}_3 \bar{U}_4 \delta_{1+2-3-4}^{(\lambda)} \,{\rm d}\lambda_{1234}\nonumber\\ &\qquad + \frac{\varepsilon^4 k_0^2}{4{\rm \pi}} \int |\lambda_1 - \lambda_3| U_1 U_2 \bar{U}_3 \bar{U}_4 \delta_{1+2-3-4}^{(\lambda)} \,{\rm d}\lambda_{1234} + O(\varepsilon^5), \end{align}$$

(5.29)\begin{align} &\int T_{NoRes}^{(2)} z_1 z_2 \bar{z}_3 \bar{z}_4 \delta_{1+2-3-4} \,{\rm d}k_{1234} = O(\varepsilon^5), \end{align}

(5.30)\begin{align} & \int T_{NoRes}^{(3)} z_1 z_2 \bar{z}_3 \bar{z}_4 \delta_{1+2-3-4} \,{\rm d}k_{1234} \nonumber\\ & \quad = \varepsilon^3 \int \left(1-\chi_{\mathcal{N}_\mu} (k_1, -k_1-k_2, k_2)\right) \left(1-\chi_{\mathcal{N}_\mu} (k_3, -k_3-k_4, k_4)\right)\nonumber\\ & \qquad \times \left(c_2^l + \varepsilon c_2^r (\lambda_2 + \lambda_3)\right) U_1 U_2 \bar{U}_3 \bar{U}_4 \delta_{1+2-3-4}^{(\lambda)} \,{\rm d}\lambda_{1234}\nonumber\\ & \qquad + \frac{\varepsilon^4 k_0^2}{4{\rm \pi}} \int \left(1-\chi_{\mathcal{N}_\mu} (k_3-k_1, k_1, -k_3)\right) \left(1-\chi_{\mathcal{N}_\mu} (k_4-k_2, -k_4, k_2)\right)\nonumber\\ & \qquad \times (1+ {\rm sgn} (\lambda_1-\lambda_3)) |\lambda_1 - \lambda_3| U_1 U_2 \bar{U}_3 \bar{U}_4 \delta_{1+2-3-4}^{(\lambda)} \,{\rm d}\lambda_{1234} + O(\varepsilon^5), \end{align}

where

(5.31) \begin{equation} \left. \begin{aligned} & c_1^l = \frac{k_0^3 \omega_0^2}{4{\rm \pi} \omega_{2k_0} (2\omega_0 + \omega_{2k_0})},\\ & c_1^r = \frac{c_1^l}{2} \left( \frac{3g-5\mathcal{P} k_0^2 + 7\mathcal{D} k_0^4}{\omega_{0}^2} + \frac{g+4\mathcal{P} k_0^2 -48 \mathcal{D} k_0^4}{\omega_{2k_0}^2}\right.\\ &\qquad - \left.\frac{g-12\mathcal{P} k_0^2 +80 \mathcal{D} k_0^4}{\omega_{2k_0}(2\omega_0 + \omega_{2k_0})} - \frac{g-3\mathcal{P} k_0^2 +5 \mathcal{D} k_0^4}{\omega_{0}(2\omega_0 + \omega_{2k_0})} \right),\\ & c_2^l ={-}\frac{k_0^3 \omega_0^2}{4{\rm \pi} \omega_{2k_0} (2\omega_0 - \omega_{2k_0})},\\ & c_2^r = \frac{c_2^l}{2} \left( \frac{3g-5\mathcal{P} k_0^2 + 7\mathcal{D} k_0^4}{\omega_{0}^2} + \frac{g+4\mathcal{P} k_0^2 -48 \mathcal{D} k_0^4}{\omega_{2k_0}^2}\right.\\ &\qquad +\left. \frac{g-12\mathcal{P} k_0^2 +80 \mathcal{D} k_0^4}{\omega_{2k_0}(2\omega_0 - \omega_{2k_0})} - \frac{g-3\mathcal{P} k_0^2 +5 \mathcal{D} k_0^4}{\omega_{0}(2\omega_0 - \omega_{2k_0})} \right). \end{aligned} \right\} \end{equation}

Lemma 5.6 Under the modulational ansatz (5.18), we have

(5.32) \begin{align} & \int T_{Res} z_1 z_2 \bar{z}_3 \bar{z}_4 \delta_{1+2-3-4} \,{\rm d}k_{1234} \nonumber\\ &\quad = \frac{\varepsilon^4 k_0^2}{8{\rm \pi}} \int \left[ \chi_{\mathcal{N}_\mu} (k_3-k_1, k_1, -k_3) + \chi_{\mathcal{N}_\mu} (k_4-k_2, -k_4, k_2)\right.\nonumber\\ &\left. \qquad -\, 2 \chi_{\mathcal{N}_\mu} (k_3-k_1, k_1, -k_3) \chi_{\mathcal{N}_\mu} (k_4-k_2, -k_4, k_2)\right]\nonumber\\ & \qquad \times (1+ {\rm sgn} (\lambda_1-\lambda_3)) |\lambda_1 - \lambda_3| U_1 U_2 \bar{U}_3 \bar{U}_4 \delta_{1+2-3-4}^{(\lambda)} \,{\rm d}\lambda_{1234}\nonumber\\ & \qquad + \frac{\varepsilon^3}{2} \int \left[ \chi_{\mathcal{N}_\mu} (k_1, -k_1-k_2, k_2) + \chi_{\mathcal{N}_\mu} (k_3, -k_3-k_4, k_4)\right.\nonumber\\ &\left. \qquad -\, 2 \chi_{\mathcal{N}_\mu} (k_1, -k_1-k_2, k_2) \chi_{\mathcal{N}_\mu} (k_3, -k_3-k_4, k_4) \right]\nonumber\\ & \qquad \times (c_2^l + \varepsilon c_2^r (\lambda_2 + \lambda_3)) U_1 U_2 \bar{U}_3 \bar{U}_4 \delta_{1+2-3-4}^{(\lambda)} \,{\rm d}\lambda_{1234}. \end{align}

Based on lemmas 5.4–5.6, the reduced Hamiltonian $H^{(4)}_+$ has the form

(5.33) $$\begin{gather} H^{(4)}_+= \varepsilon^3 \int \frac{\alpha}{4{\rm \pi}} U_1 U_2 \bar{U}_3 \bar{U}_4 \delta_{1+2-3-4}^{(\lambda)} \,{\rm d}\lambda_{1234}\nonumber\\ +\, \varepsilon^4 \int \left( \frac{\beta}{8{\rm \pi}} (\lambda_2 + \lambda_3) - \frac{k_0^2}{8{\rm \pi}} \gamma |\lambda_1 - \lambda_3| \right) U_1 U_2 \bar{U}_3 \bar{U}_4 \delta_{1+2-3-4}^{(\lambda)} \,{\rm d}\lambda_{1234} + O(\varepsilon^5), \end{gather}$$

where, using $k_j = k_0 + \varepsilon \lambda _j$ from (5.18), we denote

(5.34)\begin{gather} \frac{\alpha}{4{\rm \pi}}: = c_0^l - \frac{1}{2}c_1^l - \frac{1}{2} \left[1 - \chi_{\mathcal{N}_\mu} (k_1, -k_1-k_2, k_2) \chi_{\mathcal{N}_\mu} (k_3, -k_3-k_4, k_4) \right] c_2^l, \end{gather}

(5.35)\begin{gather} \frac{\beta}{8{\rm \pi}}:= c_0^r - \frac{1}{2}c_1^r - \frac{1}{2} \left[1 - \chi_{\mathcal{N}_\mu} (k_1, -k_1-k_2, k_2) \chi_{\mathcal{N}_\mu} (k_3, -k_3-k_4, k_4) \right] c_2^r \end{gather}

and

(5.36)\begin{equation} \gamma = 1 + \left[1-\chi_{\mathcal{N}_\mu} (k_3-k_1, k_1, -k_3) \chi_{\mathcal{N}_\mu} (k_4-k_2, -k_4, k_2)\right] (1+ {\rm sgn} (\lambda_1-\lambda_3)). \end{equation}

6.1. Explicit formulae for $\alpha, \beta$ and $\gamma$

6.1.1. Case $k_0 = 0.9$

Based on the construction of the neighbourhood around the resonant set in (3.28), we have the following identities:

(6.1)\begin{equation} \chi_{\mathcal{N}_\mu} (k_1, -k_1-k_2, k_2) = 0, \end{equation}

and

(6.2)\begin{equation} \chi_{\mathcal{N}_\mu} (k_4-k_2, -k_4, k_2) = 0 \quad \text{whenever} \ {\rm sgn} (k_4-k_2) ={+}1. \end{equation}

The validity of (6.1) follows from (5.18) and the choice $\mu = O(\varepsilon )$:

(6.3)\begin{align} & (k_1, k_2) \in \mathcal{B}_c (k_0) := \left\{ (x,y): x= k_0+O(\varepsilon), y= k_0+O(\varepsilon) \right\} \not \subseteq \mathcal{C}_\mu\nonumber\\ &\qquad \implies (k_1, -k_1-k_2, k_2) \notin \mathcal{N}_\mu. \end{align}

As for (6.2), under (5.18), we have $k_4-k_2 = \varepsilon (\lambda _4 - \lambda _2) = O(\varepsilon )$ and $k_2 = k_0+\varepsilon \lambda _2$. Then, if $k_4-k_2 >0$, the points $(k_4-k_2, k_2)$ are located far away from the set $\mathcal {C}_\mu ^+$ and so

(6.4)\begin{align} & (k_4-k_2, k_2) \in \mathcal{B}_s (k_0) := \left\{ (x,y): 0\leq x \leq O(\varepsilon), y = k_0 + O(\varepsilon) \right\} \not \subseteq \mathcal{C}^+_\mu\nonumber\\ &\qquad \implies (k_4-k_2, -k_4, k_2) \notin \mathcal{N}_\mu. \end{align}

As a result, the coefficients in (5.34)–(5.36) reduce to

(6.5a–c)\begin{equation} \frac{\alpha}{4{\rm \pi}} = c_0^l - \frac{1}{2} (c_1^l + c_2^l), \quad \frac{\beta}{8{\rm \pi}} = c_0^r - \frac{1}{2} (c_1^r + c_2^r), \quad \gamma = 2+{\rm sgn} (\lambda_1-\lambda_3). \end{equation}

When substituting $\gamma$ into the integral (5.33), the term $\textrm {sgn} (\lambda _1-\lambda _3)$ will disappear since

(6.6) \begin{align} & \int {\rm sgn} (\lambda_1-\lambda_3) |\lambda_1 - \lambda_3| U_1 U_2 \bar{U}_3 \bar{U}_4 \delta_{1+2-3-4}^{(\lambda)} \,{\rm d}\lambda_{1234}\nonumber\\ &\quad = \int (\lambda_1 - \lambda_3) U_1 U_2 \bar{U}_3 \bar{U}_4 \delta_{1+2-3-4}^{(\lambda)} \,{\rm d}\lambda_{1234}\nonumber\\ &\quad = \frac{1}{2} \int (\lambda_1+\lambda_2 - \lambda_3 - \lambda_4) U_1 U_2 \bar{U}_3 \bar{U}_4 \delta_{1+2-3-4}^{(\lambda)} \,{\rm d}\lambda_{1234} = 0, \end{align}

where we have used the index rearrangement $(\lambda _1, \lambda _2, \lambda _3, \lambda _4) \to (\lambda _2, \lambda _1, \lambda _3, \lambda _4)$ and the delta condition $\lambda _1+\lambda _2-\lambda _3-\lambda _4 = 0$. As a consequence, we will simply write $\gamma = 2$ in this case.

6.1.2. Case $k_0 = 2$

If $k_0$ is sufficiently large, we have

(6.7)\begin{equation} \chi_{\mathcal{N}_\mu} (k_1, -k_1-k_2, k_2) = 0, \end{equation}

and

(6.8)\begin{equation} \chi_{\mathcal{N}_\mu} (k_4-k_2, -k_4, k_2) = 1 \quad \text{whenever} \ {\rm sgn} (k_4-k_2) ={+}1. \end{equation}

Equation (6.7) holds due to estimates similar to (6.3). Equation (6.8) is different from (6.2) since, for a sufficiently large value of $k_0$, the points $(k_4-k_2, k_2)$ are located inside the set $\mathcal {C}_\mu ^+$,

(6.9)\begin{align} &(k_4-k_2, k_2) \in \left\{ (x,y): 0\leq x \leq O(\varepsilon), y = k_0 + O(\varepsilon) \right\} \subset \mathcal{C}^+_\mu \nonumber\\ &\quad \implies (k_4-k_2, -k_4, k_2) \in \mathcal{N}_\mu, \end{align}

and we get

(6.10a–c)\begin{equation} \frac{\alpha}{4{\rm \pi}} = c_0^l - \frac{1}{2} (c_1^l + c_2^l), \quad \frac{\beta}{8{\rm \pi}} = c_0^r - \frac{1}{2} (c_1^r + c_2^r), \quad \gamma = 1. \end{equation}

Figure 4 shows the location of points $(k_1, k_2)$ and $(k_4-k_2, k_2)$ with $k_4-k_2 >0$ for $k_0 = 0.9$ (figure 4a) and $k_0=2$ (figure 4b) relative to the neighbourhood $\mathcal {C}_\mu$. In figure 4(a) for $k_0 = 0.9$, Box 1 represents the set $\mathcal {B}_s (k_0)$ of all possible values of $(k_4-k_2, k_2)$ with $k_4-k_2 >0$ under the modulational ansatz (5.18), according to (6.4). Box 2 represents the set $\mathcal {B}_c (k_0)$ in (6.3). Similar sets are depicted in figure 4(b) for $k_0 = 2$.

Figure 4. Location of points $(k_1, k_2)$ and $(k_4-k_2, k_2)$ with $k_4-k_2 > 0$ relative to the neighbourhood $\mathcal {C}_\mu$, in the case (a) $k_0 = 0.9$ and (b) $k_0=2$ for $\mathcal {P} = 1$. Box 1 represents the set $\mathcal {B}_s (k_0)$ and Box 2 represents the set $\mathcal {B}_c (k_0)$.

6.2. Derivation of the Dysthe equation

The third-order normal form transformation has eliminated all cubic terms from the Hamiltonian $H$. In the modulational regime (5.18), the reduced Hamiltonian (now denoted by $H$) is

(6.11)\begin{equation} H = H^{(2)} + H_+^{(4)}, \end{equation}

up to fourth order. In the physical variables $(u, \bar {u})$, it reads

(6.12)\begin{align} H &= \varepsilon \int \bar{u} \, \omega (k_0 + \varepsilon D_X) \, u\, {\rm d}\kern0.07em X + \varepsilon^3 \frac{\alpha}{2} \int |u|^4 \,{\rm d}\kern0.07em X\nonumber\\ &\quad + \varepsilon^4 \frac{\beta}{2} \int |u|^2 {\rm Im} (\bar{u} \partial_X u) \,{\rm d}\kern0.07em X - \varepsilon^4 \frac{\gamma k_0^2}{4} \int |u|^2 |D_X| |u|^2 \, {\rm d}\kern0.07em X + O(\varepsilon^5), \end{align}

where $D_X = -\textrm {i} \, \partial _X$ for the slow spatial variable $X$ (so that its Fourier symbol is $\lambda$) and accordingly $|D_X|$ is the non-local operator whose Fourier symbol is $|\lambda |$. Expanding the linear dispersion relation in (6.12) around $k_0$ as

(6.13)\begin{equation} \omega (k_0+ \varepsilon D_X) = \omega_0 + \varepsilon \omega_0' D_X + \frac{1}{2}\varepsilon^2 \omega_0'' D^2_X+ \frac{1}{6}\varepsilon^3 \omega_0''' D^3_X + O(\varepsilon^4), \end{equation}

gives an alternate form of the Hamiltonian $H$ in physical variables

(6.14)$$\begin{gather} H = \int \left[ \varepsilon\, \omega_0 \, |u|^2 + \varepsilon^2 \omega_0' {\rm Im}(\bar{u} \partial_X u)+ \frac{1}{2} \varepsilon^3 \omega_0'' |\partial_X u|^2 + \frac{1}{2} \varepsilon^3 \alpha |u|^4\right.\nonumber\\ \left.- \frac{1}{6} \varepsilon^4 \omega_0''' {\rm Im} (\bar{u} \partial_X^3 u) + \frac{1}{2} \varepsilon^4 \beta |u|^2 {\rm Im} (\bar{u} \partial_X u) - \frac{1}{4} \varepsilon^4 \gamma k_0^2 |u|^2 |D_X| |u|^2 \right] \,{\rm d}\kern0.07em X + O(\varepsilon^5). \end{gather}$$

Coefficients in (6.13) have the following expressions:

(6.15)\begin{align} \left. \begin{aligned} \omega_0' & := \partial_k \omega(k_0) = \frac{g + 5 \mathcal{D} k_0^4 - 3 \mathcal{P} k_0^2}{2 \omega_0},\\ \omega_0'' & := \partial_k^2 \omega(k_0) = \frac{15 \mathcal{D}^2 k_0^8 - 22 \mathcal{D} \mathcal{P} k_0^6 + 30 g \mathcal{D} k_0^4 + 3 \mathcal{P}^2 k_0^4 - 6 g \mathcal{P} k_0^2 - g^2}{4 \omega_0^{3/2}},\\ \omega_0''' & := \partial_k^3 \omega(k_0) = 3 \left( 5 \mathcal{D}^3 k_0^{12} - 13 \mathcal{D}^2 \mathcal{P} k_0^{10} - 5 g \mathcal{D}^2 k_0^8 + 15 \mathcal{D} \mathcal{P}^2 k_0^8 \right.\\ & \quad - \left.34 g \mathcal{D} \mathcal{P} k_0^6 + \mathcal{P}^3 k_0^6 + 55 g^2 \mathcal{D} k_0^4 - 5 g \mathcal{P}^2 k_0^4 - 5 g^2 \mathcal{P} k_0^2 + g^3 \right)/\left( 8 \omega_0^{5/2} \right). \end{aligned} \right\} \end{align}

The Hamiltonian system

(6.16)\begin{equation} \partial_t \begin{pmatrix} u \\ \bar{u} \end{pmatrix} = \varepsilon^{{-}1} \begin{pmatrix} 0 & -{\rm i} \\ {\rm i} & 0 \end{pmatrix} \begin{pmatrix} \partial_u H \\ \partial_{\bar{u}} H \end{pmatrix}, \end{equation}

implies

(6.17)\begin{align} {\rm i} \, \partial_t u & = \varepsilon^{{-}1} \partial_{\bar{u}} H, \nonumber\\ & = \omega_0 u -{\rm i} \varepsilon \omega_0' \partial_X u - \frac{1}{2} \varepsilon^2 \omega_0'' \partial_X^2 u + \varepsilon^2 \alpha |u|^2 u\nonumber\\ & \quad + \frac{{\rm i}}{6} \varepsilon^3 \omega_0''' \partial_X^3 u - {\rm i}\, \varepsilon^3 \beta |u|^2 \partial_X u - \varepsilon^3 \frac{\gamma k_0^2}{2} u |D_X| |u|^2, \end{align}

which is a Hamiltonian Dysthe equation for two-dimensional hydroelastic waves on deep water. It describes modulated waves moving in the positive $x$-direction at group speed $\omega _0'$ as shown by the advection term. The non-local term $u|D_X| |u|^2$, a signature of the Dysthe equation, reflects the presence of the wave-induced mean flow as in the classical derivation using the method of multiple scales. The coefficient of this mean-flow term is similar to that in the pure gravity case ($\mathcal {D} = 0$, $\mathcal {P} = 0$) except for the dependence on $\gamma$.

The first two terms on the right-hand side of (6.17) can be eliminated via phase invariance and reduction to a moving reference frame. The latter is equivalent, in the framework of canonical transformations, to subtraction from $H$ of a multiple of the momentum (2.10) which reduces to

(6.18)\begin{equation} I = \int \eta \, (\partial_x \xi) \,{{\rm d}\kern0.07em x} = \int \big[ k_0 |u|^2 + \varepsilon \, {\rm Im}(\bar{u} \partial_{X} u) \big] \,{\rm d}\kern0.07em X, \end{equation}

while the former is equivalent to subtraction from $H$ of a multiple of the wave action

(6.19)\begin{equation} M = \varepsilon \int |u|^2 \, {\rm d}\kern0.07em X, \end{equation}

which is conserved due to the phase-invariance property of the Dysthe equation. The resulting Hamiltonian is given by

(6.20)\begin{equation} \hat{H} = H - \varepsilon \omega'_0 I - (\omega_0 - {k}_0 \omega'_0 ) M, \end{equation}

which, after introducing a new long-time scale $\tau = \varepsilon ^2 t$, leads to the following version of the Hamiltonian Dysthe equation:

(6.21)\begin{equation} {\rm i} \, \partial_\tau u ={-} \frac{1}{2} \omega_0'' \partial_X^2 u + \alpha |u|^2 u + \frac{{\rm i}}{6} \varepsilon \omega_0''' \partial_X^3 u - {\rm i}\, \varepsilon \beta |u|^2 \partial_X u - \varepsilon \frac{\gamma k_0^2}{2} u |D_X| |u|^2. \end{equation}

Remark 6.1 The Hamiltonian Dysthe equation derived by Craig et al. (Reference Craig, Guyenne and Sulem2021a) for two-dimensional deep-water surface gravity waves ($\mathcal {D} = 0$, $\mathcal {P} = 0$) reads

(6.22)\begin{equation} {\rm i} \, \partial_\tau u = \frac{g}{8\omega_0^3} \partial_X^2 u + k_0^3 |u|^2 u + {\rm i}\, \varepsilon \frac{g^3}{16 \omega_0^5} \partial_X^3 u - 3{\rm i}\, \varepsilon k_0^2 |u|^2 \partial_X u - \varepsilon k_0^2 u|D_X| |u|^2. \end{equation}

It is natural to expect that, in the limit $\mathcal {D} \to 0$ and $\mathcal {P} \to 0$, the coefficients of (6.21) converge to those of (6.22). Indeed, a direct computation verifies this convergence for the coefficients of the linear terms as well as $\alpha \to k_0^3$, $\beta \to 3k_0^2$ together with $\gamma \to 2$ for any $k_0$ because resonant triads are ruled out in this situation, hence $\chi _{\mathcal {N}_\mu } = 0$.

6.3. Envelope solitons

We take this opportunity to highlight a salient effect of ice compression when added to ice bending. As analysed by Akylas (Reference Akylas1983), a general result about the NLS equation in the focusing regime is that it admits permanent envelope solitons (i.e. solitary wave packets) provided that $c = c_g$ is realized in the wave system. If so, the wave packet has carrier wavenumber $k_{min}$ and propagates in such a way that its crests are stationary relative to its envelope.

In the flexural-gravity case (${\mathcal {P}} = 0$), it has been revealed that such envelope solitons do not exist (Milewski et al. Reference Milewski, Vanden-Broeck and Wang2011; Guyenne & Părău Reference Guyenne and Părău2012), i.e. the corresponding NLS equation is of defocusing type because of the sign of the Benjamin–Feir index (BFI) defined as the product of the dispersion coefficient $\omega_0''$ and the coefficient of the cubic term $\alpha$, namely

(6.23)\begin{equation} {\rm BFI} ={-}\omega_0'' \alpha ={-}0.006 < 0, \quad \text{at} \ k = k_{min} = 0.76, \end{equation}

according to the Benjamin–Feir (BF) criterion for modulational instability, and with $k_{min}$ defined by (3.32). On the other hand, we get

(6.24)\begin{equation} {\rm BFI} ={-}\omega_0'' \alpha = 0.2954 > 0, \quad \text{at} \ k = k_{min} = 0.88, \end{equation}

for ${\mathcal {P}} = 1$, which implies a focusing regime. Figure 5 portrays BFI for $k_{min}$ as a function of ${\mathcal {P}}$. As explained earlier, we restrict ourselves to the interval $0 \le {\mathcal {P}} < 2$ where non-conservative exponential wave growth/decay is ruled out. We find that the change from defocusing to focusing occurs around ${\mathcal {P}} = 0.39$ after which BFI quickly increases. This regime transition is significant, varying from negative BFI of low magnitude to positive BFI of high magnitude over a relatively short interval of ${\mathcal {P}}$. For example, if ${\mathcal {P}} = 1.9$, the much larger value of BFI $\simeq 55$ is beyond the range displayed in this figure.

Figure 5. The BF instability/stability criterion at $k_{min}$ for the NLS equation as a function of $\mathcal {P}$.

7.1. Stability of Stokes waves

We first give the theoretical prediction for modulational or BF instability of Stokes waves in sea ice. These are represented in (6.17) by the exact uniform solution

(7.1)\begin{equation} u_0(t) = B_0 \exp({-{\rm i} (\omega_0 + \varepsilon^2 \alpha B_0^2) t}), \end{equation}

where $B_0$ is a positive real constant. In the gravity water wave case (${\mathcal {D}} = 0$, ${\mathcal {P}} = 0$), such a solution is known to be linearly unstable with respect to sideband (i.e. long-wave) perturbations.

The formal calculation consists in linearizing (6.17) about $u_0$ by inserting a perturbation of the form

(7.2)\begin{equation} u(X,t) = u_0(t) \left[ 1 + B(X,t) \right], \end{equation}

where

(7.3)\begin{equation} B(X,t) = B_1 {\rm e}^{\varOmega t + {\rm i}\, \lambda X} + B_2 {\rm e}^{\bar{\varOmega} t - {\rm i}\, \lambda X}, \end{equation}

and $B_1, B_2$ are complex coefficients. We find that the condition $\operatorname {Re}(\varOmega ) \neq 0$ for instability implies

(7.4)\begin{equation} \alpha_1 ={-}\frac{\omega_0''}{2} \lambda^2 \left[ 2 B_0^2 \left( \alpha - \varepsilon \frac{\gamma k_0^2}{2} |\lambda| \right) + \frac{\omega_0''}{2} \lambda^2 \right] > 0. \end{equation}

This is a tedious but straightforward calculation for which we omit the details (Craig et al. Reference Craig, Guyenne and Sulem2021a; Guyenne et al. Reference Guyenne, Kairzhan and Sulem2022b).

In view of the subsequent comparison with numerical simulations of the Euler system, all variables associated with the envelope model are now rescaled to absorb $\varepsilon$ back into their definition. Figure 6 depicts the normalized growth rate

(7.5)\begin{equation} \frac{|\textrm{Re}(\varOmega)|}{\omega_0} = \frac{\sqrt{\alpha_1}}{\omega_0}, \end{equation}

delimiting the instability region as predicted by condition (7.4) for $(A_0,k_0) = (0.1,0.9)$ and $(0.01,5)$ with varying values of ${\mathcal {P}}$. These parameter regimes are representative of the two different cases $\gamma = 2$ and $1$ for the mean-flow term, as discussed in § 6.1. They correspond to initial wave steepness $\varepsilon = k_0 A_0 = 0.09$ and $0.05$, respectively. Note that the envelope amplitude $B_0$ and the surface amplitude $A_0$ are related by

(7.6)\begin{equation} B_0 = A_0 \sqrt{\frac{\omega_0}{2k_0}}, \end{equation}

according to (2.25) and (5.20).

Figure 6. Regions of BF instability according to (7.4) for (a) $(A_0,k_0) = (0.1,0.9)$ and (b) $(A_0,k_0) = (0.01,5)$. The various curves represent $\mathcal {P} = 0$ (red), $\mathcal {P} = 1$ (blue), $\mathcal {P} = 1.9$ (black).

We see in figure 6 that instability tends to be enhanced with increasing ${\mathcal {P}}$. The growth rate is especially strong for $(A_0,k_0) = (0.1,0.9)$ and ${\mathcal {P}} = 1.9$. By contrast, the variations in ${\mathcal {P}}$ are weaker for $(A_0,k_0) = (0.01,5)$. Maximum growth rate occurs around $\lambda = 0.02$ for $(A_0,k_0) = (0.1,0.9)$ and around $\lambda = 0.1$ for $(A_0,k_0) = (0.01,5)$. The wavenumber at maximum growth rate (as well as the extent of the instability region) differs by an order of magnitude between these two configurations.

The same remark can be made when comparing these results with predictions on BF instability for the gravity water wave problem with similar values of $\varepsilon$ (Craig et al. Reference Craig, Guyenne and Sulem2021a; Guyenne et al. Reference Guyenne, Kairzhan and Sulem2022b). Including elasticity leads to narrower regions of instability centred around smaller sideband wavenumbers $\lambda$, which are smaller by orders of magnitude than their counterparts in the pure gravity case. This implies that a much longer domain is needed in the present simulations to observe this instability with such long-wave modulations.

7.2. Reconstruction of the original variables

At any instant $t$, the surface elevation and velocity potential can be reconstructed from the wave envelope by inverting the normal form transformation. This is accomplished by solving the auxiliary system backward from $s = 0$ to $s = -1$, with ‘initial’ conditions given by the transformed variables

(7.7)$$\begin{gather} \eta(x,t) |_{s=0} = \frac{1}{\sqrt{2}} a^{{-}1}(D) \left[ u(x,t) {\rm e}^{{\rm i}\, k_0 x} + \bar{u}(x,t) {\rm e}^{-{\rm i}\, k_0 x} \right], \end{gather}$$

(7.8)$$\begin{gather}\xi(x,t) |_{s=0} = \frac{1}{{\rm i}\, \sqrt{2}} a(D) \left[ u(x,t) {\rm e}^{{\rm i}\, k_0 x} - \bar{u}(x,t) {\rm e}^{-{\rm i}\, k_0 x} \right], \end{gather}$$

according to (2.25) and (5.20). In these expressions, $u$ obeys (6.17) and $a^{-1}(D) = \sqrt {|D|/\omega (D)}$. The final solution at $s = -1$ represents the original variables $(\eta,\xi )$. During the evolution in $s$, higher-order harmonic contributions to $(\eta,\xi )$ are produced by the nonlinear terms in (4.20a,b) starting from the first harmonics (7.7)–(7.8) associated with carrier wavenumber $k_0$.

7.3. Simulations

We test this Hamiltonian Dysthe equation against the full equations (2.8) in the context of BF instability of Stokes waves in sea ice. Following a high-order spectral approach (Craig & Sulem Reference Craig and Sulem1993; Guyenne & Părău Reference Guyenne and Părău2012), (2.8) are discretized in space by a pseudospectral method via the fast Fourier transform (FFT). Nonlinear products are calculated in the physical space, while Fourier multipliers and spatial derivatives are evaluated in the Fourier space. The computational domain is set to $0 \le x < L$ with periodic boundary conditions and is divided into a regular mesh of $N$ collocation points. Both functions $\eta$ and $\xi$ are expressed in terms of truncated Fourier series

(7.9)\begin{equation} \begin{pmatrix} \eta(x_j,t) \\ \xi(x_j,t) \end{pmatrix} = \sum_{p ={-}N/2}^{N/2-1} \begin{pmatrix} \hat \eta_p(t) \\ \hat \xi_p(t) \end{pmatrix} {\rm e}^{{\rm i}\, \kappa_p x_j}, \quad \kappa_p = p \Delta \kappa, \quad x_j = j \Delta x, \quad j = 0, \dots, N-1, \end{equation}

where $\Delta \kappa = 2{\rm \pi} /L$ and $\Delta x = L/N$. Note that $\{ \kappa _p \}$ represents the set of discrete values associated with any wavenumber $k_\ell \in \mathbb {R}$. The DNO is approximated by its series expansion (2.16) for which a small number $M$ of terms is sufficient to achieve highly accurate results in light of its analyticity properties. The value $M = 6$ is selected based on previous extensive tests (Xu & Guyenne Reference Xu and Guyenne2009). Time integration of (2.8) is carried out in the Fourier space so that linear terms can be solved exactly by the integrating factor technique. The nonlinear terms are integrated in time by using a fourth-order Runge–Kutta scheme with constant step $\Delta t$. More details can be found in Guyenne & Părău (Reference Guyenne and Părău2012) and Xu & Guyenne (Reference Xu and Guyenne2009), with the exception that (2.8) for $\boldsymbol {v} = (\eta,\xi )^\top$ can be written as

(7.10)\begin{equation} \partial_t \boldsymbol{v} = \boldsymbol{L} \boldsymbol{v} + \boldsymbol{N}(\boldsymbol{v}), \end{equation}

where the linear part $\boldsymbol {L} \boldsymbol {v}$ is defined by

(7.11)\begin{equation} \boldsymbol{L} \boldsymbol{v} = \begin{pmatrix} 0 & G^{(0)} \\ -g - \mathcal{D} \partial_x^4 - \mathcal{P} \partial_x^2 & 0 \end{pmatrix} \begin{pmatrix} \eta \\ \xi \end{pmatrix}, \end{equation}

and the nonlinear part ${\boldsymbol N}({\boldsymbol v}) = ({\boldsymbol N}_1,{\boldsymbol N}_2)^\top$ is given by

(7.12)\begin{align} {\boldsymbol N}_1 &= G(\eta) - G^{(0)}, \end{align}

(7.13)\begin{align} {\boldsymbol N}_2 & ={-} \frac{1}{2} (\partial_x \xi)^2 + \frac{1}{2} \frac{\left[ G(\eta) \xi + (\partial_x \eta) (\partial_x \xi) \right]^2}{1 + (\partial_x \eta)^2} \nonumber\\ &\quad - \mathcal{D} \left( \partial_s^2 \kappa - \partial_x^4 \eta + \frac{1}{2} \kappa^3 \right) - \mathcal{P} (\kappa - \partial_x^2 \eta). \end{align}

In the Fourier space, the integrating factor $\varTheta _k(t)$ associated with $\boldsymbol {L} \boldsymbol {v}$ takes the form

(7.14)\begin{equation} \varTheta_k(t) = \begin{pmatrix} \cos(\omega_k t) & \sqrt{\dfrac{G^{(0)}}{g + \mathcal{D} k^4 - \mathcal{P} k^2}} \sin(\omega_k t) \\ -\sqrt{\dfrac{g + \mathcal{D} k^4 - \mathcal{P} k^2}{G^{(0)}}} \sin(\omega_k t) & \cos(\omega_k t) \end{pmatrix}, \end{equation}

for $k \neq 0$, and

(7.15)\begin{equation} \varTheta_0(t) = \begin{pmatrix} 1 & 0 \\ -g t & 1 \end{pmatrix}, \end{equation}

for $k = 0$ according to l'Hôpital's rule. The resulting fourth-order Runge–Kutta scheme reads

(7.16)$$\begin{gather} {\boldsymbol f}_1 = {\boldsymbol N}_k({\boldsymbol v}_k^n), \end{gather}$$

(7.17)$$\begin{gather}{\boldsymbol f}_2 = \varTheta_k \left( -\frac{\Delta t}{2} \right) {\boldsymbol N}_k \left[ \varTheta_k \left( \frac{\Delta t}{2} \right) \left( {\boldsymbol v}_k^n + \frac{\Delta t}{2} {\boldsymbol f}_1 \right) \right], \end{gather}$$

(7.18)$$\begin{gather}{\boldsymbol f}_3 = \varTheta_k \left( -\frac{\Delta t}{2} \right) {\boldsymbol N}_k \left[ \varTheta_k \left( \frac{\Delta t}{2} \right) \left( {\boldsymbol v}_k^n + \frac{\Delta t}{2} {\boldsymbol f}_2 \right) \right], \end{gather}$$

(7.19)$$\begin{gather}{\boldsymbol f}_4 = \varTheta_k (-\Delta t) \, {\boldsymbol N}_k \left[ \varTheta_k(\Delta t) \left( {\boldsymbol v}_k^n + \Delta t \, {\boldsymbol f}_3 \right) \right], \end{gather}$$

(7.20)$$\begin{gather}{\boldsymbol v}_k^{n+1} = \varTheta_k(\Delta t) {\boldsymbol v}_k^n + \frac{\Delta t}{6} \varTheta_k(\Delta t) \left(\,{\boldsymbol f}_1 + 2 {\boldsymbol f}_2 + 2 {\boldsymbol f}_3 + {\boldsymbol f}_4 \right), \end{gather}$$

for the solution ${\boldsymbol v}_k = (\eta _k,\xi _k)^\top$ advancing from time $t_n$ to time $t_{n+1} = t_n + \Delta t$.

The same numerical methods are used to solve the envelope equation (6.17), with the same resolutions in space and time. On the other hand, for the reconstruction procedure, the Fourier integrals in the auxiliary system are not evaluated by a pseudospectral method with the FFT because they cannot be expressed in terms of convolution integrals, unlike the case of surface gravity waves (Craig et al. Reference Craig, Guyenne and Sulem2021a; Guyenne et al. Reference Guyenne, Kairzhan and Sulem2022a). A reason is the more complicated expression of the linear dispersion relation in this problem due to ice bending and compression, which does not allow for further simplification of the interaction kernels in (4.20a,b). Instead, these Fourier integrals are computed by direct quadrature. More specifically, they are reduced to one-dimensional integrals in $k_1$ (or $k_2$) by integrating out the Dirac delta functions before applying the trapezoidal rule. The interval of integration as well as all Fourier coefficients in the integrands are confined to the truncated spectrum $-N/2 \le p \le N/2-1$ by virtue of (7.9). The coupled system (4.20a,b) is evolved in $s$ via a fourth-order Runge–Kutta scheme, with step size $\Delta s = 10 \Delta t$ or $100 \Delta t$. By construction, (4.20a,b) are purely nonlinear (i.e. quadratic in nonlinearity) and do not contain any stiff linear terms. Accordingly, the value of $\Delta s$ may be selected so that $0 < \Delta t \ll \Delta s \ll 1$, which helps speed up the $s$-marching process to mitigate the computational cost entailed by the quadrature rule. We have checked that using a smaller value of $\Delta s$ (say $\Delta s = \Delta t$) gives similar results.

While the FFT cannot be exploited for the reconstruction procedure, direct quadrature of the Fourier integrals in (4.20a,b) was not found to be a major issue in this two-dimensional setting. Because this computation is not performed at each instant $t$ (only when data are recorded) and because it is performed over a short interval $-1 \le s \le 0$, the associated cost is insignificant overall. We point out again that (4.20a,b) are an exact representation of the normal form transformation to eliminate non-resonant cubic interactions in this problem. They are solved numerically in their full form, without resorting to any asymptotic approximation. Indetermination due to $a_1^2$ or $a_3^2$ being singular at $k_1 = 0$ or $k_3 = 0$ in $P_{123}$ and $Q_{123}$ may be lifted by simply setting the corresponding contributions to zero by virtue of the zero-mass assumption. Because wavenumbers are discrete in this numerical context and to accommodate floating-point arithmetic, the characteristic function $\chi _{\mathcal {N}}(k_1,k_2,k_3)$ in (4.17) is implemented under relaxed conditions

(7.21)\begin{equation} \chi_{\mathcal{N}}(k_1, k_2, k_3) = \left\{ \begin{array}{@{}ll@{}} 1, & \text{if} \ k_1 + k_2 + k_3 = 0, \quad | \omega_1 - \omega_2 + \omega_3| < \Delta, \\ 0, & \text{otherwise}, \end{array} \right. \end{equation}

with $0 < \Delta \ll 1$. We have tested different values of $\Delta$, i.e. $\Delta = \{ 10^{-12}, 10^{-6}, 10^{-2} \}$ and obtained similar results, which suggests that such a resonant triad does not contribute substantially to the reconstruction process. This result may be attributed to the spectral discretization which limits the possibilities for resonances, or to conditions of the modulational regime being simulated for which quartic resonances are dominant.

Initial conditions of the form

(7.22)\begin{equation} u(x,0) = B_0 \left[ 1 + 0.1 \cos(\lambda x) \right], \end{equation}

are prescribed for (6.17), where a long-wave perturbation of wavenumber $\lambda \ll k_0$ is superimposed to a uniform solution of amplitude $B_0$. For the full system, the initial conditions $\eta (x,0)$ and $\xi (x,0)$ are reconstructed by solving (2.8) from transformed initial data (7.7)–(7.8) combined with (7.22). This allows for a meaningful comparison with matching initial conditions. Figure 7 illustrates the time evolution of the relative $L^2$ errors

(7.23)\begin{equation} \frac{\| \eta_f - \eta_w \|_2}{\| \eta_f \|_2}, \end{equation}

on $\eta$ between the fully ($\eta _f$) and weakly ($\eta _w$) nonlinear solutions for ${\mathcal {P}} = 1$ with $(A_0,k_0,\lambda ) = (0.1,0.9,0.02)$ and $(0.01,5,0.1)$ as considered in the previous stability analysis. The numerical parameters are set to $L = 200{\rm \pi}$, $N = 1024$ ($\Delta \kappa = 0.01$, $\Delta x = 0.61$), $\Delta t = 0.01$ for $(A_0,k_0,\lambda ) = (0.1,0.9,0.02)$, and to $L = 20{\rm \pi}$, $N = 1024$ ($\Delta \kappa = 0.1$, $\Delta x = 0.06$), $\Delta t = 0.0002$ for $(A_0,k_0,\lambda ) = (0.01,5,0.1)$. As mentioned earlier, a sufficiently long domain is specified in $x$ so that long-wave modulations can be resolved. Needless to say the time scale of these modulations is commensurate with their wavelength. For reference, errors from the cubic NLS equation (by neglecting the higher-order terms in (6.17)) are also shown in this figure. The same reconstruction procedure for $\eta$ based on (4.20a,b) is adopted in both cases. Overall, $L^2$ errors from the NLS equation are noticeably higher than those from the Dysthe equation, which confirms the superiority of the latter model in this asymptotic regime. Their values tend to gradually grow over time due to accumulation and amplification of numerical errors during the BF instability process. For $(A_0,k_0,\lambda ) = (0.1,0.9,0.02)$, these two models differ in performance by almost an order of magnitude. If we switched from $\gamma = 1$ to $\gamma = 2$ for $(A_0,k_0,\lambda ) = (0.01,5,0.1)$, these curves would be flipped with NLS errors being lower than Dysthe errors. Therefore, figure 7 may serve to validate (5.36) of $\gamma$, hence our treatment of cubic resonances.

Figure 7. Relative errors on $\eta$ between fully and weakly nonlinear solutions for $\mathcal {P} = 1$ with (a) $(A_0,k_0,\lambda ) = (0.1,0.9,0.02)$ and (b) $(A_0,k_0,\lambda ) = (0.01,5,0.1)$: blue curve, Dysthe equation; red curve, NLS equation.

A comparison of surface profiles $\eta$ obtained from these two models and the full nonlinear (Euler) system is presented in figures 8 and 9 with snapshots at various instants for ${\mathcal {P}} = 1$. Again, the two configurations $(A_0,k_0,\lambda ) = (0.1,0.9,0.02)$ and $(0.01,5,0.1)$ are examined. Due to the large aspect ratio and the presence of short oscillations, we simply plot the local maxima (i.e. crests) and local minima (i.e. troughs) of the Euler solution by using dots, without showing its entire profile to avoid cluttering these graphs. Development of the BF instability and associated near-recurrence over time are clearly observed. The snapshots at $t = 5000$ (figure 8) and at $t = 170$ (figure 9) correspond to a time of maximum wave growth with modes $\lambda = 0.02$ (two humps $|p| = 2$) and $\lambda = 0.1$ (a single hump $|p| = 1$) being amplified, respectively, which is consistent with the previous stability analysis (see figure 6).

Figure 8. Comparison of $\eta$ between fully and weakly nonlinear solutions for $(A_0,k_0,\lambda ) = (0.1,0.9,0.02)$ and $\mathcal {P} = 1$ at (a) $t = 0$, (b) $t = 3400$, (c) $t = 4000$, (d) $t = 5000$, (e) $t = 7000$, ( f) $t = 10\,000$: blue curve, Dysthe equation; red curve, NLS equation; black dots, Euler system.

Figure 9. Comparison of $\eta$ between fully and weakly nonlinear solutions for $(A_0,k_0,\lambda ) = (0.01,5,0.1)$ and $\mathcal {P} = 1$ at (a) $t = 0$, (b) $t = 170$, (c) $t = 320$, (d) $t = 620$, (e) $t = 860$, ( f) $t = 980$: blue curve, Dysthe equation; red curve, NLS equation; black dots, Euler system.

Here the waves travel from left to right but, owing to the periodic boundary conditions, they re-enter the domain from one side whenever they exit it through the other side. A striking feature of figure 8 for $(A_0,k_0,\lambda ) = (0.1,0.9,0.02)$ is the excellent match in phase, amplitude and shape between the Dysthe and Euler solutions throughout the entire simulation. By contrast, the NLS wave packet seems to move faster and the resulting shift accentuates over time. A slight left–right asymmetry is discernible near the tails of the Dysthe and Euler pulses, especially when wave modulation is strong, while such an asymmetry is absent from the NLS pulse. The presence of spatial derivatives of odd order in the Dysthe equation may explain this phenomenon as opposed to the NLS equation. For $(A_0,k_0,\lambda ) = (0.01,5,0.1)$, while very good agreement is obtained between the Dysthe and Euler solutions at early stages of their time evolution, discrepancies arise during subsequent cycles of BF instability as shown in figure 9. Notably, the Dysthe wave packet is found to be steeper (respectively, less steep) than the Euler wave packet at $t = 620$ (respectively, $t = 980$), although their phases still match well. In spite of these differences, all three solutions behave qualitatively in the same way during the process of wave modulation–demodulation. These results are in line with the $L^2$ errors reported in figure 7. Similar observations can be made for other values of $0 \le {\mathcal {P}} < 2$ and are not displayed for convenience.

We further assess the performance of our reconstruction procedure by testing it against independent predictions from the NLS equation proposed by Trichtchenko et al. (Reference Trichtchenko, Milewski, Părău and Vanden-Broeck2019) for flexural-gravity waves (${\mathcal {P}} = 0$) on deep water. In this framework, the surface variables are determined perturbatively by a Stokes expansion

(7.24a,b)\begin{equation} \eta = \eta_1 {\rm e}^{{\rm i}\, \theta} + \eta_2 {\rm e}^{2 {\rm i}\, \theta} + {\rm c.c.}, \quad Q = \omega_0 \eta_1 {\rm e}^{{\rm i} \,\theta} + 2 \omega_0 \eta_2 {\rm e}^{2 {\rm i}\, \theta} + {\rm c.c.}, \end{equation}

up to second order, where $\textrm {c.c.}$ denotes complex conjugation, $\theta = k_0 x - \omega _0 t$ and $Q = \partial _x \xi$. The surface velocity potential $\xi$ can be readily deduced from $Q$ as $\xi = \partial _x^{-1} Q$ by using the FFT combined with the zero-mass assumption (to handle indetermination at $k = 0$). The carrier wave envelope $\eta _1$ obeys

(7.25)\begin{equation} {\rm i} \, (\partial_t + \omega_0' \partial_x) \eta_1 + \frac{\omega_0''}{2} \partial_x^2 \eta_1 + \varGamma |\eta_1|^2 \eta_1 = 0, \end{equation}

with $\omega _0^2 = g k_0 + k_0^5$ and

(7.26a–c)\begin{align} \omega_0' = \frac{g + 5 k_0^4}{2 \omega_0}, \quad \omega_0'' ={-}\frac{g^2 - 30 g k_0^4 - 15 k_0^8}{4 \omega_0^3}, \quad \varGamma ={-}\frac{\omega_0 k_0^2 (4 g^2 - 27 g k_0^4 + 44 k_0^8)}{2 (g + k_0^4) (g - 14 k_0^4)}, \end{align}

while the second-harmonic component $\eta _2$ is bound to $\eta _1$ via

(7.27)\begin{equation} \eta_2 = \frac{\omega_0^2 \eta_1^2}{g - 14 k_0^4}. \end{equation}

Recall that the coefficients $\omega _0''$ and $\varGamma$ for dispersion and nonlinearity coincide with ours in this situation. The NLS equation (7.25) was derived following the method of multiple scales based on an integral reformulation of the Euler system (Trichtchenko et al. Reference Trichtchenko, Milewski, Părău and Vanden-Broeck2019). The algebraic equations (7.24a,b) to reconstruct $(\eta,\xi )$ are more explicit and more efficient than (4.20a,b), but they are perturbative and devoid of any Hamiltonian structure. On the other hand, our reconstruction procedure is non-perturbative and involves solving an auxiliary system of Hamiltonian differential equations, as motivated by the basic Hamiltonian formulation of this problem.

Restricting our attention to ${\mathcal {P}} = 0$, figure 10 compares $L^2$ errors (7.23) from our NLS equation with $(\eta,\xi )$ calculated according to (4.20a,b) and from its counterpart (7.25) with $(\eta,\xi )$ given by (7.24a,b). For each of these models, the $L^2$ error is estimated relative to the Euler solution with respective initial conditions $\eta (x,0)$ and $\xi (x,0)$. In the case where (7.25) is tested against (2.8), these initial conditions are provided by (7.24a,b) and (7.27) with

(7.28)\begin{equation} \eta_1(x,0) = \frac{A_0}{2} \left[ 1 + 0.1 \cos(\lambda x) \right], \end{equation}

noting again that $A_0$ and $B_0$ are related through (7.6) to ensure matching wave parameters between these two NLS models. We adopt the same numerical schemes as described earlier, together with the same resolutions in space and time, to solve (7.25) and evaluate (7.24a,b). The initially perturbed Stokes wave is again prescribed by $A_0 = 0.1$, $k_0 = 0.9$ and $\lambda = 0.02$. We find that the two sets of errors are indistinguishable as shown in this figure. This is not surprising given the fact that, for such mild initial conditions and for such a weakly nonlinear solution as produced by the NLS equation, high-order harmonics are not expected to contribute significantly to the wave spectrum. In light of this remarkable agreement, figure 10 helps provide a cross-validation of these two independent methods for reconstructing $(\eta,\xi )$ from the wave envelope: on one hand, our dynamical system approach based on the numerical solution of (4.20a,b) and, on the other hand, an asymptotic approach based on the second-order algebraic formulae (7.24a,b). Such a validation was not presented by Trichtchenko et al. (Reference Trichtchenko, Milewski, Părău and Vanden-Broeck2019). Finally, the conservation of energy $H$ and wave action $M$ by the Dysthe equation (6.17) is demonstrated in figure 11 where the relative errors

(7.29a,b)\begin{equation} \frac{\Delta H}{H_0} = \frac{|H - H_0|}{H_0}, \quad \frac{\Delta M}{M_0} = \frac{|M - M_0|}{M_0}, \end{equation}

are plotted as functions of time $t$ for $A_0 = 0.1$, $k_0 = 0.9$ and $\lambda = 0.02$. Integration over $x$ in (6.14) and in the $L^2$ norm (7.23) is carried out via the trapezoidal rule over the periodic domain $[0,L]$. The reference values $H_0$ and $M_0$ denote the initial values of (6.14) and (6.19) at $t = 0$. Overall, both $H$ and $M$ are very well conserved by (6.17) for various values of ${\mathcal {P}}$. The gradual loss of accuracy over time is likely attributed to accumulation of numerical errors as enhanced by the BF instability.

Figure 10. Relative errors on $\eta$ between fully and weakly nonlinear solutions for $\mathcal {P} = 0$ with $(A_0,k_0,\lambda ) = (0.1,0.9,0.02)$: blue curve, our NLS equation; red curve, NLS equation from Trichtchenko et al. (Reference Trichtchenko, Milewski, Părău and Vanden-Broeck2019).

Figure 11. Relative errors on (a) $H$ and (b) $M$ for the Dysthe equation with $(A_0,k_0,\lambda ) = (0.1,0.9,0.02)$ and $\mathcal {P} = 0$ (red), $\mathcal {P} = 1$ (blue), $\mathcal {P} = 1.9$ (black).