## Abstract

High-order asymptotic series are obtained for two- and three-dimensional gravity-capillary solitary waves. In two dimensions, the first term in the asymptotic series is the well-known sech^{2} solution of the Korteweg–de Vries equation; in three dimensions, the first term is the rational lump solution of the Kadomtsev–Petviashvili equation I. The two-dimensional series is used (with nine terms included) to investigate how small surface tension affects the height and energy of large-amplitude waves and waves close to the solitary version of Stokes’ extreme wave. In particular, for small surface tension, the solitary wave with the maximum energy is obtained. For large surface tension, the two-dimensional series is also used to study the energy of depression solitary waves. Energy considerations suggest that, for large enough surface tension, there are solitary waves that can get close to the fluid bottom. In three dimensions, analytic solutions for the high-order perturbation terms are computed numerically, and the resulting asymptotic series (to three terms) is used to obtain the speed versus maximum amplitude curve for solitary waves subject to sufficiently large surface tension. Finally, the above asymptotic method is applied to the Benney–Luke (BL) equation, and the resulting asymptotic series (to three terms) is verified to agree with the solitary-wave solution of the BL equation.

## 1. Introduction

In this paper, we develop high-order asymptotic series for capillary-gravity solitary waves in two and three dimensions. We use the two-dimensional series to investigate large-amplitude depression and elevation waves and the three-dimensional series to investigate depression waves in the presence of large surface tension.

The earliest recorded observation of a solitary water wave dates back to Russell (1844), who observed an unchanging, steadily progressing wave propagating down a shallow channel. Scott Russell’s observations were later investigated by Boussinesq (1877) and Korteweg & de Vries (1895), wherein an asymptotic theory for long waves (relative to the fluid depth) was developed. Contemporaneously, Stokes (1880) studied periodic solutions of steady water-wave equations. In particular, Stokes (1880) conjectured that steady periodic waves tend to a limiting wave—referred to here as Stokes’ extreme wave—which has sharp angles of 120^{°} at wave crests and is convex between successive crests.

A rigorous theory for small-amplitude, steady, periodic water waves was developed in Levi-Civita (1925) and Nekrasov (1921) and was later extended to large-amplitude water waves and Stokes’ extreme wave using techniques from global bifurcation theory (cf. Krasovskii 1961; Keady & Norbury 1978; Toland 1978; Plotnikov 1980, 1982; Amick *et al.* 1982; McLeod 1997; Plotnikov & Toland 2004). In Amick & Toland (1981*a*), it was shown that periodic solutions of the steady water-wave equations tend to solitary waves as the period goes to infinity, thereby establishing the existence of supercritical (*c*^{2}>*g**h*) solitary waves that are symmetric, positive and monotonically decaying away from the wave crest (see also Craig & Sternberg 1988; McLeod 1997). A comprehensive global bifurcation theory for solitary water-wave solutions can be found in Amick & Toland (1981*b*).

Asymptotic techniques have also proven fruitful in the study of large-amplitude solitary waves. Of particular relevance to this paper, in Fenton (1972), a high-order asymptotic series for the free surface was developed, where the first term in the series is the sech^{2} of the Korteweg–de Vries (KdV) equation. Although Fenton’s series is derived under the assumption that the free surface is small, it is used to predict the height, speed and energy of the maximal wave to high accuracy (Longuet-Higgins & Fenton 1974). Specifically, in Longuet-Higgins & Fenton (1974), Fenton’s series is used to predict the height of the solitary extreme wave as 0.827, which is consistent with independent numerical computations (Yamada 1957; Lenau 1966). An important observation made in Longuet-Higgins & Fenton (1974) is that the maximally energetic wave has a smaller amplitude than the extreme solitary wave.

