Theory of deep-water surface gravity waves derived from a Lagrangian

Abstract

An exact set of equations describing deep-water irrotational surface gravity waves, originally proposed by Balk (Phys. Fluids, vol. 8 (2), 1996, pp. 416–420), and studied in the case of standing waves by Longuet-Higgins (J. Fluid Mech., vol. 423, 2000, pp. 275–291) and Longuet-Higgins (Proc. R. Soc. Lond. A, vol. 457 (2006), 2001, pp. 495–510), are analytically examined and put in a form more suitable for practical applications. The utility of this approach is its simplicity. The Lagrangian is a low-order polynomial in the Fourier coefficients, leading to equations of motion that are correspondingly of low degree. The structure of these equations is examined, and the existence of solutions is considered. To gain intuition about the system of equations, a truncated model is first examined. Linear stability analysis is performed, and the evolution of the fully nonlinear system is discussed. The theory is then applied to fully resolved permanent progressive deep-water waves and a simple method for finding the eigenvalues and eigenvectors of the normal modes of this system is presented. Potential applications of these results are then discussed.

1 Introduction

Nonlinear surface wave models are a powerful tool for studying complex wave scenarios, and in particular for investigating the properties of steep and overturning waves. An improved description of the dynamics of surface water waves, and their behaviour under these nonlinear conditions, is crucial for a better understanding of air–sea interaction processes (Melville 1996). Nonlinear surface waves pose a formidable theoretical problem, and analytical studies of weakly nonlinear waves have led to a better understanding of water wave phenomena, but these results are limited by their restriction to weak nonlinearities. That is, many of these theoretical predictions are far outside of their region of validity as wave breaking is approached and during the subsequent breaking event. For this reason, numerical models of surface waves can provide insight into the physics of highly nonlinear waves processes. In this paper we derive, and put in a form more suitable for numerical applications, the model of Balk (1996), based on a Lagrangian for deep-water surface gravity waves.

The pioneering study of Longuet-Higgins & Cokelet (1976) numerically integrated the equations of motion for irrotational inviscid deep-water surface gravity waves, allowing for the systematic examination of properties of waves up to and past the point of overturning. The method employed by Longuet-Higgins & Cokelet (1976, see also Dold 1992) used a Green’s identity to solve Laplace’s equation, which was then used to time step the velocity potential and the free surface displacement. This scheme is a mixed Eulerian–Lagrangian approach, with particles on the free surface serving as the dependent variables of the system. Additionally, this method is advantageous because the points tend to cluster around the regions of highest curvature, so that one can resolve wave overturning with relatively few surface points.

The method of Balk (1996) takes a different, variational, approach to describing two dimensional surface gravity waves. At each point in time, the domain is conformally mapped to a periodic lower half-plane in the complex domain. The geometry of this domain makes solving Laplace’s equation trivial and the governing equations can then be written as a set of coupled low-order polynomial ordinary differential equations in the dependent variables, which are Fourier coefficients of the free surface displacement.

Variational models have proven to be particularly useful in fluid mechanics (Salmon 1988). This is due to several reasons, including their succinctness, the ability to find conservation laws (in particular for truncated models where conservation laws are not always straightforward to derive) and the freedom to choose coordinate systems to optimize simplicity. Furthermore, a variational approach allows for a priori knowledge of the necessarily truncated integrals of motion used as basic tests on the fidelity of the numerical model. These benefits motivated a closer inspection of Balk’s model.

Although there are currently a number of different analytical and numerical approaches to studying surface waves, Balk’s approach is examined here because of its essential simplicity. For example, a common technique is to formulate the water wave equations of motion in terms of a Hamiltonian (Zakharov 1968). However, in this method, a Taylor series around the mean water level is performed, and subsequently an infinite series is needed to describe the kinetic energy. Conformal mapping techniques help alleviate this problem, and this has been examined in detail by, for example, Fornberg (1980), Tanveer (1991) and Fokas & Nachbin (2012). However, in all of these approaches, sophisticated tools, or considerable algebraic complexity, is necessarily introduced. For practical (i.e. numerical) applications then, a straightforward method is desirable. Balk’s relatively unexploited work provides this framework. Longuet-Higgins (2000, 2001) recognized the utility of such an approach, and applied this framework to the study of standing waves. This motivated us to extend the work of Longuet-Higgins, and derive the equations of motion for the more general system.

In this paper, the scheme of Balk (1996, hereinafter referred to as Balk) is derived from Hamilton’s principle, applied to the action describing irrotational inviscid water waves. This is advantageous as one knows the form of the conserved quantities (e.g. energy, mass and momentum) when numerical approximations are necessarily applied. Furthermore, the governing equations are found to be a set of coupled lower-order polynomial ordinary differential equations, making the numerical implementation straightforward (Longuet-Higgins 2000). This also the facilities the discussion of various orders of approximation. Significantly, normal modes of steep Stokes waves are derived in this framework, facilitating analytic and numerical examination of the geometry and kinematics of steep and focusing waves (Deike, Pizzo & Melville 2017; Pizzo 2017). The scheme presented here will be used in a subsequent paper to analyse the nonlinear stability of steep Stokes waves (Longuet-Higgins & Dommermuth 1997).

The outline of this paper is as follows. First, in § 2 we provide a derivation of the model of Balk. In § 3 the equations are examined. Next, in § 4 we look at a truncated model, with only one spatial mode. In § 5, the linear stability of Stokes waves is examined. The results are then discussed in § 6.

2 Derivation of the equations of Balk

In this section we derive the governing equations for irrotational inviscid surface gravity waves in a conformally mapped reference frame. This formulation has a number of distinct advantages, including the ability to model wave overturning, a knowledge of the truncated properties (a necessary feature of any numerical model) of the waves (e.g. energy, momentum), and a relatively straightforward recipe for numerical implementation. Although the model is algebraically intensive, the tools needed to derive these results are relatively simple (cf. the boundary integral methods of Longuet-Higgins & Cokelet (1976) and Dold & Peregrine (1986)). Instead of following Balk directly, we make use of some results from Dyachenko et al. (1996) and show how this approach leads to a more natural derivation of the governing equations originally found by Balk.

To this end, we consider irrotational two-dimensional inviscid flow in a fluid of infinite depth with a free surface. Density and the gravity constant are set to 1 and the flow is assumed to be $2\unicode[STIX]{x03C0}$ periodic in the $x$ -coordinate.

The domain in the complex plane of the fluid is given by $z=x+\text{i}y$ , which can be considered to be the conformal image of the lower half-plane described by $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}+\text{i}\unicode[STIX]{x1D705}$ (these variables are time dependent, but this is not explicitly written now for clarity of presentation).

The transformation from the $(\unicode[STIX]{x1D709},\unicode[STIX]{x1D705})$ plane to the $(x,y)$ plane is given by $z(\unicode[STIX]{x1D701})=x(\unicode[STIX]{x1D709},\unicode[STIX]{x1D705})+\text{i}y(\unicode[STIX]{x1D709},\unicode[STIX]{x1D705})$ . The velocity of fluid $w=u+\text{i}v$ is given by the complex potential $f(z)=\unicode[STIX]{x1D711}(x,y)+\text{i}\unicode[STIX]{x1D6F9}(x,y)$ (here $\unicode[STIX]{x1D711}$ and $\unicode[STIX]{x1D6F9}$ represent the velocity potential and streamfunction, respectively, and we assume they satisfy the condition $v\rightarrow 0$ as $y\rightarrow -\infty$ ). That is,

(2.1) $$\begin{eqnarray}w^{\ast }=\frac{\text{d}f}{\text{d}z}=\left(\frac{\text{d}{\mathcal{C}}}{\text{d}\unicode[STIX]{x1D701}}\right)\left(\frac{\text{d}z}{\text{d}\unicode[STIX]{x1D701}}\right)^{-1},\end{eqnarray}$$

