2.1. Governing equations

We consider the evolution of a free surface $\{ y = \eta (x,t) \}$ on top of a two-dimensional fluid of infinite depth

(2.1)\begin{equation} \mathcal{S}(\eta) = \left\{ x \in \mathbb{R}, -\infty < y < \eta(x,t) \right\} , \end{equation}

under the influence of gravity. Assuming that the flow is incompressible and inviscid, the velocity field, denoted by ${{\boldsymbol {u}}}(x,y,t) = (u(x,y,t), v(x,y,t))^{\rm T}$, satisfies the Euler equations

(2.2)\begin{equation} \left.\begin{gathered} {{\partial}}_t {{\boldsymbol{u}}} + ({{\boldsymbol{u}}} \boldsymbol{\cdot}\boldsymbol{\nabla}) {{\boldsymbol{u}}} + \frac{1}{\rho}\,\boldsymbol{\nabla}P - \boldsymbol{g} = {\bf{0}}, \\ \boldsymbol{\nabla}\boldsymbol{\cdot} {{\boldsymbol{u}}} = 0 , \end{gathered}\right\} \end{equation}

where $\rho$ is the (constant) fluid density, $P(x,y,t)$ is the pressure, and $\boldsymbol {g} = (0, -g)^\textrm {T}$ is the acceleration due to gravity. Hereafter, the symbol $\boldsymbol {\nabla }$ denotes the spatial gradient $(\partial _x,\partial _y)^\textrm {T}$ when applied to functions, or the variational gradient when applied to functionals. The vorticity, defined as

(2.3)\begin{equation} \gamma = {{\partial}}_x v - {{\partial}}_y u , \end{equation}

satisfies

(2.4)\begin{equation} \partial_t \gamma +({{\boldsymbol{u}}} \boldsymbol{\cdot}\boldsymbol{\nabla}) \gamma=0 , \end{equation}

which expresses the well-known fact that in two dimensions, the vorticity is conserved along particle trajectories. In particular, if the initial vorticity is constant throughout the fluid domain, then it remains constant. The present study is devoted to flows with the property of having non-zero constant vorticity.

The boundary conditions at the free surface are composed of the dynamical condition

(2.5)\begin{equation} P = P_0 , \end{equation}

where $P_0$ is the atmospheric constant, and the kinematic condition

(2.6)\begin{equation} v= \partial_t \eta + u\,\partial_x \eta . \end{equation}

The divergence-free condition implies the existence of a stream function $\psi$ such that

(2.7a,b)\begin{equation} u = \partial_y\psi , \quad v={-}\partial_x \psi , \end{equation}

satisfying

(2.8)\begin{equation} {-\Delta \psi} = \gamma . \end{equation}

By construction, $\tilde {\psi } = \psi + \gamma y^{2}/2$ is a harmonic function. Introducing the generalized potential $\varphi$, defined as the harmonic conjugate of $\tilde {\psi }$, we have

(2.9)\begin{equation} \left.\begin{gathered} \partial_x \varphi = \partial_y \tilde{\psi} = u + \gamma y , \\ \partial_y \varphi ={-}\partial_x \tilde{\psi} = v . \end{gathered}\right\} \end{equation}

The presence of constant vorticity induces a horizontal background shear current that has a linear profile in the vertical direction. We associate $\gamma > 0$ with a co-propagating current because it contributes positively to the horizontal fluid velocity $u = \partial _x \varphi - \gamma y$ for $y < 0$ (along most of the water column), while $\gamma < 0$ is associated with a counter-propagating current.

In the variables $\varphi$ and $\tilde {\psi }$, (2.2) takes the form

(2.10)\begin{equation} \boldsymbol{\nabla}\left[ \partial_t \varphi + \tfrac{1}{2} \left( (\partial_x \varphi)^{2} + (\partial_y \varphi)^{2} \right) +\gamma \tilde{ \psi} - \gamma y\,\partial_x \varphi + P + g y \right] = \boldsymbol{0} , \end{equation}

and after integration, it reduces to

(2.11)\begin{equation} \partial_t \varphi + \tfrac{1}{2} \left( (\partial_x \varphi)^{2} + (\partial_y \varphi)^{2} \right) + \gamma \tilde{ \psi} - \gamma \eta\,\partial_x \varphi + g \eta = 0 , \end{equation}

at the free surface $y=\eta$, with $P_0 = 0$ without loss of generality. The Euler system can then be written as

(2.12)\begin{gather} \Delta \varphi = 0 \quad {\textrm{in}}\ \mathcal{S} , \end{gather}

(2.13)\begin{gather}\partial_t \eta - \partial_y \varphi + (\partial_x \varphi) (\partial_x \eta) - \gamma \eta\,\partial_x \eta = 0 ,\quad {\rm{on}} \ y=\eta(x,t) , \end{gather}

(2.14)\begin{gather}\partial_t \varphi + \tfrac{1}{2} \left( (\partial_x \varphi)^{2} + (\partial_y \varphi)^{2} \right) + \gamma \tilde{\psi} - \gamma \eta\,\partial_x \varphi + g\eta = 0 , \quad {\rm{on}} \ y=\eta(x,t) , \end{gather}

with the condition that $\varphi, \tilde {\psi } \to 0$ uniformly in $x$ as $y\to -\infty$.

2.2. Hamiltonian formulation

It is well known since the seminal paper of Zakharov (Reference Zakharov1968) that the irrotational ($\gamma =0$) water wave system has a Hamiltonian formulation in the variables $(\eta, \xi )$, where $\xi (x,t) = \varphi (x,\eta (x,t),t)$ is the trace of the velocity potential on the free surface.

Wahlén (Reference Wahlén2007) (see also Constantin et al. Reference Constantin, Ivanov and Prodanov2008) observed that in the presence of constant vorticity, the water wave system in $(\eta, \xi )$ can still be expressed as a Hamiltonian system but in non-canonical form, namely

(2.15)\begin{equation} \partial_t \begin{pmatrix} \eta \\ \xi \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & \gamma\,\partial_x^{{-}1} \end{pmatrix} \begin{pmatrix} \partial_\eta H \\ \partial_\xi H \end{pmatrix} . \end{equation}

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

(2.16)\begin{equation} H (\eta, \xi) = \frac{1}{2} \int_\mathbb{R} \left[ \xi\,G(\eta)\,\xi - \gamma \eta^{2}\,\partial_x \xi + \frac{\gamma^{2}}{3}\, \eta^{3} + g \eta^{2} \right] {{\rm d}\kern0.06em x} \,, \end{equation}

where $G(\eta )$ is the Dirichlet–Neumann operator in the fluid domain, which associates to the Dirichlet data $\xi$ on $y=\eta (x)$ the normal derivative of the harmonic function $\varphi$ with an additional normalized factor, namely

(2.17)\begin{equation} G(\eta) : \xi \longmapsto \sqrt{1+(\partial_x \eta)^{2}} \, \partial_{n} \varphi \big|_{y=\eta} . \end{equation}

Other invariants of motion are the volume (or mass)

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

and the momentum (or impulse)

(2.19)\begin{equation} I = \int_\mathbb{R} \left( \eta\,\partial_x \xi - \tfrac{1}{2} \gamma \eta^{2} \right) {{\rm d}\kern0.06em x} . \end{equation}

Wahlén (Reference Wahlén2007) found that under the change of variables

(2.20)\begin{equation} (\eta, \xi) \to \left( \eta, \zeta = \xi-\frac{\gamma}{2}\,\partial_x^{{-}1} \eta \right) , \end{equation}

system (2.15) can be transformed into canonical form

(2.21)\begin{equation} \partial_t \begin{pmatrix} \eta \\ \zeta \end{pmatrix} = J \, \boldsymbol{\nabla} \mathcal{H}(\eta,\zeta) = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} \partial_\eta \mathcal{H} \\ \partial_\zeta \mathcal{H} \end{pmatrix} , \end{equation}

where the Hamiltonian $\mathcal {H}(\eta, \zeta ) = H(\eta, \xi )$ in the new variables is

(2.22)\begin{equation} \mathcal{H} (\eta, \zeta) = \frac{1}{2} \int_\mathbb{R} \left[ \left( \zeta + \frac{\gamma}{2}\,\partial_x^{{-}1} \eta \right) G(\eta) \left( \zeta + \frac{\gamma}{2}\,\partial_x^{{-}1} \eta \right) - \gamma \eta^{2}\,\partial_x\zeta - \frac{\gamma^{2}}{6}\,\eta^{3} + g \eta^{2} \right] {{\rm d}\kern0.06em x} . \end{equation}

In the above formulas, $\partial _x^{-1} \eta$ is defined as $\partial _x^{-1} \eta (x) = \int _{-\infty }^{x} \eta (s)\,\textrm {d}s$. We assume that $\eta \to 0$ as $x\to \pm \infty$. Furthermore, we assume that $\int _{-\infty }^{\infty } \eta \,{\textrm {d}\kern0.06em x} =0$ at $t=0$, which implies that it remains so for all $t$. Indeed, integrating (2.13) in $x$ over $\mathbb {R}$, we get

(2.23)\begin{equation} \partial_t \int_\mathbb{R} \eta \,{{\rm d}\kern0.06em x} - \int_\mathbb{R} \left( G(\eta)\,\xi - \frac{ \gamma }{2}\,\partial_x \eta^{2} \right) {{\rm d}\kern0.06em x} = 0 . \end{equation}

The second term on the left-hand side rewrites as $\int _\mathbb {R} G(\eta )\,\xi \,{\textrm {d}\kern0.06em x} = \int _{\partial \mathcal {S}} \partial _n \varphi \,\textrm {d} \sigma = \int _{\mathcal {S}} \Delta \varphi \,{\textrm {d}\kern0.06em x}$ ${\textrm {d} y} = 0$, leading to $\int _\mathbb {R} \eta \,{\textrm {d}\kern0.06em x} = 0$ for all $t$.

Denoting by $\hat {f}_k = ({1}/{\sqrt {2 {\rm \pi}}}) \int _\mathbb {R} \textrm {e}^{-\textrm {i}kx}\,f(x) \,{\textrm {d}\kern0.06em x}$ the Fourier transform of a function $f(x)$, we have in particular

(2.24a,b)\begin{equation} \partial_x^{{-}1} \eta(x) ={-} \frac{{\rm i}}{\sqrt{2 {\rm \pi}}} \int_\mathbb{R} \frac{1}{k}\, \hat{ \eta}_k\,{\rm e}^{{\rm i}kx} \,{\rm d}k , \quad \partial_x \zeta(x) = \frac{{\rm i}}{\sqrt{2 {\rm \pi}}} \int_\mathbb{R} k \hat{\zeta}_k\,{\rm e}^{{\rm i}kx} \,{\rm d}k . \end{equation}

From the above assumption,

(2.25)\begin{equation} \hat{\eta}_0 = 0 , \quad \hat{\eta}_k = {O}(k)\ \text{for small} \ k . \end{equation}