More recently, a global bifurcation theory for steady capillary-gravity water waves has been developed (cf. Groves & Toland 1997). In particular, it has been shown that the (*σ*,*c*)-parameter plane (*σ* is the normalized surface tension and *c* is the normalized wave speed) is partitioned into four distinct regions. Depending on where (*σ*,*c*) falls in the parameter plane, there are either periodic waves and depression solitary waves (cf. Amick & Toland 1981*a*; Kirchgassner 1988; Sachs 1991 for analysis; cf. Hunter & Vanden-Broeck 1983; Vanden-Broeck & Dias 1992; Dias & Iooss 1993 for computational studies), periodic waves modulated by exponentially decaying envelopes (Iooss & Kirchgassner 1990), pulse-like elevation waves (called ‘generalized solitary waves’) that decay to exponentially small periodic ripples at infinity (see Beale 1991; Iooss & Kirchgassner 1992; Sun 1991; Sun & Shen 1993 for analysis; see Hunter & Vanden-Broeck 1983; Vanden-Broeck & Dias 1992; Dias & Iooss 1993 for computational studies), and pulse-like, depression waves with large troughs separated by small oscillations (Buffoni *et al.* 1996). These results show that, in particular, depression solitary waves exist when *σ*>1/3, and generalized solitary waves exist when 0<*σ*<1/3. It is still open as to whether there are true solitary water waves for 0<*σ*<1/3. However, the non-existence of true solitary water waves is suggested by both numerical computations (Turner & Vanden-Broeck 1991) and asymptotic methods (Akylas & Grimshaw 1992). For *σ* close enough to 1/3, it is shown in Sun (1999) that there are no exponentially decaying solitary waves.

Unlike their two-dimensional counterpart, there is much less known about solitary water waves in three dimensions. The existence of fully localized solitary waves is suggested by lump-type solutions to the Kadomtsev–Petviashvili equation I (KP-I) (cf. Ablowitz & Clarkson 1991) and computations on the free-surface Euler equations (cf. Vanden-Broeck & Dias 1992; Dias *et al.* 1996; Parau *et al.* 2005). Analysis of the KP-I equation indicates that sufficient surface tension is required for the existence of localized solitary water waves; this was confirmed analytically in Craig (2002), wherein it is shown that localized solitary waves cannot exist if there is no surface tension. Recently, the existence of localized, three-dimensional solitary waves in the presence of sufficient surface tension has been rigorously proven in Groves & Sun (2008). In addition to the above lump-type solutions, another class of capillary-gravity waves is found in Kim & Akylas (2005) when the surface tension is sufficiently small; these solitary waves come in the form of locally confined wavepackets and are governed, to leading order, by the Davey–Stewartson equation in finite depth and a two-dimensional nonlinear Schrodinger equation in infinite depth (see also Yang & Akylas 1997; Dias & Kharif 1999).

Here, we use a recent non-local formulation of water waves (Ablowitz *et al.* 2006)—referred to in this paper as the non-local spectral (NSP) formulation—to examine capillary-gravity solitary waves in two and three dimensions. As discussed above, there is an extensive literature on capillary-gravity solitary waves; however, as far as we know, no one has obtained high-order asymptotic series for gravity-capillary waves. Our approach is to employ the Poincare–Stokes procedure, which has the advantage of being straightforward and generalizable to three dimensions. While Fenton’s asymptotic series was used to accurately obtain the height, energy and speed of the solitary extreme wave (Longuet-Higgins & Fenton 1974), real water waves are subject to small surface tension, and it is natural to ask how small capillarity affects the energy and height of waves close to the extreme wave.

In two dimensions, we find from the NSP formulation that the asymptotic series for the wave speed *c*, free surface *η* and tangential surface velocity *U* (*X*=*c*−*x**t*) assume the form
1.1
and
1.2

