00 Paper

J. Shelton and P. H. Trinh

[*]Corresponding author<EMAIL_ADDRESS>
82023

# A model ODE for the exponential asymptotics of
nonlinear parasitic capillary ripples

Josh Shelton*0000-0001-6257-5190 and Philippe H. Trinh 0000-0003-3227-1844
Department of Mathematical Sciences, University of Bath, Bath, BA2 7AY, UK

(2023; 2023; 14)

###### Abstract

In this work, we develop a linear model ODE to study the parasitic capillary
ripples present on steep Stokes waves when a small amount of surface tension
is included in the formulation. Our methodology builds upon the exponential
asymptotic theory of Shelton & Trinh (J. Fluid. Mech., vol. 939, 2022, A17),
who demonstrated that these ripples occur beyond-all-orders of a small-
surface-tension expansion. Our model equation, a linear ODE forced by
solutions of the Stokes wave equation, forms a convenient tool to calculate
numerical and asymptotic solutions. We show analytically that the parasitic
capillary ripples that emerge in solutions to this linear model have the same
asymptotic scaling and functional behaviour as those in the fully nonlinear
problem. It is expected that this work will lead to the study of parasitic
capillary ripples that occur in more general formulations involving viscosity
or time-dependence.

###### keywords:

Exponential asymptotics; gravity-capillary waves; parasitic capillary ripples;
free-surface flows

##  1 Introduction

Consider the propagation of a steep gravity-driven water wave of permanent
form [a _Stokes wave_ , after Stokes (1847)]. It is plausible that the
addition of a small amount of surface tension perturbs the gravity wave.
Indeed, in experimental observations of high amplitude water waves, _parasitic
capillary ripples_ are observed to form due to this small capillary
effect—these appear as small-scale oscillations located ahead of the wave
crest [cf. figure 11 from Ebuchi et al. (1987)]. However, despite the
intuitive nature of such waves, the limit of small surface tension forms a
singular perturbation in the water-wave equations, and it is not clear that
such solutions are admissible in the framework of potential flow theory.
Recently, it was discovered (Shelton et al., 2021; Shelton & Trinh, 2022) that
when limited to a certain subclass of solutions, such parasitic ripples do
appear as an exponentially-small perturbation of the underlying steadily-
travelling Stokes wave. The analysis of the ripples is challenging, and our
previous presentation required significant work in applying the techniques of
exponential asymptotics for their study. In this paper, we present the
derivation and analysis of an asymptotically reduced model for these parasitic
capillary ripples. Despite its apparent simplicity (a second-order linear
forced ordinary differential equation) we shall show that the model is able to
preserve many of the key features of the full water-wave problem.

The specific formulation we consider is that of a fully nonlinear gravity-
capillary wave propagating upon an inviscid, irrotational, and incompressible
fluid of infinite depth. Solutions are characterised by an amplititude or
energy, $\mathscr{E}$, Froude number, $F$, and Bond number, $B$. As noted
above, a systematic numerical study of solutions to this formulation with
small surface tension, $B\to 0$, was recently performed by Shelton et al.
(2021). There, it was demonstrated that a discrete set of solution branches
exist in the small surface tension limit; as the branches are traversed in the
direction of $B\to 0$, more ripples appear on the surface of the underlying
Stokes wave. The authors confirmed, via a similar study to that displayed in
figure 1, that the amplitude of the parasitic ripples is exponentially small
and of $O(\mathrm{e}^{-\mathrm{const.}/B})$, and therefore their contributions
lie beyond-all-orders of a traditional asymptotic expansion in powers of $B$.
This motivated the study by Shelton & Trinh (2022), who analytically solved
for these capillary ripples using techniques in exponential asymptotics. In
the latter work, it was shown that the parasitic ripples emerge in connection
with Stokes lines in the analytic continuation of the free-surface. This
analysis was complicated partly due to the fact that the full nonlinear water-
wave model must be studied in the form of a nonlinear integro-differential
equation.

Figure 1:  Numerical solutions of the full gravity-capillary wave problem are
shown at fixed energy $\mathscr{E}=0.4$ for different values of the surface-
tension (or Bond) parameter, $B$. The main plot shows the logarithm of the
capillary-ripple amplitude as a function of $1/B$. The approximate straight-
line fit is indicative that this amplitude is of
$O(\mathrm{e}^{-\mathrm{const.}/B})$ as $B\to 0$. The insets show the physical
wave profile, $y=\eta(x)$, at two selected points of the main diagram. These
solutions are computed using the algorithm discussed in Shelton et al. (2021).
A comparison between these amplitudes, and those predicted both numerically
and asymptotically by our reduced model occurs in figure 5.

However, during the course of the research that led to the asymptotic results
of Shelton & Trinh (2022), the authors discovered that the parasitic capillary
ripples could also be approximated, in a well-defined asymptotic sense, via a
much simpler linear differential equation for a single unknown. Later, we
shall present this as the “$N=2$ model”, shown in equation (3.4a). Although
the equation is linear and therefore schematically simple, its coefficients
are specified in terms of the leading-order nonlinear Stokes wave (and its
higher-order capillary correction). This simplification was not presented in
our previous work. However, previously similar exponential-asymptotic
simplifications have used in fluid mechanical studies in different contexts
[notably water-waves driven by structures and topography by Tulin (1984); Tuck
(1991); Trinh (2016, 2017); Jamshidi & Trinh (2020); Kataoka & Akylas (2022);
and interfacial Hele-Shaw flows by Tu (1991)].

The ability to model parasitic capillary ripples on steep Stokes waves using a
second-order linear differential that is ‘forced’ by numerical data (the
Stokes wave) is a intriguing concept. The complexity of modelling this
behaviour has motivated a number of prior minimal models. For instance, we
recall the work of Longuet-Higgins (1995), where the capillary ripples were
interpreted as being generated by a surface pressure distribution. The work of
Crapper (1970) posed a variational formulation to produce a model ODE for the
capillary ripples, which was solved numerically. We note also the related
formulation of gravity-capillary solitary waves, which predicts the formation
of small-scale ripples on the surface. Models that study this include the
well-known fifth-order Kortweg-de Vries equation (Pomeau et al., 1988;
Grimshaw & Joshi, 1995), and the forced KdV equation from Akylas & Yang
(1995).

However, in contrast to previous simplified models for parasitic capillary
ripples, our model seems distinguished by its “asymptotically-exact” nature.
We show in this paper that all of the important features of the asymptotic
solutions of Shelton & Trinh (2022), such as the functional behaviour and
exponentially-small scaling, are preserved in our model equation. The main
downside to our assumptions is an under-prediction of the capillary wave
amplitude (by a factor of $\approx 0.4114$). Our work opens the door to the
development of similar minimal models in more complicated formulations of
surface waves with capillary ripples. These include time-dependent waves and
also waves with vorticity and viscosity; such extensions are reviewed near
figure 6 of our discussion.

### Outline of our paper

In §2 we formulate the free surface equations and amplitude condition used to
describe steadily travelling gravity-capillary waves. Our model equation is
then developed in §3. Numerical solutions of this are found in §4, and
comparison between the asymptotic scaling of solutions to the model and full
problem occurs in §4.3. An exponential asymptotic analysis of our model
equation occurs in §5. Future directions involving the additional physical
effects of viscosity, vorticity, and time dependence are discussed in §6.

##  2 Mathematical formulation

We consider the two-dimensional free-surface flow of an inviscid,
irrotational, and incompressible fluid of infinite depth. In including the
effects of both gravity and surface tension, we search for steadily travelling
free surface waves that are spatially periodic in the direction of
propagation. All quantities are nondimensionalised with respect to the
wavelength, $\lambda$, and wavespeed, $c$. The boundary-value problem then
consists of Laplace’s equation for the velocity potential, $\phi(x,y)$, within
the fluid domain, kinematic and dynamic boundary conditions on the unknown
free-surface, $y=\eta(x)$, decay conditions in the deep-water limit of
$y\to-\infty$, and periodicity of the solutions at $x=-1/2$ and $x=1/2$.

As the unknown free-surface, $y=\eta(x)$, is a streamline along which the
streamfunction, $\psi(x,y)$, is constant, we may map the physical
$(x,y)$-domain, $-1/2\leq x\leq 1/2$ and $-\infty<y\leq\eta(x)$, to the lower-
half potential $(\phi,\psi)$-plane. The free-surface may then be parameterised
by the velocity potential, $-1/2\leq\phi\leq 1/2$, which forms the independant
variable of our formulation. We note that the complex velocity may be written
as $u-\mathrm{i}v=q\mathrm{e}^{-\mathrm{i}\theta}$, where $q(\phi)$ and
$\theta(\phi)$ are the streamline speed and angle, respectively, along the
free-surface. We will consider the boundary-integral formulation of the
gravity-capillary wave problem, in terms of $q$ and $\theta$, which consists
of the equations

$\displaystyle
F^{2}q^{2}\frac{\mathrm{d}q}{\mathrm{d}\phi}+\sin{(\theta)}-Bq\frac{\mathrm{d}}{\mathrm{d}\phi}\bigg{(}q\frac{\mathrm{d}\theta}{\mathrm{d}\phi}\bigg{)}=0,$
(2.1a)
$\displaystyle\log{(q)}=\mathscr{H}[\theta]=\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{-1/2}^{1/2}\theta(\phi^{\prime})\cot{[\pi(\phi^{\prime}-\phi)]}\mathrm{d}\phi^{\prime}.$
(2.1b)

Together, equations (2.1a) and (2.1b) form a coupled set of nonlinear integro-
differential equations for the solutions, $q(\phi)$ and $\theta(\phi)$, with
associated constants $F$ and $B$.

In the above, Bernoulli’s equation (2.1a) contains the effects of both gravity
and surface tension, and the boundary-integral equation (2.1b), derived from
the analyticity of $u-\mathrm{i}v=\mathrm{e}^{\log{(q)}-\mathrm{i}\theta}$,
closes the system through the Hilbert transform, $\mathscr{H}[\theta]$. The
two nondimensional constants appearing in (2.1a), the Froude number, $F$, and
the Bond number, $B$, are given by

$F=\frac{c}{\sqrt{g\lambda}}\quad\text{and}\quad B=\frac{\sigma}{\rho
g\lambda^{2}},$ (2.2)

where $c$ is the wavespeed, $g$ is the acceleration due to gravity, $\lambda$
is the wavelength, $\sigma$ is the coefficient of surface tension, and $\rho$
is the fluid density. The intention of this work is to study nonlinear
solutions of (2.1) in the small surface tension limit of $B\to 0$. We will
asymptotically develop a linear second-order ODE (with coefficients and
forcing terms determined from nonlinear equations) to study solutions in this
regime.

###  2.1 The amplitude condition

Linear solutions to system (2.1) bifurcate from a Froude number given by
$F^{2}=(1+4k^{2}\pi^{2}B)/(2k\pi)$, where $k$ is the wavenumber. Thus, rather
than regard $F$ as a fixed constant for which solutions are sought, it is
convenient to introduce an amplitude parameter, $\mathscr{E}$, for which
enforcing $\mathscr{E}$ via an amplitude condition yields $F$ as an
eigenvalue. Linear solutions will then bifurcate from $\mathscr{E}=0$, and it
is intended for solutions to become increasingly nonlinear as the value of the
amplitude increases.