In the following, we will drop the hat notation. Since $\eta (x)$ and $\zeta (x)$ are real-valued functions, we have the relation $(\eta _{-k}, \zeta _{-k}) = (\bar {\eta }_k, \bar {\zeta }_k)$ where the overbar stands for complex conjugation.

2.3. Taylor expansion of the Hamiltonian near equilibrium

It can be shown that the Dirichlet–Neumann operator is analytic in $\eta$ (Coifman & Meyer Reference Coifman and Meyer1985) and admits a convergent Taylor series expansion

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

about the quiescent state $\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 = -\textrm {i} \, \partial _x$, the first three terms are

(2.27)\begin{equation} \left. \begin{array}{l@{}} G^{(0)}(\eta) = |D| , \\[3pt] G^{(1)}(\eta) = D \eta D - G^{(0)}\eta G^{(0)} , \\[3pt] 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{array} \right\} \end{equation}

In Fourier variables, substituting the expansion for $G(\eta )$ into the Hamiltonian (2.22), we get

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

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

(2.29)\begin{equation} \left.\begin{gathered} \mathcal{H}^{(2)} = \frac{1}{2} \int_\mathbb{R} \left[ |k| \left( \zeta_{k} - \frac{{\rm i} \gamma}{2k}\,\eta_{k} \right)\left( \bar{\zeta}_{k} + \frac{{\rm i} \gamma}{2k}\,\bar{\eta}_{k} \right) + g |\eta_{k}|^{2} \right] {\rm d}k , \\ \mathcal{H}^{(3)} = \frac{1}{2 \sqrt{2{\rm \pi}}} \int_{\mathbb{R}^{3}} \left[ ({-}k_1 k_3 - |k_1|\,|k_3|) \left( \zeta_1 - \frac{{\rm i} \gamma}{2k_1}\,\eta_1 \right) \eta_2 \left( \zeta_3 - \frac{{\rm i} \gamma}{2k_3}\,\eta_3 \right)\right.\\ \left. - {\rm i} \gamma k_1 \zeta_1 \eta_2 \eta_3- \frac{\gamma^{2}}{6}\,\eta_1 \eta_2 \eta_3 \right] \delta_{123} \,{\rm d}k_{123} , \\ \mathcal{H}^{(4)} ={-}\frac{1}{8 {\rm \pi}} \int_{\mathbb{R}^{4}} |k_1|\,|k_4|\,(|k_1|+|k_4|-2|k_3+k_4|)\left( \zeta_1 - \frac{{\rm i} \gamma}{2k_1}\,\eta_1 \right) \eta_2 \eta_3\\ \times \left( \zeta_4 - \frac{{\rm i} \gamma}{2k_4}\,\eta_4 \right) \delta_{1234} \,{\rm d}k_{1234}. \end{gathered}\right\} \end{equation}

In the above expressions, we use the compact notations $(\eta _j, \zeta _j)= (\eta _{k_j}, \zeta _{k_j})$, $\textrm {d}k_{1 \dots n} = \textrm {d}k_1 \dots \textrm {d}k_n$ and $\delta _{1 \dots n} = \delta (k_1+ \dots + k_n)$, where $\delta (k) = ({1}/{2{\rm \pi} })\int _\mathbb {R} \textrm {e}^{-\textrm {i}kx}\, {\textrm {d}\kern0.06em x}$ is the Dirac distribution. Hereafter, the domain of integration is omitted in integrals and is understood to be $\mathbb {R}$ for each $x_j$ or $k_j$.

2.4. Linearization near equilibrium

The linearized water wave system about still water, written in terms of $(\eta,\zeta )$, is

(2.30)\begin{equation} \partial_t \begin{pmatrix} \eta\\ \zeta \end{pmatrix} = J \, \boldsymbol{\nabla} \mathcal{H}^{(2)}(\eta,\zeta) = \begin{pmatrix} g- \dfrac{\gamma^{2}}{4}\,\partial_x^{{-}1} |D|\,\partial_x^{{-}1} & -\dfrac{\gamma}{2}\,\partial_x^{{-}1} |D| \\[6pt] \dfrac{\gamma}{2}\,\partial_x^{{-}1} |D| & |D| \end{pmatrix} \begin{pmatrix} \eta\\ \zeta \end{pmatrix} . \end{equation}

We now introduce the symplectic complex coordinate

(2.31)\begin{equation} z := \frac{1}{\sqrt 2} \left( a(D)\,\eta + {\rm i}\, a(D) ^{{-}1}\,\zeta \right) , \end{equation}

where $a^{2}(D) := \omega (D)/|D|$, and $\omega (D) := \sqrt {\gamma ^{2}/4 + g\,|D|}$. The mapping $(\eta,\zeta ) \to (z,\bar {z})$ is canonical, and in these variables, the water wave system becomes

(2.32)\begin{equation} \partial_t \begin{pmatrix} z \\ \bar{z} \end{pmatrix} = \begin{pmatrix} 0 & -{\rm i} \\ {\rm i} & 0 \end{pmatrix} \begin{pmatrix} \partial_{z} \mathcal{H} \\ \partial_{\bar{z}} \mathcal{H} \end{pmatrix} := J_1 \begin{pmatrix} \partial_{z} \mathcal{H} \\ \partial_{\bar{z}} \mathcal{H} \end{pmatrix} . \end{equation}

Equivalently, in the Fourier space,

(2.33)\begin{equation} z_k := \frac{1}{\sqrt{2}} \left( a_k \eta_k + {\rm i}\, a_k^{{-}1} \zeta_k \right) , \end{equation}

where $a_k^{2} = \omega _k /|k|$ and $\omega _k= \sqrt {\gamma ^{2}/4 + g\,|k|}$. Since the functions $\eta (x)$ and $\zeta (x)$ are real-valued in the physical space, we also have

(2.34)\begin{equation} \bar{z}_{{-}k} = \frac{1}{\sqrt{2}} \left( a_k \eta_k - {\rm i} \, a_k^{{-}1} \zeta_k \right) . \end{equation}

We can express $\eta _k$ and $\zeta _k$ in terms of $z_k$ as

(2.35a,b)\begin{equation} \eta_k = \frac{1}{\sqrt 2}\,a_k^{{-}1} (z_k + \bar{z}_{{-}k}), \quad \zeta_k = \frac{1}{{\rm i}\sqrt{2}}\,a_k (z_k - \bar{z}_{{-}k}) . \end{equation}

The quadratic term $\mathcal {H}^{(2)}$ given in (2.29) diagonalizes as

(2.36)\begin{equation} \mathcal{H}^{(2)} = \int\varOmega_k\,|z_k|^{2} \,{\rm d}k , \end{equation}

where

(2.37)\begin{equation} \varOmega_k = \varOmega(k) = \frac{\gamma}{2}\,{{\rm sgn}}(k) + \omega_k \end{equation}

is the linear dispersion relation for deep-water gravity waves with constant vorticity (Berti et al. Reference Berti, Franzoi and Maspero2021). The linearized system with $\mathcal {H}$ replaced by $\mathcal {H}^{(2)}$ in (2.32) reduces to the scalar equation $\partial _t z_k = -\textrm {i} \, \varOmega _k z_k$.

We define the Poisson bracket of two functionals $K(\eta, \zeta )$ and $H(\eta, \zeta )$ of real-valued functions $\eta$ and $\zeta$ as

(2.38)\begin{equation} \{K, H\} = \int \left( \partial_\eta H\,\partial_\zeta K - \partial_\zeta H\,\partial_\eta K \right) {{\rm d}\kern0.06em x} . \end{equation}

Assuming that $K$ and $H$ are real-valued, and using the Plancherel formula, it is written in terms of the Fourier variables as

(2.39)\begin{align} \{K, H\} &= \int \left( \partial_{\eta_k} H\,\overline{\partial_{\zeta_k} K} - \partial_{\zeta_k} H\, \overline{\partial_{\eta_k} K} \right) {\rm d}k = \int \left( \partial_{\eta_k} H\,\partial_{\bar{\zeta}_k} K - \partial_{\zeta_k} H\,\partial_{\bar{\eta}_k} K \right) {\rm d}k \nonumber\\ & = \int \left( \partial_{\eta_k} H\,\delta_{\zeta_{{-}k}} K - \partial_{\zeta_k}\,H \partial_{\eta_{{-}k}} K \right) {\rm d}k = \int \left( \partial_{\eta_{k_1}} H\,\partial_{\zeta_{k_2}} K - \partial_{\zeta_{k_1}} H\,\partial_{\eta_{k_2}} K \right) \delta_{12} \,{\rm d}k_{12} . \end{align}

In terms of the complex symplectic coordinates, it becomes

(2.40)\begin{equation} \{K, H\} = \frac{1}{{\rm i}} \int (\partial_{z_{k_1}} H\,\partial_{\bar{z}_{{-}k_2}} K - \partial_{\bar{z}_{- k_1}} H\, \partial_{z_{k_2}} K)\,\delta_{12} \,{\rm d}k_{12} . \end{equation}

3.1. Third-order term in the Hamiltonian

The cubic term $\mathcal {H}^{(3)}$ given in (2.29) can be simplified further due to symmetries. We first state a simple identity that will be useful throughout the paper.

Lemma 3.1 For any $(k_1, k_2, k_3) \in \mathbb {R}^{3}$ with $k_1+k_2+k_3 = 0$ and $k_j \ne 0$, we have

(3.1)\begin{equation} {{\rm sgn}}(k_1) \,{{\rm sgn}}(k_2) + {{\rm sgn}}(k_1) \,{{\rm sgn}}(k_3) + {{\rm sgn}}(k_2) \,{{\rm sgn}}(k_3) ={-}1 . \end{equation}

Lemma 3.2 The cubic term $\mathcal {H}^{(3)}$ in (2.29) can be written as

(3.2)\begin{equation} \mathcal{H}^{(3)} ={-}\frac{1}{2\sqrt{2 {\rm \pi}}} \int (1+{{\rm sgn}}(k_1) \,{{\rm sgn}}(k_3)) \left( |k_1|\,|k_3|\,\zeta_1 \eta_2 \zeta_3 + {\rm i}\, \frac{\gamma k_2}{2}\,\eta_1 \zeta_2 \eta_3 \right) \delta_{123} \,{\rm d}k_{123} . \end{equation}

Proof. Expanding the brackets in (2.29), we regroup the terms as follows:

(3.3)\begin{align} \mathcal{H}^{(3)} &= \frac{1}{2\sqrt{2{\rm \pi}}} \int \left[ ({-}k_1k_3\!-\!|k_1|\,|k_3|) \zeta_1 \eta_2 \zeta_3\!+\! (k_1k_3+|k_1|\, |k_3|)\,\frac{{\rm i} \gamma}{k_1}\,\eta_1 \eta_2 \zeta_3 - {\rm i} \gamma k_1 \zeta_1 \eta_2 \eta_3 \right.\nonumber\\ & \quad \left. {}+ \left( (k_1k_3+|k_1|\,|k_3|)\,\frac{\gamma^{2}}{4k_1 k_3} - \frac{\gamma^{2}}{6} \right) \eta_1 \eta_2 \eta_3 \right] \delta_{123} \,{\rm d}k_{123} . \end{align}