Here, *c*_{i},*α*, *d*_{i−1,j} and *f*(*i*−1,*j*) depend on the water-wave parameters, including surface tension; the coefficients are explicitly given in Haut (2008) for *N*=9. Equation (1.1) is a generalization of Fenton’s (1972) series to fluids with non-zero surface tension. We choose our non-dimensional scaling so that *ε*=*η*_{max}/*h*, where *h* is the depth of the undisturbed fluid and *η*_{max} is the maximum (dimensional) wave height.

Using equations (1.1) and (1.2) gives a high-order representation for the energy of capillary-gravity solitary waves. Using this, we study how small capillarity affects large-amplitude waves and waves close to the solitary extreme wave. As we mentioned above, Fenton’s series has been used to show (Longuet-Higgins & Fenton 1974) that the solitary wave with the greatest energy has a smaller height than the extreme wave. Physically, it is natural to expect that waves with energy greater than the maximally energetic wave are unstable and are, therefore, unlikely to be observed. Using the zero surface tension case as a guide, we use equations (1.1) and (1.2) to examine how the height, energy and speed of the maximally energetic wave are perturbed in the presence of small surface tension. We find that the height and energy of the maximally energetic wave both decrease with increasing surface tension.

It is important to note that, as previously discussed, there may not be true solitary water waves when small surface tension is included. In particular, it was established (Beale 1991; Sun 1991; Iooss & Kirchgassner 1992; Sun & Shen 1993) that there are travelling waves *η* that satisfy, for 0<*σ*<1/3,
where
is the solitary-wave solution of the KdV equation and the remainder term *R*(*x*,*ε*) has exponentially small oscillations as , for *ε*≪1. That is, if one were to evolve the water-wave equations with the series (1.1) and (1.2) as the initial conditions, extremely small transients might develop in the tails. We are assuming here that the exponentially small oscillations at infinity can be neglected in the computation of the energy and height of the (near) extreme wave. For zero surface tension, equations (1.1) and (1.2) allow one to predict the height of the solitary extreme wave as *ε*=0.827 (Longuet-Higgins & Fenton 1974), which is consistent with independent numerical calculations (cf. Yamada 1957; Lenau 1966).

Solitary depression waves, for 0<*σ*<1/3, are also computed from equations (1.1) and (1.2). We find that, for large enough surface tension, the asymptotic series for *η* (with *N*=9) monotonically increases to zero away from its minimum at the origin, for all 0<*ε*<1. For smaller surface tension, the series begins to develop oscillations as *ε* increases to one. Since we choose our non-dimensionalization so that , this suggests that, for large enough surface tension, there are solitary waves that get close to the fluid bottom. Additionally, equations (1.1) and (1.2) are used to compute how the energy and speed of depression solitary waves depend on surface tension. We find that the energy of the maximally energetic solitary wave increases with larger surface tension. Moreover, for fixed *σ*≥2, the energy increases, as a function of the non-dimensional wave depth *ε*, until its maximum value is obtained just slightly below *ε*=1. This again indicates that, for large enough surface tension, there are solitary waves that can get close to the fluid bottom.

In the second part of this paper, we employ the Poincare–Stokes method on the (2+1)-dimensional NSP formulation of water waves. The first term in the resulting series satisfies the KP-I equation and has an explicit rational solution (cf. Ablowitz & Clarkson 1991). Although analytic solutions for the higher order terms are not available, we solve for the next two perturbations numerically, and use the resulting series for the free surface to obtain how the wave speed depends on the maximum wave amplitude.

As a check on the validity of the above asymptotic procedure, we employ the Poincare–Stokes method on the Benney–Luke (BL) equation (Benney & Luke 1964), which is known to have lump-type solutions (cf. Pego & Quintero 1999; Berger & Milewski 2000). The benefit of using the BL equation is that we can compute a solitary-wave solution to this equation in a straightforward manner (see Ablowitz & Musslimani 2005). As in the full water-wave case, the first term in the resulting asymptotic series is the rational solution to the KP-I equation, and the next two perturbation terms are obtained numerically. We verify that the computed asymptotic series gives an accurate approximation to the exact (numerically computed) solution of the BL equation.

## 2. Equations of motions

We use the following reformulation of water waves (Ablowitz *et al.* 2006), which is expressed in terms of the free surface *η* and *q*=*φ*(*x*,*η*,*t*), the velocity potential evaluated on the free surface
2.1
and
2.2
In equations (2.1) and (2.2), *x*=(*x*_{1},*x*_{2}), *k*=(*k*_{1},*k*_{2}), *k**x*≡*k*⋅*x* and the constants *g*, *ρ* and *σ* denote gravity, density and surface tension, respectively. Equation (2.2) is the well-known Bernoulli equation, formulated in terms of *η* and *q*. Equation (2.1) is a non-local equation that depends on a free spectral parameter *k*. Equations (2.1) and (2.2) will be referred to in this paper as the NSP formulation of water waves, and equation (2.1) will be referred to as the NSP equation. An advantage of the above formulation is that the vertical component of the velocity is eliminated, thus reducing the dimensionality and fixing the domain in which the equations are posed; for this reason, the NSP formulation is convenient for the subsequent asymptotic calculations.

We now non-dimensionalize equations (2.1) and (2.2). To do so, we define the non-dimensional variables *x*′, *t*′, *η*′ and *q*′ by
where *l* is the characteristic wave length, *a* is the characteristic wave height and *γ* is a non-dimensional parameter. Then, equations (2.1) and (2.2) become, upon dropping primes,
2.3
and
2.4
In equations (2.3) and (2.4), *ε*≡*a*/*h* is a measure of the nonlinearity and *μ*≡*h*/*l* is a measure of dispersion.

Finally, since we are interested in studying steady solutions of equations (2.3) and (2.4), we change coordinates by taking *X*=*x*−*c*_{X}*t* and *Y* =*y*−*c*_{Y}*t*. Equations (2.3) and (2.4) become, in *X* and *Y* variables,
2.5
and
2.6

## 3. Derivation of asymptotic series in two dimensions

In this section, we develop high-order, two-dimensional asymptotic series for the free surface, tangential surface velocity and wave speed. Although we do not give all terms of order *O*(*ε*^{j}), *j*=1,…,9, our calculations are carried out to nine terms using Mathematica (see Haut 2008).

We start with the two-dimensional reduction of equations (2.5) and (2.6)
3.1
and
3.2
In equations (3.1) and (3.2), *U*=∂*η*/∂*X*, and we assume that
3.3