where ${\mathcal{C}}(\unicode[STIX]{x1D701})=a(\unicode[STIX]{x1D709},\unicode[STIX]{x1D705})+\text{i}b(\unicode[STIX]{x1D709},\unicode[STIX]{x1D705})$ is the complex potential in the $\unicode[STIX]{x1D701}$ plane and $^{\ast }$ denotes the complex conjugate. The surface of the fluid is going to be the image of the real axis in the $\unicode[STIX]{x1D701}$ -plane, i.e.

(2.2a,b ) $$\begin{eqnarray}X(\unicode[STIX]{x1D709})\equiv x(\unicode[STIX]{x1D709},0),\quad Y(\unicode[STIX]{x1D709})\equiv y(\unicode[STIX]{x1D709},0),\end{eqnarray}$$

and the surface values of the velocity potential and streamfunction are

(2.3a,b ) $$\begin{eqnarray}A(\unicode[STIX]{x1D709})\equiv a(\unicode[STIX]{x1D709},0),\quad B(\unicode[STIX]{x1D709})\equiv b(\unicode[STIX]{x1D709},0).\end{eqnarray}$$

We note that $\{X,Y,A,B\}$ completely determine the state of the fluid. Furthermore, Balk notes that knowledge of $\{Y(\unicode[STIX]{x1D709}),B(\unicode[STIX]{x1D709})\}$ are sufficient to completely characterize the fluid, since the boundary values at $\unicode[STIX]{x1D705}=0$ of the real parts of analytic functions in the lower half-plane can be found from the imaginary parts via the Hilbert transform.

Recall, in the $x$ – $y$ plane, we have from Zakharov (1968, see also Miles 1977) that the action for water waves can be written as

(2.4) $$\begin{eqnarray}S=\int \unicode[STIX]{x1D713}\unicode[STIX]{x1D702}_{t}-\left(\int _{-\infty }^{\unicode[STIX]{x1D702}}\frac{1}{2}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719})^{2}\,\text{d}y+\frac{1}{2}\unicode[STIX]{x1D702}^{2}\right)\,\text{d}x\,\text{d}t,\end{eqnarray}$$

where $\unicode[STIX]{x1D713}=\unicode[STIX]{x1D719}(x,\unicode[STIX]{x1D702},t)$ .

Now, from the chain rule the time derivative of the free surface $\unicode[STIX]{x1D702}$ can be rewritten in terms of the conformal variables as (cf. Dyachenko et al. 1996, equation (2.14))

(2.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}_{t}\,\text{d}x & = & \displaystyle \left(\frac{\text{d}Y}{\text{d}t}\frac{\text{d}X}{\text{d}\unicode[STIX]{x1D709}}-\frac{\text{d}X}{\text{d}t}\frac{\text{d}Y}{\text{d}\unicode[STIX]{x1D709}}\right)\text{d}\unicode[STIX]{x1D709}\nonumber\\ \displaystyle & = & \displaystyle {\dot{Y}}X^{\prime }\,\text{d}\unicode[STIX]{x1D709}-{\dot{X}}Y^{\prime }\,\text{d}\unicode[STIX]{x1D709},\end{eqnarray}$$

where dots and primes over the summation symbol represent differentiation with respect to time and $\unicode[STIX]{x1D709}$ , respectively. Therefore, this term in the action becomes

(2.6) $$\begin{eqnarray}\int \unicode[STIX]{x1D713}\unicode[STIX]{x1D702}_{t}\,\text{d}x=\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}A({\dot{Y}}X^{\prime }-{\dot{X}}Y^{\prime })\,\text{d}\unicode[STIX]{x1D709}.\end{eqnarray}$$

Next, the kinetic and potential energy densities, per unit wavelength, can be written as

(2.7) $$\begin{eqnarray}V=\frac{1}{2}\int \unicode[STIX]{x1D702}^{2}\,\text{d}x=\frac{1}{2}\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}Y^{2}X^{\prime }\,\text{d}\unicode[STIX]{x1D709},\end{eqnarray}$$

while the kinetic energy density is given by (via a Green integral identity)

(2.8) $$\begin{eqnarray}T=\frac{1}{2}\int ((\unicode[STIX]{x1D6FC}_{x})^{2}+(\unicode[STIX]{x1D6FC}_{y})^{2})\,\text{d}y\,\text{d}x=-\frac{1}{2}\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}AB^{\prime }\,\text{d}\unicode[STIX]{x1D709}.\end{eqnarray}$$

Finally, exploiting the analyticity of ${\mathcal{C}}$ , we can use the Cauchy–Riemann equations to rewrite (2.4) as

(2.9) $$\begin{eqnarray}S=\int _{t_{a}}^{t_{b}}\left(\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}A({\dot{Y}}X^{\prime }-{\dot{X}}Y^{\prime })-\frac{1}{2}AB^{\prime }-\frac{1}{2}Y^{2}X^{\prime }\,\text{d}\unicode[STIX]{x1D709}\right)\,\text{d}t,\end{eqnarray}$$

with the time interval $(t_{a},t_{b})$ chosen so that the dependent variables are fixed at these times.

Now, variations of $S$ with respect to $A$ give

(2.10) $$\begin{eqnarray}B^{\prime }=({\dot{X}}Y^{\prime }-{\dot{Y}}X^{\prime }),\end{eqnarray}$$

which is precisely the kinematic boundary condition, in the conformally mapped coordinates (equation (2) in the paper of Balk). The constraint that mass is conserved is the condition

(2.11) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\int \unicode[STIX]{x1D702}\,\text{d}x=\frac{\text{d}}{\text{d}t}\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}YX^{\prime }\,\text{d}\unicode[STIX]{x1D709}=0,\end{eqnarray}$$

so that taking the constant of integration to be 0, corresponding to the mean water level being at $y=0$ , we find

(2.12) $$\begin{eqnarray}\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}YX^{\prime }\,\text{d}\unicode[STIX]{x1D709}=0,\end{eqnarray}$$

which is equation (3) in Balk.

Substituting (2.10) into (2.9) we find that the action becomes

(2.13) $$\begin{eqnarray}S=\int _{t_{a}}^{t_{b}}T-V\,\text{d}t,\end{eqnarray}$$

with the constraint given in (2.12). In the derivation given above, the constraints on the system enter naturally by a change of variables in the action, whereas in Balk’s derivation, he simply imposes these conditions.

Note, a priori there is no reason that the Lagrangian should be the difference of kinetic and potential energies (e.g. see (2.4)) in an Eulerian framework, which has led to significant confusion in the literature (see, for example, the discussion in Luke (1967) and Salmon (1988)). However, both Balk’s original formulation and equation (8) in Luke (1967) constrain the system to obey conservation of mass, the kinematic boundary condition, and the condition of no flow at the bottom. When these conditions are met, then the Lagrangian for water waves may be set equal to the difference of kinetic and potential energies (see the discussion in Luke (1967)).

Conservation of mass can explicitly be built into the Lagrangian by exploiting the periodicity of the flow in the $x$ -direction. This property implies the following Fourier expansions:

(2.14) $$\begin{eqnarray}X(\unicode[STIX]{x1D709},t)+\text{i}Y(\unicode[STIX]{x1D709},t)=\unicode[STIX]{x1D709}+\mathop{\sum }_{k=-\infty }^{\infty }(X_{k}(t)+\text{i}Y_{k}(t))\text{e}^{-\text{i}k\unicode[STIX]{x1D709}},\end{eqnarray}$$

and

(2.15) $$\begin{eqnarray}A(\unicode[STIX]{x1D709},t)+\text{i}B(\unicode[STIX]{x1D709},t)=\mathop{\sum }_{k=-\infty }^{\infty }(A_{k}(t)+\text{i}B_{k}(t))\text{e}^{-\text{i}k\unicode[STIX]{x1D709}},\end{eqnarray}$$