The term on the second line of (3.3) corresponding to $\eta _1 \eta _2 \eta _3$ vanishes because of the symmetry in $\eta _1 \eta _2 \eta _3$ under index rearrangements, and identity (3.1). Indeed, we have

(3.4)\begin{align} \int (k_1k_3+|k_1|\,|k_3|)\,\frac{\gamma^{2}}{4k_1 k_3}\,\eta_1 \eta_2 \eta_3\,\delta_{123} \,{\rm d}k_{123} & = \frac{\gamma^{2}}{4} \int (1+{{\rm sgn}}(k_1) \,{{\rm sgn}}(k_3)) \eta_1 \eta_2 \eta_3\,\delta_{123} \,{\rm d}k_{123} \nonumber\\ & = \frac{\gamma^{2}}{6} \int \eta_1 \eta_2 \eta_3\,\delta_{123} \,{\rm d}k_{123} . \end{align}

As a result, $\mathcal {H}^{(3)}$ simplifies to

(3.5)\begin{align} \mathcal{H}^{(3)} &={-} \frac{1}{2\sqrt{2{\rm \pi}}} \int \left[ (1+{{\rm sgn}}(k_1) \,{{\rm sgn}}(k_3))\,|k_1|\,|k_3|\,\zeta_1 \eta_2 \zeta_3\right.\nonumber\\ &\quad \left. - {\rm i} \gamma\,|k_1| \,{{\rm sgn}}(k_3)\,\zeta_1 \eta_2 \eta_3 \right] \delta_{123} \,{\rm d}k_{123} .\end{align}

The first term in (3.5) identifies to the first term in (3.2), while its second term can be transformed using index rearrangements and identity (3.1) as follows:

(3.6)\begin{align} \int |k_1| \,{{\rm sgn}}(k_3)\,\zeta_1 \eta_2 \eta_3\,\delta_{123} \,{\rm d}k_{123} &= \int |k_2| \,{{\rm sgn}}(k_3)\,\eta_1 \zeta_2 \eta_3\,\delta_{123} \,{\rm d}k_{123} \nonumber\\ &= \int \frac{|k_2|}{2}\,({{\rm sgn}}(k_1) + {{\rm sgn}}(k_3)) \eta_1 \zeta_2 \eta_3\,\delta_{123} \,{\rm d}k_{123}\nonumber\\ &={-} \int \frac{k_2}{2}\left(1 + {{\rm sgn}}(k_1) \,{{\rm sgn}}(k_3) \right) \eta_1 \zeta_2 \eta_3\,\delta_{123} \,{\rm d}k_{123}. \end{align}

In terms of the complex symplectic coordinates, $\mathcal {H}^{(3)}$ is a linear combination of third-order monomials:

(3.7)\begin{align} \mathcal{H}^{(3)} &= \frac{1}{8 \sqrt{\rm \pi}} \int \frac{1+{{\rm sgn}}(k_1) \,{{\rm sgn}}(k_3)}{a_1 a_2 a_3} \left( \omega_1 \omega_3 - \frac{\gamma \omega_2}{2}\,{{\rm sgn}}(k_2) \right) \nonumber\\ &\quad \times \left( z_1 z_2 z_3 + \bar{z}_1 \bar{z}_2 \bar{z}_3 - 2( \bar{z}_{{-}1} \bar{z}_{{-}2} z_3 + z_{{-}1} z_{{-}2} \bar{z}_3 ) + \bar{z}_{{-}1} z_2 \bar{z}_{{-}3} + z_{{-}1} \bar{z}_2 z_{{-}3} \right) \delta_{123} \,{\rm d}k_{123} . \end{align}

3.2. Canonical transformations

We are looking for a canonical transformation of the physical variables

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

defined in a neighbourhood of the origin, such that the transformed Hamiltonian $\mathcal {H}'$satisfies

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

and reduces to