We now expand equation (3.1) in *ε*, and take the inverse Fourier transform of the resulting expression
3.4
Using equation (3.4), we solve (asymptotically) for *U* in terms of *η*,
3.5
Similarly, expanding equation (3.2) in *ε* yields
3.6
We use equation (3.5) in equation (3.6),
3.7

To obtain an analytic solution to equation (3.7), we assume that *η* has the form
3.8
where *α* is an arbitrary constant; we have also taken *c*_{1}=−2/3(3*α*^{2}*σ*−*α*^{2}) so that equation (3.7) is satisfied to leading order. Putting the above ansatz for *η* into equation (3.7), using and setting powers of sech^{2}(*X**α*) to zero yield a linear system of equations for *d*_{1,1}, *d*_{1,2} and *c*_{2}, with solution
Therefore, *η* is given by
3.9
Finally, using equation (3.9) in equation (3.5) gives us the following expression for *U*:
3.10

Until now, the parameters *α* and *ε* have been arbitrary. We now choose *α*=*α*(*ε*) in equation (3.9) so that
3.11
For example, setting *X*=0 in equation (3.9) and considering the case when *σ*>1/3,
Assuming that *α* can be expanded in terms of *ε* as *α*^{2}=*α*_{0}+*α*_{1}*ε*+*O*(*ε*^{2}), the previous equation implies that
In a similar way, we find that, when 0<*σ*<1/3,

If the maximum amplitude of *η* occurs at the origin, then *ε* is the (non-dimensional) maximum amplitude.

Importantly, we can also obtain a representation of *c* in terms of *ε* by an asymptotic argument. A similar argument will be used in three dimensions, where there are no analytic solutions available for the terms in the series. Namely, we assume that *η*(*X*)∼e^{−2αX} as . Therefore, *η*^{(m)}(*X*)∼(−2*α*)^{m}e^{−2α} as , and we deduce from equation (3.7) that
3.12
Equating the coefficients of *ε* in equation (3.12) to zero, we obtain, as before, expressions for *c*_{1} and *c*_{2} in terms of *α*.

## 4. Maximum energy and height of solitary waves with surface tension

We now investigate how the maximally energetic solitary wave depends on small surface tension. The rationale for looking at the maximally energetic wave is that it is natural to expect that larger waves are unstable and are, therefore, unlikely to be observed.

We begin by substituting the series (1.1) and (1.2), with *N*=9, into the following formula for the total energy
where KE and PE denote the kinetic and potential energies, respectively.

Doing so results in a series for *E*
4.1

We now express the previous series for the energy in terms of a new coordinate, ω, defined by ω=1−(*ε**U*(0)−*c*)^{2}. The significance of ω is that ω=1 corresponds to a stagnation point at the wave crest and was used in Longuet-Higgins & Fenton (1974) to accurately calculate the height, speed and energy of the maximal wave for zero surface tension. From our analytic expressions for *c* and *U*, it follows that
We invert the previous series for *ε* in terms of ω, neglecting terms of order *O*(ω^{9})
4.2
Finally, we use equation (4.2) in equation (4.1) to obtain an expression for the energy in terms of ω
4.3