where we also must have $Y_{-k}=Y_{k}^{\ast },X_{-k}=X_{k}^{\ast },A_{-k}=A_{k}^{\ast },B_{-k}=B_{k}^{\ast }$ for $\{X,Y,A,B\}$ to be real. The Hilbert transform then implies ( $k\neq 0$ )

(2.16a,b ) $$\begin{eqnarray}X_{k}=\text{i}\unicode[STIX]{x1D70E}_{k}Y_{k},\quad A_{k}=\text{i}\unicode[STIX]{x1D70E}_{k}B_{k},\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{k}=(1,0,-1)$ for $(k>0,k=0,k<0)$ , respectively. This gives us relationships between the real and imaginary parts of the conjugate pairs. We note explicitly

(2.17) $$\begin{eqnarray}\displaystyle & \displaystyle X=\unicode[STIX]{x1D709}+\mathop{\sum }_{1}^{\infty }2\text{Re}[Y_{k}]\sin k\unicode[STIX]{x1D709}-2\text{Im}[Y_{k}]\cos k\unicode[STIX]{x1D709}, & \displaystyle\end{eqnarray}$$

(2.18) $$\begin{eqnarray}\displaystyle & \displaystyle Y=Y_{0}+\mathop{\sum }_{1}^{\infty }2\text{Re}[Y_{k}]\cos k\unicode[STIX]{x1D709}+2\text{Im}[Y_{k}]\sin k\unicode[STIX]{x1D709}, & \displaystyle\end{eqnarray}$$

(2.19) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D709}\in (0,2\unicode[STIX]{x03C0}), & \displaystyle\end{eqnarray}$$

where $\text{Re}/\text{Im}$ stand for the real and imaginary parts, respectively.

Using the kinematic boundary condition, i.e. equation (2.10), we find (for $n\neq 0$ )

(2.20) $$\begin{eqnarray}B_{n}=\frac{\text{i}}{n}\left(-{\dot{Y}}_{n}+\mathop{\sum }_{j+k=n}{\dot{Y}}_{j}Y_{k}k(\unicode[STIX]{x1D70E}_{j}-\unicode[STIX]{x1D70E}_{k})\right).\end{eqnarray}$$

The complex potential is defined up to an arbitrary constant, so that without loss of generality, we may set $A_{0}=B_{0}=0$ .

Substituting the Fourier expansions of $\{A,B\}$ into the equation for the kinetic energy, i.e. equation (2.8), we find

(2.21) $$\begin{eqnarray}T=\frac{1}{2}\mathop{\sum }_{k=-\infty }^{\infty }|k|B_{k}B_{-k}.\end{eqnarray}$$

Similarly, with these expansions the potential energy is found directly from (2.7) and (2.14) (see also Longuet-Higgins 2000)

(2.22) $$\begin{eqnarray}V=\frac{1}{2}\mathop{\sum }_{k=-\infty }^{\infty }Y_{k}Y_{-k}+\frac{1}{2}\mathop{\sum }_{k+j+l=0}|l|Y_{k}Y_{j}Y_{l}.\end{eqnarray}$$

The Lagrangian, $L=T-V$ , can now be written as

(2.23) $$\begin{eqnarray}L=\frac{1}{2}\mathop{\sum }_{k=-\infty }^{\infty }|k|B_{k}B_{-k}-\frac{1}{2}\mathop{\sum }_{k=-\infty }^{\infty }Y_{k}Y_{-k}-\frac{1}{2}\mathop{\sum }_{k+j+l=0}|l|Y_{k}Y_{j}Y_{l},\end{eqnarray}$$

where $B_{k}$ is given in (2.20) and we choose $Y_{0}$ to satisfy the constraint given in (2.12), namely,

(2.24) $$\begin{eqnarray}Y_{0}=-\mathop{\sum }_{k=-\infty }^{\infty }|k|Y_{k}Y_{-k}.\end{eqnarray}$$

The two most obvious conserved quantities associated with this Lagrangian are the $x$ -momentum (associated with phase shift invariance), defined as

(2.25) $$\begin{eqnarray}{\mathcal{P}}=\int \unicode[STIX]{x1D6FD}_{y}\,\text{d}x\,\text{d}y=\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}BX^{\prime }\,\text{d}\unicode[STIX]{x1D709}=\mathop{\sum }_{k=-\infty }^{\infty }|k|B_{k}Y_{-k},\end{eqnarray}$$

and energy (associated with time shift invariance) given by

(2.26) $$\begin{eqnarray}H\equiv T+V=\frac{1}{2}\mathop{\sum }_{k=-\infty }^{\infty }|k|B_{k}B_{-k}+\frac{1}{2}\mathop{\sum }_{k=-\infty }^{\infty }Y_{k}Y_{-k}+\frac{1}{2}\mathop{\sum }_{k+j+l=0}|l|Y_{k}Y_{j}Y_{l}.\end{eqnarray}$$

Next, the governing equations are found by varying the Lagrangian (2.23) with respect to $Y_{k}$ . This yields the Euler–Lagrange equations for each mode $k$

(2.27) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\left(\frac{\unicode[STIX]{x2202}L}{\unicode[STIX]{x2202}{\dot{Y}}_{k}}\right)=\frac{\unicode[STIX]{x2202}L}{\unicode[STIX]{x2202}Y_{k}}.\end{eqnarray}$$

For partial wave solutions, the infinite sums used in the Fourier expansion of the dependent variables will be truncated at a finite $N$ , so that we will have a truncated Lagrangian $L_{N}=T_{N}-V_{N}$ . This is advantageous since we know the conserved quantities associated with this truncated Lagrangian before deriving the dynamical evolution equations (see Salmon (1988) for a thorough discussion).

3 Equations of motion

In order to better understand the governing equations, we extend the work of Longuet-Higgins (2000, 2001) from the special case of standing waves to the more general equations of motion.

First, from (2.21) we have for $n>0$ ,

(3.1) $$\begin{eqnarray}\displaystyle \text{i}nB_{n} & = & \displaystyle \cdots +2|n+2|Y_{n+2}{\dot{Y}}_{-2}+2|n+1|Y_{n+1}{\dot{Y}}_{-1}\nonumber\\ \displaystyle & & \displaystyle +\,|n|Y_{n}{\dot{Y}}_{0}+{\dot{Y}}_{n}+2|1|Y_{-1}{\dot{Y}}_{n+1}+2|2|Y_{-2}{\dot{Y}}_{n+2}+\cdots \,,\end{eqnarray}$$

and similarly for $n<0$ ,

(3.2) $$\begin{eqnarray}\text{i}nB_{n}=\cdots +2|2|Y_{2}{\dot{Y}}_{n-2}+2Y_{1}{\dot{Y}}_{n-1}+{\dot{Y}}_{n}+|n|Y_{n}{\dot{Y}}_{0}+2|n-1|Y_{n-1}{\dot{Y}}_{1}+2|n-2|Y_{n-2}{\dot{Y}}_{2}+\cdots \,.\end{eqnarray}$$

By substituting these two expressions into (2.21) we find that we can write the kinetic energy as

(3.3) $$\begin{eqnarray}T=\mathop{\sum }_{m=1}^{N}\left(\mathop{\sum }_{k=-N}^{N}P_{mk}{\dot{Y}}_{k}\right)\left(\mathop{\sum }_{l=-N}^{N}P_{-ml}{\dot{Y}}_{l}\right),\end{eqnarray}$$

where

(3.4) $$\begin{eqnarray}P_{mn}=\frac{2}{\sqrt{[m]}}(f(m-n)|m-n|Y_{m-n}-|mn|Y_{m}Y_{-n}).\end{eqnarray}$$

The operator $[\ast ]$ is defined as