(3.10)\begin{equation} \mathcal{H}'(w') = \mathcal{H}^{(2)}(w') + Z^{(4)} + Z^{(5 )} + \cdots , \end{equation}

where each term $Z^{(m)}$ is of degree $m$, and all cubic terms are eliminated. We construct the transformation $\tau$ by the Lie transform method as a Hamiltonian flow $\phi$ from ‘time’ $s=-1$ to ‘time’ $s=0$ governed by

(3.11a–c)\begin{equation} \partial_s \phi = J \, \boldsymbol{\nabla} K(\phi) , \quad \phi(w')|_{s=0} = w' , \quad \phi(w')|_{s={-}1} = w , \end{equation}

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

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

Abusing notation, we further use $w = (\eta,\zeta )^\textrm {T}$ to denote the new variable $w'$. Terms in this expansion can be expressed using Poisson brackets as

(3.13)\begin{equation} \mathcal{H}(\phi(w))|_{s=0} = \mathcal{H}(w) , \end{equation}

(3.14)\begin{align} \frac{{\rm d} \mathcal{H}}{{\rm d}s}(\phi(w))|_{s=0} & = \int \left( \partial_\eta \mathcal{H}\,\partial_s \eta + \partial_\zeta \mathcal{H}\, \partial_s \zeta \right) {{\rm d}\kern0.06em x} \nonumber\\ & = \int \left( \partial_\eta \mathcal{H}\,\partial_\zeta K - \partial_\zeta \mathcal{H}\,\partial_\eta K \right) {{\rm d}\kern0.06em x} = \{K, \mathcal{H}\}(w) , \end{align}

and similarly for the remaining terms. The Taylor expansion of $\mathcal {H}'$ around $s=0$ now has the form

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

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

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

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

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

which would eliminate all cubic terms from the transformed Hamiltonian $\mathcal {H}'$. In the following, we will show that it is possible to do this, and that there are no resonant cubic terms in the Hamiltonian.

3.3. Third-order Birkhoff normal form

To find the auxiliary Hamiltonian $K^{(3)}$ from (3.17), we use the following diagonal property of the coadjoint operator $\textrm {coad}_{\mathcal {H}^{(2)}} := \{\cdot, \mathcal {H}^{(2)}\}$ when applied to monomial terms (Craig & Sulem Reference Craig and Sulem2016). For example, taking $\mathcal {I} := \int z_1 z_2 \bar {z}_{-3}\,\delta _{123} \,\textrm {d}k_{123},$ we have

(3.18)\begin{equation} \{ \mathcal{I}, \mathcal{H}^{(2)}\} = {\rm i} \int (\varOmega_1 + \varOmega_2 - \varOmega_{{-}3}) z_1 z_2 \bar{z}_{{-}3}\, \delta_{123} \,{\rm d}k_{123} , \end{equation}

where $\varOmega _{\pm j} := \varOmega _{\pm k_j}$.

Proposition 3.1 The cohomological equation (3.17) has a unique solution $K^{(3)}$ that in complex symplectic coordinates is

(3.19)\begin{align} K^{(3)} &= \frac{1}{8{\rm i} \sqrt{\rm \pi}} \int \frac{1+{{\rm sgn}}(k_1)\, {{\rm sgn}}(k_3)}{a_1 a_2 a_3} \left( \omega_1 \omega_3 - \frac{\gamma \omega_2}{2}\,{{\rm sgn}}(k_2) \right) \delta_{123} \nonumber\\ &\quad \times\left( \frac{z_1 z_2 z_3 - \bar{z}_1 \bar{z}_2 \bar{z}_3}{\varOmega_1 + \varOmega_2 + \varOmega_3} + 2\,\frac{\bar{z}_{{-}1} \bar{z}_{{-}2} z_3 - z_{{-}1} z_{{-}2} \bar{z}_3 }{\varOmega_{{-}1} + \varOmega_{{-}2} - \varOmega_3} - \frac{\bar{z}_{{-}1} z_2 \bar{z}_{{-}3} - z_{{-}1} \bar{z}_2 z_{{-}3}}{\varOmega_{{-}1} - \varOmega_2 + \varOmega_{{-}3}}\right) {\rm d}k_{123} . \end{align}

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

(3.20)\begin{align} K^{(3)} &= \frac{1}{4 {\rm i} \sqrt{2 {\rm \pi}}} \int (1+{{\rm sgn}}(k_1)\, {{\rm sgn}}(k_3))\,\delta_{123} \nonumber\\ &\quad \times \left[-\frac{\dfrac{\gamma^{4}}{8} + \dfrac{\gamma^{2}}{2}\,g\,|k_2| + g^{2}\,|k_1|\,|k_3|}{g^{2}\,|k_1|\, |k_3|}\right.\nonumber\\&\quad \left.\left( -\frac{\gamma}{2}\,{{\rm sgn}}(k_2)\,\eta_1 \eta_2 \eta_3 + {\rm i}\,|k_2|\,\eta_1 \zeta_2 \eta_3 - 2{\rm i}\,|k_3|\,\eta_1 \eta_2 \zeta_3 \right) \right.\nonumber\\ &\quad + \frac{\gamma \, {{\rm sgn}}(k_2) }{g^{2}\,|k_1|\,|k_3|} \left( 2 \omega_3^{2}\,|k_1|\,|k_2|\,\zeta_1 \zeta_2 \eta_3 - \omega_2^{2}\,|k_1|\,|k_3|\,\zeta_1 \eta_2 \zeta_3 \right.\nonumber\\ &\quad\left.\left. {}+ {\rm i} \frac{\gamma \, {{\rm sgn}}(k_2) }{2}\,|k_1|\,|k_2|\,|k_3|\,\zeta_1 \zeta_2 \zeta_3 \right) \vphantom{\frac{\dfrac{\gamma^{4}}{8} + \dfrac{\gamma^{2}}{2} g|k_2| + g^{2} |k_1| |k_3|}{g^{2}|k_1| |k_3|}}\right] {\rm d}k_{123} . \end{align}

Below, we write an alternative representation of the auxiliary Hamiltonian $K^{(3)}$ that will be useful later when we compute the Poisson bracket $\{ \mathcal {H}^{(3)}, K^{(3)} \}$. Introduce the coefficient functions

(3.21a,b)\begin{equation} S_{123} := \frac{1\!+\!{{\rm sgn}} (k_1)\, {{\rm sgn}} (k_3)}{a_1 a_2 a_3} \left( k_1 k_3 a_1^{2} a_3^{2} - \frac{\gamma}{2}\,k_2 a_2^{2} \right) , \quad A_{123} := \frac{1}{8\sqrt{\rm \pi}}\,(S_{123} \!+\! S_{312} - S_{231}) . \end{equation}

Note that $S_{123} = S_{321}$.

Lemma 3.3 We have

(3.22)\begin{equation} \left.\begin{gathered} \mathcal{H}^{(3)} = \int A_{123} (z_1 z_2 z_3 + \bar{z}_1 \bar{z}_2 \bar{z}_3 - z_{{-}1} z_{{-}2} \bar{z}_3 - \bar{z}_{{-}1} \bar{z}_{{-}2} z_3)\,\delta_{123} \,{\rm d}k_{123} , \\ K^{(3)} = \frac{1}{{\rm i}} \int \left[ \frac{A_{123}(z_1 z_2 z_3 - \bar{z}_{1} \bar{z}_{2} \bar{z}_{3})}{\varOmega_1 + \varOmega_2 + \varOmega_3} - \frac{A_{123}(z_{{-}1} z_{{-}2} \bar{z}_{3} - \bar{z}_{{-}1} \bar{z}_{{-}2} z_3)}{\varOmega_{{-}1} + \varOmega_{{-}2} - \varOmega_3} \right] \delta_{123} \,{\rm d}k_{123} . \end{gathered}\right\} \end{equation}

Proof. From (3.7), using that $k_1 k_3 = |k_1|\,|k_3|$ whenever $1+{\textrm {sgn}}(k_1)\, {\textrm {sgn}}(k_3) \neq 0$, we have

(3.23)\begin{align} \mathcal{H}^{(3)} &= \frac{1}{8 \sqrt{\rm \pi}} \int S_{123} \left( (z_1 z_2 z_3 + \bar{z}_1 \bar{z}_2 \bar{z}_3)- 2 (z_{{-}1} z_{{-}2} \bar{z}_3 + \bar{z}_{{-}1} \bar{z}_{{-}2} z_3 )\right. \nonumber\\ &\quad\left. {}+ (\bar{z}_{{-}1} z_2 \bar{z}_{{-}3} + z_{{-}1} \bar{z}_2 z_{{-}3}) \right) \delta_{123} \,{\rm d}k_{123} . \end{align}

By symmetry, the first term on the right-hand side in (3.23) becomes

(3.24)\begin{equation} \frac{1}{8 \sqrt{\rm \pi}} \int S_{123} (z_1 z_2 z_3 + \bar{z}_1 \bar{z}_2 \bar{z}_3)\,\delta_{123} \,{\rm d}k_{123} = \int A_{123} (z_1 z_2 z_3 + \bar{z}_1 \bar{z}_2 \bar{z}_3)\,\delta_{123} \,{\rm d}k_{123} . \end{equation}

The second term in (3.23) has two copies of $(z_{-1} z_{-2} \bar {z}_3 + \bar {z}_{-1} \bar {z}_{-2} z_3 )$. We keep one copy as it is, and apply two index rearrangements to the second copy. First, $(k_1,k_2,k_3)\to (k_3,k_2,k_1)$ gives

(3.25)\begin{equation} \int S_{123} (z_{{-}1} z_{{-}2} \bar{z}_3 + \bar{z}_{{-}1} \bar{z}_{{-}2} z_3)\,\delta_{123}\,{\rm d}k_{123} = \int S_{123}(\bar{z}_{1} z_{{-}2} z_{{-}3} + z_{1} \bar{z}_{{-}2} \bar{z}_{{-}3})\,\delta_{123}\,{\rm d}k_{123} , \end{equation}

where we use that $S_{123} = S_{321}$. Second, $(1,2,3) \to (3,1,2)$ implies

(3.26)\begin{equation} \int S_{123} (\bar{z}_{1} z_{{-}2} z_{{-}3} + z_{1} \bar{ z}_{{-}2} \bar{z}_{{-}3})\,\delta_{123}\,{\rm d}k_{123} = \int S_{312} ( z_{{-}1} z_{{-}2} \bar{z}_{3} + \bar{ z}_{{-}1} \bar{z}_{{-}2} z_3)\,\delta_{123}\,{\rm d}k_{123} . \end{equation}

For the third term in (3.23), we apply $(k_1,k_2,k_3)\to (k_2,k_3,k_1)$ and leave details to the reader. Combining all the transformations above, we obtain

(3.27)\begin{align} &\frac{1}{8 \sqrt{\rm \pi}} \int S_{123} \left(- 2(z_{{-}1} z_{{-}2} \bar{z}_3 + \bar{z}_{{-}1} \bar{z}_{{-}2} z_3 ) + \bar{z}_{{-}1} z_2 \bar{z}_{{-}3} + z_{{-}1} \bar{z}_2 z_{{-}3} \right) \delta_{123} \,{\rm d}k_{123} \nonumber\\ &\quad = \int A_{123}({-}z_{{-}1} z_{{-}2} \bar{z}_3 - \bar{z}_{{-}1} \bar{z}_{{-}2} z_3 )\,\delta_{123} \,{\rm d}k_{123} , \end{align}

which together with (3.24) implies the first relation of (3.22). The second equation in (3.22) is derived from (3.17) by using the diagonal properties (3.18).

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

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

with initial conditions at $s= -1$ being the original variables.

To write $K^{(3)}$ in the physical space, it is convenient to introduce $\tilde {\eta } = {{\rm H}} \eta$, $\tilde {\zeta } = {{\rm H}} \zeta$, where ${{\rm H}} = - \textrm {i}\, {\textrm {sgn}}(D)$ is the Hilbert transform. In the following, we will use the identities ${{\rm H}}\,|D| = -\partial _x$, ${{\rm H}}\,|D|^{-1} = \partial _x^{-1}$, and similar ones.

Proposition 3.2 The auxiliary Hamiltonian $K^{(3)}$ in (3.20) can be written in the physical space as

(3.29)\begin{align} K^{(3)} &= \int \left[ \frac{1}{2}\,\tilde{\eta}^{2}\,\partial_x \tilde{\zeta} - \frac{\gamma}{4g}\,(g \eta^{2} \tilde{\eta} - \zeta^{2}\,\partial_x \eta +2 \zeta \tilde{\eta}\,\partial_x \tilde{\zeta}) \right.\nonumber\\ &\quad -\frac{\gamma^{2}}{4g^{2}} \left(g (\partial_x^{{-}1} \eta) (\partial_x \eta) \zeta + g\eta \tilde{\eta}\tilde{\zeta} - g \tilde{\eta} (\partial_x \tilde{\zeta})\,\partial_x^{{-}1} \eta - \frac{1}{2}\,\zeta^{2}\,\partial_x \tilde{\zeta}\right)\nonumber\\ &\quad - \frac{\gamma^{3}}{16g^{2}} \left( \zeta^{2} \tilde{\eta} +2 \zeta(\partial_x \tilde{\zeta}) \partial_x^{{-}1}\,\eta -2g (\partial_x^{{-}1} \eta) (\partial_x \tilde{\eta}) (\partial_x^{{-}1}\tilde{\eta}) \right) \nonumber\\ &\quad \left. - \frac{\gamma^{4}}{16g^{2}}\,(\eta \tilde{\zeta} -\tilde{\eta} \zeta)\,\partial_x^{{-}1} \eta - \frac{\gamma^{5}}{64g^{2}}\,\tilde{\eta} (\partial_x^{{-}1} \eta)^{2} \right] {{\rm d}\kern0.06em x} . \end{align}

Proof. The main idea is to expand the brackets in the expression (3.20) for $K^{(3)}$, identify terms and combine them appropriately. We decompose $K^{(3)}$ as

(3.30)\begin{equation} K^{(3)} = {\rm I} + {\rm II} + {\rm III} + {\rm IV} , \end{equation}

where $\textrm {I}$ is the part of $K^{(3)}$ associated with $\eta \eta \eta$-type terms, namely

(3.31)\begin{align} {\rm I} &= \frac{1}{4 {\rm i} \sqrt{2 {\rm \pi}}} \int (1+{{\rm sgn}}(k_1)\, {{\rm sgn}}(k_3)) \left( \frac{\dfrac{\gamma^{4}}{8} + \dfrac{\gamma^{2}}{2}\,g\,|k_2| + g^{2}\,|k_1|\,|k_3|}{g^{2}\,|k_1|\,|k_3|} \right) \frac{\gamma}{2}\nonumber\\ &\quad\times {{\rm sgn}}(k_2)\,\eta_1 \eta_2 \eta_3\,\delta_{123} \,{\rm d}k_{123} . \end{align}

The part $\textrm {II}$ is associated with $\eta \eta \zeta$-type terms, $\textrm {III}$ is associated with $\eta \zeta \zeta$-type terms, and $\textrm {IV}$ is associated with $\zeta \zeta \zeta$-type terms. Below, we give the computations for $\textrm {I}$. The remaining terms can be computed in a similar way. We further decompose $\textrm {I}$ into a sum of three terms $\textrm {I} = \textrm {I}_1 + \textrm {I}_2 + \textrm {I}_3$ based on the power of $\gamma$ involved:

(3.32)\begin{equation} \left.\begin{gathered} {\rm I}_1 = \frac{\gamma^{5}}{64{\rm i} g^{2} \sqrt{2{\rm \pi}}}\int (1+{{\rm sgn}}(k_1)\, {{\rm sgn}}(k_3))\, \frac{1}{|k_1|\,|k_3|}\, {{\rm sgn}}(k_2)\,\eta_1 \eta_2 \eta_3\,\delta_{123} \,{\rm d}k_{123} , \\ {\rm I}_2 = \frac{\gamma^{3}}{16 {\rm i} g \sqrt{2{\rm \pi}}}\int (1+{{\rm sgn}}(k_1)\, {{\rm sgn}}(k_3))\, \frac{|k_2|}{|k_1|\,|k_3|}\, {{\rm sgn}}(k_2)\,\eta_1 \eta_2 \eta_3\,\delta_{123} \,{\rm d}k_{123} , \\ {\rm I}_3 = \frac{\gamma}{8{\rm i} \sqrt{2{\rm \pi}}}\int (1+{{\rm sgn}}(k_1)\, {{\rm sgn}}(k_3))\, {{\rm sgn}}(k_2)\,\eta_1 \eta_2 \eta_3\,\delta_{123} \,{\rm d}k_{123} . \end{gathered}\right\} \end{equation}

Writing ${\textrm {sgn}}(k_2) = k_2/|k_2|$ in $\textrm {I}_1$, we get

(3.33)\begin{align} {\rm I}_1& = \frac{\gamma^{5}}{64{\rm i}g^{2} \sqrt{2{\rm \pi}}}\int \frac{k_2}{|k_1|\,|k_2|\,|k_3|}\,\eta_1 \eta_2 \eta_3\,\delta_{123} \,{\rm d}k_{123} \nonumber\\ & \quad + \frac{\gamma^{5}}{64{\rm i}g^{2} \sqrt{2{\rm \pi}}}\int \frac{{{\rm sgn}}(k_1)\,{{\rm sgn}}(k_2)\, {{\rm sgn}}(k_3)}{|k_1|\,|k_3|}\, \eta_1 \eta_2 \eta_3\,\delta_{123} \,{\rm d}k_{123} = {\rm I}_1^{(1)} + {\rm I}_1^{(2)} . \end{align}

Using the index rearrangements $(k_1,k_2,k_3) \to (k_2,k_1,k_3)$ and $(k_1,k_2,k_3)\to (k_1,k_3,k_2)$, we see that $\textrm {I}_1^{(1)} = 0$ since the integration is over $k_1+k_2+k_3=0$. The part $\textrm {I}_1$ thus reduces to

(3.34)\begin{equation} {\rm I}_1 = \frac{\gamma^{5}}{64{\rm i} g^{2} \sqrt{2{\rm \pi}}} \int \frac{{{\rm sgn}}(k_2)}{k_1 k_3}\, \eta_1 \eta_2 \eta_3\,\delta_{123} \,{\rm d}k_{123} ={-} \frac{\gamma^{5}}{64 g^{2}} \int (\partial_x^{{-}1} \eta)^{2} \tilde{\eta} \, {{\rm d}\kern0.06em x} . \end{equation}

Similar steps imply

(3.35a,b)\begin{equation} {\rm I}_2 = \frac{\gamma^{3}}{8g} \int (\partial_x^{{-}1} \eta) (\partial_x \tilde{\eta}) (\partial_x^{{-}1} \eta) \, {{\rm d}\kern0.06em x} , \quad {\rm I}_3 ={-}\frac{\gamma}{4} \int \eta^{2} \tilde{\eta}\, {{\rm d}\kern0.06em x} . \end{equation}

Combining (3.34)–(3.35a,b),

(3.36)\begin{equation} {\rm I} ={-}\frac{\gamma}{4} \int \eta^{2} \tilde{\eta} \, {{\rm d}\kern0.06em x} + \frac{\gamma^{3}}{8g} \int (\partial_x^{{-}1} \eta) (\partial_x\tilde{\eta}) (\partial_x^{{-}1} \tilde{\eta}) \, {{\rm d}\kern0.06em x} - \frac{\gamma^{5}}{64 g^{2}} \int (\partial_x^{{-}1} \eta)^{2} \tilde{\eta}\, {{\rm d}\kern0.06em x} , \end{equation}

which are the second, eleventh and fourteenth terms in (3.29). Terms $\textrm {II}$, $\textrm {III}$ and $\textrm {IV}$ are calculated in a similar fashion and identify to the remaining terms in (3.29).

We now apply the variational derivatives $\partial _\eta$ and $\partial _\zeta$ on $K^{(3)}$ to obtain the third-order normal form transformation defining the Hamiltonian flow (3.28).

Proposition 3.3 The Hamiltonian system that defines the third-order normal form transformation has the form of a system of two partial differential equations:

(3.37)\begin{align} \partial_s \eta & = \partial_{\zeta} K^{(3)} = \frac{1}{2}\,{{\rm H}}\,\partial_x \tilde{\eta}^{2} + \frac{\gamma}{2g} \left( \zeta\,\partial_x \eta - \tilde{\eta}\,\partial_x \tilde{\zeta} - |D|\,(\zeta \tilde{\eta}) \right) \nonumber\\ & \quad +\frac{\gamma^{2}}{4g^{2}} \left( \zeta\,\partial_x \tilde{\zeta} + \frac{1}{2}\,|D|\,\zeta^{2} \right)+ \frac{\gamma^{2}}{4g} \left(-(\partial_x^{{-}1} \eta) (\partial_x \eta)+ {{\rm H}}( \eta \tilde{\eta} )+ |D|\,(\tilde{\eta}\, \partial_x^{{-}1} \eta)\right) \nonumber\\ &\quad- \frac{\gamma^{3}}{8g^{2}} \left( \zeta \tilde{\eta} + (\partial_x \tilde{\zeta}) (\partial_x^{{-}1} \eta) + |D|\,(\zeta\, \partial_x^{{-}1} \eta)\right) + \frac{\gamma^{4}}{16g^{2}} \left( {{\rm H}} ( \eta\,\partial_x^{{-}1} \eta) + \tilde{\eta}\,\partial_x^{{-}1} \eta \right), \end{align}

(3.38)\begin{align} \partial_s \zeta & ={-}\partial_{\eta} K^{(3)} ={{\rm H}} ( \tilde{\eta}\,\partial_x \tilde{\zeta} ) + \frac{\gamma}{2} \left( \eta \tilde{\eta} - \frac{1}{2}\,{{\rm H}} (\eta^{2}) \right)+ \frac{\gamma}{2g} \left( \zeta\,\partial_x \zeta - {{\rm H}} (\zeta\,\partial_x \widetilde \zeta) \right)\nonumber\\ & \quad - \frac{\gamma^{2}}{4g}\left( \partial_x (\zeta\,\partial_x^{{-}1} \eta) + \partial_x^{{-}1} (\zeta\,\partial_x \eta) - \tilde{\zeta} \tilde{\eta} + {{\rm H}}( \eta \tilde{\zeta})\right.\nonumber\\ &\quad\left. - {{\rm H}} \left( (\partial_x \tilde{\zeta}) (\partial_x^{{-}1} \eta) \right) - \partial_x^{{-}1} \left( (\partial_x \tilde{\zeta}) \tilde{\eta} \right) \right) \nonumber\\ &\quad + \frac{\gamma^{3}}{8g} \left( \partial_x^{{-}1} \left( (\partial_x \tilde{\eta}) ( \partial_x^{{-}1} \tilde{\eta}) \right) - |D| \left( (\partial_x^{{-}1} \eta) (\partial_x^{{-}1} \tilde{\eta}) \right) + |D|^{{-}1} \left( (\partial_x \tilde{\eta}) (\partial_x^{{-}1} \eta) \right) \right) \nonumber\\ & \quad - \frac{\gamma^{3}}{16g^{2}} \left({{\rm H}} (\zeta^{2}) + 2\,\partial_x^{{-}1} (\zeta\,\partial_x \tilde{\zeta})\right)\nonumber\\ &\quad + \frac{\gamma^{4}}{16g^{2}} \left( \tilde{\zeta}\,\partial_x^{{-}1} \eta - \partial_x^{{-}1} (\eta\tilde{\zeta}) + {{\rm H}} (\zeta\, \partial_x^{{-}1} \eta) + \partial_x^{{-}1} (\zeta \tilde{\eta}) \right)\nonumber\\ &\quad -\frac{\gamma^{5}}{64g^{2}} \left( {{\rm H}} \left( (\partial_x^{{-}1} \eta)^{2} \right) + 2\,\partial_x^{{-}1} \left( \tilde{\eta}\, \partial_x^{{-}1} \eta \right) \right) .\end{align}

In the absence of vorticity ($\gamma = 0$), the equation for $\eta$ simplifies to $\partial _s \eta = \frac {1}{2} {{\rm H}}\,\partial _x ( \tilde {\eta })^{2}$, or equivalently, to the inviscid Burgers equation for $\tilde {\eta }$,

(3.39)\begin{equation} \partial_s \tilde{\eta} + \tilde{\eta}\,\partial_x \tilde{\eta}=0 , \end{equation}

as obtained by Craig & Sulem (Reference Craig and Sulem2016), while $\tilde {\zeta }$ satisfies

(3.40)\begin{equation} \partial_s \tilde{\zeta} + \tilde{\eta}\,\partial_x \tilde{\zeta}=0 , \end{equation}

which is its linearization along the Burgers flow. These equations were tested in Craig et al. (Reference Craig, Guyenne and Sulem2021a) and Guyenne et al. (Reference Guyenne, Kairzhan, Sulem and Xu2021, Reference Guyenne, Kairzhan and Sulem2022) in the context of irrotational gravity waves on deep water.

Due to the complexity of the formulas for $\partial _{\eta } K^{(3)}$ and $\partial _{\zeta } K^{(3)}$, we have checked a posteriori that the cohomology equation (3.17) is indeed satisfied.

6.1. Hamiltonian in the physical variables

Lemma 6.1 In the physical variables $(u, \bar {u})$, the Hamiltonian $\mathcal {H}$ in (6.1) reads

(6.2)\begin{align} \mathcal{H} &= \varepsilon \int \bar{u}\,\varOmega (k_0 + \varepsilon D_X) u \, {\rm d}X + \varepsilon^{3}\,\frac{k_0^{3} (\omega_0 - \gamma)(\gamma^{2} + 4 \omega_0^{2})}{4 \omega_0 \varOmega_0 (2 \omega_0 - \gamma)}\int|u|^{4} \,{\rm d}X\nonumber\\ &\quad + \varepsilon^{4}\, \frac{\beta}{2} \int |u|^{2} \,{\rm Im} (\bar{u} \partial_X u) \,{\rm d}X - \varepsilon^{4}\,\frac{k_0^{2} \omega_0^{2}}{2 \varOmega_0^{2}} \int |u|^{2}\,|D_X|\,|u|^{2} \,{\rm d}X + {{{O}}} (\varepsilon^{5}) , \end{align}

where $D_X = -\mathrm{i} \partial _X$ in the slow spatial variable $X$.

Proof. The first term in (6.2) comes out by applying the change of variables (5.2a,b) to the quadratic Hamiltonian

(6.3)\begin{equation} \mathcal{H}^{(2)} = \int \varOmega(k_0+\varepsilon \lambda)\,|z(k_0 +\varepsilon \lambda)|^{2} \,{\rm d}k= \varepsilon \int \bar{u} \, \varOmega(k_0+\varepsilon D_X) u \, {\rm d}X , \end{equation}

where we use that $u(X) = ({1}/{\sqrt {2{\rm \pi} }} )\int \textrm {e}^{\textrm {i}\lambda X}\,U(\lambda ) \,\textrm {d}\lambda$ is the inverse Fourier transform of $U(\lambda )$. Furthermore,

(6.4)\begin{equation} \int U_1 U_2 \bar{U}_3 \bar{U}_4\,\delta_{1+2-3-4} \,{\rm d}\lambda_{1234} = 2 {\rm \pi}\int |u|^{4} \,{\rm d} X . \end{equation}

Similarly,

(6.5)\begin{equation} \int (\lambda_2 + \lambda_3) U_1 U_2 \bar{U}_3 \bar{U}_4\,\delta_{1+2-3-4} \,{\rm d}\lambda_{1234} = 4 {\rm \pi}\int |u|^{2} \,{\rm Im}(\bar{u}\partial_X u) \, {\rm d}X , \end{equation}

while the remaining term of $\mathcal {H}_+^{(4)}$ identifies to

(6.6)\begin{equation} \int |\lambda_1-\lambda_3|\,U_1 U_2 \bar{U}_3 \bar{U}_4\,\delta_{1+2-3-4} \,{\rm d}\lambda_{1234} = 2 {\rm \pi}\int |u|^{2}\,|D_X|\,|u|^{2} \, {\rm d}X . \end{equation}

Expanding the linear dispersion relation in (6.2) around $k_0$ as

(6.7)\begin{equation} \varOmega (k_0+ \varepsilon D_X) = \varOmega_0 + \frac{g D_X}{2 \omega_0}\,\varepsilon - \frac{g^{2} D_X^{2}}{8 \omega_0^{3}}\,\varepsilon^{2} + \frac{g^{3} D_X^{3}}{16 \omega_0^{5}}\,\varepsilon^{3} + {{{O}}}(\varepsilon^{4}) \end{equation}

gives an alternative form of the Hamiltonian $\mathcal {H}$ in the physical variables:

(6.8)\begin{align} \mathcal{H} &= \int \left[\varepsilon\, \varOmega_0 \, |u|^{2} + \varepsilon^{2}\,\frac{g}{2\omega_0}\,{\rm Im}(\overline u \partial_X u) \right.\nonumber\\ & \quad - \varepsilon^{3}\,\frac{g^{2}}{8 \omega_0^{3}}\,|\partial_X u|^{2} + \varepsilon^{3}\,\frac{k_0^{3} (\omega_0 - \gamma)(\gamma^{2} + 4 \omega_0^{2})}{4 \omega_0 \varOmega_0 (2 \omega_0 - \gamma)}\,|u|^{4} - \varepsilon^{4}\,\frac{g^{3}}{16 \omega_0^{5}}\,{\rm Im}(\overline u \partial_X^{3} u) \nonumber\\ & \quad\left. {}+ \varepsilon^{4}\, \frac{\beta}{2}\,|u|^{2}\, {\rm Im} (\bar{u} \partial_X u) - \varepsilon^{4}\,\frac{k_0^{2} \omega_0^{2}}{2 \varOmega_0^{2}}\,|u|^{2} |D_X|\,|u|^{2} \right] {\rm d}X + {{{O}}}(\varepsilon^{5}) . \end{align}

6.2. Derivation of the Dysthe equation

Using the relation (5.3) rewritten as

(6.9)\begin{equation} \begin{pmatrix} u \\ \bar{u} \end{pmatrix} = P_2 \begin{pmatrix} z \\ \bar{z} \end{pmatrix} = \varepsilon^{{-}1} \begin{pmatrix} {\rm e}^{-{\rm i} k_0 x} & 0 \\ 0 & {\rm e}^{{\rm i} k_0 x} \end{pmatrix} \begin{pmatrix} z \\ \bar{z} \end{pmatrix} , \end{equation}

the Hamiltonian system (2.32) takes the form

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

where $J_2 = \varepsilon P_2 J_1 P_2^{*}$ (Craig, Guyenne & Sulem Reference Craig, Guyenne and Sulem2010). The additional factor $\varepsilon$ in the definition of $J_2$ reflects the change in symplectic structure associated with the spatial rescaling $X = \varepsilon x$.

Substituting the reduced Hamiltonian (6.8) into (6.10), we get

(6.11)\begin{align} {\rm i}\,\partial_t u & = \varepsilon^{{-}1} \partial_{\bar{u}} \mathcal{H} \nonumber\\ & = \varOmega_0 u -{\rm i}\varepsilon\,\frac{g}{2 \omega_0}\,\partial_X u + \varepsilon^{2} \frac{g^{2}}{8 \omega_0^{3}}\,\partial_X^{2} u + \varepsilon^{2}\,\frac{k_0^{3} (\omega_0 - \gamma)(\gamma^{2} + 4 \omega_0^{2})}{2 \omega_0 \varOmega_0 (2 \omega_0 - \gamma)}\,|u|^{2} u \nonumber\\ & \quad + {\rm i}\varepsilon^{3}\,\frac{g^{3}}{16 \omega_0^{5}}\,\partial_X^{3} u - {\rm i} \varepsilon^{3} \beta\,|u|^{2}\, \partial_X u - \varepsilon^{3}\,\frac{k_0^{2} \omega_0^{2}}{\varOmega_0^{2}}\,u\,|D_X|\,|u|^{2} , \end{align}

which is a Hamiltonian Dysthe equation for two-dimensional gravity waves on deep water with constant vorticity. It describes modulated waves moving in the positive $x$-direction at group velocity $\varOmega _0' = \partial _{k} \varOmega (k_0) = g/(2 \omega _0)$ as indicated by the advection term. The non-local term $u\,|D_X|\,|u|^{2}$ is a signature of the Dysthe equation, which reflects the presence of the wave-induced mean flow. The coefficient $\beta$ is given by (5.18). The coefficient of the nonlinear term $|u|^{2} u$ above agrees with that in the NLS equation of Thomas et al. (Reference Thomas, Kharif and Manna2012) and in the Dysthe equation of Curtis et al. (Reference Curtis, Carter and Kalisch2018) (see (7.15) below) up to scaling factors consistent with the difference in definition for the wave envelope described in the various models.

The first two terms on the right-hand side of (6.11) 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 ${{\mathcal {H}}}$ of a multiple of the momentum (2.19), which reduces to

(6.12)\begin{equation} I = \int \eta\,\partial_x\zeta \,{{\rm d}\kern0.06em x} = \varepsilon \int \left[ k_0\,|u|^{2} + \varepsilon \operatorname{Im}(\bar{u}\,\partial_{X} u) \right] {\rm d}X, \end{equation}

while the former is equivalent to subtraction from ${{\mathcal {H}}}$ of a multiple of the wave action

(6.13)\begin{equation} M = \varepsilon \int |u|^{2} \, {\rm d}X , \end{equation}

which is conserved due to the phase-invariance property of the Dysthe equation. Because $I$ and $M$ Poisson commute with ${{\mathcal {H}}}$, this transformation preserves the symplectic structure $J_2$ (Craig et al. Reference Craig, Guyenne and Sulem2021b). The resulting Hamiltonian is given by

(6.14)\begin{equation} \hat{ {{\mathcal{H}}}} = {{\mathcal{H}}} - \varOmega_0' I - \left( \varOmega_0 - {k}_0 \varOmega_0' \right) 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.15)\begin{align} {\rm i} \, \partial_\tau u & = \frac{g^{2}}{8 \omega_0^{3}}\,\partial_X^{2} u + \frac{k_0^{3} (\omega_0 - \gamma)(\gamma^{2} + 4 \omega_0^{2})}{2 \omega_0 \varOmega_0 (2 \omega_0 - \gamma)}\,|u|^{2} u + {\rm i}\varepsilon\,\frac{g^{3}}{16 \omega_0^{5}}\, \partial_X^{3} u \nonumber\\ & \quad - {\rm i} \varepsilon \beta\,|u|^{2}\,\partial_X u - \varepsilon\,\frac{k_0^{2} \omega_0^{2}}{\varOmega_0^{2}}\,u\,|D_X|\,|u|^{2} . \end{align}