We approximate the resulting energy (4.3) and wave height *ε* (4.2) by Pade(3,3) rational functions. As a check, we re-derive the wave height of the extreme wave (Longuet-Higgins & Fenton 1974) as *ε*=0.827 by setting ω=1 in the Pade(3,3) function obtained from equation (4.2).

Figures 1 and 2 show the energy *E* as a function of ω and *ε*, respectively, for different values of *σ*. We see from figure 1 that the maximal energy in each case peaks before the stagnation point is reached (recall that this corresponds to ω=1). Figures 3 and 4 display, as a function of surface tension, the energy and height of the maximally energetic wave; we see from these figures that both the height and energy decrease as a function of the surface tension.

## 5. Solitary waves of depression

In this section, we use equations (1.1) and (1.2), for *σ*>1/3, to study depression waves.

Figure 5 shows plots of *η* for fixed *σ*=2/3 and *ε*=0.1, 0.2, 0.3 and 0.4. As we can see from these graphs, the series for *η* has oscillations in *X* for large enough *ε*. In particular, when *ε*≈0.4, the graph of *η* assumes positive maxima; therefore, the assumed scaling that was used to obtain *α* as a function of *ε* is not relevant.

In contrast, figure 6 shows plots of *η* for *σ*=1, 1.1, 1.2 and 1.3 and *ε*=0.95. These plots highlight that, when , the series for *η* is observed to monotonically increase to zero away from its minimum at the origin, for every value of 0<*ε*<1. This suggests that, for large enough surface tension, there can be solitary waves that get close to the fluid bottom.

As in §4, we solve for the energy *E* in terms of the wave depth *ε*
5.1
Our discussion below of the energy *E* is based on the Pade(3,3) approximation of equation (5.1).

Figure 7 shows plots of the energy *E* as a function of the (non-dimensional) wave depth *ε*, for *σ*=2,5,10 and 100. It can be observed from these plots that the maximal energy, for fixed *σ*≥2, is obtained by solitary waves whose depth is just slightly less than the fluid bottom depth of *ε*=1. Indeed, for *σ*≥2, we find numerically that the maximal energy is always obtained with 0.96≤*ε*≤0.99. This also indicates that, for large enough surface tension, there are solitary waves that get close to the fluid bottom.

Finally, figure 8 displays that, for all values of surface tension examined, the speed *c* monotonically decreases as a function of maximal wave depth *ε*. To generate this figure, we used a Pade(3,3) rational approximation of the series for *c*.

## 6. Derivation of the asymptotic series in three dimensions

In this section, we employ the Poincare–Stokes method on equations (2.5) and (2.6) in three dimensions. That is, we assume that
We find that *q*_{0} satisfies the KP-I equation,
6.1
and *q*_{1} and *q*_{2} satisfy
6.2
6.3
In equations (6.2) and (6.3), *L* is the differential operator
*F*_{1} is an explicit expression depending on *q*_{0} and derivatives of *q*_{0} and *F*_{2} is an expression depending on *q*_{0}, *q*_{1} and derivatives of *q*_{0} and *q*_{1}. Because of the size of these expressions, they are not explicitly displayed here.

Equation (6.1) has the rational solution (cf. Ablowitz & Clarkson 1991)
6.4
Although analytic solutions to equations (6.2) and (6.3) are not obtained, we compute *q*_{1} and *q*_{2} numerically and use the resulting asymptotic series for *q* to obtain the speed versus maximum amplitude relationship for *η*.

To illustrate this method, we derive the first perturbation to the KP-I equation from equations (2.5) and (2.6). To do so, we first expand equation (2.6) in *ε*
6.5
We (asymptotically) solve for *η* in terms of *q* to get
6.6

We now expand equation (2.5) in *ε* and take the inverse Fourier transform of the resulting expression
6.7
Using equation (6.6) to express *η* in terms of *q* in equation (6.7) yields
6.8

Finally, taking *q*=*q*_{0}+*ε**q*_{1}+*ε*^{2}*q*_{2}+··· in equation (6.8) and setting the coefficients in the resulting powers of *ε* to zero gives equations (6.1)–(6.3).