$$\begin{eqnarray}\displaystyle [j]=\left\{\begin{array}{@{}ll@{}}|\,j|: & j\neq 0,\\ 1: & j=0,\end{array}\right. & & \displaystyle \nonumber\end{eqnarray}$$

and

$$\begin{eqnarray}\displaystyle f(m-k)=\left\{\begin{array}{@{}ll@{}}0: & \text{Sign}(k)=\text{Sign}(m-k),\\ 1: & \text{otherwise},\end{array}\right. & & \displaystyle \nonumber\end{eqnarray}$$

and $Y_{k}\equiv 0$ when $k>N$ for $N$ the number of Fourier modes in our model, i.e. the resolution. Note, by construction, the first term in parentheses in (3.3) is the complex conjugate of the second term.

By using these definitions, and reversing the order of summation, the kinetic energy can be rewritten as

(3.5) $$\begin{eqnarray}T=\mathop{\sum }_{k=-N}^{N}\mathop{\sum }_{l=-N}^{N}Q_{kl}{\dot{Y}}_{k}{\dot{Y}}_{l},\end{eqnarray}$$

where

(3.6) $$\begin{eqnarray}Q_{kl}=\mathop{\sum }_{m=1}^{N}P_{mk}P_{-ml}.\end{eqnarray}$$

That is, the matrix $Q_{kl}$ is related to the column by column multiplication of two $P$ matrices.

We now use these definitions to solve for the equations of motion. Substituting the kinetic energy into the first term of equation (2.27) we find

(3.7) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}\left(\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}{\dot{Y}}_{n}}\right) & = & \displaystyle \frac{\text{d}}{\text{d}t}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}{\dot{Y}}_{n}}\right)\mathop{\sum }_{k}\mathop{\sum }_{l}Q_{kl}{\dot{Y}}_{k}{\dot{Y}}_{l}=\frac{\text{d}}{\text{d}t}\left(\mathop{\sum }_{k}Q_{kn}{\dot{Y}}_{k}+\mathop{\sum }_{l}Q_{nl}{\dot{Y}}_{l}\right)\nonumber\\ \displaystyle & = & \displaystyle \mathop{\sum }_{k}Q_{nk}\ddot{Y} _{k}+\mathop{\sum }_{l}Q_{nl}\ddot{Y} _{l}+\mathop{\sum }_{k}\frac{\text{d}Q_{kn}}{\text{d}t}{\dot{Y}}_{k}+\mathop{\sum }_{l}\frac{\text{d}Q_{nl}}{\text{d}t}{\dot{Y}}_{l}.\end{eqnarray}$$

Next, we note that

(3.8) $$\begin{eqnarray}\frac{\text{d}Q_{kn}}{\text{d}t}=\mathop{\sum }_{j}\frac{\unicode[STIX]{x2202}Q_{kn}}{\unicode[STIX]{x2202}Y_{j}}{\dot{Y}}_{j},\end{eqnarray}$$

which implies that (3.7) becomes

(3.9) $$\begin{eqnarray}\mathop{\sum }_{k}Q_{nk}\ddot{Y} _{k}+\mathop{\sum }_{l}Q_{nl}\ddot{Y} _{l}+\mathop{\sum }_{k}\mathop{\sum }_{l}\frac{\unicode[STIX]{x2202}Q_{kn}}{\unicode[STIX]{x2202}Y_{l}}{\dot{Y}}_{l}{\dot{Y}}_{k}+\frac{\unicode[STIX]{x2202}Q_{nl}}{\unicode[STIX]{x2202}Y_{k}}{\dot{Y}}_{l}{\dot{Y}}_{k}.\end{eqnarray}$$

Finally, the second term in the Euler–Lagrange equation, i.e. equation (2.27), is given by

(3.10) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}Y_{n}}=\mathop{\sum }_{k}\mathop{\sum }_{l}\frac{\unicode[STIX]{x2202}Q_{kl}}{\unicode[STIX]{x2202}Y_{n}}{\dot{Y}}_{k}{\dot{Y}}_{l}.\end{eqnarray}$$

Putting this together, we see that our governing equations are now of the form

(3.11) $$\begin{eqnarray}\mathop{\sum }_{k}Q_{kn}\ddot{Y} _{k}+\mathop{\sum }_{l}Q_{nl}\ddot{Y} _{l}+\mathop{\sum }_{k}\mathop{\sum }_{l}\left(\frac{\unicode[STIX]{x2202}Q_{kn}}{\unicode[STIX]{x2202}Y_{l}}{\dot{Y}}_{l}{\dot{Y}}_{k}+\frac{\unicode[STIX]{x2202}Q_{nl}}{\unicode[STIX]{x2202}Y_{k}}{\dot{Y}}_{l}{\dot{Y}}_{k}-\frac{\unicode[STIX]{x2202}Q_{kl}}{\unicode[STIX]{x2202}Y_{n}}{\dot{Y}}_{k}{\dot{Y}}_{l}\right)-\frac{\unicode[STIX]{x2202}V}{\unicode[STIX]{x2202}Y_{n}}=0,\end{eqnarray}$$

for each $n=(\pm 1,\pm 2,\ldots \pm N)$ .

Further simplifications of these equations can be made by realizing that

(3.12) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}Q_{kl}}{\unicode[STIX]{x2202}Y_{n}}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}Y_{n}}\mathop{\sum }_{m}P_{mk}P_{-ml}=\mathop{\sum }_{m}\frac{\unicode[STIX]{x2202}P_{mk}}{\unicode[STIX]{x2202}Y_{n}}P_{-ml}+P_{mk}\frac{\unicode[STIX]{x2202}P_{-ml}}{\unicode[STIX]{x2202}Y_{n}}.\end{eqnarray}$$

From the definition of $P_{ij}$ , i.e. equation (3.4), we find that