This governs the long-time evolution of the envelope of modulated waves in a reference frame moving in the positive horizontal direction at group velocity $\varOmega _0'$. The non-local operator $|D_X|$ is the Fourier multiplier with symbol $|\lambda |$. The associated Hamiltonian reads more explicitly (after multiplying by $\varepsilon ^{-3}$ and dropping the hat)

(6.16)\begin{align} \mathcal{H} &= \int \left[ - \frac{g^{2}}{8 \omega_0^{3}}\,|\partial_X u|^{2} + \frac{k_0^{3} (\omega_0 - \gamma)(\gamma^{2} + 4 \omega_0^{2})}{4 \omega_0 \varOmega_0 (2 \omega_0 - \gamma)}\, |u|^{4} - \varepsilon\,\frac{g^{3}}{16 \omega_0^{5}}\,{\rm Im} (\bar{u}\,\partial_X^{3} u)\right. \nonumber\\ & \quad \left.{}+ \varepsilon\, \frac{\beta}{2}\,|u|^{2} \,{\rm Im} (\bar{u}\,\partial_X u) - \varepsilon\,\frac{k_0^{2} \omega_0^{2}}{2 \varOmega_0^{2}}\, |u|^{2}\,|D_X|\,|u|^{2} \right] {\rm d}X . \end{align}