Using equation (6.8), we now determine the unknown speed perturbations *c*_{X1}, *c*_{Y1} and *c*_{X2}, and *c*_{Y2} in terms of the two ‘width’ parameters *α*_{1} and *α*_{2} occurring in the lump solution (6.4). In the one-dimensional case, the free surface *η* satisfies
6.9
where *η*_{0} is the sech^{2} solution to the KdV equation and depends on a free parameter *α*. We found, in this case, that *c*_{1},*c*_{2},…, are uniquely determined in terms of *α* if we assume that *η*(*X*) ∼ *η*_{0}(*X*) ∼ e^{−2α|X|} as .

Analogously, in the three-dimensional case, *q* satisfies
6.10
where *q*_{0} is the rational solution to the KP-I equation (6.1). It can be seen from equation (6.4) that *q*_{0} depends on the two free parameters *α*_{1} and *α*_{2}. Using the one-dimensional case as a guide, we assume that the full solution *q* decays like the leading-order lump solution *q*_{0} at infinity. Specifically, using the assumption
6.11
in equation (6.8) yields an equation of the form *P*(*X*,*Y*)/*Q*(*X*,*Y*)∼0, where *Q* and *P* are polynomials that depend on *ε*, *α*_{1}, *α*_{2}, *c*_{X1}, *c*_{Y1}, *c*_{X2} and *c*_{Y2}. More specifically, *P*(*X*,*Y*) is a polynomial of the form *A**X*+*B**Y* +*C**X*^{3}+*D**Y*^{3}+*E**X*^{2}*Y* +*F**X**Y*^{2}, where *A*, *B*, *C*, *D*, *E* and *F* depend on *α*_{1}, *α*_{2}, *c*_{X1}, *c*_{Y1}, *c*_{X2} and *c*_{Y2}. Neglecting the *X* and *Y* terms (since we are interested only in leading-order decay) and setting the coefficients *C*,…,*F* to zero, we obtain four equations (two of which are degenerate) relating *α*_{1}, *α*_{2}, *c*_{X1}, *c*_{Y1}, *c*_{X2} and *c*_{Y2}. Equating coefficients of *ε*^{0} and *ε*^{1} to zero and solving the resulting system of equations, we find that
6.12
6.13
6.14
6.15

### (a) Numerical solution of the reduced equations

In this section, we describe how to numerically solve for *q*_{1} and *q*_{2}. As *q*_{0} is given by equation (6.4) and *q* and *η* satisfy
once we have *q*_{1} and *q*_{2}, we can obtain an approximation to *q* and *η* to order *O*(*ε*^{2}). By solving for *q* and *η* using a range of values for *α*_{1} and *α*_{2}, we compute, in §6*b*, the speed versus amplitude relationship for *η*.

We first change variables by letting
Then, in the new variables, the equations for *q*_{1}′(*X*′,*Y* ′) and *q*_{2}′(*X*′,*Y* ′) become, after dropping primes,
6.16
and
6.17
where *q*_{0}(*X*,*Y*)=*X*/(*X*^{2}+*Y*^{2}+1). Notice that, with this change of variables, the linear operator occuring on the left-hand side of equations (6.16) and (6.17) is independent of the parameters *α*_{1} and *α*_{2}, and therefore only needs to be inverted once.

We solve equations (6.16) and (6.17) numerically, employing two different procedures.

(i) In our first numerical scheme, we solve equations (6.16) and (6.17) for *q*_{1} and *q*_{2} by employing a sixth-order finite-difference scheme with periodic boundary conditions. Equations (6.16) and (6.17) can then be written as the matrix equations
where *q*_{jmn} ≈ *q*_{j}(*X*_{m},*Y*_{n}), *X*_{m} = −*a* + 2*a*/*N**m* and *Y*_{n} = −*a* + 2*a*/*N**n*, and the vector *F*_{j} approximates the function *F*_{j}(*q*_{0},…,*q*_{j−1})(*X*,*Y*) at the grid points *X*=*X*_{m} and *Y* =*Y*_{m}. Since the expression *F*_{1}(*q*_{0}) only depends on *q*_{0}, this can be evaluated analytically and the resulting vector *F*_{1} computed in a straightforward way. However, since the expression *F*_{2}(*q*_{0},*q*_{1}) depends on *q*_{0}, *q*_{1} and derivatives of *q*_{0} and *q*_{1}, some scheme is needed to approximate derivatives of *q*_{1} from the approximation *q*_{1}.