$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}P_{mk}}{\unicode[STIX]{x2202}Y_{n}}=\frac{2}{\sqrt{[m]}}\left\{\begin{array}{@{}l@{}}f(m-k)|m-k|:\quad m-k-n=0,\\ -|mk|Y_{m}:\quad k+n=0,\\ -|mk|Y_{-k}:\quad m-n=0.\end{array}\right. & & \displaystyle \nonumber\end{eqnarray}$$

Next, we consider the potential energy term, which, from equation (2.7), can be written as

(3.13) $$\begin{eqnarray}\displaystyle 2V & = & \displaystyle -4\left(\mathop{\sum }_{n=1}^{\infty }|n|Y_{n}Y_{-n}\right)^{2}+2\mathop{\sum }_{n=1}^{\infty }Y_{n}Y_{-n}\nonumber\\ \displaystyle & & \displaystyle +\,8\mathop{\sum }_{i=-\infty }^{\infty }\mathop{\sum }_{j=-\infty }^{\infty }\mathop{\sum }_{k=-\infty }^{\infty }F(i,j,k)\frac{Y_{i}Y_{j}Y_{k}}{8[j][k]},\end{eqnarray}$$

where

$$\begin{eqnarray}\displaystyle F(i,j,k)=\left\{\begin{array}{@{}ll@{}}|ijk|: & i+j+k=0\;\text{and}\;i,j,k\neq 0,\\ 0: & \text{otherwise}.\end{array}\right. & & \displaystyle \nonumber\end{eqnarray}$$

Finally, we define

(3.14) $$\begin{eqnarray}S_{nkl}=\left(\frac{\unicode[STIX]{x2202}Q_{lk}}{\unicode[STIX]{x2202}Y_{n}}-\frac{\unicode[STIX]{x2202}Q_{nl}}{\unicode[STIX]{x2202}Y_{k}}-\frac{\unicode[STIX]{x2202}Q_{ln}}{\unicode[STIX]{x2202}Y_{k}}\right),\end{eqnarray}$$

so that the governing equation becomes

(3.15) $$\begin{eqnarray}\mathop{\sum }_{l}(Q_{nl}+Q_{ln})\ddot{Y} _{l}=\mathop{\sum }_{l}\mathop{\sum }_{k}S_{nkl}{\dot{Y}}_{k}{\dot{Y}}_{l}+\frac{\unicode[STIX]{x2202}V}{\unicode[STIX]{x2202}Y_{n}}\quad \quad (n=\pm 1,\pm 2,\ldots ).\end{eqnarray}$$

Note, this can be written more concisely as

(3.16) $$\begin{eqnarray}\mathbb{Q}_{nl}\ddot{Y} _{l}-S_{nkl}{\dot{Y}}_{k}{\dot{Y}}_{l}-\frac{\unicode[STIX]{x2202}V}{\unicode[STIX]{x2202}Y_{n}}=0,\quad (n=\pm 1,\pm 2,\ldots )=0,\end{eqnarray}$$

where

(3.17) $$\begin{eqnarray}\mathbb{Q}_{nl}=Q_{nl}+Q_{ln},\end{eqnarray}$$

and Einstein summation is assumed. Together with initial conditions at time $t=0$ for $Y_{n}$ and ${\dot{Y}}_{n}$ , equation (3.16) are the complete equations for surface gravity waves written in a compact form.

Analytically, we may also check that our equations reduce to those of Longuet-Higgins (2000), for the case where all $Y_{k}$ are real, which corresponds to standing waves. It is a simple matter to show that under the condition that $Y_{-k}=Y_{k}$ , and by mapping into the coordinates used in that work, namely $a_{n}=2|n|Y_{n}$ with $a_{0}=1$ , then the entries of $P$ become

(3.18) $$\begin{eqnarray}P_{mn}\rightarrow \frac{1}{2}\frac{1}{\sqrt{|m|}|n|}(a_{n-m}+a_{m+n}-a_{m}a_{n}).\end{eqnarray}$$

With this identification, the rest of our variables can easily be shown to be equivalent to those of Longuet-Higgins (2000).

4 Example: $N=1$

In this section we consider a simple example where only one mode is present. We perform linear stability analysis of permanent progressive waves, as this system illustrates many properties of the large $N$ Stokes waves to be considered in the next section. Furthermore, we examine properties of general solutions to this system of equations.

4.1 Governing equations

To begin, we let $\unicode[STIX]{x1D6FC}=Y_{1}$ and $\unicode[STIX]{x1D6FD}=Y_{-1}$ , to find

(4.1) $$\begin{eqnarray}T=[(1-2\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD})\dot{\unicode[STIX]{x1D6FC}}-2\unicode[STIX]{x1D6FC}^{2}\dot{\unicode[STIX]{x1D6FD}}][(1-2\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD})\dot{\unicode[STIX]{x1D6FD}}-2\unicode[STIX]{x1D6FD}^{2}\dot{\unicode[STIX]{x1D6FC}}],\end{eqnarray}$$

and

(4.2) $$\begin{eqnarray}V=\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}-2(\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD})^{2}.\end{eqnarray}$$

The Euler–Lagrange equations for $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ then imply

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle L\ddot{\unicode[STIX]{x1D6FC}}+M\ddot{\unicode[STIX]{x1D6FD}}=-2[L_{\unicode[STIX]{x1D6FC}}\dot{\unicode[STIX]{x1D6FC}}^{2}+2L_{\unicode[STIX]{x1D6FD}}\dot{\unicode[STIX]{x1D6FC}}\dot{\unicode[STIX]{x1D6FD}}+(2M_{\unicode[STIX]{x1D6FD}}-N_{\unicode[STIX]{x1D6FC}})\dot{\unicode[STIX]{x1D6FD}}^{2}]-C_{1}, & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle M\ddot{\unicode[STIX]{x1D6FC}}+N\ddot{\unicode[STIX]{x1D6FD}}=-2[(2M_{\unicode[STIX]{x1D6FC}}-L_{\unicode[STIX]{x1D6FD}})\dot{\unicode[STIX]{x1D6FC}}^{2}+2N_{\unicode[STIX]{x1D6FC}}\dot{\unicode[STIX]{x1D6FC}}\dot{\unicode[STIX]{x1D6FD}}+N_{\unicode[STIX]{x1D6FD}}\dot{\unicode[STIX]{x1D6FD}}^{2}]-C_{2}, & \displaystyle\end{eqnarray}$$

where

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle L=4\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}^{3}-2\unicode[STIX]{x1D6FD}^{2}, & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle M={\textstyle \frac{1}{2}}+4\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D6FD}^{2}-2\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}, & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle N=4\unicode[STIX]{x1D6FC}^{3}\unicode[STIX]{x1D6FD}-2\unicode[STIX]{x1D6FC}^{2}, & \displaystyle\end{eqnarray}$$

and

(4.8a,b ) $$\begin{eqnarray}C_{1}=4\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}^{2}-\unicode[STIX]{x1D6FD},\quad C_{2}=4\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D6FD}-\unicode[STIX]{x1D6FC}.\end{eqnarray}$$

Note, these governing equations can be rewritten as

(4.9) $$\begin{eqnarray}\left[\begin{array}{@{}cc@{}}L & M\\ M & N\end{array}\right]\left[\begin{array}{@{}c@{}}\ddot{\unicode[STIX]{x1D6FC}}\\ \ddot{\unicode[STIX]{x1D6FD}}\end{array}\right]=\left[\begin{array}{@{}c@{}}F_{1}\\ F_{2}\end{array}\right],\end{eqnarray}$$

where $F_{1},F_{2}$ are the right-hand sides of equations (4.3) and (4.4), respectively. We define the coefficient matrix $\mathbb{Q}$ as

(4.10) $$\begin{eqnarray}\mathbb{Q}=\left[\begin{array}{@{}cc@{}}L & M\\ M & N\end{array}\right].\end{eqnarray}$$

Now, equation (4.9) has solutions provided that $\mathbb{Q}$ is invertible. Therefore, there is some interest in looking at how the wave evolution breaks down as $\text{det}(\mathbb{Q})\rightarrow 0$ . To this end, the determinant of $\mathbb{Q}$ is

(4.11) $$\begin{eqnarray}\text{det}(\mathbb{Q})=LN-M^{2}=0\;\Longrightarrow \;\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}=1/4.\end{eqnarray}$$

The reality condition on the Fourier components mandates that $\unicode[STIX]{x1D6FC}^{\ast }=\unicode[STIX]{x1D6FD}$ , so that letting $\unicode[STIX]{x1D6FC}={\mathcal{A}}+i{\mathcal{B}}=\unicode[STIX]{x1D6FD}^{\ast }$ , this condition implies that $({\mathcal{A}},{\mathcal{B}})$ must lie within a circle of radius $1/2$ , i.e.

(4.12) $$\begin{eqnarray}{\mathcal{A}}^{2}+{\mathcal{B}}^{2}<1/4.\end{eqnarray}$$

In the case ${\mathcal{B}}=0$ , we have the result of the $N=1$ case for standing waves, as discussed by Longuet-Higgins (2000).

4.2 Permanent progressive waves

A simplified class of solutions to equation (4.9) are those of the permanent progressive type. That is, we choose

(4.13a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}_{0}\text{e}^{\text{i}ct},\quad \unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FC}_{0}\text{e}^{-\text{i}ct},\end{eqnarray}$$

where $c$ is a to be determined phase velocity of the waves. When this ansatz is made, we find from (4.9) the following relationship between $c$ and $\unicode[STIX]{x1D6FC}_{0}$ :

(4.14) $$\begin{eqnarray}c^{2}=1-4\unicode[STIX]{x1D6FC}_{0}^{2}.\end{eqnarray}$$