For comparison with the numerical and asymptotic results by Shelton et al.
(2021) and Shelton & Trinh (2022), we specify the amplitude parameter to be
the wave energy, which is an integral condition comprising of kinetic,
capillary, and gravitational potential energy components. In the time-
dependant problem this is a conserved quantity, and thus has also been used as
an amplitude condition in the computation of nonlinear standing waves in the
small surface tension limit by Shelton et al. (2023b). When expressed in terms
of $q$ and $\theta$, and the integrand is rearranged into powers of $B$, the
energy is given by

$\displaystyle\mathscr{E}=\frac{1}{E_{hw}}\int_{-1/2}^{1/2}\bigg{[}$
$\displaystyle\frac{F^{4}}{8q}(1-q^{2})\Big{(}(3-q^{2})\cos{\theta}-2q\Big{)}$
(2.3)
$\displaystyle\qquad+B\bigg{(}\frac{(1-\cos{\theta})}{q}+\frac{F^{2}\theta_{\phi}}{2}\Big{(}(2-q^{2})\cos{\theta}-q\Big{)}\bigg{)}+B^{2}\frac{q^{2}\theta_{\phi}^{2}\cos{\theta}}{2}\bigg{]}\mathrm{d}\phi.$

In the above, we have normalised with respect to the energy of the highest
Stokes wave, $E_{hw}\approx 0.00184$.

Other choices of amplitude could be used, such as the wave displacement,
$\eta(0)-\eta(1/2)$, from Schwartz & Vanden-Broeck (1979), individual Fourier
coefficients from Chen & Saffman (1980), or the parameter
$1-q(0)^{2}q(1/2)^{2}$ from Longuet-Higgins (1975) and Cokelet (1977). These
amplitude conditions were used either to solve for weakly nonlinear solutions
through expanding in powers of the small amplitude parameter, or to calculate
fully nonlinear solutions numerically. We note that, inspired from the work by
Shimizu & Shōji (2012), there have been many recent studies which calculate
asymmetric solutions by fixing individual Fourier coefficients (corresponding
to asymmetric modes). While each of these amplitude conditions only affects
the ordering of solutions within the three-dimensional Bond, Froude, and
Amplitude bifurcation space, there are a number of significant differences
that arise when each condition is fixed.

  1. (i)

Bifurcation structure.
It was shown by Shelton et al. (2021) that when the wave energy (2.3) is fixed
as an amplitude parameter, a self-similar bifurcation structure emerges under
the small surface tension limit of $B\to 0$. Portions of analogous bifurcation
structures, determined with other fixed amplitude measures, have been
calculated by e.g. Schwartz & Vanden-Broeck (1979), and Champneys et al.
(2002) for finite depth. The choice of amplitude is seen to have a significant
effect on the resultant bifurcation structure.

  2. (ii)

Asymptotic solutions for small surface tension.
In order to satisfy each order of the amplitude condition when considering
asymptotic solutions for $q$ and $\theta$ in the limit of $B\to 0$, it is also
necessary to asymptotically expand the eigenvalue, $F$. The $n$th order of
each asymptotic expansion, $q_{n}$, $\theta_{n}$, and $F_{n}$, is determined
by the $n$th order Bernoulli’s equation, the boundary-integral equation, and
the amplitude condition. There then exists a value of each different amplitude
parameter such that the same gravity-wave solutions, $q_{0}$, $\theta_{0}$,
and $F_{0}$, are found at leading-order. However, for this respective fixed
value, the lower-order solutions, $q_{1}$, $\theta_{1}$, and $F_{1}$ will
differ between each choice of amplitude. These lower-order components feed
into the solution for the exponentially-small parasitic ripples specified
later in §5.2. Thus, we conclude that different choices of amplitude parameter
lead to subtly different behaviour of the parasitic capillary ripples as $B\to
0$.

###  2.2 Analytic continuation

A key concept in the exponential asymptotics of singular perturbation problems
is that singularities of an asymptotic solution generate a divergent
expansion. This divergence is then responsible for the Stokes phenomenon that
occurs on exponentially small components of the solution. The parasitic
capillary ripples our model problem intends to capture are known to be
exponentially-small in magnitude as $B\to 0$, and thus it is necessary to
consider singularities of the leading-order solution for $B=0$ (a Stokes
wave). However, for solutions of less than maximal amplitude, no singular
behaviour is present when $\phi\in\mathbb{R}$. It is only when the complex-
valued domain, $\phi\mapsto f\in\mathbb{C}$, is considered that these singular
points are identified.

To find solutions in the complex-valued domain, we must analytically continue
the governing equations by taking $\phi\mapsto f$, where $f\in\mathbb{C}$.
Bernoulli’s equation (2.1a) may be analytically continued by simply
relabelling $\phi$ by $f$. Analytic continuation of the boundary-integral
equation (2.1b) is more complicated due to a singularity in the integrand of
the principal-valued Hilbert transform, which is lost when $\text{Im}[f]\neq
0$. Thus, when we analytically continue the Hilbert transform, we produce a
residue contribution $-a\mathrm{i}\theta$, which yields

$\underbrace{\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{-1/2}^{1/2}\theta(\phi^{\prime})\cot{[\pi(\phi^{\prime}-\phi)]}\mathrm{d}\phi^{\prime}}_{\textstyle\vphantom{\hat{\mathscr{H}}[\theta](f)}\mathscr{H}[\theta](\phi)}=-a\mathrm{i}\theta(f)+\underbrace{\int_{-1/2}^{1/2}\theta(\phi^{\prime})\cot{[\pi(\phi^{\prime}-f)]}\mathrm{d}\phi^{\prime}}_{\textstyle\hat{\mathscr{H}}[\theta](f)}.$
(2.4)

Here, $a=\pm 1$ denotes the direction of analytic continuation into either the
upper- or lower-half $f$-plane by