To approximate *q*_{1}^{(r,s)}(*X*_{m},*Y*_{n}), we start with the following pair of equations:
6.18
and
6.19
Then, with *X*_{m}=−*a*+2*a*/*N**m* and *Y*_{n}=−*a*+2*a*/*N**n*, we compute that
Therefore, setting *θ*_{k}=*π*/*a*(−*N*/2+*k*) and ϕ_{l}=*π*/*a*(−*N*/2+*l*), *k*,*l*=0,…,*N*−1, yields
Discretizing equation (6.18), it can be shown that the function *q*_{1}^{(r,s)}(*X*,*Y*) is approximated by the vector
where
The above expressions can be evaluated rapidly using a fast Fourier transform.

In our finite-difference scheme, we took *a*=15 and a step size of *h*=0.1.

(ii) In our second numerical scheme, we solve for *q*_{1} and *q*_{2} using a pseudo-spectral method (see Boyd (2001) for the one-dimensional version). Specifically, we employ the change of variables and , where 0<*θ*_{1},*θ*_{2}<*π*. In these new variables, we have
We look for a solution by approximating
Changing variables in equations (6.16) and (6.17) to *θ*_{1} and *θ*_{2}, putting in the above ansatz for and forcing equality at the points *θ*_{1}_{m}=(*m**π*)/(2*N*+1) and *θ*_{2}_{n}=(*n**π*)/(2*N*+1) for *m*,*n*=1,…,2*N*, we obtain the resulting linear system for the coefficients *a*_{mn} and *b*_{mn}. Approximating the derivatives in the forcing terms *F*_{j}(*q*_{0},…,*q*_{j−1}) is straightforward using this method.

We note that , where *S*_{n}(*X*) is a family of rational functions defined by the three-term recurrence relation
This relationship allows us to pass easily between and *q*_{j}(*X*,*Y*).

Choosing *N*=20 and *L*_{x}=*L*_{y}=2, we solved for *q*_{1} and *q*_{2}.

### (b) Computation of the speed versus amplitude relationship

Once we numerically obtain *q*, for a given choice of the parameters *α*_{1} and *α*_{2}, we have, from equation (6.6), a numerical approximation for *η*. We write *η* symbolically as
6.20
where *α*_{1} and *α*_{2} have been specified in order to solve for *q*. Furthermore, each choice of *α*_{1} and *α*_{2} corresponds to speeds *c*_{X} and *c*_{Y}, using equations (6.13) and (6.14) to within an error of order *O*(*ε*^{3}). Therefore, we can write the expansion (6.20) as
6.21

Now, suppose that, in our scaling of the non-dimensional variables, we take
where *a* is the characteristic height of the solitary wave and the primes denote dimensional quantities. Then, since *a**η*(*X*,*Y*)=*η*′(*X*′,*Y* ′),
where *η*, *X* and *Y* are non-dimensional. Furthermore, from equation (6.21), we conclude that the following holds to within an error of order *O*(*ε*^{3}):
By numerically solving for *ε* in the above equation, we obtain, for each choice of speed , the (non-dimensional) maximum depth *ε*. Since
when we pass back to dimensional variables, we obtain the dimensional speed versus maximum depth relationship.

Figure 9 shows the resulting speed/amplitude curve, which was generated by choosing, for *α*_{1}=*α*_{2}, 11 equally spaced points between 0.148 and 0.175 and connecting the resulting speed/amplitude pairs with an interpolating function (recall that the speed is determined once *α*_{1} and *α*_{2} are fixed). The dashed curve in figure 9 is the speed/amplitude curve generated using the finite-difference scheme; the solid curve in figure 9 is the speed/amplitude curve generated using the pseudo-spectral scheme. The relative difference between these two curves is about 2 per cent.

In figure 10, we plot the *X* and *Y* cross sections of the solution *η* corresponding to *α*_{1}=*α*_{2}=0.15, *ε*=0.56 and *c*^{2}−1=0.47; notice that *ε* was chosen so that the maximum depth corresponds to unity. Qualitatively similar *X* and *Y* cross sections for *η* are obtained in Parau *et al.* (2005) by direct numerical simulation of the full water-wave equations.