Figure 1. Permanent progressive solutions to the $N=1$ equations of motion. For larger values of $\unicode[STIX]{x1D6FC}_{0}$ , the waves become peakier at the trough and flatter at the crests. For $\unicode[STIX]{x1D6FC}_{0}\rightarrow 1/2$ , the solution forms a cusp so that there is a limiting wave in this system.

For $\unicode[STIX]{x1D6FC}_{0}$ small we approach the infinitesimal sinusoidal wave solution we expect from linear theory for deep-water surface gravity waves (Phillips 1977). As $\unicode[STIX]{x1D6FC}_{0}\rightarrow 1/2$ , the free surface forms a cusp (note this is also when the coefficient matrix $\mathbb{Q}$ becomes singular), taking the form of a cycloid (see figure 1). This may also be seen as the point where $X_{\unicode[STIX]{x1D709}}(\unicode[STIX]{x1D709}=\unicode[STIX]{x03C0})=Y_{\unicode[STIX]{x1D709}}(\unicode[STIX]{x1D709}=\unicode[STIX]{x03C0})=0$ .

The total energy of these waves $T+V$ is conserved, and takes the form

(4.15) $$\begin{eqnarray}T+V=2\unicode[STIX]{x1D6FC}_{0}^{2}(1-3\unicode[STIX]{x1D6FC}_{0}^{2}).\end{eqnarray}$$

Similar to Stokes waves (Schwartz 1974; Longuet-Higgins 1975), the energy does not monotonically increase with $\unicode[STIX]{x1D6FC}_{0}$ . The energy maximum occurs when $\unicode[STIX]{x1D6FC}_{0}=1/\sqrt{6}$ . This has implications for the stability of the system, as is discussed below.

4.3 Normal mode stability analysis

The linear stability analysis of this system is performed by letting

(4.16a,b ) $$\begin{eqnarray}Y_{1}=(\unicode[STIX]{x1D6FC}_{1}+\unicode[STIX]{x1D716}{\mathcal{A}}(t))\text{e}^{\text{i}ct},\quad \quad Y_{-1}=(\unicode[STIX]{x1D6FC}_{1}+\unicode[STIX]{x1D716}{\mathcal{A}}^{\ast }(t))\text{e}^{-\text{i}ct},\end{eqnarray}$$

where $\unicode[STIX]{x1D716}\ll 1$ is an ordering parameter. Applying Hamilton’s principle to the Lagrangian at second order in $\unicode[STIX]{x1D716}$ , i.e. $O(\unicode[STIX]{x1D716}^{2})$ , we find

(4.17) $$\begin{eqnarray}\mathbb{M}\ddot{\boldsymbol{A}}+2\text{i}c\mathbb{K}\dot{\boldsymbol{A}}+\unicode[STIX]{x1D708}\mathbb{C}\boldsymbol{A}=\mathbf{0},\end{eqnarray}$$

where $(\mathbb{M},\mathbb{K},\mathbb{C})$ are $2\times 2$ matrices and $\boldsymbol{A}=(A,A^{\ast })$ . Defining $\unicode[STIX]{x1D708}=-4\unicode[STIX]{x1D6FC}_{0}^{2}$ and $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D708}(1-2\unicode[STIX]{x1D6FC}_{0}^{2})$ , the matrices are given by

(4.18a-c ) $$\begin{eqnarray}\mathbb{M}=\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D707} & 1+\unicode[STIX]{x1D707}\\ 1+\unicode[STIX]{x1D707} & \unicode[STIX]{x1D707}\end{array}\right],\quad \mathbb{K}=\left[\begin{array}{@{}cc@{}}0 & 1\\ -1 & 0\end{array}\right],\quad \mathbb{C}=\left[\begin{array}{@{}cc@{}}1 & 1\\ 1 & 1\end{array}\right].\end{eqnarray}$$

Note, $\mathbb{M}$ and $\mathbb{C}$ are symmetric while $\mathbb{K}$ is skew–symmetric.

We now let $A=A_{0}\text{e}^{\unicode[STIX]{x1D706}t}$ , where $\unicode[STIX]{x1D706}$ is the normal mode eigenvalue and $(A_{0},A_{0}^{\ast })$ the corresponding eigenvector. The eigenvalues are found by requiring that the determinant of the left-hand side of equation (4.17) vanishes. That is

(4.19) $$\begin{eqnarray}\text{det}\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D707}\unicode[STIX]{x1D706}^{2}+\unicode[STIX]{x1D708} & \unicode[STIX]{x1D706}^{2}(1+\unicode[STIX]{x1D707})+\unicode[STIX]{x1D708}+2c\unicode[STIX]{x1D706}\\ \unicode[STIX]{x1D706}^{2}(1+\unicode[STIX]{x1D707})+\unicode[STIX]{x1D708}-2c\unicode[STIX]{x1D706} & \unicode[STIX]{x1D707}\unicode[STIX]{x1D706}^{2}+\unicode[STIX]{x1D708}\end{array}\right]=0.\end{eqnarray}$$

This implies that the eigenvalues are given by

(4.20) $$\begin{eqnarray}-\unicode[STIX]{x1D706}^{2}=\frac{4(1-6\unicode[STIX]{x1D6FC}_{0}^{2})}{(1-4\unicode[STIX]{x1D6FC}_{0}^{2})^{2}}=\frac{4(1-6\unicode[STIX]{x1D6FC}_{0}^{2})}{c^{4}}.\end{eqnarray}$$

We see immediately that the eigenvalues are real when

(4.21) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{0}>\sqrt{1/6},\end{eqnarray}$$

which is equivalent to the point at which the energy becomes a maximum. This is in agreement with existing studies on Stokes waves (Saffman 1985; MacKay & Saffman 1986).

The corresponding eigenvectors are found by solving equation (4.17) with the eigenvalues substituted in. These are unique up to a rescaling by a complex constant.

4.4 General solutions

In this subsection we consider more general solutions to (4.3) and (4.4). To this end, we consider two integrals of motion: the total energy $H$ , i.e. equation (2.26) and the momentum ${\mathcal{P}}$ , i.e. equation (2.25). These take relatively simple form when we let

(4.22) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D70C}(t)\text{e}^{\text{i}\unicode[STIX]{x1D703}(t)}=\unicode[STIX]{x1D6FD}^{\ast }.\end{eqnarray}$$

Then, we have

(4.23) $$\begin{eqnarray}(1-4\unicode[STIX]{x1D70C}^{2})^{2}\dot{\unicode[STIX]{x1D70C}}^{2}+\unicode[STIX]{x1D70C}^{2}(1+\dot{\unicode[STIX]{x1D703}}^{2}-2\unicode[STIX]{x1D70C}^{2})=H,\end{eqnarray}$$

and

(4.24) $$\begin{eqnarray}2\unicode[STIX]{x1D70C}^{2}\dot{\unicode[STIX]{x1D703}}={\mathcal{P}}.\end{eqnarray}$$

From the last equation we find $\dot{\unicode[STIX]{x1D703}}$ as a function of $\unicode[STIX]{x1D70C}$ , which may be substituted into (4.23) to give

(4.25) $$\begin{eqnarray}(1-4\unicode[STIX]{x1D70C}^{2})^{2}\dot{\unicode[STIX]{x1D70C}}^{2}+\frac{{\mathcal{P}}^{2}}{4\unicode[STIX]{x1D70C}^{2}}+\unicode[STIX]{x1D70C}^{2}(1-2\unicode[STIX]{x1D70C}^{2})=H,\end{eqnarray}$$

which may be rearranged to find

(4.26) $$\begin{eqnarray}{\textstyle \frac{1}{2}}\dot{\unicode[STIX]{x1D70C}}^{2}+\unicode[STIX]{x1D6F1}(\unicode[STIX]{x1D70C})=0,\end{eqnarray}$$

where