$a=\left\\{\begin{aligned} +1\qquad\text{for $\text{Im}[f]>0$},\\\
-1\qquad\text{for $\text{Im}[f]<0$}.\end{aligned}\right.$ (2.5)

We have introduced in (2.4) the notation $\hat{\mathscr{H}}$ for the complex-
valued Hilbert transform, which is a function of $f\in\mathbb{C}$, but
involves integration over $\phi^{\prime}\in\mathbb{R}$. The fact that the
complex-valued Hilbert transform only requires knowledge of
$\theta(\phi^{\prime})$ for $\phi^{\prime}\in\mathbb{R}$ along the real axis
will motivate our reduced model equation, in which $\hat{\mathscr{H}}$ is
subdominant for later components of the asymptotic expansion of $\theta$.

The analytically continued equations are then given by

$\displaystyle
F^{2}q^{2}\frac{\mathrm{d}q}{\mathrm{d}f}+\sin{(\theta)}-Bq\frac{\mathrm{d}}{\mathrm{d}f}\bigg{(}q\frac{\mathrm{d}\theta}{\mathrm{d}f}\bigg{)}=0,$
(2.6a) $\displaystyle\log{(q)}+a\mathrm{i}\theta=\hat{\mathscr{H}}[\theta],$
(2.6b)

for the solutions $q(f)$ and $\theta(f)$, with $f\in\mathbb{C}$. We note that
the energy condition (2.3) is enforced along the real-valued free surface,
upon which $\text{Im}[f]=0$. There are two ways to interpret system (2.6).

  1. (i)

First, we may treat (2.6) as an initial value problem for which the solution
along the real-axis, $q(\phi)$, $\theta(\phi)$, $B$, and $F$, is known from
the real-valued system (2.1). This known solution feeds into the complex-
valued Hilbert transform on the right-hand side of (2.6b), and we may then
solve (2.6) with $B$ and $F$ specified to find $q(f)$ and $\theta(f)$ for $f$
along paths in $\mathbb{C}$. This method was implemented numerically with
$B=0$ by Crew & Trinh (2016) to find branch points in the analytic
continuation of gravity waves, and with $B\neq 0$ by Shelton & Trinh (2022) to
observe the parasitic capillary ripple amplitude increasing in the complex
domain.

  2. (ii)

Second, we may solve (2.6) in conjunction with energy condition (2.3), for the
solutions $q$, $\theta$, and $F$. Numerically, this requires for the
implementation of a coupled scheme, in which $q\rvert_{a=1}$ and
$\theta\rvert_{a=1}$ are solutions of (2.6) with $a=1$, and $q\rvert_{a=-1}$
and $\theta\rvert_{a=-1}$ are solutions with $a=-1$ (these are coupled as they
both involve the same eigenvalue $F$). The energy condition (2.3) is then
evaluated along the real-axis with $q=(q\rvert_{a=1}+q\rvert_{a=-1})/2$ and
$\theta=(\theta\rvert_{a=1}+\theta\rvert_{a=-1})/2$. This is the method we use
in §4 to numerically validate solutions of our model equation.

##  3 Development of the model equation

In the small surface tension limit of $B\to 0$, the parasitic capillary
ripples we intend to capture are exponentially-small in magnitude. Thus, they
will not be present in a perturbative expansion in algebraic powers of $B$. To
detect this behaviour of the solution of equations (2.6a) and (2.6b), we now
consider truncated asymptotic expansions of the form

$q(f)=\sum_{n=0}^{N-1}B^{n}q_{n}(f)+\bar{q}(f),\qquad\theta(f)=\sum_{n=0}^{N-1}B^{n}\theta_{n}(f)+\bar{\theta}(f),\qquad
F=\sum_{n=0}^{N-1}B^{n}F_{n}+\bar{F}.$ (3.1–c)

In the above, the expansions are truncated at $O(B^{N-1})$, and the remainders
$\bar{q}$, $\bar{\theta}$, and $\bar{F}$ will be of $O(B^{N})$. It was shown
by Shelton & Trinh (2022), following exponential asymptotic techniques
developed by Berry (1989) and Olde Daalhuis et al. (1995), that when this
truncation point is chosen optimally to be $N=O(1/B)$, for which $N\to\infty$
as $B\to 0$, the remainders in (3.1) are exponentially-small.

Equations for the remainder functions $\bar{q}$, $\bar{\theta}$, and
$\bar{F}$, may be found by substituting expansions (3.1) into Bernoulli’s
equation (2.6a) and the boundary-integral equation (2.6b). These may be
simplified to give a single equation for $\bar{q}$ by eliminating
$\bar{\theta}$, an expression for which is found by rearranging (2.6b) and
substituting this into (2.6a). This yields

$\displaystyle\Big{[}a\mathrm{i}Bq_{0}+a\mathrm{i}B^{2}q_{1}+\cdots\Big{]}\bar{q}^{\prime\prime}+\Big{[}-F_{0}^{2}q_{0}^{2}+B\big{(}q_{0}\theta_{0}^{\prime}-a\mathrm{i}q_{0}^{\prime}-2F_{0}^{2}q_{0}q_{1}-2F_{0}F_{1}q_{0}^{2}\big{)}+\cdots\Big{]}\bar{q}^{\prime}$
(3.2a)
$\displaystyle-\Big{[}2F_{0}^{2}q_{0}q_{0}^{\prime}+a\mathrm{i}q_{0}^{-1}\cos{(\theta_{0})}+\cdots\Big{]}\bar{q}-\Big{[}2F_{0}q_{0}^{2}q_{0}^{\prime}+\cdots\Big{]}\bar{F}=\mathcal{R}(q_{\text{reg}},\theta_{\text{reg}},F_{\text{reg}})+O(\bar{q}^{2};\hat{\mathscr{H}}[\bar{\theta}]),$
which is a linear second-order differential equation for $\bar{q}$ with
eigenvalue $\bar{F}$. The forcing term on the right-hand side of (3.2a)
depends on the regular asymptotic series of $q$, $\theta$, and $F$ introduced
in (3.1), which we have denoted by
$q_{\text{reg}}=q_{0}+Bq_{1}+\cdots+B^{N-1}q_{N-1}$ for instance. In defining
Bernoulli’s equation (2.6a) and the boundary integral equation (2.6b)
evaluated on these regular expansions as
$\xi_{\text{bern}}=F_{\text{reg}}^{2}q_{\text{reg}}^{2}q_{\text{reg}}^{\prime}+\sin{(\theta_{\text{reg}})}-Bq_{\text{reg}}(q_{\text{reg}}\theta_{\text{reg}}^{\prime})^{\prime}$
and
$\xi_{\text{int}}=\log{(q_{\text{reg}})}+a\mathrm{i}\theta_{\text{reg}}-\hat{\mathscr{H}}[\theta_{\text{reg}}]$,
the forcing term $\mathcal{R}$ is given by
$\mathcal{R}=\xi_{\text{bern}}+a\mathrm{i}\cos{(\theta_{\text{reg}})}\xi_{\text{int}}-a\mathrm{i}Bq_{\text{reg}}q_{\text{reg}}^{\prime}\xi_{\text{int}}^{\prime}-a\mathrm{i}Bq_{\text{reg}}^{2}\xi_{\text{int}}^{\prime\prime}.$
(3.2b)

When solving equation (3.2a), we regard each order of $q_{\text{reg}}$,
$\theta_{\text{reg}}$, and $F_{\text{reg}}$ as known (either numerically or
analytically) from the corresponding order of Bernoulli’s equation (2.6a), the
boundary-integral equation (2.6b), and the energy constraint (2.3).
Consequently, each order of $\xi_{\text{bern}}$ and $\xi_{\text{int}}$ in
(3.2b) is identically zero up to and including $O(B^{N-1})$. Thus, the forcing
term $\mathcal{R}$ is of $O(B^{N})$. The exponentially-small components of the
asymptotic solutions were derived by Shelton & Trinh (2022), who solved (3.2a)
with $N=O(1/B)$, for which the forcing term $\mathcal{R}$ was shown to be
exponentially-small in $B$ due to divergence of the base expansions. We now
review this beyond-all-orders theory in §3.1 to motivate our reduced model of
§3.2, for which expansions (3.1) are truncated sub-optimally with $N=2$.

###  3.1 The asymptotic theory of Shelton & Trinh (2022)

Using the methodology of exponential asymptotics, Shelton & Trinh (2022)
solved equation (3.2a) under the limit of $B\to 0$. The resultant asymptotic
solutions had an exponentially-small magnitude and an algebraically small
phase, corresponding to the manifestation of highly oscillatory parasitic
ripples under this limit. The steps required in this beyond-all-order analysis
are to:

  1. (i)

Solve for the divergence of $q_{n}$ and $\theta_{n}$ in the limit of
$n\to\infty$.
This is performed by approximating the $O(B^{n})$ terms in equations (2.6a)
and (2.6b) in the limit of $n\to\infty$, for which a differential-difference
equation emerges. The singular perturbative nature of Bernoulli’s equation
(2.6a) results in an $O(B^{n})$ equation for which $q_{n}$ is determined from
$q_{n-1}^{\prime}$, and this leads to the solution diverging as $n\to\infty$
in the factorial-over-power manner of

$q_{n}(f)\sim\Lambda
q_{0}^{2}\frac{\exp{\Big{(}a\mathrm{i}\theta_{0}+a\mathrm{i}\int_{0}^{f}\Big{[}\tfrac{\cos{(\theta_{0})}}{F_{0}^{2}q_{0}^{3}}-F_{0}^{2}q_{1}-2F_{0}F_{1}q_{0}\Big{]}\mathrm{d}f^{\prime}\Big{)}}\Gamma(n+4/5)}{\Big{(}a\mathrm{i}F_{0}^{2}\int_{af^{*}}^{f}q_{0}(f^{\prime})\mathrm{d}f^{\prime}\Big{)}^{n+4/5}}.$
(3.3a) Here, $\Lambda$ is a constant of integration and $f^{*}$ is the
location of the closest branch point of $q_{0}$ to the real-axis. The
divergence of the $\theta$ expansion may then be found from $\theta_{n}\sim
a\mathrm{i}q_{0}^{-1}q_{n}$.

  2. (ii)

Truncate optimally when $N=O(1/B)$.
The base asymptotic series reorders in the limit of $n\to\infty$ when
$B^{n}q_{n}\sim B^{n+1}q_{n+1}$, due to the divergence of $q_{n}$ in (3.3a).
This yields the optimal truncation point $N=\rho+\lvert\chi\rvert/B$, where
$0\leq\rho<1$ ensures that $N$ is an integer, and the singulant
$\chi=a\mathrm{i}F_{0}^{2}\int_{af^{*}}^{f}q_{0}\mathrm{d}f^{\prime}$ is the
function in the denominator of (3.3a). The forcing term, $\mathcal{R}$, from
(3.2b) may then be simplified further by noting that: as $B\to 0$,
$\mathcal{R}$ is of $O(B^{N})$; as $N\to\infty$, only one component,
$-q_{0}^{2}\theta_{N-1}^{\prime\prime}$, of this order is dominant due to the
solution divergence; and when $N\sim\lvert\chi\rvert/B$ this term may be
expanded to find $\mathcal{R}=O(\mathrm{e}^{-\lvert\chi\rvert/B})$. This
yields

$\mathcal{R}\sim-
q_{0}^{2}\theta_{N-1}^{\prime\prime}B^{N}=O(\mathrm{e}^{-\lvert\chi\rvert/B}).$
(3.3b)

  3. (iii)

Capture the Stokes phenomenon by localising in a boundary layer about the
Stokes lines, where $\text{Im}[\chi]=0$ and $\text{Re}[\chi]\geq 0$.
When the divergent solution expansions are truncated optimally, the forcing
term in equation (3.2a) is exponentially small. WKBJ solutions of this
equation display the Stokes phenomenon, in which their magnitude rapidly
varies across Stokes lines. The Stokes lines of interest are those that cross
the real-valued domain. These lie along the imaginary $f$-axis, and intersect
at the wave crest, $\phi=0$. Imposing periodicity then yields the asymptotic
solutions for parasitic capillary waves from Shelton & Trinh (2022).

The authors verified these asymptotic solutions through comparison to
numerical solutions of the original nonlinear equations (2.1a) and (2.1b). In
§3.2 we introduce a linear model, with $N=2$, for capturing the same small
surface tension phenomenon. This reduced model will be seen to predict the
correct exponentially-small scaling and functional behaviour of the parasitic
capillary ripples, and only fails by producing an incorrect constant
magnitude.

###  3.2 The model $N=2$ equation

We now present a linear second-order model ODE to capture the parasitic
capillary ripples present in gravity-capillary waves for small surface
tension. We begin with equation (3.2a), and make three main simplifications.
First, the expansions are truncated at $N=2$. Second, in the coefficients of
$\bar{q}^{\prime\prime}$, $\bar{q}^{\prime}$, $\bar{q}$, and $\bar{F}$, terms
of orders $B^{3}$, $B^{2}$, $B$, and $B$, respectively, are neglected. Third,
only one term in the $O(B^{2})$ component of the forcing term is retained.
This yields our model $N=2$ equation

$\displaystyle\Big{[}a\mathrm{i}Bq_{0}+a\mathrm{i}B^{2}q_{1}\Big{]}\bar{q}_{a}^{\prime\prime}+\Big{[}-F_{0}^{2}q_{0}^{2}+B\big{(}q_{0}\theta_{0}^{\prime}-a\mathrm{i}q_{0}^{\prime}-2F_{0}^{2}q_{0}q_{1}-2F_{0}F_{1}q_{0}^{2}\big{)}\Big{]}\bar{q}^{\prime}_{a}$
(3.4a)
$\displaystyle-\Big{[}2F_{0}^{2}q_{0}q_{0}^{\prime}+a\mathrm{i}q_{0}^{-1}\cos{(\theta_{0})}\Big{]}\bar{q}_{a}-\Big{[}2F_{0}q_{0}^{2}q_{0}^{\prime}\Big{]}\bar{F}=-q_{0}^{2}\theta_{1}^{\prime\prime}B^{2},$
with solution $\bar{q}_{a}(f)$ and eigenvalue $\bar{F}$. To emphasise the
dependence of $\bar{q}$ on the direction of analytic continuation $a=\pm 1$,
we have introduced above the notation $\bar{q}_{a}$. We note that while (3.4a)
is linear, the coefficients contain $q_{0}$, $\theta_{0}$, $F_{0}$, $q_{1}$,
$\theta_{1}$, and $F_{1}$, which are solutions to the nonlinear equations
found at $O(1)$ and $O(B)$ of (2.6a) and (2.6b). The advantage of this model
is that these forcing terms only need to be determined once, numerically for
instance, and when known, the effect of small surface tension may be
investigated by repeatedly solving the linear equation (3.4a).

Furthermore, in equation (3.4a) we consider $\bar{F}$ to be an eigenvalue.
This is determined by enforcing an amplitude condition for $\bar{q}$, the
remainder energy equation, which is given by

$\displaystyle\int_{-1/2}^{1/2}\bigg{[}\bar{F}\frac{F_{0}(q_{0}^{2}-1)\big{(}2q_{0}+[q_{0}^{2}-3]\cos{\theta_{0}}\big{)}}{q_{0}}+\big{(}\bar{q}_{1}+\bar{q}_{-1}\big{)}\frac{F_{0}^{2}\big{(}4q_{0}^{3}+[3q_{0}^{4}-4q_{0}^{2}-3]\cos{\theta_{0}}\big{)}}{4q_{0}^{2}}$
(3.4b)
$\displaystyle+\big{(}\bar{q}_{1}-\bar{q}_{-1}\big{)}\frac{\mathrm{i}F_{0}^{2}(1-q_{0}^{2})(q_{0}^{2}-3)\sin{\theta_{0}}}{4q_{0}^{2}}+\big{(}\bar{q}^{\prime}_{1}-\bar{q}^{\prime}_{-1}\big{)}\frac{\mathrm{i}B\big{(}q_{0}+[q_{0}^{2}-2]\cos{\theta_{0}}\big{)}}{q_{0}}\bigg{]}\mathrm{d}\phi$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad=B^{2}\int_{-1/2}^{1/2}\Big{(}q_{0}+[q_{0}^{2}-2]\cos{\theta_{0}}\Big{)}\theta_{1}^{\prime}\mathrm{d}\phi.$

Condition (3.4b) is derived in Appendix A through consideration of the full
energy equation (2.3). We note that while equation (3.4a) yields solutions in
either the upper- or lower-half $f$-plane, the amplitude condition (3.4b) is
evaluated along the real $\phi$-axis where $\text{Im}[f]=0$. Thus, the unknown
functions in (3.4b), $\bar{q}_{1}$ and $\bar{q}_{-1}$, which denote solutions
of (3.4a) with $a=1$ or $a=-1$, are evaluated upon $\text{Im}[f]=0$.

##  4 Numerical results

We now calculate numerical solutions to equation (3.4a) subject to amplitude
condition (3.4b), both of which were developed in §3.2. Solutions to this
model system, which was justified through asymptotic arguments, are intended
to contain the parasitic capillary ripples found in solutions to the fully
nonlinear gravity-capillary wave problem. The purpose of this section is to
calculate the exponentially-small scaling of the ripple amplitude predicted by
our model in order to validate our later exponential asymptotic work in §5,
and to display the limitations of the model $N=2$ equation. The main disparity
between our model linear system (3.4) and the nonlinear gravity-capillary
equations is shown to be that only solutions corresponding to a typical
asymptotic expansion as $B\to 0$ are determined correctly. Our model $N=2$
equation also produces solutions which violate our initial assumption of
$\bar{q}\ll 1$; these are multiple-scales solutions which do not agree with
those of the fully nonlinear equations.

###  4.1 Numerical methodology

We employ a spectral method to evaluate equations (3.4a) and (3.4b) on the
real-valued periodic domain $-1/2\leq\phi\leq 1/2$, and solutions are found by
Newton iteration. The computer code used to solve this problem is provided in
Appendix C. Firstly, the leading- and first-order solutions, $q_{0}$,
$\theta_{0}$, $F_{0}$, $q_{1}$, $\theta_{1}$, and $F_{1}$ must be found for a
specified energy, $\mathscr{E}$. These are solutions to nonlinear integro-
differential equations, which may be evaluated in Fourier space by using the
Fourier multipliers for differentiation, $2\pi\mathrm{i}k$, and the Hilbert
transform, $\mathrm{i}\cdot\text{sgn}(k)$, where $k$ is the wavenumber. The
solutions found in this section have an energy of $\mathscr{E}=0.4$.

We then find $\bar{q}_{1}$, $\bar{q}_{-1}$, and $\bar{F}$ as solutions to
system (3.4), which includes our model ODE (3.4a) and the energy condition
(3.4b). These are solved numerically for a given value of $B$. The domain is
discretised into $N$ collocation points given by $\phi_{j}=-1/2+(j-1)/N$ for
$j=1,\ldots,N$, upon which we evaluate the spectral representation of the
solution to find values for $\bar{q}_{a}$, $\bar{q}^{\prime}_{a}$, and
$\bar{q}^{\prime\prime}_{a}$. For specified initial guesses of $\bar{q}_{1}$,
$\bar{q}_{-1}$, and $\bar{F}$ (each either zero, or a previous numerical
solution found with a different value of $B$), we calculate the energy (3.4b)
and evaluate equation (3.4a) at $\phi_{j}$, both for $a=1$ and $a=-1$. This
yields $2N+1$ unknowns for $2N+1$ equations, specified by

Figure 2:  Branches of solutions to the $N=2$ model equation (3.4a) are shown
at fixed energy $\mathscr{E}=0.4$ in $(a)$ the $B$ vs. $\bar{F}$ plane; and
$(b)$ the $B$ vs. $F$ plane. In subplot $(b)$, we have taken
$F=F_{0}+BF_{1}+\bar{F}$ and the branches have been plotted for
$-0.0015<\bar{F}<0.0015$ only.

$\displaystyle\text{$2N+1$ Solutions}:$
$\displaystyle\quad[{\bar{q}_{1}(\phi_{1}),\ldots,\bar{q}_{1}(\phi_{N})},{\bar{q}_{-1}(\phi_{1}),\ldots,\bar{q}_{-1}(\phi_{N})},\bar{F}],$
(4.1) $\displaystyle\text{$2N+1$ Equations}:$
$\displaystyle\quad[\underbrace{\text{\eqref{eq:remaindereqN2} with
$a=1$}}_{\text{at $\phi_{j}$ for
$j=1,\ldots,N$}},\underbrace{\text{\eqref{eq:remaindereqN2} with
$a=-1$}}_{\text{at $\phi_{j}$ for
$j=1,\ldots,N$}},\eqref{eq:remainderEnergyeq}].$

The residual of the $2N+1$ equations is then minimised by Newton iteration.

###  4.2 Numerical solutions

We now find numerical solutions to system (3.4) using the method detailed in
§4.1. For a given value of the energy, $\mathscr{E}=0.4$, we calculate the
forcing terms $q_{0}$, $\theta_{0}$, $F_{0}$, $q_{1}$, $\theta_{1}$, and
$F_{1}$ with the code from figure 9. Then, we solve for $\bar{q}_{1}$,
$\bar{q}_{-1}$, and $\bar{F}$ with a specified Bond number, $B$. The code in
figure 8 does this for a range of values of $B$. This yields an initial coarse
search of the solution space. We then use each of these as an initial guess to
explore the solution space through numerical continuation.

The branches of solutions to the model $N=2$ system (3.4) found by numerical
continuation are shown in figure 2. These are shown in the $(B,\bar{F})$-plane
in figure 2$(a)$, and the $(B,F_{0}+BF_{1}+\bar{F})$-plane in figure 2$(b)$.
We note that in the fully nonlinear gravity-capillary wave problem, adjacent
branches connect to one another in a multiple-scales regime. This allows for
solutions with parasitic capillary ripples of different wavenumbers to be
connected to one another in the bifurcation diagram. However, adjacent
branches of the $N=2$ model in figure 2 do not connect. This is not
surprising, because it is only the solutions in the middle of each solution
branch (where $\bar{q}$ and $\bar{F}$ are small), for which the assumptions
made for the $N=2$ model are valid.

Figure 3:  A single solution branch from figure 2$(a)$ with $\mathscr{E}=0.4$
is shown in the $(B,\bar{F})$-plane. Insets $(a)$–$(f)$ show individual
solutions across this branch. These contain oscillatory parasitic ripples, for
which the wavenumber is seen to increase by one as the branch is travelled,
from $(f)$ to $(a)$, towards smaller $B$. Our exponential asymptotic results
in §5 will be valid across the center of this branch, where the ripple
amplitude is small.

For the solutions presented in this section, we first calculate
$(\bar{q}_{1}+\bar{q}_{-1})/2$, which produces a real-valued contribution
along the free surface, and then plot
$q=q_{0}+Bq_{1}+(\bar{q}_{1}+\bar{q}_{-1})/2$. Solutions across one of the
solution branches are shown in figure 3. As the surface tension decreases
across this branch, the parasitic capillary ripple gains an extra wavelength.
Solutions $(b)$ to $(e)$ in the middle of this branch are those whose ripple
amplitude is exponentially small as $B\to 0$, and are therefore correctly
predicted by the model equation. In the derivation of the $N=2$ equation, we
neglected terms of $O(\bar{q}^{2})$. Thus, solutions $(a)$ and $(f)$ in figure
3 are artifacts of our model and do not correspond to solutions found in the
original water-wave problem (2.1).

The capillary ripple amplitude, $\bar{q}(1/2)$, may be measured for solutions
along the branch in figure 3, and the minimum is seen to occur in the center
between solutions $(c)$ and $(d)$. This profile, which has $B=0.001330$, is
shown in figure 4$(b)$ alongside analogous solutions from other solution
branches. As the surface tension decreases ($B\to 0$), we see that the
wavenumber of the ripples increases and their amplitude decreases. This
numerical ripple amplitude is shown in more detail in the
$\log{[\bar{q}(1/2)]}$ vs $1/B$ plot of figure 4$(a)$.

Figure 4:  Subfigure $(a)$ shows the capillary ripple amplitude found for
solutions across the branches shown in figure 2. Fifteen numerical solutions
have been selected across each branch. Each of these branches has a solution
for which the ripple amplitude, $\bar{q}(1/2)$ is minimal, and these solutions
are shown in subfigure $(b)$.

In this plot, a straight line with gradient $\Delta$ would correspond to the
exponential scaling of $\mathrm{e}^{\Delta/B}$ as $B\to 0$. This value may be
estimated from the slope in figure 4$(a)$ to be $\Delta\approx-0.0066$, which
one may compare to the later exponential-asymptotic prediction of
$\Delta\approx-0.007718$ from equation (5.9). Furthermore, we see in figure
4$(a)$ that along each branch, there is a significant variation in the
amplitude. We show in §4.3 how this is predicted by the exponential asymptotic
theory of Shelton & Trinh (2022), and is due to a constant prefactor in the
exponentially-small solution, determined by enforcing periodicity, growing as
the location between adjacent branches is approached.

###  4.3 Comparing solutions of the $N=2$ model to the full problem

We now compare numerical solutions of the $N=2$ model equation (3.4a) to those
found from the original equations (2.1). To facilitate comparison we first
take numerical solutions of the full system, such as those in figure 1, for
which the physical values of the free surface are parameterised in terms of
the velocity potential by $x(\phi)$ and $y(\phi)$, and convert these into the
surface speed, $q(\phi)$, by the relation

$q(\phi)=\bigg{[}\bigg{(}\frac{\mathrm{d}x}{\mathrm{d}\phi}\bigg{)}^{2}+\bigg{(}\frac{\mathrm{d}y}{\mathrm{d}\phi}\bigg{)}^{2}\bigg{]}^{-1/2}.$
(4.2)

Comparison between the solution profiles, $q(\phi)$, of the full problem and
the $N=2$ equation occurs in figure 5$(a,b)$. In figure 5, we also plot the
capillary ripple amplitude, $\bar{q}(1/2)$, for numerical and asymptotic $N=2$
solutions, and numerical solutions of the full problem. To determine
$\bar{q}(1/2)$ for the numerical solutions of the full problem, such as the
grey profiles in figure 5$(a,b)$, we have evaluated $q-q_{0}-Bq_{1}$ at
$\phi=1/2$.

Figure 5:  Comparison is shown between numerical solutions of the full problem
(2.1), and numerical and analytical solutions to the $N=2$ equation (3.4a).
The main figure plots the capillary ripple amplitude, $\bar{q}(1/2)$, for
different values of the surface tension, $B$. The black dots show numerical
values for the $N=2$ equation from figure 4. Grey dots show the equivalent
amplitude for numerical solutions of the original system (2.1). The black line
is the exponential asymptotic prediction from §5 for the $N=2$ equation.
Numerical solutions to the model (black) and full system (grey) are shown in
subfigures $(a)$ and $(b)$. The $N=2$ equation is seen to underpredict the
amplitude of the parasitic capillary ripples.

The two main insights gained from figure 5 are in the comparison between
numerical solutions of the full and reduced problems, and the comparison
between numerical and exponential asymptotic predictions for the $N=2$ reduced
problem. Regarding the numerical comparisons shown in figure 5$(a,b)$, we see
that while the functional behaviour of the parasitic capillary ripples looks
similar in both systems, the $N=2$ model has substantially underpredicted
their amplitude. We derive this underprediction through a beyond-all-orders
analysis in §5, in which this constant is specified in equation (5.6) to be
smaller than that in the full problem by a factor of approximately $0.4118$.
This is the main drawback of our linear model equation. We note that since
$\log{(0.4118)}\approx-0.8872$, this amplitude disparity is visible in the
semilog plot of figure 5, and corresponds to the vertical shift between the
ripple amplitudes in the $N=2$ solutions (black) and the full solutions
(grey).

##  5 Exponential asymptotics

In this section, we perform a beyond-all-orders study of the $N=2$ model
equation (3.4a) in the limit of $B\to 0$. This serves two purposes. Firstly,
in §5.1 and §5.2 we demonstrate precisely how the parasitic capillary ripples
produced by our model differ from those of the full problem by only a
multiplicative constant. In §5.3 we investigate the additional effect of
eigenvalue divergence (in a Froude number expansion), and the associated
eigenfunction divergence that this induces in the solution, $\bar{q}(f)$. This
study of divergent eigenvalues is more tractable in the model $N=2$ equation
than in the full problem. While the eigenvalue divergence itself is of little
importance in this specific problem, the value of this analysis is in the
corresponding eigenfunction divergence which is often required to satisfy
boundary conditions at later orders of the asymptotic series.

###  5.1 An initial asymptotic expansion

We now calculate the exponentially-small component of solutions to the $N=2$
equation (3.4a), which corresponds to the parasitic ripples observed in the
numerical solutions of §4. To determine these, one must first asymptotically
expand the solutions, $\bar{q}$ and $\bar{F}$, as $B\to 0$, approximate later
orders of these expansions, and then resolve the associated Stokes phenomenon.

In substituting expansions of the form

$\bar{q}(f)=\sum_{n=2}^{\infty}B^{n}q_{n}(f)\qquad\text{and}\qquad\bar{F}=\sum_{n=2}^{\infty}B^{n}F_{n}$
(5.1,b)

into the $N=2$ equation (3.4a), we find at $O(B^{n})$ for $n=2$ and $n\geq 4$,
the equations

$F_{0}^{2}q_{0}^{2}q_{2}^{\prime}+\big{(}2F_{0}^{2}q_{0}q_{0}^{\prime}+a\mathrm{i}q_{0}^{-1}\cos{\theta_{0}}\big{)}q_{2}=q_{0}^{2}\theta_{1}^{\prime\prime}-2F_{0}q_{0}^{2}q_{0}^{\prime}F_{2},$
(5.2a) $\displaystyle
a\mathrm{i}q_{0}q_{n-1}^{\prime\prime}+a\mathrm{i}q_{1}q_{n-2}^{\prime\prime}-F_{0}^{2}q_{0}^{2}q_{n}^{\prime}+\big{(}q_{0}\theta_{0}^{\prime}-a\mathrm{i}q_{0}^{\prime}$
$\displaystyle-2F_{0}^{2}q_{0}q_{1}-2F_{0}F_{1}q_{0}^{2}\big{)}q_{n-1}^{\prime}$
(5.2b)
$\displaystyle-\big{(}2F_{0}^{2}q_{0}q_{0}^{\prime}+a\mathrm{i}q_{0}^{-1}\cos{\theta_{0}}\big{)}q_{n}=2F_{0}q_{0}^{2}q_{0}^{\prime}F_{n}.$

In (5.2a) above, we determine the solution $q_{2}$ and eigenvalue $F_{2}$ by
also enforcing the $O(B^{2})$ component of the $N=2$ energy condition (3.4b).
The $O(B^{3})$ equation has not been given above as only knowledge of $q_{n}$
as $n\to\infty$ and the singular behaviour of $q_{2}$ will be required. Rather
than solve (5.2b) with a fixed value of $n$, for which $q_{n-1}$ and $q_{n-2}$
are assumed known and $q_{n}$ unknown, we solve under the limit of
$n\to\infty$ for which a differential-difference equation emerges.

####  5.1.1 Late-order divergence

Due to singular behaviour of $q_{2}$ at $f=af^{*}$, and the fact that later
orders of the asymptotic solution are determined from repeated differentiation
of previous orders, the power of the singularity of $q_{n}$ will increase as
$n\to\infty$. Furthermore, factorial divergence will emerge due to this
differentiation of increasingly singularity behaviour. We thus posit a
factorial-over-power ansatz of the form

$q_{n}(f)\sim Q(f)\frac{\Gamma(n+\gamma)}{\chi(f)^{n+\gamma}}$ (5.3)

for this divergence. Substitution of (5.3) into the $O(\epsilon^{n})$ equation
(5.2b) then yields at leading order as $n\to\infty$ an equation for the
singulant $\chi(f)$, and at the next order an equation for the amplitude
function $Q(f)$. The constants of integration from these two equations, and
$\gamma$, will then be determined by matching this outer solution for $q_{n}$
to an inner solution at the singular points, $f=af^{*}$. As these equations
(by design of our model equation) are identical to those found by Shelton &
Trinh (2022), only the solutions will be given here. These solutions are

$\chi_{a}(f)=a\mathrm{i}F_{0}^{2}\int_{af^{*}}^{f}q_{0}\mathrm{d}f^{\prime}\quad\text{and}\quad
Q_{a}(f)=\Lambda_{a}q_{0}^{2}\mathrm{e}^{a\mathrm{i}\theta_{0}+a\mathrm{i}\int_{0}^{f}\Big{[}\frac{\cos{\theta_{0}}}{F_{0}^{2}q_{0}^{3}}-F_{0}^{2}q_{1}-2F_{0}F_{1}q_{0}\Big{]}\mathrm{d}f^{\prime}},$
(5.4,b)

where the constant of integration for $\chi$ has been specified by the
matching condition $\chi_{a}(af^{*})=0$. In (5.4) above, we have introduced
the notation $\chi_{a}$ and $Q_{a}$ to emphasise the dependence of these
quantities on the direction of analytic continuation, $a=\pm 1$.

The constants $\gamma$ and $\Lambda_{a}$ are found by matching with an inner
solution at $f=af^{*}$, which is performed in Appendix B. This yields

$\gamma=\frac{4}{5}\qquad\text{and}\qquad\Lambda_{a}=\frac{9\mathrm{i}f^{*}}{5F_{0}^{2}c_{a}^{4}}\mathrm{e}^{-\mathcal{P}(af^{*})}\bigg{(}\frac{4a\mathrm{i}F_{0}^{2}c_{a}}{5}\bigg{)}^{\gamma}\lim_{n\to\infty}\frac{\hat{q}_{n}}{\Gamma(n+\gamma+2)}.$
(5.5,b)

Here, $c_{a}$ is a constant from the singular behaviour $q_{0}\sim
c_{a}(f-af^{*})^{1/4}$, and $\mathcal{P}(af^{*})$, which is defined in
equation (A7) of Shelton & Trinh (2022), arises from taking the inner limit of
$Q_{a}(f)$ in (5.4). Additionally, $\hat{q}_{n}$ is determined from recurrence
relation (B.4). Evaluation of each component in (5.5) is difficult, so we
focus here on the difference between the value of $\Lambda_{a}$ in the $N=2$
model equation and the fully nonlinear problem. This factor, which will be
seen to correspond to the incorrect parasitic capillary ripple amplitude
predicted by our model equation, is given by

$\frac{\Lambda_{a}^{(N=2)}}{\Lambda_{a}^{(\text{full})}}=-\frac{9}{10}\lim_{n\to\infty}\frac{\hat{q}_{n-2}^{(N=2)}}{\hat{q}_{n}^{(\text{full})}}\approx
0.4118.$ (5.6)

In the above, we have taken the ratio of $\Lambda_{a}^{(N=2)}$ from (5.5) to
$\Lambda_{a}^{(\text{full})}$ from equation (6.13) of Shelton & Trinh (2022).
The superscript notation $\rvert^{(N=2)}$ and $\rvert^{(\text{full})}$ has
been introduced to distinguish similar variables between these works. The
quantity $\hat{q}_{n-2}^{(N=2)}$ is determined by evaluating recurrence
relation (B.4) numerically, and $\hat{q}_{n}^{(\text{full})}$ is found from
recurrence relation (B14) of Shelton & Trinh (2022). Since
$\lim_{n\to\infty}\hat{q}_{n-2}^{(N=2)}/\Gamma(n+\gamma)\approx 0.7393$, and
$\lim_{n\to\infty}\hat{q}_{n}^{(\text{full})}/\Gamma(n+\gamma)\approx-1.616$,
we have demonstrated in (5.6) that our model $N=2$ equation underpredicts the
amplitude of the solution divergence by a factor of $0.4118$. This is
significant, because the exponentially-small parasitic ripples determined next
in §5.2 are closely related to this divergence. We will show that this factor
also corresponds to the incorrect capillary ripple amplitude predicted by the
$N=2$ model that was observed in figure 5.

Essentially, this incorrect amplitude factor is the result of a second-order
truncation at $N=2$ in the asymptotic expansions. In order to obtain the
correct beyond-all-orders behaviour one must truncate optimally with $N\sim
O(1/B)$. One would expect this disparity to improve as more terms are retained
in the model, but this is inconvenient from a numerical perspective.

###  5.2 The exponentially-small solution

Now that the functional behaviour of the late-term divergence of the
asymptotic solutions is known, the parasitic capillary ripples may be
determined. These are found by considering a remainder to an optimally
truncated expansion, which is exponentially-small as $B\to 0$. The main
difficulty in determining this is the Stokes phenomenon, which causes the
exponentially-small component of the solution to rapidly change across Stokes
lines in the complex plane. Since this analysis is generic in the study of
beyond-all-order asymptotics, and the full details for this specific problem
were given in §7.2 of Shelton & Trinh (2022), only the result will be given
here.

The parasitic ripples are given analytically by the formula

$q_{\text{ripples}}(\phi)\sim\frac{1}{B^{\gamma}}\bigg{[}-\frac{\pi\cos{(G)}}{\sin{(G)}}-\mathrm{i}\pi+\sqrt{2\pi}\mathrm{i}\int_{-\infty}^{\frac{\arg{(\chi}_{1})\lvert\chi_{1}\rvert^{1/2}}{B}}\mathrm{e}^{-t^{2}/2}\mathrm{d}t\bigg{]}Q_{1}(\phi)\mathrm{e}^{-\chi_{1}(\phi)/B}+{c.c}.$
(5.7)

The last component in (5.7) above arises from the lower-half plane
contribution with $a=-1$, which is the complex-conjugate of the first
component when $f$ is specified to take real-values along the free surface.
The constant $G$ is given by

$G=\int_{0}^{1/2}\bigg{[}\frac{\cos{(\theta_{0})}}{F_{0}^{2}q_{0}^{3}}-F_{0}^{2}q_{1}-2F_{0}F_{1}q_{0}-\frac{F_{0}^{2}q_{0}}{B}\bigg{]}\mathrm{d}\phi.$
(5.8)

Analytically, the Stokes phenomenon is captured by the integral in (5.7). The
behaviour of this depends on the sign on $\arg{(\chi_{1})}$ in the upper limit
of integration; if this is negative, the integral yields zero, and when
positive, the integral gives $2\pi\mathrm{i}$. There is a boundary layer with
width of $O(B^{1/2})$ at $\arg{(\chi_{1})}=0$ [the Stokes line], across which
this change smoothly occurs. The constants
$-\pi\cos{(G)}/\sin{(G)}-\mathrm{i}\pi$ in (5.7) were found by enforcing
periodicity on this solution. We note that for certain values of $B$, the
constant $G$ in (5.8) is zero. This results in $q_{\text{ripples}}$ growing
without bound as these values of $B$ are approached, due to the factor of
$1/\sin{(G)}$ in the solution. The perturbation solution is invalid near these
points, and it was shown in figure 9 of Shelton & Trinh (2022) how these
correspond to the values of the Bond number between adjacent branches of
solutions in the $(B,F)$-plane.

There are a few complications when evaluating (5.7) numerically, as a means to
compare the asymptotic prediction with the numerical results of the $N=2$
equation. Since it is assumed that one already knows $q_{0}$, $\theta_{0}$,
$F_{0}$, $q_{1}$, and $F_{1}$ along the real-axis as part of this model (the
code in figure 9 calculates these), the constant $G$, and the function
$Q_{1}(\phi)$, may be evaluated straightforwardly. The difficulty is in
calculating $\chi_{1}(\phi)$, as the expression in (5.4) involves integration
of $q_{0}$ through the analytically continued domain. To evaluate this term,
we split the path of integration into a component along the imaginary axis
between $f^{*}$ and $0$, and another along the real axis between $0$ and
$\phi$ by writing

$\chi_{1}(\phi)=\underbrace{\mathrm{i}F_{0}^{2}\int_{f^{*}}^{0}q_{0}(f)\mathrm{d}f}_{\text{$\text{Re}[\chi_{1}]$;
real
constant}}+\underbrace{\mathrm{i}F_{0}^{2}\int_{0}^{\phi}q_{0}(\phi^{\prime})\mathrm{d}\phi^{\prime}}_{\text{$\text{Im}[\chi_{1}]$;
imaginary function}}.$ (5.9)

The first component of (5.9), a constant found by integrating $q_{0}$ through
the analytically continued domain, corresponds to the exponentially-small
scaling of the parasitic ripples. When $\mathscr{E}=0.4$, we may evaluate the
integral in (5.9) with $f^{*}\approx 0.07411\mathrm{i}$, which yields
$\text{Re}[\chi_{1}]\approx 0.007718$. This is the asymptotic prediction for
the gradient shown in the $\log{(\text{ripple amplitude})}$ vs $1/B$ plots in
figure 1, figure 4$(a)$, and figure 5. We have used the numerical method of
Crew & Trinh (2016) to determine values of the analytically continued Stokes
wave, $q_{0}$. Comparison between the current exponential asymptotic theory
and numerical solutions to the $N=2$ problem is shown in figure 5, where we
have used equation (5.7) to find $q_{\text{ripples}}(1/2)$.

###  5.3 The divergent eigenvalue expansion

In our formulation, we have imposed an amplitude condition (3.4b), for which
the Froude number, $\bar{F}$, is obtained as an eigenvalue of the problem.
When considering an asymptotic expansion of $\bar{q}$ for small Bond number,
$B$, it is only possible to satisfy the energetic condition (3.4b) at each
order if the Froude number is also similarly expanded. This had motivated our
earlier expansion of $\bar{F}=\sum_{n=2}^{\infty}B^{n}F_{n}$ in equation
(5.1). It is likely (and indeed turns out) that this expansion also diverges
and has an associated exponentially-small component—this is due to the
coupling with the divergent eigenfunction. In our beyond-all-order analysis of
§5.1.1 and §5.2, this feature has been neglected for the reason that the
particular solution induced by each component of the eigenvalue is subdominant
near the singularities, $f=af^{*}$, and their corresponding Stokes lines.
Below, we present the methodology required to determine the eigenvalue
divergence, and the corresponding correction to the eigenfunction. A similar
methodology was used by Shelton et al. (2023a) for the study of equatorially
travelling waves with exponentially-small (imaginary) eigenvalues, which
corresponded to the growth rate of a temporal instability.

In the context of our current water-wave problem, a divergent Froude number of
the form

$F_{n}\sim\delta\frac{\Gamma(n+\alpha)}{\Delta^{n+\alpha}},$ (5.10)

will induce an exponentially-small component that lies beyond-all-orders of
the expansion. Therefore, we expect that, in addition to asymptotic
corrections of the Froude number in powers of $B$, there is a contribution of
$F_{\text{exp}}=O(B^{-\alpha}\mathrm{e}^{-\Delta/B})$ that should be
considered. We note that the constants $\delta$ and $\Delta$ can take complex
values; analogous to (5.10) there will be multiple exponentially-small
contributions that sum so that the final (optimally truncated) $F$ is real-
valued. Our focus here will be on deriving the components of the
eigenfunction, $\bar{q}$, produced by the inclusion of the divergent
eigenvalue at $O(B^{n})$. In the late-term equation (5.2b), $F_{n}$ appears as
a forcing term, and the particular solution associated with this may be found
by substituting the divergent solutions for $q_{n}$ and $F_{n}$ from (5.3) and
(5.10). At leading-order as $n\to\infty$, one finds $\chi=\Delta$. At the next
order in $n$, the equation
$F_{0}^{2}q_{0}^{2}Q^{\prime}+(2F_{0}^{2}q_{0}q_{0}^{\prime}+a\mathrm{i}q_{0}^{-1}\cos{\theta_{0}})Q=-2\delta
F_{0}q_{0}^{2}q_{0}^{\prime}$ is obtained for the amplitude function, $Q(f)$.
Solving this equation yields the contribution to the late-order solution as

$q_{n}(f)\sim-\frac{2\delta}{F_{0}q_{0}^{2}}\exp{\bigg{(}-\int_{0}^{f}\frac{a\mathrm{i}\cos{\theta_{0}}}{F_{0}^{2}q_{0}^{3}}\mathrm{d}t}\bigg{)}\int_{0}^{f}\bigg{[}q_{0}^{2}q_{0}^{\prime}\exp{\bigg{(}-\int_{0}^{p}\frac{a\mathrm{i}\cos{\theta_{0}}}{F_{0}^{2}q_{0}^{3}}\mathrm{d}t}\bigg{)}\bigg{]}\mathrm{d}p\frac{\Gamma(n+\alpha)}{\Delta^{n+\alpha}}.$
(5.11)

Furthermore, it should be noted that the exponentially-small component of the
eigenvalue will contribute to the exponentially-small solution a component
similar to (5.11), of $O(B^{-\alpha}\mathrm{e}^{-\Delta/B})$. Since $\Delta$
here is constant, this particular solution will not be oscillatory across the
real-valued domain, $\phi$; therefore it does not form part of the spatially-
varying parasitic capillary ripple solution.

##  6 Discussion

In this work, we have developed a linear model equation to capture the
parasitic capillary ripples present in an inviscid travelling gravity-
capillary wave. The distinction between our reduced formulation (the $N=2$
equation (3.4a)) and previous idealised models by _e.g._ Longuet-Higgins
(1995); Crapper (1957) for gravity-capillary waves is that, in ours, the
ripples are predicted correctly when compared to solutions of the fully
nonlinear problem. This is by design; in essence we have considered the
exponential asymptotic analysis of the full problem and ensured that key
contributions are included in the $N=2$ model equation. As noted in the
introduction, similar model equations have previously been developed by
others, though for other contexts in exponential asymptotics and fluid
mechanics (Tulin, 1984; Tuck, 1991; Tu, 1991; Trinh, 2016, 2017; Jamshidi &
Trinh, 2020; Kataoka & Akylas, 2022).

To validate our methodology, we have considered one of the simplest settings
for parasitic ripples: a steadily travelling inviscid, irrotational, and
incompressible water-wave on infinite depth.

Figure 6:  Numerical solutions of different mathematical surface wave
formulations are shown. Each of these surface profiles contains oscillatory
ripples which are induced by the inclusion of a small amount of surface
tension. The formulation considered in this paper, shown in subfigure $(a)$,
is an inviscid, irrotational, and incompressible water wave on infinite depth
which includes the effects of gravity and surface tension.

However, parasitic capillary waves are not merely limited to our current
framework; they are ubiquitous in surface-wave formulations when a small
amount of interfacial tension is included (cf. the review by Perlin & Shultz
2000). We have illustrated a number of these alternative formulations in
figure 6. Examples include the viscous formulations by Mui & Dommermuth
(1995), Fedorov & Melville (1998), and Xu & Perlin (2023), where the parasitic
ripples are predominately located on the forward face of the travelling wave
[fig. 6$(b)$]; time-evolution of viscous surface waves by Hung & Tsai (2009)
[fig. 6$(c)$]; time-dependant standing gravity-capillary waves by Shelton et
al. (2023b) [fig. 6$(d)$]; and solitary waves on finite depth by Champneys et
al. (2002) [fig. 6$(g)$]. Similar phenomena have been numerically observed in
exotic surface profiles where there is a large volume of fluid near the crest,
such as in the context of the gravity-capillary solitary waves with constant
vorticity studied by Kang & Vanden-Broeck (2000) [fig. 6$(e)$]; and large-
amplitude gravity-capillary waves with trapped bubbles by Debiane & Kharif
(1996) [fig. 6$(f)$].

We are excited of the prospect that this present work will inspire the
development of further minimal models for alternative formulations that
involve some of the physically-relevant effects mentioned above. Generically,
it would be expected that, in many of these systems, the ripples will continue
to be exponentially suppressed in the Bond number and hence governed by
exponential asymptotics.

## Acknowledgments

The authors thank Prof. Paul Milewski (Bath, Penn State) for many useful
discussions concerning this problem. We also thank Epsie Trinh (Bath) for many
productive and illuminating discussions.

## Funding

JS and PHT acknowledge support by the Engineering and Physical Sciences
Research Council (EPSRC grant no. EP/V012479/1). JS is additionally supported
by the Engineering and Physical Sciences Research Council (EPSRC grant no.
EP/W522491/1).

## References

  * Akylas & Yang (1995) Akylas, T. R. & Yang, T.-S. 1995 On short-scale oscillatory tails of long-wave disturbances. Stud. Appl. Math. 94 (1), 1–20.
  * Berry (1989) Berry, M. V. 1989 Uniform asymptotic smoothing of Stokes discontinuities. Proc. Roy. Soc. London A 422, 7–21.
  * Champneys et al. (2002) Champneys, A. R., Vanden-Broeck, J-M. & Lord, G. J. 2002 Do true elevation gravity-capillary solitary waves exist? A numerical investigation. J. Fluid Mech. 454, 403–417.
  * Chen & Saffman (1980) Chen, B. & Saffman, P. G. 1980 Steady gravity-capillary waves on deep water—II. Numerical results for finite amplitude. Stud. Appl. Math. 62 (2), 95–111.
  * Cokelet (1977) Cokelet, E. D. 1977 Steep gravity waves in water of arbitrary uniform depth. Philos. Trans. R. Soc. Lond. A. 286 (1335), 183–230.
  * Crapper (1957) Crapper, G. D. 1957 An exact solution for progressive capillary waves of arbitrary amplitude. J. Fluid Mech. 2 (6), 532–540.
  * Crapper (1970) Crapper, G. D. 1970 Non-linear capillary waves generated by steep gravity waves. J. Fluid Mech. 40 (1), 149–159.
  * Crew & Trinh (2016) Crew, S. C. & Trinh, P. H. 2016 New singularities for Stokes waves. J. Fluid Mech. 798, 256–283.
  * Debiane & Kharif (1996) Debiane, M. & Kharif, C. 1996 A new limiting form for steady periodic gravity waves with surface tension on deep water. Phys. Fluids 8 (10), 2780–2782.
  * Ebuchi et al. (1987) Ebuchi, N., Kawamura, H. & Toba, Y. 1987 Fine structure of laboratory wind-wave surfaces studied under an optical method. Boundary-Layer Met. 39, 133–151.
  * Fedorov & Melville (1998) Fedorov, A. V. & Melville, W. K. 1998 Nonlinear gravity–capillary waves with forcing and dissipation. J. Fluid Mech. 354, 1–42.
  * Grimshaw & Joshi (1995) Grimshaw, R. & Joshi, N. 1995 Weakly nonlocal solitary waves in a singularly perturbed Korteweg-de Vries equation. SIAM J. Appl. Math. 55, 124–135.
  * Hung & Tsai (2009) Hung, L.-P. & Tsai, W.-T. 2009 The formation of parasitic capillary ripples on gravity–capillary waves and the underlying vortical structures. J. Phys. Ocean. 39 (2), 263–289.
  * Jamshidi & Trinh (2020) Jamshidi, S. & Trinh, P. H. 2020 Gravity–capillary waves in reduced models for wave–structure interactions. J. Fluid Mech. 890, A18.
  * Jervis (1996) Jervis, M. T. 1996 Some effects of surface tension on water waves and water waves at a wall. PhD thesis, University of Bristol.
  * Kang & Vanden-Broeck (2000) Kang, Y. & Vanden-Broeck, J.-M. 2000 Gravity-capillary waves in the presence of constant vorticity. Euro. J. Mech. 19 (2), 253–268.
  * Kataoka & Akylas (2022) Kataoka, T. & Akylas, T.R. 2022 Nonlinear effects in steady radiating waves: An exponential asymptotics approach. Physica D 435, 133272.
  * Longuet-Higgins (1975) Longuet-Higgins, M. S. 1975 Integral properties of periodic gravity waves of finite amplitude. Proc. R. Soc. Lond. A. 342 (1629), 157–174.
  * Longuet-Higgins (1995) Longuet-Higgins, M. S. 1995 Parasitic capillary waves: a direct calculation. J. Fluid Mech. 301, 79–107.
  * Mui & Dommermuth (1995) Mui, R. C. Y. & Dommermuth, D. G. 1995 The vortical structure of parasitic capillary waves. J. Fluids Eng. 117 (3), 355–361.
  * Olde Daalhuis et al. (1995) Olde Daalhuis, A. B., Chapman, S. J., King, J. R., Ockendon, J. R. & Tew, R. H. 1995 Stokes phenomenon and matched asymptotic expansions. SIAM J. Appl. Math. 55(6), 1469–1483.
  * Perlin & Shultz (2000) Perlin, M. & Shultz, W. W. 2000 Capillary effects of surface waves. Annu. Rev. Fluid Mech. 32, 241–274.
  * Pomeau et al. (1988) Pomeau, Y., Ramani, A. & Grammaticos, B. 1988 Structural stability of the Korteweg-de Vries solitons under a singular perturbation. Physica D 31, 127–134.
  * Schultz et al. (1998) Schultz, W. W., Vanden-Broeck, J.-M., Jiang, L. & Perlin, M. 1998 Highly nonlinear standing water waves with small capillary effect. J. Fluid Mech. 369, 253–272.
  * Schwartz & Vanden-Broeck (1979) Schwartz, L. W. & Vanden-Broeck, J.-M. 1979 Numerical solution of the exact equations for capillary-gravity waves. J. Fluid. Mech. 95, 119–139.
  * Shelton et al. (2023a) Shelton, J., Chapman, S. J. & Trinh, P. H. 2023a Pathological exponential asymptotics for a model problem of an equatorially trapped rossby wave. Submitted .
  * Shelton et al. (2021) Shelton, J., Milewski, P. & Trinh, P. H. 2021 On the structure of steady parasitic gravity-capillary waves in the small surface tension limit. J. Fluid Mech. 922.
  * Shelton et al. (2023b) Shelton, J., Milewski, P. & Trinh, P. H. 2023b On the structure of parasitic gravity-capillary standing waves in the small surface tension limit. Submitted .
  * Shelton & Trinh (2022) Shelton, J. & Trinh, P. H. 2022 Exponential asymptotics for steady parasitic capillary ripples on steep gravity waves. J. Fluid Mech. 939, A17.
  * Shimizu & Shōji (2012) Shimizu, C. & Shōji, M. 2012 Appearance and disappearance of non-symmetric progressive capillary–gravity waves of deep water. Jpn. J. Ind. Appl. Math. 29, 331–353.
  * Stokes (1847) Stokes, G. G. 1847 On the theory of oscillatory waves. Trans. Camb. Philos. Soc. 8, 441–455.
  * Trinh (2016) Trinh, P. H. 2016 A topological study of gravity free-surface waves generated by bluff bodies using the method of steepest descents. Proc. Roy. Soc. A 472 (2191), 20150833.
  * Trinh (2017) Trinh, P. H. 2017 On reduced models for gravity waves generated by moving bodies. J. Fluid Mech. 813, 824–859.
  * Tu (1991) Tu, Y. 1991 Saffman-Taylor problem in sector geometry: Solution and selection. Phys. Rev. A 44 (2), 1203.
  * Tuck (1991) Tuck, E. O. 1991 Ship-hydrodynamic free-surface problems without waves. J. Ship Res. 35 (4), 277–287.
  * Tulin (1984) Tulin, M. P. 1984 Surface waves from the ray point of view. In Proc. 14th Symp. Naval. Hydr., pp. 9–19. Hamburg, Germany: National Academy Press.
  * Xu & Perlin (2023) Xu, C. & Perlin, M. 2023 Parasitic waves and micro-breaking on highly nonlinear gravity–capillary waves in a convergent channel. J. Fluid Mech. 962, A46.

##  Appendix A Reduced energy condition

We now motivate the energy condition (3.4b) applied to solutions of our $N=2$
model equation. Firstly we find the full energy condition for the asymptotic
remainders, $\bar{q}$, $\bar{\theta}$, and $\bar{F}$, by substituting the
truncated expansions (3.1) into the energy constraint (2.3), which yields

$\displaystyle\int_{-1/2}^{1/2}$
$\displaystyle\bigg{[}\bigg{(}\frac{4F_{0}^{3}(q_{0}^{2}-1)(2q_{0}+[q_{0}^{2}-3]\cos{\theta_{0}})}{8q_{0}}+\cdots\bigg{)}\bar{F}+\bigg{(}\frac{F_{0}^{4}(4q_{0}^{3}+[3q_{0}^{4}-4q_{0}^{2}-3]\cos{\theta_{0}})}{8q_{0}^{2}}+\cdots\bigg{)}\bar{q}$
(A.1)
$\displaystyle+\bigg{(}\frac{F_{0}^{4}(1-q_{0}^{2})(q_{0}^{2}-3)\sin{\theta_{0}}}{8q_{0}}+\cdots\bigg{)}\bar{\theta}-\bigg{(}\frac{BF_{0}^{2}(q_{0}+[q_{0}^{2}-2]\cos{\theta_{0}})}{2}+\cdots\bigg{)}\bar{\theta}^{\prime}\bigg{]}\mathrm{d}\phi=-\xi_{\text{en}}.$

In the above, we have neglected terms of $O(\bar{q}^{2})$ and
$O(\hat{\mathscr{H}}[\bar{\theta}])$ as they are subdominant as $N\to\infty$.
Additionally, we have defined $\xi_{\text{en}}$ as the energy condition (2.3)
evaluated with the regular expansions $q_{r}$, $\theta_{r}$, and $F_{r}$,
where for instance $q_{r}=q_{0}+Bq_{1}+\cdots+B^{N-1}q_{N-1}$. Since this is
identically zero up to and including $O(B^{N-1})$, we have
$\xi_{\text{en}}=O(B^{N})$. Furthermore, due to the divergence of the late-
terms of the regular expansions as $N\to\infty$, only one component,

$\xi_{\text{en}}\sim\frac{F_{0}^{2}B^{N}}{2}\int_{-1/2}^{1/2}\big{(}-q_{0}+[2-q_{0}^{2}]\cos{\theta_{0}}\big{)}\theta_{N-1}^{\prime}\mathrm{d}\phi,$
(A.2)

of this order is dominant. The condition (3.4b) we use for our model equation
is obtained by retaining only the leading-order terms as $B\to 0$ in each
component of the integrand of (A.1) and setting $N=2$ in (A.2).

We note that to evaluate the energy condition (A.2), each solution is required
along the real $\phi$-axis (the free-surface). However, due to the analytic
continuation of the Hilbert transform, the governing equation (3.4a) for
$\bar{q}_{a}$ yields solutions in the complex $f$-plane. The analytically
continued solution in the upper-half plane, $\text{Im}[f]>0$, is found with
$a=1$, and that in the lower-half plane, $\text{Im}[f]<0$, for $a=-1$. The
functions $\bar{q}(\phi)$, $\bar{\theta}(\phi)$, and
$\bar{\theta}^{\prime}(\phi)$ required to evaluate (A.1) are then given by

$\bar{q}(\phi)=\frac{\bar{q}_{1}+\bar{q}_{-1}}{2},\qquad\bar{\theta}(\phi)=\frac{\mathrm{i}\big{(}\bar{q}_{1}-\bar{q}_{-1}\big{)}}{2q_{0}},\qquad\bar{\theta}^{\prime}(\phi)=\frac{\mathrm{i}\big{(}\bar{q}^{\prime}_{1}-\bar{q}^{\prime}_{-1}\big{)}}{2q_{0}},$
(A.3–c)

where each term on the right-hand side is evaluated on $\text{Im}[f]=0$. In
(A.3) and (A.3), we have applied the boundary-integral relation
$\bar{\theta}_{a}\sim a\mathrm{i}q_{0}^{-1}\bar{q}_{a}$ to
$\bar{\theta}(\phi)=(\bar{\theta}_{1}+\bar{\theta}_{-1})/2$ in order to
express the energy condition on $\bar{q}$ alone.

##  Appendix B Inner solution

Each order of the asymptotic solution to the $N=2$ equation is singular at the
crest singularities $f=f^{*}$ and $f=-f^{*}$. This asymptotic expansion
reorders as these points are approached, and thus an inner problem must be
introduced to resolve the solution here. We now solve this inner problem in
order to determine the two unknown constants, $\gamma$ and $\Lambda_{a}$, of
the outer divergent solution. Since $\Lambda_{a}$ corresponds to the constant
amplitude of the divergent solution and the parasitic capillary ripples, the
purpose of determining this precisely is to demonstrate how the $N=2$ model
incorrectly predicts this constant.

We now introduce the inner variable, $z$, and inner solution,
$\bar{q}_{\text{inner}}(z)$, by defining

$z=\frac{4a\mathrm{i}F_{0}^{2}c_{a}}{5B}(f-af^{*})^{5/4}\qquad\text{and}\qquad\bar{q}_{\text{outer}}=\frac{45B^{2}}{32F_{0}^{4}c_{a}}(f-af^{*})^{-9/4}\bar{q}_{\text{inner}}.$
(B.1,b)

Here, the right-hand side of (B.1) has been found by taking the inner limit of
the singulant, $\chi$, from (5.4) and dividing this by $B$. This has required
the singular behaviour $q_{0}\sim c_{a}(f-af^{*})^{1/4}$. The inner solution
in (B.1) has been defined by the singular behaviour of $q_{2}$, which is found
from the outer $O(B^{2})$ equation (5.2a) to be $q_{2}\sim
45(32F_{0}^{4}c_{a})^{-1}(f-af^{*})^{-9/4}$. This has required the singular
behaviour of
$\cos{(\theta_{0})}\sim-a\mathrm{i}F_{0}^{2}c_{a}^{3}(f-af^{*})^{-1/4}/4$,
$q_{1}\sim 3a\mathrm{i}(4F_{0}^{2})^{-1}(f-af^{*})^{-1}$, and
$\theta_{1}\sim-3(4F_{0}^{2}c_{a})^{-1}(f-af^{*})^{-5/4}$, and specifies that
$\bar{q}_{\text{inner}}\sim 1$ as $B\to 0$ and $z\to\infty$ for the inner
solution.

The inner equation is then found by substituting definitions (B.1) into the
$N=2$ equation (3.4a), and then retaining only the leading-order components as
$B\to 0$. This yields

$\bigg{(}\frac{1}{2z}-\frac{5}{6}\bigg{)}\frac{\mathrm{d^{2}}\bar{q}_{\text{inner}}}{\mathrm{d}z^{2}}+\bigg{(}-\frac{17}{10z^{2}}+\frac{23}{6z}-\frac{5}{6}\bigg{)}\frac{\mathrm{d}\bar{q}_{\text{inner}}}{\mathrm{d}z}+\bigg{(}\frac{117}{50z^{3}}-\frac{57}{10z^{2}}+\frac{1}{z}\bigg{)}\bar{q}_{\text{inner}}=\frac{1}{z},$
(B.2)

to which we consider solutions under the outer limit of $z\to\infty$ by
specifying

$\bar{q}_{\text{inner}}(z)=\sum_{n=0}^{\infty}\frac{\hat{q}_{n}}{z^{n}}.$
(B.3)

Substitution of (B.3) into the inner equation (B.2) gives the solutions
$\hat{q}_{0}=1$, $\hat{q}_{1}=171/55$, and for $n\geq 2$ the recurrence
relation

$\frac{(5n+6)}{6}\hat{q}_{n}=\frac{(5n+14)(5n+4)}{30}\hat{q}_{n-1}-\frac{(5n+3)(5n-1)}{50}\hat{q}_{n-2}.$
(B.4)

One may then determine the outer limit of the inner solution by writing (B.3)
in outer variables, which yields

$\bar{q}\sim-\frac{9c_{a}}{10}(f-af^{*})^{1/4}\sum_{n=2}^{\infty}B^{n}\frac{\hat{q}_{n-2}}{\Big{[}\frac{4a\mathrm{i}F_{0}^{2}c_{a}}{5}(f-af^{*})^{5/4}\Big{]}^{n}}.$
(B.5)

We can then match the $O(B^{n})$ component of (B.5) with the inner limit of
the outer divergent solution (5.3) to find the unknown constants $\gamma$ and
$\Lambda_{a}$. An explicit expression for $\Lambda_{a}$ is given in equation
(5.5), which is primarily used to find the difference between this constant
(which corresponds to a scaling factor of the parasitic capillary ripples) in
the $N=2$ model equation, and the corresponding constant from the full problem
determined by Shelton & Trinh (2022). This is performed in §5.1.1.

##  Appendix C Numerical codes

In this section, the MATLAB code used to find numerical solutions to the $N=2$
model equation (3.4a) is provided. The main code that implements this is given
in figure 7. First, the forcing terms must be found. These are solutions to
the $O(1)$ and $O(B)$ components of the gravity-capillary wave equations
(2.1), and are found numerically by Newton iteration with the code in figure
9. For the amplitude used here, $\mathscr{E}=0.4$, it is sufficient to begin
with an initial guess from linear theory. For values of the energy closer to
unity, one must instead take the initial guess to be a previously calculated
numerical solution with a smaller value of $\mathscr{E}$. The main code then
determines the solution $\bar{q}$ for different values of the surface tension,
$B$, using the code in figure 8. In this example, we began with an initial
guess of $\bar{q}_{a}(\phi)=0$, but to produce the branches in figure 2 a
continuation routine was used instead.

    L = 0.5; n = 1024; dphi = 2*L/n; phi = -L:dphi:L-dphi;
    h = 1i*[0 ones(1,n/2-1) 0 -ones(1,n/2-1)]; k = 1i*[0:n/2-1 0 -n/2+1:-1]*pi/L;
    % Calculates O(1) solution.
    th0 = -0.001*sin(phi); F0 = sqrt(1/2/pi);           %Guess from linear theory.
    [q0,th0,F0,err0] = N2_q0(th0,F0,0.4,n,h,k,dphi);    %Finds theta0, q0, F0.
    % Calculates O(B) solution.
    th1=zeros(1,n); F1=0;                               %Initial guess
    [q1,th1,F1,err1] = N2_q1(th1,F1,th0,F0,n,h,k,dphi); %Finds theta1, q1, F1.
    %% Finds qbar with combined method. (coarse search) %%
    Bcont = 0; Fcont = 0; Bvec = 0.001:0.000002:0.002;
    for i = 1:length(Bvec)
        B = Bvec(i); qBar1=zeros(1,n); qBar2=zeros(1,n); FBar=0;
        [qBar1,qBar2,Fbar,errb] = N2_qBar(qBar1,qBar2,FBar,q0,q1,th0,th1,B,F0,F1,n,k,dphi);
        qBar = (qBar1 + qBar2)/2; qBarR = real(qBar); Fbar = real(Fbar); Fvec(i) = F0+B*F1+Fbar;
        subplot(2,1,1);    plot(phi,q0+B*q1+qBarR)
        subplot(2,1,2);    scatter(Bvec(1:i),Fvec(1:i)); drawnow
    end

Figure 7:  MATLAB code used to calculate $\bar{q}(\phi)$ and $\bar{F}$ from
the functions in figure 8, subject to known forcing terms determined by the
functions in figure 9.

    function [qBar1,qBar2,FBar,err] = N2_qBar(qBar1,qBar2,FBar,q0,q1,th0,th1,B,F0,F1,n,k,dphi)
    [sol,ber] = fsolve(@(qBarF)N2_Disc(qBarF,q0,q1,th0,th1,F0,F1,B,n,k,dphi),...
                [qBar1 qBar2 FBar],optimoptions(’fsolve’,’Display’,’iter’));
    qBar1 = sol(1:n); qBar2 = sol(n+1:2*n); FBar = real(sol(2*n+1)); err = norm(ber,inf)^2;
    end

    function out = N2_Disc(in,q0,q1,th0,th1,F0,F1,B,n,k,dphi)
    qB1 = in(1:n); qB2 = in(n+1:2*n); FB = in(2*n+1);
    qB1phi = ifft(k.*fft(qB1)); qB1phi2 = ifft(k.^2.*fft(qB1));
    qB2phi = ifft(k.*fft(qB2)); qB2phi2 = ifft(k.^2.*fft(qB2));
    th0phi = ifft(k.*fft(th0)); q0phi = ifft(k.*fft(q0));
    th1phi = ifft(k.*fft(th1)); th1phi2=ifft(k.^2.*fft(th1));
    %% Main equations. %%
    for i=1:2; avec=[1 -1]; a=avec(i);
        T1 = -a*1i*B^2*q1-a*1i*B*q0; RHS = q0.^2.*th1phi2*B^2;
        T2 = B*((2*F0*F1*q0.^2)+(2*F0^2*q0.*q1)-(q0.*th0phi)+(a*1i*q0phi))+F0^2*q0.^2;
        T3 = (2*F0^2*q0.*q0phi)+(a*1i*cos(th0)./q0); T4 = 2*F0*q0.^2.*q0phi;
        if i==1;      bern1 = T1.*qB1phi2+T2.*qB1phi+T3.*qB1+T4*FB-RHS;
        elseif i==2;  bern2 = T1.*qB2phi2+T2.*qB2phi+T3.*qB2+T4*FB-RHS;
        end
    end
    %% Energy. %%
    qB = (qB1+qB2)/2; thBar = 1i*(qB1-qB2)./q0/2; thBphi = ifft(k.*fft(thBar));
    En1 = F0^3/2*(1-q0.^2)./q0.*((3-q0.^2).*cos(th0)-2*q0)*FB;
    En2 = F0^4/8./q0.^(2).*(4*q0.^3+(3*q0.^4-4*q0.^2-3).*cos(th0)).*qB;
    En3 = F0^4/8*(1-q0.^2)./q0.*(((q0.^2-3).*sin(th0))).*thBar;
    En4 = B*F0^2/2*((2-q0.^2).*cos(th0)-q0).*thBphi;
    EnRHS = B^2*F0^2*th1phi.*((q0.^2-2).*cos(th0)+q0);
    energy = dphi*sum(En1+En2+En3+En4-EnRHS);
    out = [bern1, bern2, energy];
    end

Figure 8:  MATLAB functions used to calculate the solution $\bar{q}(\phi)$ and
eigenvalue $\bar{F}$.

    function [q0,th0,F0,err0] = N2_q0(th0,F0,E0,n,h,k,dphi)%Finds theta0,q0
    options = optimoptions(’fsolve’,’Display’,’iter’);
    [sol,ber] = fsolve(@(initial0)N2_q0Disc(initial0,E0,h,k,n,dphi),[th0 F0],options);
    th0 = sol(1:n); F0 = real(sol(n+1)); q0 = exp(real(ifft(h.*fft(th0)))); err0 = norm(ber,inf);
    end

    function out = N2_q0Disc(in,E0,h,k,n,dphi)
    th0 = in(1:n); F0 = real(in(n+1));
    q0 = exp(real(ifft(h.*fft(th0)))); q0phi = ifft(k.*fft(q0));
    q0phi(1) = 0; q0phi(n/2+1) = 0;
    bern = F0^2*q0.^2.*q0phi + sin(th0);
    En = F0^4*(1-q0.^2)./(8*q0).*((3*cos(th0))-(2*q0)-(q0.^2.*cos(th0)));
    Energy0 = dphi*sum(En)/0.00184;
    out = [bern Energy0-E0];
    end

    function [q1,th1,F1,err1] = N2_q1(th1,F1,th0,F0,n,h,k,dphi)%Finds theta1,q1
    options = optimoptions(’fsolve’,’Display’,’iter’);
    [sol,ber] = fsolve(@(initial1)N2_q1Disc(initial1,th0,F0,h,k,n,dphi),[th1 F1],options);
    q0 = exp(real(ifft(h.*fft(th0))));
    th1 = sol(1:n); F1 = real(sol(n+1)); q1 = q0.*real(ifft(h.*fft(th1))); err1 = norm(ber,inf);
    end

    function out = N2_q1Disc(in,th0,F0,h,k,n,dphi)
    th1 = in(1:n); F1 = real(in(n+1));
    th0phi = ifft(k.*fft(th0));
    q0 = exp(real(ifft(h.*fft(th0)))); q0phi = ifft(k.*fft(q0));
    q1 = q0.*ifft(h.*fft(th1));        q1phi = ifft(k.*fft(q1));
    q1phi(1) = 0; q1phi(n/2+1) = 0;
    bern = (2*F0*F1*q0.^2.*q0phi+F0^2*(q0.^2.*q1phi+2*q0.*q1.*q0phi))+...
           th1.*cos(th0)-q0.*ifft(k.*fft(q0.*th0phi));
    en1 = 4*F0^3*F1*(1-q0.^2)./(8*q0).*((3*cos(th0))-2*q0-(q0.^2.*cos(th0)));
    en2 = -F0^4*q1/8.*(1+q0.^(-2)).*((3*cos(th0))-2*q0-(q0.^2.*cos(th0)));
    en3 = F0^4*(1-q0.^2)./(8*q0).*((-3*th1.*sin(th0))-2*q1-(2*q0.*q1.*cos(th0))+...
          (q0.^2.*th1.*sin(th0)));
    en4 = ((1-cos(th0))./q0)+((.5*F0^2.*th0phi).*((2*cos(th0))-q0-(q0.^2.*cos(th0))));
    energy1 = dphi*sum(en1+en2+en3+en4);
    out = [bern energy1];
    end

Figure 9:  MATLAB functions used to calculate the forcing terms,
$q_{0}(\phi)$, $\theta_{0}(\phi)$, $F_{0}$, $q_{1}(\phi)$, $\theta_{1}(\phi)$,
and $F_{1}$.