As suggested by Trulsen et al. (Reference Trulsen, Kliakhandler, Dysthe and Velarde2000), retaining the exact linear dispersion relation, rather than expanding it in powers of $\varepsilon$, may provide an overall better approximation of the wave envelope. On a related note, Obrecht & Saut (Reference Obrecht and Saut2015) proposed a full-dispersion Davey–Stewartson system and compared its analytical properties to those of the classical version. In the present context, the Dysthe equation with full linear dispersion takes the form

(6.17)\begin{align} {\rm i}\,\partial_t u &= \varOmega(k_0+\varepsilon D_X) u + \varepsilon^{2}\,\frac{k_0^{3} (\omega_0 - \gamma)(\gamma^{2} + 4 \omega_0^{2})}{2 \omega_0 \varOmega_0 (2 \omega_0 - \gamma)}\,|u|^{2} u\nonumber\\ &\quad - {\rm i} \varepsilon^{3} \beta\,|u|^{2}\,\partial_X u - \varepsilon^{3}\,\frac{k_0^{2} \omega_0^{2}}{\varOmega_0^{2}}\,u\,|D_X|\,|u|^{2} , \end{align}

and the corresponding Hamiltonian is given by (6.2).

7.1. Stability of Stokes waves

We first give the theoretical prediction for modulational or Benjamin–Feir (BF) instability of Stokes waves. These are represented by the exact uniform solution

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

for (6.11), where $B_0$ is a positive real constant, and

(7.2)\begin{equation} \beta_0 = \frac{k_0^{3} (\omega_0 - \gamma)(\gamma^{2} + 4 \omega_0^{2})}{2 \omega_0 \varOmega_0 (2 \omega_0 - \gamma)} . \end{equation}

In the irrotational case ($\gamma = 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.11) about $u_0$ by inserting a perturbation of the form

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

where

(7.4)\begin{equation} B(X,t) = B_1 \exp({\sigma t + {\rm i} \lambda X}) + B_2 \exp({\bar{\sigma} t - {\rm i} \lambda X} ), \end{equation}

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

(7.5)\begin{equation} \alpha = \frac{g^{2}}{8 \omega_0^{3}}\,\lambda^{2} \varGamma > 0 , \end{equation}

with

(7.6a,b)\begin{equation} \varGamma = 2 B_0^{2} \left( \beta_0 - \varepsilon \beta_3\,|\lambda| \right) - \frac{g^{2}}{8 \omega_0^{3}}\, \lambda^{2} , \quad \beta_3 = \frac{k_0^{2} \omega_0^{2}}{\varOmega_0^{2}} . \end{equation}

This is a tedious but straightforward calculation, for which we skip the details. Similar calculations can be found in Curtis et al. (Reference Curtis, Carter and Kalisch2018), Dysthe (Reference Dysthe1979) and Gramstad & Trulsen (Reference Gramstad and Trulsen2011).