(4.27) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}(\unicode[STIX]{x1D70C})=\frac{1}{2}\frac{1}{(1-4\unicode[STIX]{x1D70C}^{2})^{2}}\left(\frac{{\mathcal{P}}^{2}}{4\unicode[STIX]{x1D70C}^{2}}-H+\unicode[STIX]{x1D70C}^{2}(1-2\unicode[STIX]{x1D70C}^{2})\right).\end{eqnarray}$$

Equation (4.26) describes a particle trapped in a potential well, immediately yielding insightful qualitative information about the evolution of the amplitudes. Examples of potentials are shown in figure 2.

Figure 2. Potential $\unicode[STIX]{x1D6F1}(\unicode[STIX]{x1D70C})$ when $\unicode[STIX]{x1D6EC}>0$ (black line) and $\unicode[STIX]{x1D6EC}<0$ (red line), where $\unicode[STIX]{x1D6EC}$ is defined in (4.28) and $\unicode[STIX]{x1D70C}=|\unicode[STIX]{x1D6FC}|$ . For the case of $\unicode[STIX]{x1D6EC}>0$ , the solutions are stable for all times. When $\unicode[STIX]{x1D6EC}<0$ , solutions may be unstable, leading to the formation of a cusp in the free surface as $\unicode[STIX]{x1D70C}\rightarrow 1/2$ .

For $\unicode[STIX]{x1D70C}$ small we see the potential goes to positive infinity, as ${\mathcal{P}}^{2}\geqslant 0$ . For $\unicode[STIX]{x1D70C}\rightarrow 1/2$ , the sign of the potential is given by the sign of

(4.28) $$\begin{eqnarray}\unicode[STIX]{x1D6EC}={\mathcal{P}}^{2}-H+1/8.\end{eqnarray}$$

When $\unicode[STIX]{x1D6EC}<0$ , solutions might lead to cusp formation as the value of $\unicode[STIX]{x1D70C}$ can approach $1/2$ .

5 Stokes waves

The normal mode linear stability of Stokes waves has been studied in detail by, for example, Longuet-Higgins (1978a ), McLean et al. (1981), Tanaka (1983) and MacKay & Saffman (1986). Here, we present an alternative derivation of the linear stability analysis with the added benefit that the system of equations takes a particularly simple form and is general. That is, once the Stokes coefficients (this may, for instance, be for class $M$ Stokes waves (Saffman 1980)) are determined, the stability analysis follows readily by solving an eigenvalue problem with prescribed matrix entries which are low-order polynomials in the dependent variables, making this approach particularly suitable for numerical implementation.

5.1 Stokes waves

We now consider permanent progressive solutions to the system where $N$ is taken to large enough to properly resolve the physical scenario in question.

To this end, we note that Stokes wave solutions to our system of equations take the simple form (Balk 1996)

(5.1) $$\begin{eqnarray}Y_{k}=\unicode[STIX]{x1D6FC}_{k}\text{e}^{\text{i}ckt},\end{eqnarray}$$

so that $B_{k}=cY_{k}$ . The Lagrangian can then be written as

(5.2) $$\begin{eqnarray}L=-\frac{c^{2}Y_{0}}{2}+V,\end{eqnarray}$$

with $V$ given by equation (2.7). Therefore, the system governing the coefficients $\unicode[STIX]{x1D6FC}_{k}$ is (Longuet-Higgins (1985), see also Longuet-Higgins (1978c ))

(5.3) $$\begin{eqnarray}F_{n}=\unicode[STIX]{x1D6FC}_{n}(1+n\unicode[STIX]{x1D6FC}_{0})+\frac{1}{4}\sum (|n|+2|\,j|)\unicode[STIX]{x1D6FC}_{j}\unicode[STIX]{x1D6FC}_{k},\end{eqnarray}$$

for $F_{n}=0$ . There are $N+1$ equations here, so we must prescribe a parameter to close the system. Generally, the phase speed $c$ has been used. However, $c$ does not increase monotonically with increasing slope and contains critical points, which leads to loss in numerical accuracy (see § 5 of Longuet-Higgins (1985)). Alternatively, a variety of other parameters which do increase monotonically with increasing $ak$ have been proposed. Following Longuet-Higgins (1985), we take the parameter $Q=1-q_{crest}^{2}$ , where $q_{crest}$ is the particle speed at the crest. From Bernoulli’s equation, we have $Q=1+Y(0)=1+\frac{1}{2}\unicode[STIX]{x1D6FC}_{0}+\unicode[STIX]{x1D6FC}_{1}+\unicode[STIX]{x1D6FC}_{2}+\cdots \,.$ To solve the system of equations, i.e. (5.3), we employ a Newton–Raphson-type method.

5.2 Linear stability analysis

To find the form of the linear normal modes, we next let

(5.4) $$\begin{eqnarray}Y_{k}=(\unicode[STIX]{x1D6FC}_{k}+\unicode[STIX]{x1D716}A_{k}\text{e}^{\unicode[STIX]{x1D706}t})\text{e}^{\text{i}kct},\end{eqnarray}$$

where $\unicode[STIX]{x1D716}$ is a small parameter, $\unicode[STIX]{x1D706}$ is the normal mode eigenvalue and $A_{k}$ are the coefficients of the eigenvector.

The Lagrangian to $O(\unicode[STIX]{x1D716}^{2})$ becomes

(5.5) $$\begin{eqnarray}L=\frac{1}{2}\mathop{\sum }_{-N}^{N}|k|B_{k}B_{-k}+\frac{1}{2}Y_{0}^{2}-\mathop{\sum }_{1}^{N}Y_{k}Y_{-k}-\frac{1}{2}\mathop{\sum }_{j+k+l=0}^{\prime }|l|Y_{j}Y_{k}Y_{l},\end{eqnarray}$$

where

(5.6) $$\begin{eqnarray}B_{k}=\frac{\text{i}}{k}\left(-{\dot{Y}}_{k}-{\dot{Y}}_{0}Y_{k}|k|+\sum ^{\prime }{\dot{Y}}_{j}Y_{l}l(\unicode[STIX]{x1D70E}_{j}-\unicode[STIX]{x1D70E}_{l})\right),\end{eqnarray}$$

and primes represent sums over all indices except 0. Explicitly, the mean surface displacement is given by

(5.7) $$\begin{eqnarray}Y_{0}=-2\mathop{\sum }_{1}^{N}|k|\unicode[STIX]{x1D6FC}_{k}^{2}-\unicode[STIX]{x1D716}\mathop{\sum }_{-N}^{N}|k|\unicode[STIX]{x1D6FC}_{k}A_{-k}-\unicode[STIX]{x1D716}^{2}\mathop{\sum }_{-N}^{N}|k|A_{k}A_{-k},\end{eqnarray}$$

so that its time derivative is

(5.8) $$\begin{eqnarray}{\dot{Y}}_{0}=-\unicode[STIX]{x1D716}\mathop{\sum }_{-N}^{N}|k|\unicode[STIX]{x1D6FC}_{k}{\dot{A}}_{-k}-2\unicode[STIX]{x1D716}^{2}\mathop{\sum }_{-N}^{N}|k|{\dot{A}}_{k}A_{-k}.\end{eqnarray}$$

Substituting these relationships into the $B_{k}$ term, we find

(5.9) $$\begin{eqnarray}\displaystyle B_{k} & = & \displaystyle c\unicode[STIX]{x1D6FC}_{k}+\frac{\text{i}\unicode[STIX]{x1D716}}{k}\left(-({\dot{A}}_{k}+\text{i}ckA_{k})+\unicode[STIX]{x1D6FC}_{k}|k|\mathop{\sum }_{-N}^{N}|l|\unicode[STIX]{x1D6FC}_{l}{\dot{A}}_{-l}\right.\end{eqnarray}$$

(5.10) $$\begin{eqnarray}\displaystyle & & \displaystyle +\left.\sum [2\text{i}c(\unicode[STIX]{x1D6FC}_{j}A_{l}+\unicode[STIX]{x1D6FC}_{l}A_{j})+(\unicode[STIX]{x1D6FC}_{j}{\dot{A}}_{l}+\unicode[STIX]{x1D6FC}_{l}{\dot{A}}_{j})]\,jl(\unicode[STIX]{x1D70E}_{j}-\unicode[STIX]{x1D70E}_{l})\right)\end{eqnarray}$$