## 7. Asymptotics on the Benney–Luke equation

In this section, we verify the accuracy of the Poincare–Stokes method on the BL equation. The advantage of using the BL equation is that solitary-wave solutions of the BL equation can readily be obtained numerically (cf. Berger & Milewski 2000) and compared with solutions obtained from the asymptotic procedure.

Specifically, we consider
7.1
Now, let *X*=*x*−*c*_{X}*t* and *Y* =*y*−*c*_{Y}*t*. Then, equation (7.1) becomes, after taking ,
7.2
One can recover from equation (7.2) the Kadomtsev–Petviashvili (KP) equation; therefore, as in the asymptotic calculation given for full water waves, we assume that *q* decays like the lump solution to KP,
as . Using this in equation (7.2), we get four equations relating *c*_{X}, *c*_{Y}, *α*_{1} and *α*_{2}, any two of which are independent. Two such equations are
7.3
and
7.4
Using equations (7.3) and (7.4), *c*_{X} and *c*_{Y} are determined once *α*_{1} and *α*_{2} are specified.

We now proceed, as in the full water-wave case, to derive an asymptotic approximation to equation (7.2). To do so, assume that *c*_{X}∼*c*_{X1}+*ε**c*_{X2}+··· ,*c*_{Y}∼*c*_{Y1}+*ε**c*_{Y2}+··· and *q*∼*q*_{0}+*ε**q*_{1}+··· .

Using that *c*_{X}∼*c*_{X1}+*ε**c*_{X2}+··· and *c*_{Y}∼*c*_{Y1}+*ε**c*_{Y2}+··· in equations (7.3) and (7.4), we get that
7.5
We then use *q*∼*q*_{0}+*ε**q*_{1}+··· in equation (7.2) to get the following equations for *q*_{i}:
7.6
and
7.7

To compare how well our asymptotic procedure works, we proceed as follows. First, for a given *ε*, *α*_{1} and *α*_{2}, we solve equations (7.3) and (7.4) for the corresponding speeds *c*_{X} and *c*_{Y}. With these values, we use an iterative spectral scheme (see Ablowitz & Musslimani 2005) to numerically solve equation (7.2) for *q*.

With the same values for *ε*, *α*_{1} and *α*_{2}, we then solve equation (7.7) for *q*_{1} and *q*_{2}, where the speeds *c*_{X} and *c*_{Y} are computed from equations (7.5). Our numerical method is the same as that employed for full water waves.

Finally, to compare how well the Poincare–Stokes method compares with the exact solution, we examine the residuals
Figures 11–13 show the residuals *q*−*q*_{0}, *q*−(*q*_{0}+*ε**q*_{1}) and *q*−(*q*_{0}+*ε**q*_{1}+*ε*^{2}*q*_{2}), respectively, for *ε*=0.2,0.3,…,0.7. The solid lines in figures 11–13 are plots of the functions 0.741*ε*, 0.377*ε*^{2} and 0.533*ε*^{3}, where the coefficients in front of *ε*^{j} are chosen to minimize the least squares error.

## 8. Conclusion

In this paper, we examined gravity-capillary solitary water waves in two and three dimensions, using a recent NSP formulation of interfacial fluids with a free surface. Employing the Poincare–Stokes procedure, we developed a high-order asymptotic series for the free interface, where the first term in the series is the well-known solitary-wave solution of the KdV equation. From this series, we obtained how the energy and height of the maximally energetic wave depend on the surface tension. We also used this series to investigate two-dimensional solitary waves of depression and found that there are likely to be solitary waves that get close to the fluid bottom. We then adopted this asymptotic approach to three dimensions, where the first term in the corresponding asymptotic series is the rational solution to the KP-I equation; however, unlike the situation in two dimensions, the terms in the asymptotic series, beyond the first term, were computed numerically. The three-dimensional series (with three terms) was used to compute the associated speed versus amplitude curve. Finally, the Poincare–Stokes method was applied to the BL equation; the resulting asymptotic series (to three terms) was shown to agree favourably with the (numerically computed) solution of the BL equation.

## Acknowledgements

This work was partially supported by NSF grant DMS0604151. We thank H. Segur for his valuable input.

## Footnotes

- Received February 27, 2009.
- Accepted May 21, 2009.

- © 2009 The Royal Society