Figure 1 depicts the normalized growth rate

(7.7)\begin{equation} \frac{|\textrm{Re}(\sigma)|}{\omega_0} = \frac{\sqrt{\alpha}}{\omega_0} , \end{equation}

delimiting the instability region as predicted by condition (7.5) for $(B_0,k_0) = (0.002,10)$ and various values of $\gamma$. The growth rate (and instability region) for $\gamma = 0$ is also included as a reference. Hereafter, all the variables are rescaled to absorb $\varepsilon$ back into their definition, and all the equations are non-dimensionalized by using $1/k_0$ and $1/\sqrt {g k_0}$ as characteristic length and time scales, respectively, so that $g = 1$. For convenience, we retain the same notations for all the dimensionless quantities. We set $\varepsilon = k_0 A_0$ (surface wave steepness), noting that the envelope amplitude $B_0$ and the surface amplitude $A_0$ are related via

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

according to (2.31) and (5.3). The graphs in figure 1 correspond to a wave steepness of about $\varepsilon = 0.05$. Clearly, the vorticity $\gamma$ (both its magnitude and sign) has an influence on (7.5). We see that $\gamma < 0$ tends to enhance the instability by amplifying the growth rate and enlarging the instability region to higher sideband wavenumbers. On the other hand, $\gamma > 0$ tends to diminish it. Figure 1 even suggests that for sufficiently large $\gamma > 0$, instability no longer occurs. This is confirmed by figure 2, which shows that the factor $\varGamma$ in (7.5) is no longer positive at any wavenumber $\lambda$ when $\gamma > 3.5$ for $(B_0,k_0) = (0.002,10)$. A positive vorticity (co-propagating current) therefore has a stabilizing effect on the dynamics of Stokes waves.

Figure 1. Regions of modulational instability according to (7.5) for $(B_0,k_0) = (0.002,10)$. The blue curves correspond to $\gamma \neq 0$, while the red curve corresponds to $\gamma = 0$: (a) $\gamma = \{ -0.5, -1, -1.5, -2, -3.5 \}$ (expanding curves with decreasing $\gamma$); (b) $\gamma = \{ +0.5, +1, +1.5, +2, +3.5 \}$ (shrinking curves with increasing $\gamma$).

Figure 2. Plots of $\varGamma$ versus $\lambda$ for $(B_0,k_0) = (0.002,10)$ and $\gamma = \{ 0, +1, +2, +3, +4, +5 \}$ (falling curves with increasing $\gamma$). The red curve corresponds to $\gamma = 0$.

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 (3.37)–(3.38) backward from $s = 0$ to $s = -1$, with ‘initial’ conditions given by the transformed variables

(7.9)$$\begin{gather} \left. \eta(x,t) \right|_{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.10)$$\begin{gather}\left.\zeta(x,t) \right|_{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.35{a},{b}) and (5.3). In these expressions, $u$ obeys (6.11) and $a^{-1}(D) = \sqrt {|D|/\omega (D)}$. The final solution at $s = -1$ represents the original variables $(\eta,\zeta )$. Starting from the first harmonics (with carrier wavenumber $k_0$) in the initial conditions (7.9)–(7.10), the evolutionary process in $s$ will generate automatically the next-order contributions from lower and higher harmonics via nonlinear interactions according to (3.37)–(3.38). Recall also that the non-canonical velocity potential $\xi$ can be recovered from the canonical one $\zeta$ via the direct relation (2.20).

7.3. Simulations and comparisons

For the comparison, the full nonlinear system (2.15) is solved numerically following a high-order spectral approach (Craig & Sulem Reference Craig and Sulem1993). The corresponding equations read more explicitly

(7.11)\begin{gather} \partial_t \eta = G(\eta)\,\xi + \gamma \eta\,\partial_x \eta , \end{gather}

(7.12)\begin{gather}\partial_t \xi ={-}g \eta - \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}} + \gamma \eta\,\partial_x \xi + \gamma\,\partial_x^{{-}1} G(\eta)\,\xi . \end{gather}

These are discretized in space by a pseudo-spectral method based on the fast Fourier transform (FFT). The computational domain is taken to be $0 \leqslant x \leqslant 2{\rm \pi}$ with periodic boundary conditions, and is divided into a regular mesh of $N$ collocation points. The Dirichlet–Neumann operator is computed via its series expansion (2.26) for which a small number $m$ of terms is sufficient to achieve highly accurate results by virtue of its analyticity properties. The value $m = 6$ is selected based on previous extensive tests (Xu & Guyenne Reference Xu and Guyenne2009). Time integration of (7.11) and (7.12) 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 (Reference Guyenne2017, Reference Guyenne2018).

The same numerical methods are applied to the envelope equation (6.11), as well as to the reconstruction procedure, with the same resolutions in space and time. In particular, the auxiliary system (3.37)–(3.38) is integrated in $s$ by using the same step size $\Delta s = \Delta t$. While this system of equations may look complicated, its numerical treatment is straightforward and efficient via the FFT. Moreover, because this computation is not performed at each instant $t$ (only when data on $\eta$ are required), and because it is performed over a short interval $-1 \leqslant s \leqslant 0$, the associated cost is insignificant. Note that by virtue of the zero-mass assumption (2.25), indetermination at $k = 0$ in the evaluation of any quantity involving a Fourier multiplier such as $\partial _x^{-1}$ or $|D|^{-1}$ may be lifted by simply setting its zeroth-mode component to zero.

To examine the stability of Stokes waves in the presence of a shear current, initial conditions of the form

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

are specified for (6.11), where $\lambda$ denotes the wavenumber of some long-wave perturbation. For the purpose of comparing with the full system (7.11)–(7.12), initial conditions $\eta (x,0)$ and $\xi (x,0)$ are reconstructed by solving (3.37)–(3.38) from transformed initial data (7.9)–(7.10) given in terms of (7.13).

The following tests focus on the case $(B_0,k_0,\lambda ) = (0.002,10,1)$ as considered in the previous stability analysis. The spatial and temporal resolutions are set to $\Delta x = 0.012$ ($N = 512$) and $\Delta t = 0.005$. Figure 3 shows the time evolution of the relative $L^{2}$ error

(7.14)\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 various values of $\gamma \lessgtr 0$. We see that the errors remain under unity (i.e. under 100 %) over the time interval $0 \leqslant t \leqslant 1000$, noting that the validity of the Dysthe equation deteriorates faster as $\gamma$ is decreased. This is expected in view of the stability analysis because the solution tends to become more unstable (and thus more nonlinear) with decreasing $\gamma$. Development of the BF instability is especially apparent for $\gamma = -1$ and $-2$, as indicated by a hump in their error plots.

Figure 3. Relative $L^{2}$ errors on $\eta$ between fully and weakly nonlinear solutions for $(B_0,k_0,\lambda ) = (0.002,10,1)$ and $\gamma = \{ -2, -1, 0, +1, +2 \}$. The red curve corresponds to $\gamma = 0$.

Comparison of surface elevations $\eta$ predicted from the weakly nonlinear equation (6.11) and the full nonlinear system (7.11)–(7.12) is presented in figure 4 for the same set of values of $\gamma$ at their respective times of maximum wave growth. The perturbed Stokes wave at $t = 0$ (which is the same initial condition for all cases considered) is depicted in figure 4(a). These results are consistent with our previous observations from figures 1 and 3. Excitation and growth of the most unstable sideband mode $\lambda = 1$ (according to the stability analysis) are clearly revealed in these plots. A more negative $\gamma$ promotes the BF instability (by making it happen sooner with a stronger wave amplification), while a more positive $\gamma$ tends to reduce and even offset it. In all these cases, the Dysthe model is found to provide a very good approximation up to at least $t = 1000$. As expected, for $\gamma = -2$, discrepancies are more pronounced due to the higher nonlinearity reached; a slight phase lag and drop in wave amplitude can be discerned around the main peak at $t = 500$.

Figure 4. Comparison of surface elevations $\eta$ between fully and weakly nonlinear solutions for $(B_0,k_0,\lambda ) = (0.002,10,1)$ with (a) $\gamma = 0$ ($t = 0$), (b) $\gamma = 0$ ($t = 940$), (c) $\gamma = -1$ ($t = 680$), (d) $\gamma = -2$ ($t = 500$), (e) $\gamma = +1$ ($t = 1000$), (f) $\gamma = +2$ ($t = 1000$). The blue curve represents the Hamiltonian Dysthe equation, while the red curve represents the full nonlinear system.

It is suitable to compare our Hamiltonian Dysthe equation (6.11) with another related model that has recently been derived by Curtis et al. (Reference Curtis, Carter and Kalisch2018) and Curtis & Murphy (Reference Curtis and Murphy2020) in the same physical setting. Note that these authors additionally considered surface tension, but we will only examine the gravity wave version of their model. Moreover, because they expressed their model in a form that contains a first derivative in time as well as a mixed derivative in space and time (see (2.37) in Curtis et al. Reference Curtis, Carter and Kalisch2018), we find it more appropriate to rewrite it in a more standard form with a single time derivative (as it is typically so for the Dysthe equation; Dysthe Reference Dysthe1979) to allow for a fairer comparison. We also take into account the fact that vorticity in the mathematical formulation used by Curtis et al. (Reference Curtis, Carter and Kalisch2018) is defined as the opposite of ours. The resulting model for the first-harmonic envelope $\eta _1$ is given by

(7.15)\begin{align} \partial_t \eta_1 & ={-}c_g\,\partial_x \eta_1 - {\rm i}\,\frac{c_g^{2}}{d_0}\,\partial_x^{2} \eta_1+ 2\, \frac{c_g^{3}}{d_0^{2}}\,\partial_x^{3} \eta_1 - {\rm i}\,\frac{\tilde{\alpha}_0}{d_0}\,|\eta_1|^{2} \eta_1 \nonumber\\ &\quad - \frac{\tilde{\alpha}_1}{d_0}\,|\eta_1|^{2}\,\partial_x \eta_1 - \frac{\tilde{\alpha}_2}{d_0}\,\eta_1^{2}\,\partial_x \bar{\eta}_1 + {\rm i}\,\frac{\tilde{\alpha}_3}{d_0}\,\eta_1 {{\rm H}}\,\partial_x |\eta_1|^{2} \,, \end{align}

where

(7.16)\begin{equation} \left.\begin{gathered} d_0 = 2 \varOmega_0 - \gamma , \quad c_g = \frac{1}{d_0} , \quad d_1 = \frac{d_0}{1 + \gamma c_g} , \quad \tilde{\alpha}_0 = \alpha_0 - \gamma^{2} k_0 d_0 d_1 ,\\ \tilde{\alpha}_1 = \alpha_1 - \gamma^{2} k_0 d_1 (c_g + \alpha_d \gamma d_1) - \gamma^{2} d_1 (2 d_0 + 3 k_0 c_g) + 4 c_g (\gamma^{2} k_0 d_1 - \alpha_0 c_g) ,\\ \tilde{\alpha}_2 = \alpha_2 + \gamma^{2} k_0 d_1 (c_g + \alpha_d \gamma d_1) - \gamma^{2} d_1 (d_0 + k_0 c_g) + 2 c_g (\gamma^{2} k_0 d_1 - \alpha_0 c_g) ,\\ \tilde{\alpha}_3 = \alpha_3 - \gamma k_0 c_g d_1 (d_0 + d_1) . \end{gathered}\right\} \end{equation}