(5.11) $$\begin{eqnarray}\displaystyle & & \displaystyle +\,\frac{\text{i}\unicode[STIX]{x1D716}^{2}}{k}\left(-{\dot{Y}}_{0}^{(2)}|k|\unicode[STIX]{x1D6FC}_{k}-{\dot{Y}}_{0}^{(1)}|k|A_{k}+\sum ^{\prime }{\dot{A}}_{j}A_{l}l(\unicode[STIX]{x1D70E}_{j}-\unicode[STIX]{x1D70E}_{l})\right),\end{eqnarray}$$

where superscripts imply we are taking terms of that order from the expansion.

This is sufficient to define the Lagrangian. Applying the Euler–Lagrange equations, we find that the governing equation takes the form (cf. the discussion for the $N=1$ case in § 4) of a quadratic eigenvalue problem (Tisseur & Meerbergen 2001). Namely, we have

(5.12) $$\begin{eqnarray}(\unicode[STIX]{x1D706}^{2}\mathbb{M}+\unicode[STIX]{x1D706}\mathbb{C}+\mathbb{K})\boldsymbol{A}=\mathbf{0},\end{eqnarray}$$

where

(5.13) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{n,m}=\frac{\text{i}}{n}\unicode[STIX]{x1D6FF}_{-n,m}+\frac{\text{i}}{n}\unicode[STIX]{x1D6FC}_{-n-m}(n+m)(\unicode[STIX]{x1D70E}_{m}-\unicode[STIX]{x1D70E}_{-n-m})-2\text{i}\unicode[STIX]{x1D70E}_{n}\unicode[STIX]{x1D6FC}_{n}|m|\unicode[STIX]{x1D6FC}_{m},\end{eqnarray}$$

for $\unicode[STIX]{x1D6FF}_{m,n}$ the Kronecker delta function, while

(5.14) $$\begin{eqnarray}\displaystyle \mathbb{M}_{n,m}=-\text{i}\unicode[STIX]{x1D70E}_{n}\unicode[STIX]{x1D6FE}_{n,m}+2\text{i}|n|\unicode[STIX]{x1D6FC}_{n}\left(\mathop{\sum }_{-N}^{N}l\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D6FE}_{l,m}\right)+\text{i}\mathop{\sum }_{-N}^{N}\unicode[STIX]{x1D70E}_{l}\unicode[STIX]{x1D6FC}_{l-n}(l-n)\unicode[STIX]{x1D6FE}_{l,m}(\unicode[STIX]{x1D70E}_{n}-\unicode[STIX]{x1D70E}_{n-l}), & & \displaystyle\end{eqnarray}$$

(5.15) $$\begin{eqnarray}\displaystyle \mathbb{C}_{n,m} & = & \displaystyle 2\text{i}c\unicode[STIX]{x1D6FC}_{n+m}\unicode[STIX]{x1D70E}_{n+m}(\unicode[STIX]{x1D70E}_{m}-\unicode[STIX]{x1D70E}_{n})G(n,m)+2\text{i}cn|m|\unicode[STIX]{x1D6FC}_{n}\unicode[STIX]{x1D6FC}_{m}+c|n|\unicode[STIX]{x1D6FE}_{n,m}\nonumber\\ \displaystyle & & \displaystyle +\,\text{i}c\unicode[STIX]{x1D70E}_{n}\unicode[STIX]{x1D6FF}_{-n,m}+\text{i}c\unicode[STIX]{x1D70E}_{-m}(n+m)\unicode[STIX]{x1D6FC}_{n+m}(\unicode[STIX]{x1D70E}_{n}-\unicode[STIX]{x1D70E}_{-n-m})G(n,m),\end{eqnarray}$$

(5.16) $$\begin{eqnarray}\displaystyle \mathbb{K}_{n,m} & = & \displaystyle (1-c^{2}|n|)\unicode[STIX]{x1D6FF}_{-n,m}-4\left(\mathop{\sum }_{1}l\unicode[STIX]{x1D6FC}_{l}^{2}\right)|m|\unicode[STIX]{x1D6FF}_{n,-m}\nonumber\\ \displaystyle & & \displaystyle +\,-4|nm|\unicode[STIX]{x1D6FC}_{n}\unicode[STIX]{x1D6FC}_{m}+(|n|+|m|+G(n,m)|n+m|)\unicode[STIX]{x1D6FC}_{n+m},\end{eqnarray}$$

where $G(n,m)=1$ if $-N-m\leqslant n\leqslant N-m$ and 0 otherwise. The entries of these matrices are low-order polynomials in the Stokes coefficients, making them particularly suitable for numerical implementation. Furthermore, these results are general in the sense that one may examine subharmonic instabilities (Longuet-Higgins 1978b ) by considering class $M$ , for $M$ a positive integer, Stokes waves (Saffman 1980; Zufiria 1987).

To validate these algebraically complex computations, the eigenvalues corresponding to superharmonic perturbations of steep Stokes waves are computed from this method and are shown on top of the previous computations due to Tanaka (1983) in figure 3. Agreement between the predictions of equation (5.12) and Tanaka’s results are found. Furthermore, the change in (minus) total energy versus slope, $ak$ , is shown in the plot, showing that the change in stability corresponds to the critical point in the energy.

Figure 3. A comparison of the eigenvalues squared of Tanaka (1983) (red dashed-dot line) versus the prediction given by Balk’s system, i.e. equation (5.12), shown by the black dots. The change in (minus) energy, $-E$ , with respect to $ak$ is shown by the blue dashed line.

6 Conclusion

In this paper we have transformed the equations of Balk into a form more suitable for theoretical and numerical examination. Furthermore, we have derived the equations of motion that result from Balk’s Lagrangian. This represents a generalization of the work done by Longuet-Higgins (2000, 2001) for standing waves. The system of equations are found for independent coefficients $Y_{n}(t),n=1,2\ldots ,$ with partial wave solutions found by truncating the system at $n=N$ . Solutions exist as long as the matrix $Q$ , i.e. equation (3.17), is invertible.

A one mode example was considered, and the stability of permanent progressive waves was examined. Linear stability conditions were derived, and the evolution of the fully nonlinear system was considered. The linear stability of Stokes waves (here the number of modes $N$ was increased until convergence was found) was then inspected. A general system of equations to derive the eigenvalues was presented, and comparison was made to the values found by Tanaka (1983).

The main advantage of Balk’s approach is its simplicity, allowing for easy numerical implementation. The variational method allows for a priori knowledge of the necessarily truncated integrals of motion present during numerical implementation, providing additional numerical benchmarks to test the fidelity of the integration.

The spectral formulation considered here suffers from the usual drawback of slow series convergence when singularities in the upper half-complex plane approach the free surface (Fornberg 1980; Baker & Xie 2011). However, this may be overcome to a certain extent through conformal mapping techniques (e.g. Tanaka 1983), which increase the resolution in these regions of high curvature, effectively moving the singularities further away from the free surface, accelerating the rate at which the Fourier series converge. More formally, following Tanveer (1991), one may study the singularities in the unphysical plane to improve numerical efficiency. Related problems, such as the examination of the Rayleigh–Taylor and Richtmyer–Meshkov instability, may be studied in a similar fashion (Baker, Meiron & Orszag 1980; Tanveer 1991).

An application of these equations to understanding the nonlinear evolution of the normal modes of steep Stokes waves (Longuet-Higgins & Dommermuth 1997; Bridges 2004) will be discussed in a future paper.