1. Introduction
In this paper we consider the two-dimensional formulation of a travelling gravity-capillary wave on a fluid of infinite depth. When posed in a travelling frame, the steady non-dimensionalised problem is to determine a velocity potential, 
$\phi (x, y)$, which is harmonic in a periodic domain, 
$-\tfrac {1}{2} \leqslant x \leqslant \tfrac {1}{2}$ and 
$-\infty < y \leqslant \eta (x)$. On the unknown free surface, 
$y = \eta (x)$, Bernoulli's condition requires that
\begin{equation} \frac{F^2}{2} |\boldsymbol{\nabla} \phi|^2 + \eta - B \kappa = \text{const.} \end{equation}
Here, the Froude number, 
$F$, characterises the balance between inertia and gravity, and the inverse Bond number, 
$B$, characterises the balance between gravity and surface tension; with this latter effect depending on the surface curvature, 
$\kappa$. These non-dimensional constants are given by
\begin{equation} F=c/\sqrt{gL_\lambda} \quad \text{and} \quad B=\sigma/\rho g L_\lambda^2, \end{equation}
where 
$c$ is the wave speed, 
$g$ is the gravitational constant, 
$L_{\lambda }$ is the wavelength, 
$\rho$ is the fluid density and 
$\sigma$ is the coefficient of surface tension.
Typically, a wave energy or amplitude parameter, 
$\mathcal {E}$, is fixed and prescribes the degree of nonlinearity. Solutions are then characterised by bifurcation curves in 
$(B,F)$ or 
$(B,F,\mathcal {E})$-solution space. The small surface tension limit corresponds to 
$B \to 0$.
Extensive results are known for the case with 
$B=0$ when surface tension is neglected, and this originates from the seminal work of Stokes (Reference Stokes1847); cf. the reviews by Okamoto & Shõji (Reference Okamoto and Shõji2001) and Toland (Reference Toland1996). Intuitively, we might expect that the inclusion of a small amount of surface tension results in a small change in the profile of the pure gravity wave. However, since the limit of 
$B \to 0$ is singularly perturbed, this is not necessarily the case, and it is known that the introduction of surface tension has a significant impact on the existence and uniqueness of solutions, their bifurcations and their profiles.
The goal of this paper is to present a numerical study of nonlinear solutions in the singular limit of 
$B \to 0$, for which we know one solution to be the Stokes wave. We demonstrate the numerical existence of a cohesive structure of branches of solutions existing under this limit. Importantly this suggests that, with fixed wave energy, 
$\mathcal {E}$, only one of a family of solutions approaches the classical Stokes wave as 
$B \to 0$.
We firstly discuss the analytical and numerical difficulties of the 
$B \to 0$ limit.
1.1. Longuet-Higgins and parasitic ripples
It is well-known observationally that under the action of both gravity and surface tension, ripples of small wavelength form on the forward face of a propagating wave. As shown by the experimental results of Cox (Reference Cox1958) and Ebuchi et al. (Reference Ebuchi, Kawamura and Toba1987) for instance, the amplitude of these parasitic capillary ripples increases when the overall amplitude of the wave (measured by crest to trough displacement) increases. An example of such parasitic ripples, as photographed by Ebuchi et al. (Reference Ebuchi, Kawamura and Toba1987), is shown in figure 1, where it is seen that these ripples are asymmetric about the wave crest and unsteady in the frame of the propagating wave. Note that in this paper we shall refer to solutions exhibiting parasitic capillary ripples as those where a short wavelength and small amplitude wave is present on what appears to be a gravity-dominated water wave.
Figure 1. Experimental picture showing parasitic ripples located near the crests of steep gravity waves from Ebuchi, Kawamura & Toba (Reference Ebuchi, Kawamura and Toba1987) (reproduced with permission).
A seminal advance in the modelling of these parasitic capillary waves arose from the methodology of Longuet-Higgins (Reference Longuet-Higgins1963), who predicted that for small surface tension these parasitic ripples would be exponentially small in both amplitude and wavelength. In this simplest steady framework one assumes that the parasitic ripples are fixed to the same travelling frame of reference as the underlying gravity wave. However, it was noted by Perlin, Lin & Ting (Reference Perlin, Lin and Ting1993) that these analytical predictions produced poor agreement with both experimental wave profiles and numerical solutions of the steady nonlinear equations, a result of asymptotic inconsistencies in his method. We shall provide a preliminary discussion of these issues in § 7, in which we comment that the Longuet-Higgins approximation can be improved through the use of exponential asymptotics.
Nevertheless, it remains an open question as to whether parasitic capillary ripples similar to those shown in figure 1 may be found as either symmetric or asymmetric solutions of the steady framework of (1.1). In this work we present clear numerical evidence that steady symmetric parasitic ripples do exist within the solution space of the classical potential framework in the 
$B \to 0$ limit.
1.2. Schwartz & Vanden-Broeck and the complexity of 
$(B,F)$-space
In their seminal work Schwartz & Vanden-Broeck (Reference Schwartz and Vanden-Broeck1979) developed a numerical scheme using a series truncation method to compute periodic gravity-capillary waves of the exact nonlinear equations. Imposing symmetry at 
$x=0$ and an amplitude condition on the crest-to-trough displacement, they presented a preliminary classification of solutions in 
$(B,F)$-space of types 1, 2, 3 and 4. Each type number was associated with a distinct branch of solutions, and corresponded to the number of observed ‘dimples’ or inflexion points on a (half) wave profile.
A reproduction of their original bifurcation diagram, which is computed at fixed crest-to-trough amplitude, is shown in figure 2(b). Our intention in reproducing this figure is to convince the reader that indeed the bifurcation space of the gravity-capillary problem is certainly non-trivial, and it is difficult to observe any clear structure. We also show the computed (Schwartz & Vanden-Broeck Reference Schwartz and Vanden-Broeck1979) bifurcation curves in figure 2 alongside our solutions of fixed energy.
Figure 2. Our solutions of fixed energy from § 5 are shown in (a). Panel (b) shows a reproduction of the fixed amplitude results previously published in Schwartz & Vanden-Broeck (Reference Schwartz and Vanden-Broeck1979) (reproduced with permission). These branches of fixed amplitude are shown in the 
$(\kappa ,\mu )$-plane, where the boxed type number indicates the number of observed ‘dimples’ or inflexion points on a (half) wave profile. In (c) we show the type-11 fixed amplitude branch.
One of their solutions, Schwartz & Vanden-Broeck (Reference Schwartz and Vanden-Broeck1979, figure 10), is of particular interest in the context of parasitic ripples. This profile, similar to that shown in figure 3, appears to contain small-scale capillary ripples as a perturbation to the main Stokes wave. This is one of the types of solution that we will be expanding upon in this work. Note in addition that the type 1 to 4 branches, as shown in their figure, have a non-trivial and unstructured shape in the bifurcation diagram; it is not obvious if a more consistent pattern emerges upon increasing the type number, or whether these solution curves can be taken as 
$B \to 0$. We shall explain the reason for these issues in this work.
Later, in seeking to compare new experimental data with the previous analytical approximations of Longuet-Higgins (Reference Longuet-Higgins1963) and numerical solutions of Schwartz & Vanden-Broeck (Reference Schwartz and Vanden-Broeck1979), Perlin et al. (Reference Perlin, Lin and Ting1993) made extensive remarks on the challenges of navigating the solution space of the full nonlinear problem, noting that ‘there is no known method for determining the number of solutions to the numerical formulation…’ (p. 618). Indeed, they state that (p. 598),
Surprisingly little information is available on these waves of disparate scales, presumably due to the analytical/numerical, as well as experimental, difficulties involved. Perlin et al. (Reference Perlin, Lin and Ting1993)
1.3. Chen & Saffman and the impossibility of the $B \to 0$ limit
Nearly in parallel with the work by Schwartz & Vanden-Broeck (Reference Schwartz and Vanden-Broeck1979), Chen & Saffman (Reference Chen and Saffman1979, Reference Chen and Saffman1980a,Reference Chen and Saffmanb) produced a series of works where they examined the Stokes wave problem, largely from the perspective of weakly nonlinear theory (and its numerical consequences on the full nonlinear problem).
In Chen & Saffman (Reference Chen and Saffman1979) they considered weakly nonlinear solutions of (1.1) in powers of a small wave amplitude, 
$\epsilon$. Expressing the solution, 
$y=\eta (x)$, as a Fourier series, this permits analytical solutions for the Fourier coefficients, 
$A_n$. They discovered that in fixing the point of symmetry of the wave profile to be at 
$x=0$, the branches of solutions in the 
$(\kappa ,A_n)$-bifurcation space (where 
$\kappa =4{\rm \pi} ^2 B$) are discontinuous either side of the point 
$\kappa =1/n$. Due to this discontinuity and the analytical criterion, they commented (p. 204):
The gravity wave (
$\kappa =0$) is therefore a singular limit which cannot be reached smoothly by applying the limit 
$\kappa \to 0$ to a gravity capillary wave. Chen & Saffman (Reference Chen and Saffman1979)
In our work we will note how this statement is misleading since their non-smoothness is a consequence of their initial assumption of a fixed point of symmetry at 
$x=0$. We demonstrate this in § 6.4, noting that this is due to the presence of a symmetry shifting bifurcation. Thus, if their enforced symmetry at 
$x=0$ were to be relaxed, the branches of solutions in 
$(\kappa , A_n)$-space either side of the point 
$\kappa =1/n$ would be continuous.
In a second work by Chen & Saffman (Reference Chen and Saffman1980b), the numerical solution space was explored for waves of finite amplitude. Their choice of amplitude was a linear combination of Fourier coefficients, typically chosen to be that of the fundamental mode, 
$A_1$, or the 
$n$th mode 
$A_n$. Nonlinear solutions were computed. However, the resultant branches of solutions in their bifurcation diagram did not connect, from which they concluded:
These results confirm the impossibility of going continuously from a pure capillary- gravity wave to a gravity wave by letting 
$\kappa \to 0$. Chen & Saffman (Reference Chen and Saffman1980b)
We later note in § 6.2 that if the wave energy instead is fixed as an amplitude parameter then the continuous set of solutions as 
$B \to 0$, discovered within this paper, bifurcate from solutions with fundamental mode 
$A_1=0$. Hence, if 
$A_1$ is fixed to be a non-zero constant, as in the numerical work of Chen & Saffman (Reference Chen and Saffman1980b), this bifurcation point would remain undiscovered. We shall conclude that in order to achieve a continuous transition as 
$B \to 0$, the first Fourier coefficient should not be fixed.
1.4. Outline of the paper
In this work we shall consider the numerical behaviour of steady symmetric parasitic ripples for small values of the Bond number, 
$B$. Starting in § 2, we introduce the governing equations for the gravity-capillary wave problem, which we transform to depend on the velocity potential, 
$\phi$, alone. In § 3 the well-known linear solutions are derived. These form a starting point for our numerical method of § 4. Solutions are presented in § 5, which we use to demonstrate that as 
$B \to 0$ for fixed energy, 
$\mathcal {E}$, this bifurcation structure appears to form a countably infinite number of connecting branches of solutions in the 
$(B,F)$-plane. Each branch forms a ‘finger’ in the solution space, which is connected continuously to the proceeding branch. These branches then accumulate in the limit of 
$B \to 0$ such that solutions are conjectured to exist in an 
$O(1)$ interval in the Froude number, 
$F$. These branches of solutions are connected at the point where they bifurcate from a wave with smaller fundamental wavelength, resulting in numerical evidence for the 
$B \to 0$ limit of the steep gravity-capillary wave problem having a continuous set of solutions.
This allows us to show why previous authors have failed to reveal this underlying structure, which we comment upon in §§ 6 and 6.5. Lastly, in § 7 we comment upon the asymptotic properties as 
$B \to 0$ of our discovered solutions, and how this uncovered structure is likely to be present in other numerical problems for small values of the Bond number, 
$B$.
2. Mathematical formulation
Consider a two-dimensional free surface flow of an inviscid, irrotational and incompressible fluid of infinite depth. The velocity potential, 
$\phi$, is defined by 
$\boldsymbol {u} = \boldsymbol {\nabla } \phi$. We assume that, in the lab frame, the free surface, 
$y = \eta (x, t)$, is periodic in 
$x$ with wavelength 
$L_\lambda$, and moves to the right with wave speed 
$c$. We non-dimensionalise with unit length, 
$L_\lambda$, and velocity, 
$c$. We consider steady travelling-wave solutions by introducing a subflow of unit horizontal velocity in the opposite direction of wave propagation. This negates the movement of the free surface. Then 
$\eta _t = 0 = \phi _t$ yields the steady governing equations (compare with Vanden-Broeck Reference Vanden-Broeck2010, (2.48)–(2.55) for example),
$$\begin{gather} \phi_{xx}+\phi_{yy}=0 \quad \text{for} \ y \leqslant \eta, \end{gather}$$
$$\begin{gather}\phi_{y}=\eta_{x}\phi_{x} \quad \text{at} \ y=\eta, \end{gather}$$
$$\begin{gather}\frac{F^2}{2}(\phi^{2}_{x}+\phi^{2}_{y})+y-B\frac{\eta_{xx}}{(1+\eta_x^2)^{{3/2}}}=0 \quad \text{at} \ y=\eta, \end{gather}$$
$$\begin{gather}\phi_{y} \to 0 \quad \text{and} \quad \phi_{x} \to -1 \quad \text{as} \ y \to -\infty, \end{gather}$$
for the travelling wave now in 
$x \in [-\tfrac {1}{2}, \tfrac {1}{2})$. Thus, the system is governed by Laplace's equation (2.1a) within the fluid, kinematic and dynamic conditions on the free surface (2.1b) and (2.1c), and uniform flow conditions in the deep-water limit (2.1d). The horizontal velocity condition (2.1d), our subflow, indicates a uniform flow moving towards the left. The spatial subscripts in (2.1) correspond to partial differentiation.
Remark on terminology: note that in the mathematical formulation above, we have non-dimensionalised lengths by a fixed physical wavelength, 
$L_{\lambda }$, and, hence, we shall seek solutions that are 
$1$-periodic in the non-dimensional travelling frame. However, these solutions may have a smaller wavelength which is less than unity. We thus define 
$\lambda$ to be the non-dimensional fundamental wavelength (the smallest such wavelength). Moreover, in this work we shall refer to a wave with fundamental wavenumber 
$k = 1/\lambda$ as a pure 
$k$-wave. Thus, a pure 
$k$-wave has a dimensional wavelength of 
$\lambda L_\lambda$.
2.1. The conformal mapping to the $(\phi ,\psi )$-plane
We now formulate the governing equations (2.1) in the potential 
$(\phi ,\psi )$-plane, as shown in figure 4. We assume that the free surface is located along 
$\psi = 0$, and introduce the notation of 
$X$ and 
$Y$ for the fluid quantities evaluated on the free surface. Thus,
\begin{equation} X(\phi)\equiv x(\phi,0) \quad \text{and} \quad Y(\phi)\equiv \eta(x(\phi,0)). \end{equation}We may now obtain expressions for the surface derivative and curvature by differentiating (2.2a,b). This yields
\begin{equation} \eta_{x}=\frac{Y_{\phi}}{X_{\phi}} \quad \text{and} \quad \eta_{x x}=\frac{X_{\phi}Y_{\phi \phi}-Y_{\phi}X_{\phi \phi}}{X_{\phi}^3}. \end{equation}
We now seek to rewrite the kinematic (2.1b) and dynamic (2.1c) boundary conditions on the surface in terms of the conformal variables 
$X$ and 
$Y$. First, the velocities, 
$\phi _x$ and 
$\phi _y$, may be inverted, yielding
\begin{equation} \phi_x = \frac{x_{\phi}}{x_{\phi}^2+y_{\phi}^2} \quad \text{and} \quad \phi_y = \frac{y_{\phi}}{x_{\phi}^2+y_{\phi}^2}. \end{equation}
Figure 4. The conformal mapping from 
$(x,y)$ to 
$(\phi ,\psi )$ is shown.
Finally, substitution of (2.4a,b) for 
$\phi _x$ and 
$\phi _y$, (2.3a,b) for 
$\eta _x$ and 
$\eta _{xx}$, and 
$Y(\phi )=\eta (x(\phi ,0))$ into Bernoulli's equation (2.1c) and setting 
$\psi =0$ yields our governing equation (compare with (6.12) of Vanden-Broeck Reference Vanden-Broeck2010)
\begin{equation} \frac{F^2}{2J}+Y+B\frac{(Y_{\phi}X_{\phi \phi}-X_{\phi}Y_{\phi \phi})}{J^{{3/2}}}=0. \end{equation}
Above we have introduced the surface Jacobian, 
$J$, via
\begin{equation} J(\phi) = X_{\phi}^2+Y_{\phi}^2. \end{equation}
Note that in the conformal formulation, the kinematic condition (2.1b) can be verified to be satisfied identically once (2.3a,b) and (2.4a,b) are used on the streamline 
$\psi = 0$.
In addition to Bernoulli's equation (2.5), in order to close the system, we require a harmonic relationship between 
$X$ and 
$Y$. Note that within the fluid, 
$y(\phi , \psi )$ can be written as a Fourier series of the form
\begin{equation} y(\phi,\psi)=\psi+A_0+\sum_{n=1}^{\infty}\textrm{e}^{2n{\rm \pi} \psi}\left[A_n\cos{(2n{\rm \pi} \phi)}+B_n\sin{(2n{\rm \pi} \phi)} \right], \end{equation}
where 
$A_n$ and 
$B_n$ are real-valued for all 
$n$. Indeed, the above ansatz satisfies 
$y_{\phi \phi } + y_{\psi \psi } = 0$ along with the depth condition 
$y \sim \psi$ as 
$\psi \to -\infty$.
We define the Hilbert transform on 
$Y$ by
\begin{equation} \mathscr{H}[Y](\phi^{\prime})=\int\kern-8pt -_{-\infty}^{\infty} \frac{Y(\phi)}{\phi-\phi^{'}} \, \textrm{{d}}{\phi}, \end{equation}where the integral is of principal-value type. Then by the assumed periodicity of the solution, this implies that
\begin{equation} \mathscr{H}[Y](\phi^{\prime}) = \int\kern-8pt -_{{-}1/2}^{1/2} Y(\phi)\cot{[{\rm \pi}(\phi-\phi^{\prime})]} \,\textrm{{d}}{\phi}. \end{equation} We can then verify that the individual Fourier modes can be related using 
$\mathscr {H}[\sin (2n{\rm \pi} \phi )]=\cos (2n {\rm \pi}\phi )$ and 
$\mathscr {H}[\cos (2n {\rm \pi}\phi )]=-\sin (2n {\rm \pi}\phi )$. From using the Cauchy–Riemann relations of 
$x_{\phi }=y_{\psi }$ and 
$x_{\psi }=-y_{\phi }$, we obtain the harmonic relationships between 
$X$ and 
$Y$ on the free surface via
\begin{equation} X_{\phi}(\phi)=1-\mathscr{H}[Y_{\phi}(\phi)] \quad \text{and} \quad Y_{\phi}(\phi) = \mathscr{H}[X_{\phi}(\phi) - 1]. \end{equation}
A choice of any one of the relations in (2.10a,b), combined with Bernoulli's equation (2.5) allows 
$X$ and 
$Y$ to be solved.
2.2. The energy constraint
In order to fully close the formulation, we shall impose an energy constraint on the solution, which can be viewed as equivalent to a measurement of the wave amplitude. We define the wave energy, 
$E$, to be
\begin{equation} E=\frac{F^2}{2}\int_{-{1/2}}^{{1/2}}Y(X_{\phi}-1) \,\textrm{{d}} \phi +B \int_{-{1/2}}^{{1/2}}(\sqrt{J}-X_{\phi})\,\textrm{{d}} \phi +\frac{1}{2}\int_{-{1/2}}^{{1/2}}Y^2 X_{\phi} \,\textrm{{d}} \phi, \end{equation}where the first integral on the right-hand side corresponds to the kinetic energy, the second to the capillary potential energy and the third to the gravitational potential energy. The derivation of (2.11) from the bulk energy is given in Appendix A.
For comparison purposes, it will be convenient for us to rescale the energy in (2.11) by the energy of the highest (fundamental) gravity wave, 
$E_{hw}$. Thus, we write
\begin{equation} \mathcal{E} = \frac{E}{E_{{hw}}}, \end{equation}
where 
$E_{hw} \approx 0.00184$ (to 5 decimal places) is calculated using the numerical scheme of § 4 applied to the pure gravity wave using 
$n=4096$ Fourier coefficients.
The choice of how to define an amplitude or energy condition for the wave is a subtle one. In this paper we shall comment on the following three choices of amplitude:
\begin{equation}
\mathscr{A}=\begin{cases}
\mathcal{E} & \text{[energy definition from (2.11)]},\\
A_1 & \text{[first Fourier coefficient from (2.7)]}, \\
Y(0)-Y(1/2) & \text{[crest-to-trough displacement].}
\end{cases}
\end{equation}
The second choice of 
$A_{1}$, as used in Chen & Saffman (Reference Chen and Saffman1980b), designates the amplitude to be the first Fourier coefficient, while the third choice of 
$Y(0)-Y(1/2)$, as used by Schwartz & Vanden-Broeck (Reference Schwartz and Vanden-Broeck1979), is a sensible choice to measure the physical wave height of the fundamental Stokes wave.
Note that both definitions of amplitude, 
$A_{\text {1}}$ and 
$Y(0)-Y(1/2)$, have the problem that strongly nonlinear waves (as measured by a lack of decay in the Fourier coefficients) can occur, even at small amplitude values. This is particularly affected by the fact that gravity-capillary waves may take a variety of shapes beyond the simple fundamental wave considered by Stokes (Reference Stokes1847). Similar difficulties were encountered by Chen & Saffman (Reference Chen and Saffman1979, Reference Chen and Saffman1980b), who chose 
$\mathscr {A} = A_1$ but commented that:
We found from experience that none of these parameters were universally useful for describing the bifurcation phenomenon to be described in this work, and in fact we have been unable to construct a parameter which characterized the magnitude of the wave for all the phenomena in a satisfactory way. Chen & Saffman (Reference Chen and Saffman1979)
It may be that using the energetic definition of amplitude with 
$\mathscr {A}=\mathcal {E}$ is the modification required to fix these issues; indeed within the context of our numerical investigations this does seem to be the case in the small surface tension limit.
3. Linear theory, Wilton ripples and type $(n,m)$-waves
It will be useful for us to review linear solutions in the notation of § 2.1. The results of linear theory are found from the first two terms of a Stokes expansion in powers of a small amplitude parameter, 
$\epsilon$ (see, e.g. Vanden-Broeck Reference Vanden-Broeck2010, § 2.4.2). Thus, we shall consider equations (4.1a) and (4.1b) and take 
$X\sim X_0+ \epsilon X_1$ and 
$Y \sim Y_0 + \epsilon Y_1$. Solving the resultant equations yields 
$X_0=\phi$ and 
$Y_0=0$ at leading order. At 
$O(\epsilon )$, we write 
$X_1$ and 
$Y_1$ as Fourier series and assume that the two solutions are respectively odd and even about 
$\phi =0$. This yields the necessary equation that
\begin{equation} \sum_{k=1}^{\infty}\left[ F^2(2k{\rm \pi})-1-(2k{\rm \pi})^2B \right] a_k \cos{(2k{\rm \pi} \phi)}=0. \end{equation}In order to obtain non-trivial solutions, we require the linear dispersion relation of
\begin{equation} 2k{\rm \pi} F^2-1-4k^2{\rm \pi}^2B=0, \end{equation}
and obtain 
$X_1 = a_k \sin (2k {\rm \pi}\phi )$ and 
$Y_1=a_k \cos {(2k {\rm \pi}\phi )}$. Thus, the linear solution, a pure-
$k$ wave, is approximated by
\begin{equation} X \sim \phi + \epsilon\left[ a_k \sin(2k {\rm \pi}\phi) \right] \quad \text{and} \quad Y \sim \epsilon \left[ a_k \cos(2k {\rm \pi}\phi) \right], \end{equation}to the first two orders. Substitution into the energy expression (2.12) yields
\begin{equation} \mathcal{E} \sim \epsilon^2 \frac{ 2 K^2 {\rm \pi}^2 B a_k^2}{E_{{hw}}}. \end{equation} The linear solution (3.3a,b) was assumed to satisfy the single dispersion relation (3.2) for the 
$k$th Fourier mode only. Note that other solutions may be constructed that satisfy the dispersion relation for more than one mode. For instance, if the modes with 
$k=1$ and 
$k=n$ are assumed to be non-degenerate, then we require that both 
$2{\rm \pi} F^2-1-4{\rm \pi} ^2B=0$ and 
$2n{\rm \pi} F^2-1-4n^2{\rm \pi} ^2B=0$. This yields the so-called Wilton ripples predicted by Wilton (Reference Wilton1915), located wherever
\begin{equation} B_\text{wilton} = \frac{1}{4{\rm \pi}^2 n} \quad \text{and} \quad F_\text{wilton}^2 = \frac{(1+n)}{2{\rm \pi} n}, \end{equation}
with 
$n \in \mathbb {Z}^+$. The Wilton ripples are then given by
\begin{equation} X_1 = a_1 \sin(2 {\rm \pi}\phi)+a_n \sin(2n {\rm \pi}\phi) \quad \text{and} \quad Y_1 = a_1 \cos(2 {\rm \pi}\phi) +a_n \cos(2n {\rm \pi}\phi) . \end{equation}
In the numerics of § 4, we shall often initialise the numerical continuation method with the linear solution (3.3a,b) using a small value of 
$\epsilon a_k$. Crucially, since this linear solution is invalid near points (3.5a,b), we must ensure that our initial choice lies away from the critical numbers of 
$B_\text {wilton}$ and/or 
$F_\text {wilton}$.
As introduced by Chen & Saffman (Reference Chen and Saffman1979), linear solutions that consist of a combination of pure 
$n$- and 
$m$-waves, and with fundamental wavelengths of 
$\lambda =1/n$ and 
$1/m$, respectively, are described as a type 
$({n},{m})$-wave. Thus, under this terminology, Wilton's solution in (3.5a,b) is an example of a type 
$(1,n)$-wave. Our numerical results presented in § 5 will contain solutions that are the nonlinear analogue of a type 
$(1,n)$-wave.
4. The numerical method
In this section we describe the numerical procedure for solving Bernoulli's equation (2.5) and the harmonic relationship (2.10a,b) for 
$X(\phi )$ and 
$Y(\phi )$ subject to a given value of the energy, 
$\mathcal {E}$, from (2.11). Thus,
$$\begin{gather} \frac{F^2}{2J}+Y+B\frac{(Y_{\phi}X_{\phi \phi}-X_{\phi}Y_{\phi \phi})}{J^{{3/2}}}=0, \end{gather}$$
$$\begin{gather}X_{\phi}(\phi)=1-\mathscr{H}[Y_{\phi}(\phi)], \end{gather}$$
$$\begin{gather}\mathcal{E} = \frac{F^2}{2 E_{{hw}}}\int_{-{1/2}}^{{1/2}}Y(X_{\phi}-1) \,\textrm{{d}} \phi + \frac{B}{E_{{hw}}} \int_{-{1/2}}^{{1/2}}(\sqrt{J}-X_{\phi}) \,\textrm{{d}} \phi + \frac{1}{2 E_{{hw}}}\int_{-{1/2}}^{{1/2}}Y^2 X_{\phi} \,\textrm{{d}} \phi. \end{gather}$$
Recall we define 
$E_{hw}$ as the energy of the highest Stokes wave (
$E_{hw} \approx 0.00184$).
Solutions of the above problem are regarded as lying within 
$(B, F, \mathcal {E})$-space. We solve these equations using Newton iteration on a truncated Fourier series using the fast Fourier transform. The procedure is as follows.
(i) An initial guess for
$Y(\phi )$ is carefully chosen using either linear theory (3.3a,b) or from a previously computed solution (cf. § 5 for specific details).(ii) Of the triplet
$(B, F,\mathcal {E})$, we choose to fix two parameters and treat the last parameter as an unknown.(iii) The collocation variable,
$\phi$, is discretised using $N$ grid points, with $\phi _k = -1/2 + k \Delta \phi$ for $k = 0, \ldots , N-1$ and $\Delta \phi = 1/N$. We define the solution $Y(\phi _k) = Y_k$ at each of these points. Note that once $Y_k$ is known for all $k$, $X_k$ can be calculated using the harmonic relation (4.1b).(iv) Combined with the unknown parameter (either
$B$, $F$ or $\mathcal {E}$), this yields $N+1$ unknowns. Bernoulli's equation (4.1a), evaluated at $\phi _k$, provides $N$ equations and the system is closed with the additional energy constraint (4.1c). Newton iteration is then used to solve the nonlinear system of equations until a certain tolerance (typically $10^{-11}$) on the norm of the residual is met.
In our numerical scheme we leverage the Fourier transform for efficient manipulation of the solutions. In particular, note that the Hilbert transform, 
$\mathscr {H}$, needed for the harmonic relation (4.1b), can be evaluated via 
$\mathscr {H}[Y]=\mathscr {F}^{-1}[\mathrm {i} \cdot \text {sgn}(k)\mathscr {F}[Y]]$, where 
$\mathscr {F}$ denotes the Fourier transform and sgn is the signum function. Both the Fourier and inverse-Fourier transforms are calculated with the fast Fourier transform algorithm. The derivatives of 
$Y$ are also computed in Fourier space using the relationship 
$Y^{(n)}(\phi )=\mathscr {F}^{-1}[(2 {\rm \pi}\mathrm {i} k)^n \mathscr {F}[Y]]$. In order to obtain the numerical results presented in § 5, we find it sufficient to use 
$N = 1024$ mesh points. The computations are performed using a desktop computer and individual solutions are typically computed in under a second.
In essence, our goal will be to study the 
$(B, F, \mathcal {E})$-solution space, particularly as 
$B \to 0$. We start from a low-energy solution, and increase the parameter 
$\mathcal {E}$ until the desired value is reached. In order to initialise this continuation procedure at small values of 
$\mathcal {E}$, we select an initial Bond number which is chosen away from the Wilton ripples value of 
$B_\text {wilton}$ in (3.5a,b). Then the Froude number is approximated by the linear dispersion relation (3.2) with 
$k = 1$, and we use the linear approximations of 
$X$ and 
$Y$ from (3.3a,b) with a small arbitrary choice of 
$\epsilon a_k$ (typically 
$10^{-5}$). For this linear solution, 
$\mathcal {E}$ is then calculated; the above serves as the initialisation procedure for the Newton scheme which solves for values of 
$Y_i$ and 
$F$.
Once solutions are found at desired values of 
$\mathcal {E}$, we establish the 
$(B, F)$-bifurcation space by continuation from a previously calculated solution. Note that in some cases, it will be necessary to fix 
$B$ or 
$F$ and solve for the other value, depending on the gradient of the bifurcation curves.
5. Numerical results for fixed energy, $\mathcal {E}$
The numerical results we now present suggest that at a fixed value of 
$\mathcal {E}$, certain solutions in the 
$(B, F)$-bifurcation space can be classified according to ‘finger’-type structures and ‘side-branch’-type structures. An example of this structure, as drawn in the 
$(B,F)$-plane, is shown in figure 5.
Figure 5. A typical component of the bifurcation diagram illustrated in 
$(B, F)$-space consisting of a single finger, 
$G_{n \to n+1}$, (shown bolded) and two side curves 
$S_n$ and 
$S_{n+1}$. This is considered at a fixed value of the energy, 
$\mathcal {E}$.
First, let us first define the side branch, 
$S_n$, as
\begin{equation} S_n = \{
\text{Bifurcation curve of solutions analogous to type (0, n)-waves}\}.
\end{equation}
Thus, 
$S_n$ corresponds to those points in 
$(B,F)$-space associated with a certain type of solution. These solutions are pure 
$n$-waves (
$1/n$-periodic solutions in the interval); they are the nonlinear analogue of the linear type 
$(0, n)$-waves introduced in § 3, i.e. a sine or cosine wave with wavenumber 
$n$ about a constant mean value.
In addition, adjacent side branches are connected by fingers, say 
$G_{n \to n+1}$. We define such a structure as
\begin{equation} G_{n \to n+1} = \{\text{Bifurcation curve of solutions connecting}\ S_{n}\ {\rm to}\ S_{n+1}\}. \end{equation}
The finger can be interpreted as follows. Along 
$S_n$, solutions are pure 
$(n)$-waves; following this set of solutions, there exists a bifurcation point where the 
$1$-mode grows. Following this new branch, which is labelled 
$G_{n \to n+1}$, yields a solution analogous to a type 
$(1, n)$-wave. Continuing along 
$G_{n \to n+1}$, the solution transitions to type 
$(1, n+1)$ and then finally to a pure-
$(n+1)$ wave where it connects to 
$S_{n+1}$. An illustration of these classifications is shown in figure 5.
In the following sections we present solutions along the side branches, 
$S_n$, and fingers, 
$G_{n \to n+1}$, for waves that are approximately half the height of the highest fundamental gravity wave. For our choice of energy in (4.1c), this occurs at 
$\mathcal {E}=0.3804$. Starting in § 5.1, we describe the structure of solutions across a prototypical finger, 
$G_{13 \to 14}$, and then in § 5.2 we demonstrate how this finger bifurcates from side branches 
$S_{13}$ and 
$S_{14}$. Multiple fingers are then shown in § 5.3 for 
$n=7$ to 
$28$, demonstrating their behaviour as 
$B \to 0$.
5.1. Analysis of a single finger, $G_{n \to n+1}$
The prototypical finger 
$G_{13 \to 14}$ is shown in figure 6 for a value of 
$\mathcal {E} = 0.3804$. Note that solutions near the ‘tip’ of the finger seem to correspond to the phenomena of parasitic ripples discussed in § 1; that is, we observe a series of small-scale capillary-dominated ripples riding on the surface of a steep gravity wave. This is shown in insets 
$(c)$, 
$(d)$ and 
$(e)$ in figure 6. Below, we will continue to refer to solutions as being separated into capillary ripples and an underlying gravity wave, even though this classification may be ambiguous.
Figure 6. A single finger of solutions, 
$G_{13 \to 14}$, is displayed in the 
$(B,F)$-plane for a value of 
$\mathcal {E}=0.3804$. The dashed profiles in 
$(a)$ and 
$(g)$ are solutions from branches 
$S_{13}$ and 
$S_{14}$ near the bifurcation point. The (b)–(g) axis limits are the same as inset (a).
As we move down either side of 
$G_{13 \to 14}$ by decreasing the Froude number, the amplitude of the ripples increases while the amplitude of the underlying gravity wave decreases. This is shown in figure 6 via the transitions 
$(c)\to (b)\to (a)$ and 
$(e)\to (\,f)\to (g)$. It becomes extremely challenging to numerically compute solutions below 
$(a)$ and 
$(g)$.
Finally, as we travel from right to left across the finger, the wavelength of the ripples decreases as an extra ripple is formed. This can be seen by comparing solutions in insets 
$(g)$ and 
$(a)$, where 
$(g)$ has 13 maxima and 
$(a)$ has 14 maxima. The increase in the number of ripples can be observed as occurring near the tip of the finger between insets 
$(c)$ and 
$(e)$. We will discuss the structure of this process in § 7.
5.2. Analysis of side branches $S_n$ and $S_{n+1}$
We now discuss the side branches. In the case of the prototypical finger, 
$G_{13 \to 14}$, displayed in figure 6, we observe that this finger connects two side branches 
$S_{13}$ and 
$S_{14}$, as shown in figure 7. The branch 
$S_{13}$ contains pure 
$(13)$-waves, which have a fundamental wavelength of 
$\lambda =1/{13}$. The branch 
$S_{14}$ consists of pure 
$14$-solutions, which have 
$\lambda =1/{14}$.
Figure 7. The finger 
$G_{13 \to 14}$ is shown against the two side branches 
$S_{13}$ and 
$S_{14}$. The two side branches terminate at points 
$a$ and 
$c$ (black circle) through the trapping of bubbles. Circles represent the locations where solutions of 
$S_{n}$ for 
$n \geqslant 15$, shown for every fifth point, become self-intersecting, found from the numerical predictions of Appendix B. The (b)–(d) axis limits are the same as inset (c).
We next observe that at fixed energy, 
$\mathcal {E}=0.3804$, the solutions in 
$S_{13}$ and 
$S_{14}$ reach a limiting configuration through the trapping of bubbles, shown by solutions 
$(a)$ and 
$(c)$. These branches of solutions are the large amplitude analogue of those predicted by linear theory in § 3, given by
\begin{equation} F^2 \sim {1/(2n{\rm \pi})+(2n{\rm \pi})B} \quad \text{and} \quad F^2 \sim {1/(2{\rm \pi}(n+1))+2{\rm \pi}(n+1)B}. \end{equation}
These were obtained by taking the values of 
$k=n$ and 
$k=n+1$ in the linear dispersion relation (3.2).
In order to compute these branches numerically, an initial pure–
$n$ solution was taken from linear theory with the dispersion relation (3.2) satisfied for 
$k=n$. This gives a cosine profile with 
$n$ peaks across the periodic domain. Slowly increasing the energy of this solution across multiple runs yields a single solution for each branch at 
$\mathcal {E}=0.3804$, from which these branches were calculated by continuation at fixed 
$\mathcal {E}$.
The location along the branch for which solutions reach a limiting configuration through a trapped bubble can be numerically predicted by the results of Appendix B. These points are shown in figure 7 for 
$n \geqslant 15$.
We see that as the value of 
$n$ for these limiting solutions increases, the value of 
$F$ at these points increases beyond that of the original finger. Thus, below a certain value of the Bond number, we expect that each finger will instead bifurcate from self-intersecting solutions. As we then proceed to increase the Froude number and transverse the side of each of these fingers for 
$B< B_{{crit}}$, we anticipate that the solutions will turn physical. This would result in the tip of each finger consisting of purely physical solutions.
5.3. The unveiled structure for $B \to 0$
This process of generating an individual finger may be repeated across different values of the Bond number, resulting in a remarkable structure that holds in the limit of 
${B \to 0}$. Many of these fingers are shown in figure 8 from 
$n=7$ to 
$n=28$; for clarity, the side branches have been omitted from this figure. As the Bond number decreases over each finger, the wavelength of the ripples decreases from 
$1/n$ to 
$1/(n+1)$, resulting in the formation of an additional crest. Consecutive fingers are connected at the point from which they bifurcate from the side branches of pure 
$n$-waves, demonstrated previously in § 5.2 and shown by solutions 
$(d1)$ and 
$(d2)$ in figure 8. The solutions at these bifurcation points display a phase shift of 
$1/n$ between them. Due to this phase shift, the 
$n$th Fourier coefficient changes sign between these solutions at this bifurcation point. It is this phase shift that led Chen & Saffman (Reference Chen and Saffman1979) to misleadingly state (on p. 204) that the weakly nonlinear solutions are discontinuous with respect to the 
$n$th Fourier coefficient at this point.
Figure 8. Multiple fingers 
$G_{7 \to 8}$ to 
$G_{28 \to 29}$ are shown to form one connected branch in the 
$(B,F)$-plane at fixed 
$\mathcal {E}=0.3804$. The (b)–(c) axes limits are the same as inset (a1) and the (e)–(d) limits are the same as inset (f).
From solutions 
$(a1)$, 
$(b)$ and 
$(c)$, labelled at the top of the fingers in figure 8, we observe that as 
$B \to 0$, the amplitude of the ripples decreases and the overall solution appears to tend towards the fundamental Stokes wave with energy 
$\mathcal {E}$. Although the profile in figure 8(a1) seems to indicate a pure gravity wave, the capillary ripples can be detected under closer inspection. In order to quantify this, we isolate the pure gravity wave solution, 
$y_0$, and plot 
$y - y_0$ in figure 8(a2). This shows that the ripples are still present in the solution, but with a very small amplitude. Moreover, one can verify that the profile norm, 
$\lvert y-y_0 \rvert$, is of 
$O(B)$ by repeating this procedure for multiple solutions along the top of the fingers in figure 8. We shall comment on this algebraic error and the exponentially small ripples in § 7.
We note that the presence of this bifurcation from the side branches 
$S_n$ to the fingers 
$G_{n \to n+1}$ can be observed by a change in sign of the Jacobian along 
$S_n$ as the bifurcation point is passed. Solutions close to this bifurcation point are shown by 
$(a)$ and 
$(g)$ in figure 6 for 
$S_n$ (dashed) and 
$G_{n \to n+1}$ (solid). A further change of sign in the Jacobian occurs at the top of each of the fingers. Since our numerical scheme permits asymmetric solutions through the retention of the asymmetric coefficients in the Fourier series expression (2.7), it may be that this additional bifurcation involves symmetry breaking. A brief overview of gravity-capillary works containing asymmetry is provided in § 7. However, no such solutions were found during our investigation.
Furthermore, the range of 
$F$ between the tip of each finger and the bottom remains of 
$O(1)$ as 
$B \to 0$ for the solutions calculated in figure 8. Consider, for instance, the range between solutions 
$(a1)$ and 
$(\,f)$. This suggests the existence of an interval of solutions holding under the 
$B \to 0$ limit. The solution with the largest value of 
$F$ is expected to be the fundamental Stokes wave with 
$B=0$ and 
$\mathcal {E}=0.3804$, shown by the point 
$y_0$ in figure 8. We predict that, as 
$B \to 0$, the solutions with the smallest Froude number in this interval will contain a self-intersecting free surface. This is because, for 
$B< B_{{crit}}$, the pure 
$n$-solutions near the bifurcation point on the side branches are also anticipated to self-intersect. The interval would then contain a range of solutions, which transition from unphysical to physical as the Froude number increases. A further detail of the structure of these solutions may be seen from figure 9, which shows the exchange between the three components of kinetic, capillary and gravitational potential energies. For the parasitic solutions at the top of each of the fingers, which resemble a perturbation about a gravity-dominated wave, the capillary energy is small and appears to decrease to zero as 
$B \to 0$. Conversely, for the highly oscillatory solutions close to the bifurcation point between adjacent fingers, the capillary and kinetic energies are seen to tend to an 
$O(1)$ constant while the gravitational potential energy tends to zero. Thus, these highly oscillatory solutions of 
$G_{n \to n+1}$, as well as those on the side branch 
$S_n$, appear to tend towards a pure-capillary solution as 
$B \to 0$. The asymptotic properties of the solutions on 
$G_{n \to n+1}$ and 
$S_n$ will be discussed in § 7 for the limit of 
$B \to 0$.
Figure 9. 
$(a)$ The corresponding kinetic (dotted), gravitational (dash–dotted) and capillary (solid) energies for the solutions from figure 8. In the two lower figures, the kinetic, capillary and gravitational energies are also shown for 
$(b)$ 
$B$ vs 
$\mathcal {E}$ and 
$(c)$ branch arclength vs 
$\mathcal {E}$ for the solutions corresponding to the single finger 
$G_{28 \to 29}$.
6. Relation to previous numerical attempts
A key challenge is to understand the relationship between our solutions of fixed energy, 
$\mathcal {E}$, and those of previous authors with a different amplitude condition, say 
$\mathscr {A}$. In this section we demonstrate that a key limitation of previous choices of amplitude is the existence of highly energetic (and subsequently nonlinear) solutions at small values of 
$\mathscr {A}$. Thus, somewhat surprisingly, alternative choices of the amplitude measure may admit nonlinear solutions in the naive linear limit of 
$\mathscr {A} \to 0$ – this occurs due to the singular nature of 
$B \to 0$ and, in particular, the nature of the solutions between adjacent fingers.
6.1. Solutions at different values of the energy $\mathcal {E}$
In the previous section we demonstrated the structure of the bifurcation diagram and associated solutions at fixed energy 
$\mathcal {E}=0.3804$. In fact, this bifurcation structure is only perturbed in a regular fashion as the energy changes near this value. Thus, the full structure of solutions, which holds as 
$B \to 0$, can be computed for different values of 
$\mathcal {E}$ in a straightforward manner.
We show an example of this in figure 10, where we display the finger 
$G_{11 \to 12}$ and the side branches 
$S_{11}$ and 
$S_{12}$ for three different values of 
$\mathcal {E}$. In the figure the value of 
$\mathcal {E}$ decreases from 
$\mathcal {E}=0.67$ in 
$(a)$ to 
$\mathcal {E}=0.3804$ in 
$(b)$ to 
$\mathcal {E}=0.046$ in 
$(c)$. The following three changes to either the solution or branch structure are noticeable as 
$\mathcal {E}$ decreases:
(i) the amplitude of the ripples decreases;
(ii) the range of
$F$ between the top and bottom of the finger decreases; and(iii) the finger becomes more rectangular.
Figure 10. Shown in the left subplots are the fingers 
$G_{11 \to 12}$ and the side branches 
$S_{11}$ and 
$S_{12}$, plotted in the 
$(B,F)$-plane, for 
$(a)$ 
$\mathcal {E} = 0.67$; 
$(b)$ 
$\mathcal {E} = 0.3804$; 
$(c)$ 
$\mathcal {E} = 0.046$. Example solutions, near the tops of the fingers, are shown in the corresponding right subplots labelled 
$(a1)$, 
$(b1)$ and 
$(c1)$.
In 
$(c)$ the amplitude of the ripples has decreased to the point at which they are no longer observable visually.
6.2. Choice of amplitude parameters in previous works
We now revisit the alternative choices of the amplitude or energy parameter in (2.13). A few of the solutions displayed in figure 8 are similar to those previously calculated by Schwartz & Vanden-Broeck (Reference Schwartz and Vanden-Broeck1979), who plotted remnants of this figure at larger values of 
$B$ for a different amplitude parameter. Since their choice of amplitude,
\begin{equation} \mathscr{A} = A \equiv [y(0)-y({\rm \pi})]/2{\rm \pi}, \end{equation}relies on local values at the centre and edge of the periodic domain, they found these branches to behave somewhat differently than how we have described them in our § 5.
Notice that according to their choice of norm (6.1), waves with an even number of crests that are equally spaced throughout the domain will have 
$y(0)=y({\rm \pi} )$ and, consequently, 
$A=0$. This corresponds to every other branch of solutions with a fundamental wavelength smaller than the periodic domain, 
$S_{n}$ with 
$n$ even. Hence, for the branch of solutions 
$G_{n \to n+1}$ at fixed 
$\mathcal {E}$, 
$A$ grows smaller tending towards solutions near the bifurcation points of 
$S_{n+1}$ and 
$S_{n}$ – despite the high nonlinearity of these solutions. This is demonstrated for 
$n=13$ in figure 6 with solution 
$(a)$, which approaches 
$S_{14}$, and solution 
$(g)$, which approaches 
$S_{13}$. Thus, the bifurcation from 
$S_{14}$ connecting finger 
$G_{14 \to 15}$ to 
$G_{13 \to 14}$ will occur from an amplitude value of 
$A=0$. This is one reason why the full structure of solutions was not revealed through smooth continuation at fixed 
$A$ by the investigations of Schwartz & Vanden-Broeck (Reference Schwartz and Vanden-Broeck1979).
Next, let us turn to the numerical investigation of the 
$B \to 0$ limit performed by Chen & Saffman (Reference Chen and Saffman1980b), who fixed the first Fourier coefficient,
\begin{equation} \mathscr{A}=A_1, \end{equation}
as an amplitude parameter. We now know that, since the bifurcation between distinct fingers in the 
$(B,F)$-plane occurs via the side branches 
$S_n$, which have a first Fourier coefficient of zero for 
$n \geqslant 2$, it is impossible to recover the structure shown in figure 8 with a fixed value of 
$A_1$. Chen & Saffman (Reference Chen and Saffman1980b) had indicated the impossibility of a continuous deformation to the pure Stokes gravity wave as 
$B \to 0$, but we now see that this occurred as a direct consequence of their chosen amplitude parameter.
6.3. An insufficient number of Fourier coefficients
As we have noted, it is crucial to select the right continuation parameter in order to recover the 
$B \to 0$ limit. There are other possible reasons why others may have struggled to reproduce an accurate structure of the parasitic ripple phenomena. In particular, a large number of Fourier modes are required in order to capture the regions between adjacent fingers, and this is primarily due to the bifurcation occurring from the side branches, 
$S_{n}$, which contains solutions that approach pure 
$n$-waves. Thus, solutions within the finger 
$G_{n \to n+1}$, which are located near to side branches are then dominated by the 
$n$th Fourier coefficient. If in our numerical scheme we consider a series truncation at the 
$N$th Fourier coefficient, then the main coefficients contributing to the capillary-dominated ripples will be a multiple of 
$n$. Hence, an effective number of 
$N/n$ Fourier coefficients will describe the behaviour of the wave near to this bifurcation point.
For the computation of the gravity-capillary wave with the parasitic ripples, Schwartz & Vanden-Broeck (Reference Schwartz and Vanden-Broeck1979) (their figure 10) used 
$N = 40$ in order to capture a wave with 
$n=11$ ripples. Thus, in order to investigate the side branch bifurcation associated with this solution, their Fourier expansions would have contained an effective number of 
$N/n \approx 4$ Fourier coefficients – which is insufficient. Within this work, we have been using 
$1024$ Fourier coefficients, which corresponds to 
$35$ effective coefficients for solutions near the bifurcation point of the finger with the smallest Bond number in figure 8.
6.4. The symmetry shifting bifurcation
In addition to the importance of selecting an appropriate amplitude measure, let us discuss the relationship between the bifurcation structure presented earlier (e.g. in our figure 8) with the constraint on the symmetry in the travelling-wave frame. For the solutions displayed in figure 8, each finger is computed beginning with an initial solution that lies on the finger and then, with a fixed 
$\mathcal {E}$, solutions are obtained by continuation to either side of the starting point until the entire finger is computed. As a result of this continuation scheme, the solutions at the bottom of adjacent fingers are out of phase with one another; this can be seen in solutions 
$(d1)$ and 
$(d2)$ in figure 8. This method of continuation is depicted more clearly in figure 11
$(a)$, where 
$(a1)$ and 
$(a2)$ are two starting solutions, while 
$(a3)$ and 
$(a4)$ are out of phase. Note also that this phase shift is observable between the profiles 
$(g)$ and 
$(e)$ in figures 6 and 8, respectively, as these are solutions either side of the same bifurcation point.
Figure 11. Two methods for numerical continuation are depicted. The starting location for continuation is denoted by a cross, and the arrows indicate the direction travelled by continuation.
Alternatively, we could formulate a continuation scheme where the solutions in each finger are connected to those in adjacent fingers in a continuous fashion, depicted in figure 11(b). Thus, for example, the scheme is started with a single initial point, 
$(b1)$, shown in figure 11. This finger is then found via the typical continuation method. Having located a solution, 
$(b3)$, at the bifurcation point, the adjacent finger is completed by using 
$(b3)$ as a starting solution for continuation. This alternative method shown in figure 11(b) results in solutions 
$(b3)$ and 
$(b4)$ at the bottom of consecutive fingers with no phase shift. The result of this approach is a continuous set of solutions as 
$B \to 0$.
In using this alternative method, solutions at the top of consecutive fingers have a shifted point of symmetry, as demonstrated by comparing solutions 
$(b1)$ and 
$(b2)$ in figure 11. This point of symmetry has been moved from 
$x=0$ to 
$x=-1/n$ for all solutions on the new finger. We denote this to be a symmetry shifting bifurcation, which is unable to be captured if the point of symmetry of the wave is prespecified. This assumption of a fixed point of symmetry is often used in the two following methods.
(i) Numerical procedures that solve for the half-domain
$x=[0,1/2]$ and enforce a turning point at $x=0$, such as that by Schwartz & Vanden-Broeck (Reference Schwartz and Vanden-Broeck1979).(ii) The analytical work of Chen & Saffman (Reference Chen and Saffman1979), who posit a weakly nonlinear solution with assumed symmetry at
$x=0$.
Both of these methods will be unable to capture this 
$B \to 0$ limit with continuous solutions at the bifurcation point.
The relaxation of the fixed point of symmetry, in contrast to the assumption in Chen & Saffman (Reference Chen and Saffman1979), is the modification required to correct their earlier statement on the validity of the 
$B \to 0$ limit for gravity-capillary waves.
6.5. Conclusions
We have studied the bifurcation structures derived in the previous works by Schwartz & Vanden-Broeck (Reference Schwartz and Vanden-Broeck1979) and Chen & Saffman (Reference Chen and Saffman1979, Reference Chen and Saffman1980b), and have highlighted the following three core issues that are important for unifying and extending their diagrams to the small Bond number regime.
(i) The chosen amplitude parameters, relying either on local values of the wave height or specific Fourier coefficients (cf. § 6.2).
(ii) A small number of Fourier coefficients retained in the numerical schemes, the issues for which become more prominent near the bifurcation points (cf. § 6.3).
(iii) Assumptions made on the point of symmetry of the wave profile (usually fixed to be at
$x=0$), due to the symmetry shifting bifurcation connecting adjacent fingers (cf. § 6.4).
In being aware of these, we have introduced alternative methods, either by solving or mitigating the issues. For instance, we have used the wave energy, 
$\mathcal {E}$, as an amplitude parameter. This has allowed us to be able to find a number of different types of solutions to the steep gravity-capillary wave problem existing under the limit of 
$B \to 0$. One of these, the steady symmetric parasitic ripple, is similar to the asymmetric parasitic waves encountered physically and will be the focus of our forthcoming analytical work.
7. Discussion
In § 5.3 we described the two types of solutions that are found as 
$B \to 0$. The first of these are found along the sides of the fingers and the side branches; they can be described by a multiple-scales type expansion that captures the rapid oscillations about a slowly varying mean. As they become increasingly oscillatory (with diminishing 
$F$), they reach an unphysical configuration through the trapping of bubbles for 
$B< B_{{crit}}$ (cf. end of § 5.2). These highly oscillatory solutions near the bifurcation point between adjacent fingers are very interesting. This is because a multiple-scales ansatz yields Crapper's pure-capillary equation, that is, Bernoulli's equation (2.1c) in the absence of gravity, at leading order for the small-scale ripples. Since an exact solution for this is well known by the work of Crapper (Reference Crapper1970), it may be possible that once the boundary conditions of periodicity and energy have been applied to the solvability condition (obtained at the next order) our snaking structure of the fingers can be found analytically. This multiple-scales approach will be the focus of future analytical work by the current authors. We note that this snaking structure is very similar to that found by Chapman & Kozyreff (Reference Chapman and Kozyreff2009) for a version of the Swift–Hohenberg equation appearing in nonlinear optics.
Our focus has been more on the second of these types of solutions – those which correspond to waves with parasitic ripples lying on a gravity wave; these solutions are found near the tops of each finger. In a forthcoming work, Shelton & Trinh (Reference Shelton and Trinh2021), we shall present an asymptotic theory for the description of these parasitic ripples.
The essential details are as follows. For those solutions that correspond to parasitic ripples on gravity waves, we may expand their form as a naive expansion in powers of the Bond number,
\begin{equation} y(x) = \sum_{n=0}^{\infty}B^n y_n(x), \end{equation}
such that in the limit of 
$B \to 0$ we recover the pure gravity wave, 
$y_0$. However, as it turns out, the magnitude of the short wavelength parasitic ripples is exponentially small in the Bond number, 
$B$. Thus,
\begin{equation} \lvert y_{{ripples}} \rvert \sim \exp\left({-\frac{\text{const.}}{B}}\right). \end{equation}
The above can be validated based on our numerical computations in the following way. For each finger, 
$G_{{n}\to {n+1}}$, a solution profile is calculated at the tip of the finger (i.e. the vertex in the 
$(B,F)$-plane). At this point, the magnitude of the parasitic ripple is approximated by examining 
$y - y_0$ and measuring the crest-to-trough amplitude for the oscillation nearest to the edge of the domain, 
$x = 1/2$ (a typical profile is shown in 
$(a2)$ of figure 8). The result of this numerical experiment is shown in figure 12; indeed, the ripple amplitude lies approximately on a straight line in the semilog plot, confirming the exponential smallness of the ripples as 
$B \to 0$.
Figure 12. The exponential scaling in (7.2) is shown (circles) for our numerical solutions by plotting 
$1/B$ vs 
$\log (y_{{ripples}})$. One numerical solution is chosen from each finger, corresponding to that of minimal ripple amplitude. From our forthcoming work, the analytical prediction (line), which depends on the value of 
$\mathcal {E}$, predicts a gradient of approximately 
$-8.2 \times 10^{-3}$.
Thus, in light of (7.2) these parasitic ripples will fail to be captured by (7.1) and must be found beyond-all-orders of the naive asymptotic expansion. Consequently, the use of specialised tools in asymptotic analysis, known as exponential asymptotics, are required (see, e.g. Chapman, King & Adams Reference Chapman, King and Adams1998; Chapman & Vanden-Broeck Reference Chapman and Vanden-Broeck2006; Trinh & Chapman Reference Trinh and Chapman2013). Here, the necessary theory for prediction of the parasitic ripples is analogous, in spirit, to theories for the prediction of generalised solitary waves (Boyd Reference Boyd1998), but there are a number of additional challenges due to the more involved boundary-integral framework and the lack of a closed-form leading-order Stokes solution, 
$y_0$. Similar bifurcation structures have also been noted in the context of wave-structure interactions of gravity waves, as seen in the works of , for example, Dias & Vanden-Broeck (Reference Dias and Vanden-Broeck2004), Binder, Dias & Vanden-Broeck (Reference Binder, Dias and Vanden-Broeck2008), Holmes et al. (Reference Holmes, Hocking, Forbes and Baillard2013), Hocking, Holmes & Forbes (Reference Hocking, Holmes and Forbes2013); we would expect that our work here with freely propagating waves can be related to more complex problems where the bottom topography has an appreciable effect.
We remark, in addition, that the idea of exponentially small parasitic capillary ripples in the classic Stokes wave problem is not a new one. Indeed, as we discussed in § 1.1, Longuet-Higgins (Reference Longuet-Higgins1963) had proposed an analytical methodology for the derivation of parasitic ripples (followed by a similar approach in Longuet-Higgins Reference Longuet-Higgins1995). However, both of these approaches are ad hoc in nature, and as noted by Perlin et al. (Reference Perlin, Lin and Ting1993), fail to predict the correct magnitude of the ripples. Other theories, such as the averaged Lagrangian methodology of Crapper (Reference Crapper1970), exhibit similar difficulties in providing rigorous comparison to numerical results – the latter author notes that ‘… [the theory] is not very accurate, but at the timing of writing is probably the best available’ (p. 154). Consequently, the point we emphasise is that the systematic 
$B \to 0$ results we have presented in this paper are crucial for validation of the small surface tension limit. The presentation of a complete exponential asymptotics treatment of 
$B \to 0$ and the importance of prior approaches in inspiring such a methodology will be the focus of our forthcoming work.
Finally, we have only found symmetric solutions of the parasitic ripples problem in this work. For general values of 
$B$ and 
$F$, and not necessarily only for the regime of small 
$B$, we note that there are extensive efforts to search for steady asymmetric solutions. See, for instance, the works on gravity-capillary waves by Zufiria (Reference Zufiria1987b) for finite depth, Shimizu & Shōji (Reference Shimizu and Shōji2012) for infinite depth, and Zufiria (Reference Zufiria1987a) for pure gravity on infinite depth. Indeed, many of the asymmetric profiles in the works of, for example, Shimizu & Shōji (Reference Shimizu and Shōji2012) exhibit similarities to the profiles shown in our work. Thus, it seems likely that the bifurcation structures we have presented in this work form a subset of a much more complicated structure that includes the potential for asymmetry. It remains to be seen if this asymmetry of the steady system would account for that observed in the experimental results, or whether it is necessary to consider unsteady flows, such as the time-dependent Navier–Stokes formulation considered numerically by Mui & Dommermuth (Reference Mui and Dommermuth1995) and Hung & Tsai (Reference Hung and Tsai2009). We further note that asymmetric wave profiles have also been found with the unsteady potential flow formulation considered by Moreira & Peregrine (Reference Moreira and Peregrine2010) over a submerged cylinder.
Acknowledgements
We are grateful to Professor N. Ebuchi for generously providing the photo displayed in figure 1, captured during the experimental work of Ebuchi et al. (Reference Ebuchi, Kawamura and Toba1987). We also thank the anonymous reviewers for their suggestions to improve the clarify of this manuscript.
Declaration of interests
The authors report no conflict of interest.
Appendix A. The wave energy
The non-dimensionalised bulk energy in the physical domain is given by
\begin{align} \bar{E}=\frac{F^2}{2}\int_{-{1/2}}^{{1/2}}\int_{-\infty}^{\eta}(\phi_{x}^2+\phi_{y}^2) \,\textrm{{d}}{y} \,\textrm{{d}}{x} + B \int_{-{1/2}}^{{1/2}}\left([1+\eta_{x}^2]^{{1/2}}-1\right) \textrm{{d}}{x} +\int_{-{1/2}}^{{1/2}}\int_{-\infty}^{\eta} y \, \textrm{{d}}{y} \,\textrm{{d}}{x}. \end{align}On the right-hand side, the three groups correspond to the kinetic, capillary potential and gravitational potential energies.
Note that due to our subflow (where 
$\phi _x \to -1$ as 
$y \to -\infty$), the first and third integrals on the right-hand side of (A1) will be unbounded. We thus define the energy, 
$E$, to be the difference between (A1) and ‘no-flow’, 
$\eta =0$, energy, which yields a finite value. This is then transformed to act on the free surface, 
$\psi =0$, only by the method of Longuet-Higgins (Reference Longuet-Higgins1989), with which we change variables from 
$(x,y)$ to find the wave energy under the 
$(\phi ,\psi )$ mapping,
\begin{equation} E = \frac{F^2}{2}\int_{-{1/2}}^{{1/2}}Y(X_{\phi}-1) \, \textrm{{d}}{\phi} +B \int_{-{1/2}}^{{1/2}}(\sqrt{J}-X_{\phi}) \, \textrm{{d}}{\phi} + \frac{1}{2}\int_{-{1/2}}^{{1/2}}Y^2 X_{\phi} \,\textrm{{d}}{\phi}. \end{equation}
With this choice of amplitude parameter, 
$E_{{hw}} \approx 0.00184$ corresponds to the fundamental Stokes wave of maximum height. We rescale the energy by this value to obtain the amplitude parameter, 
$\mathcal {E}$, used within this report, given by 
$\mathcal {E} = {E}/{E_{{hw}}}$, in (4.1c).
Appendix B. Limiting solutions of smaller fundamental wavelength
If one solution is known to the gravity-capillary wave problem with fundamental wavelength 
$\lambda =1$, another can be constructed with 
$\lambda =1/ \alpha$, where 
$\alpha$ is a positive integer. This is visualised in figure 13. Suppose we have a solution to Bernoulli's equation (4.1a) and the harmonic relation (4.1b). In rescaling 
$Y=\alpha \hat {Y}$, 
$X=\alpha \hat {X}$ and 
$\phi =\alpha \hat {\phi }$, we repeat the first solution 
$\alpha$ times to map the original domain from 
$\phi \in [-{\alpha }/{2}, {\alpha }/{2})$ to the new domain 
$\hat {\phi } \in [-{1}/{2}, {1}/{2})$. This new solution, 
$\hat {X}$ and 
$\hat {Y}$, also satisfies the two governing equations with rescaled Froude and Bond numbers 
$\hat {F}$ and 
$\hat {B}$, given by
\begin{equation} \hat{F}=\frac{F}{\sqrt{\alpha}} \quad \text{and} \quad \hat{B}=\frac{B}{\alpha^2}. \end{equation}
The energy of this new solution can be found by substituting the rescaled variables 
$\hat {X}$, 
$\hat {Y}$, 
$\hat {F}$ and 
$\hat {B}$ into (4.1c), yielding
\begin{equation} \hat{\mathcal{E}}=\frac{ \mathcal{E}}{\alpha^2}. \end{equation}
Since 
$\alpha >1$, the energy of the new pure 
$\alpha$-wave, 
$\hat {\mathcal {E}},$ will always be smaller than that of the original wave, 
$\mathcal {E}$.
Figure 13. The rescaling used to produce another solution with a smaller fundamental wavelength is shown.
This ability to construct new solutions with a fundamental wavelength shorter than the periodic domain allows us to numerically predict the point at which solutions in the side branches 
$S_{\alpha }$ begin to trap a bubble through a self-intersecting free surface. These locations are shown in the bifurcation diagram of figure 7. As all of these side branch solutions have the same energy, 
$\hat {\mathcal {E}}=0.3804$, but different values of 
$\alpha$, the energy of the original wave is given by 
$\mathcal {E}=\alpha ^2 \hat {\mathcal {E}}$.
The procedure to find the location at which solutions in 
$S_{\alpha }$ become self-intersecting is as follows.
(i) First, we numerically calculate the
$(B,F)$-solution space of pure $1$-waves with a single trapped bubble amplitude condition. The energy of these solutions will vary.(ii) Second, we select the profile with
$\mathcal {E}=\alpha ^2 \hat {\mathcal {E}}$ and obtain values for $B$ and $F$.(iii) Third, we rescale these by using (B1a) to find
$\hat {F}$ and $\hat {B}$. This yields the location at which solutions within the side branch, $S_{\alpha }$, become self-intersecting.
Repeating this process for multiple values of 
$\alpha$ yields the predictions displayed with the circles in figure 7. With this method we are able to calculate solutions with a large value of 
$\alpha$ while keeping the number of Fourier coefficients used during Newton iteration fixed, and, thus, do not encounter the issue discussed in § 6.3.
It would also be possible to use this method to compute all of the solutions along the side branch 
$S_{\alpha }$. However, this requires the entire sheet of pure 
$1$-wave solutions to be found in the three-dimensional 
$(B,F,\mathcal {E})$-solution space, which we consider to be prohibitively expensive computationally. By restricting only to profiles displaying a single trapped bubble, this solution space simplifies to a single branch throughout the 
$(B,F,\mathcal {E})$-bifurcation diagram, which we projected to the 
$(B,F)$-plane for simplicity.
B.1. Limiting pure-$1$ solution space
This branch of limiting solutions is displayed in figure 14(a). These solutions were found from the same numerical procedure as in § 4. An initial limiting solution, displaying one trapped bubble, is found by increasing 
$\mathcal {E}$. The energy constraint is then replaced by a trapped bubble condition, which forces the second turning point of 
$X(\phi )$ to a value of 
$-0.5$. We then explore the 
$(B,F)$-solution space by continuation.
Figure 14. 
$(a)$ The branch of pure-
$1$ solutions displaying an enclosed bubble is shown in the 
$(B,F)$-bifurcation diagram. Note that we have displayed the domain 
$x \in [0,1]$ to demonstrate this limiting behaviour.
We note that as 
$B \to \infty$ along this branch, the wave profile approaches the limiting pure-capillary solution found by Crapper (Reference Crapper1957) (see their figure 1) with an amplitude of 
$0.730$. Solution 
$(b)$, with 
$B=4.050$, is an example of this. As 
$B \to 0$ along the same branch, the solution approaches a depressive solitary wave, demonstrated by solution 
$(d)$. This is since the small 
$B$ limit is related to the solitary wave limit of 
$L_{\lambda } \to \infty$. The solutions calculated by Schwartz & Vanden-Broeck (Reference Schwartz and Vanden-Broeck1979) (see their figure 2) form the intermediate range between these two limits.