The reader is directed to Curtis et al. (Reference Curtis, Carter and Kalisch2018), where the expressions for $\alpha _d$ and $\alpha _j$ ($j=0,\dots,3$) can be found. Note that ${{\rm H}}\,\partial _x = |D|$ for the non-local term in (7.15). For the purpose of comparing with the full system (7.11)–(7.12), we have also re-expressed the Curtis et al. model in a fixed reference frame, hence the additional advection term in (7.15) as compared to (2.37) in Curtis et al. (Reference Curtis, Carter and Kalisch2018). In the following discussion, we will refer to (7.15) as the ‘classical’ Dysthe equation for this problem, because it is not Hamiltonian and has the same typical form as in the irrotational case. Furthermore, its derivation is based on a perturbative Stokes-type expansion for the dependent variables $\eta$ and $\xi$, which is similar to the classical derivation by the method of multiple scales (Dysthe Reference Dysthe1979). Indeed, following Curtis et al. (Reference Curtis, Carter and Kalisch2018), the surface elevation and velocity potential at any instant $t$ can be reconstructed perturbatively from $\eta _1$ as

(7.17)\begin{equation} \left.\begin{gathered} \eta(x,t) = \gamma d_1\,|\eta_1|^{2} + 2 \, {\rm Re} \left( \eta_1\,{\rm e}^{{\rm i} \theta} + \ell_0 \eta_1^{2}\,{\rm e}^{2 {\rm i}\theta} + \cdots \right) ,\\ \xi(x,t) = {{\rm H}} \left( \frac{\varOmega_0}{k_0}\,\eta - \frac{\varOmega_0}{k_0}\,\eta {{\rm H}}\,\partial_x \eta- \frac{\gamma}{2}\,\eta^{2} + \frac{\varOmega_0}{2 k_0}\,{{\rm H}}\,\partial_x ( \eta^{2} ) + \cdots \right) , \end{gathered}\right\} \end{equation}

for which the expression of $\ell _0$ can be found in Curtis et al. (Reference Curtis, Carter and Kalisch2018), and only contributions from up to the second harmonics are included here because Curtis et al. (Reference Curtis, Carter and Kalisch2018) did not provide expressions for contributions from higher harmonics. The phase function is given by $\theta = k_0 x - \varOmega _0 t$. The ‘classical’ reconstruction procedure based on (7.17) clearly differs from the present approach. It is more explicit and thus computationally more efficient, but it is perturbative. Contributions at each order up to the desired one need to be derived, and their expressions become increasingly complicated. On the other hand, our Hamiltonian procedure requires solving an auxiliary system of PDEs to reconstruct $\eta$ and $\zeta$ (or $\xi$) from $u$, but it is non-perturbative. Indeed, (3.37) and (3.38) constitute an exact representation of the Birkhoff normal form transformation that eliminates non-resonant triads in this problem.

As an illustration, figure 5 compares the $L^{2}$ errors (7.14) on $\eta$ from the classical and Hamiltonian Dysthe equations in the large-vorticity cases $\gamma = \pm 2$. For each of these models, the error is calculated relative to the fully nonlinear solution with respective initial conditions $\eta (x,0)$ and $\xi (x,0)$. These are provided by (7.17) with

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

when (7.15) is tested against (7.11)–(7.12). Recall that $A_0$ and $B_0$ are related through (7.8). We use the same numerical methods as described earlier (and specify the same resolutions in space and time) to solve (7.15) and evaluate (7.17). For $\gamma = -2$, both Dysthe solutions are found to perform similarly, with the error from the Hamiltonian model being slightly lower than that from the classical model. The relatively quick loss of accuracy in this case, which is common to both models (with errors reaching near 50 % at $t \simeq 500$) should be attributed to deterioration of the Dysthe approximation during development of the BF instability, rather than to the reconstruction procedure. By contrast, for $\gamma = +2$, the errors remain small and do not vary much over the time interval $0 \leqslant t \leqslant 1000$, which is expected considering that the solution is more stable in this case. We see, however, that the present approach outperforms the classical one by about an order of magnitude. In all these error plots, the seemingly sharp values near $t = 0$ are already an indication of the level of approximation associated with the different equations, as they represent adjustment of the full system (7.11)–(7.12) to the imposed initial conditions during early stages of the simulation.

Figure 5. Relative $L^{2}$ errors on $\eta$ between fully and weakly nonlinear solutions for $(B_0,k_0,\lambda ) = (0.002,10,1)$. The blue curve represents the Hamiltonian Dysthe equation, while the red curve represents the classical Dysthe equation: (a) $\gamma = -2$, (b) $\gamma = +2$.

The corresponding surface profiles are depicted in figure 6 for the unstable case $\gamma = -2$, with predictions from each Dysthe model being compared to the fully nonlinear solution. Snapshots of $\eta$ at $t = 390$ (early stage of BF instability), $t = 500$ (around the time of maximum wave growth) and $t = 1000$ (near the end of the quasi-recurrent cycle of modulation–demodulation) are presented. The satisfactory performance of both Dysthe solutions as indicated in this figure is consistent with the error plots in figure 5. A noticeable discrepancy between the weakly and fully nonlinear predictions is a phase lag that tends to develop over time. Otherwise, salient features of the wave dynamics (including the shape of the steep wave at $t = 500$) seem to be well captured, even in this highly focusing regime. Regarding the comparison of surface profiles for $\gamma = +2$, these look indistinguishable from figure 4(f) at the graphical scale and thus are not displayed for convenience.

Figure 6. Comparison on $\eta$ between fully and weakly nonlinear solutions for $(B_0,k_0,\lambda ) = (0.002,10,1)$ and $\gamma = -2$ at $t = 390$, $500$, $1000$ (from left to right). (a–c) Hamiltonian Dysthe equation in blue; (d–f) classical Dysthe equation in red. The black curve represents the full nonlinear system.

We point out in passing that the main purpose of these tests is not to show whether one modulational approach is better than the other. In particular, regarding the reconstruction procedure for the classical Dysthe equation, we understand that adding contributions from higher harmonics to formulas (7.17) would likely improve their accuracy and lead to closer agreement with the full system. Rather, a goal here is to validate our new Hamiltonian approach against other existing formulations. As a byproduct of this comparison, given the overall positive assessment based on figures 5 and 6, we in turn provide an independent validation of the Curtis et al. model. Such a validation was not conducted in their earlier study (Curtis et al. Reference Curtis, Carter and Kalisch2018; Curtis & Murphy Reference Curtis and Murphy2020).

It is comforting to see that the solution of (3.37)–(3.38) helps to achieve an accurate computation of the free surface in our Hamiltonian framework, which was not obvious considering the rather lengthy expressions of (3.37)–(3.38). The good agreement found also confirms the validity of the zero-mass assumption (2.25) since it is used to evaluate non-local terms in (3.37)–(3.38). To further demonstrate the effectiveness of this reconstruction scheme (which we will refer to as full reconstruction by solving (3.37)–(3.38)), we now test the Hamiltonian Dysthe equation (6.11) by simply using (7.9)–(7.10) to recover $\eta$ and $\zeta$ from $u$ at any instant $t$ (which we will refer to as partial reconstruction). This simplified procedure is equivalent to retaining only contributions from the first harmonics in the representation of $\eta$ and $\zeta$.

The $L^{2}$ errors (7.14) associated with these two versions of our Hamiltonian approach are illustrated in figure 7 for $\gamma = \pm 2$. We have made sure again that suitable initial conditions are specified for the full system (7.11)–(7.12) when comparing it to each version. These results confirm that the decline in performance (for partial versus full reconstruction of $\eta$) can be considerable. The difference is found to be by about an order of magnitude for $\gamma = -2$, and by more than two orders of magnitude for $\gamma = +2$. In both cases, the errors grow quickly to exceed 100 % at some point during the time interval $0 \leqslant t \leqslant 1000$.

Figure 7. Relative $L^{2}$ errors on $\eta$ between fully and weakly nonlinear solutions for $(B_0,k_0,\lambda ) = (0.002,10,1)$. The blue curve represents the Hamiltonian Dysthe equation with full reconstruction, while the red curve represents the Hamiltonian Dysthe equation with partial reconstruction: (a) $\gamma = -2$, (b) $\gamma = +2$.

Examination of the surface profiles obtained from partial reconstruction as compared to the fully nonlinear solution is provided in figure 8 for $\gamma = \pm 2$. Consistent with the error plots in figure 7, we see that the discrepancies in wave amplitude and phase tend to develop faster. The phase lag is clearly noticeable and affects the entire wave train, even in the stabilizing case $\gamma = +2$. It is so severe for $\gamma = -2$ that the weakly nonlinear solution appears completely out of phase at $t = 1000$ during the near-recurrent stage.

Figure 8. Comparison on $\eta$ between fully and weakly nonlinear solutions for $(B_0,k_0,\lambda ) = (0.002,10,1)$ at $t = 390$, $500$, $1000$ (from left to right). The red curve represents the Hamiltonian Dysthe equation with partial reconstruction, while the black curve represents the full nonlinear system: (a–c) $\gamma = -2$, (d–f) $\gamma = +2$.

Finally, the time evolution of the relative error

(7.19)\begin{equation} \frac{\Delta {{\mathcal{H}}}}{{{\mathcal{H}}}_0} = \frac{|{{\mathcal{H}}} - {{\mathcal{H}}}_0|}{{{\mathcal{H}}}_0} \end{equation}

on energy (6.8) and on wave action (6.13) associated with the Hamiltonian Dysthe equation (6.11) is shown in figure 9 for various values of $\gamma$. Integrals in (6.8) and in the $L^{2}$ norm (7.14) are computed via the trapezoidal rule over the periodic cell $[0,2{\rm \pi} ]$. The reference value ${{\mathcal {H}}}_0$ denotes the initial value of (6.8) at $t = 0$. Overall, ${{\mathcal {H}}}$ is very well conserved in all these cases. The gradual loss of accuracy over time, which becomes more pronounced as $\gamma$ is decreased, is likely due to amplification of numerical errors triggered by the BF instability.

Figure 9. Relative errors on (a) $M$ and (b) ${{\mathcal {H}}}$ for the Hamiltonian Dysthe equation with $(B_0,k_0,\lambda ) = (0.002,10,1)$ and $\gamma = \{ -2, -1, 0, +1, +2 \}$. The red curve corresponds to $\gamma = 0$.