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

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.


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, ℰ, Froude number, , and Bond number, .As noted above, a systematic numerical study of solutions to this formulation with small surface tension,  → 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  → 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  (e −const./), and The approximate straight-line fit is indicative that this amplitude is of  (e −const./ ) as  → 0. The insets show the physical wave profile,  =  (  ), 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.
therefore their contributions lie beyond-all-orders of a traditional asymptotic expansion in powers of .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.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 " = 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 (2016Trinh ( , 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 ≈ 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.

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, , and wavespeed, .The boundary-value problem then consists of Laplace's equation for the velocity potential, (, ), within the fluid domain, kinematic and dynamic boundary conditions on the unknown free-surface,  = (), decay conditions in the deep-water limit of  → −∞, and periodicity of the solutions at  = −1/2 and  = 1/2.
As the unknown free-surface,  = (), is a streamline along which the streamfunction, (, ), is constant, we may map the physical (, )-domain, −1/2 ≤  ≤ 1/2 and −∞ <  ≤ (), to the lower-half potential (, )-plane.The free-surface may then be parameterised by the velocity potential, −1/2 ≤  ≤ 1/2, which forms the independant variable of our formulation.We note that the complex velocity may be written as  − i = e −i , where () and  () 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  and , which consists of the equations Together, equations (2.1a) and (2.1b) form a coupled set of nonlinear integro-differential equations for the solutions, () and  (), with associated constants  and .
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  − i = e log () −i , closes the system through the Hilbert transform, ℋ [].The two nondimensional constants appearing in (2.1a), the Froude number, , and the Bond number, , are given by where  is the wavespeed,  is the acceleration due to gravity,  is the wavelength,  is the coefficient of surface tension, and  is the fluid density.The intention of this work is to study nonlinear solutions of (2.1) in the small surface tension limit of  → 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.

The amplitude condition
Linear solutions to system (2.1) bifurcate from a Froude number given by  2 = (1 + 4 2  2 )/(2 ), where  is the wavenumber.Thus, rather than regard  as a fixed constant for which solutions are sought, it is convenient to introduce an amplitude parameter, ℰ, for which enforcing ℰ via an amplitude condition yields  as an eigenvalue.Linear solutions will then bifurcate from ℰ = 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  and , and the integrand is rearranged into powers of , the energy is given by (2.3) In the above, we have normalised with respect to the energy of the highest Stokes wave,  ℎ ≈ 0.00184.Other choices of amplitude could be used, such as the wave displacement, (0) − (1/2), from Schwartz & Vanden-Broeck (1979), individual Fourier coefficients from Chen & Saffman (1980), or the parameter 1 − (0) 2 (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.(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  → 0. Portions of analogous bifurcation structures, determined with other fixed amplitude measures, have been calculated by e.g.Schwartz &Vanden-Broeck (1979), andChampneys et al. (2002) for finite depth.The choice of amplitude is seen to have a significant effect on the resultant bifurcation structure.(ii) Asymptotic solutions for small surface tension.
In order to satisfy each order of the amplitude condition when considering asymptotic solutions for  and  in the limit of  → 0, it is also necessary to asymptotically expand the eigenvalue, .The th order of each asymptotic expansion,   ,   , and   , is determined by the 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,  0 ,  0 , and  0 , are found at leading-order.However, for this respective fixed value, the lower-order solutions,  1 ,  1 , and  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  → 0.

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  → 0, and thus it is necessary to consider singularities of the leading-order solution for  = 0 (a Stokes wave).However, for solutions of less than maximal amplitude, no singular behaviour is present when  ∈ R. It is only when the complex-valued domain,  ↦ →  ∈ 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  ↦ →  , where  ∈ C. Bernoulli's equation (2.1a) may be analytically continued by simply relabelling  by  .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 Im[  ] ≠ 0. Thus, when we analytically continue the Hilbert transform, we produce a residue contribution −i, which yields . (2.4) Here,  = ±1 denotes the direction of analytic continuation into either the upper-or lower-half  -plane by We have introduced in (2.4) the notation Ĥ for the complex-valued Hilbert transform, which is a function of  ∈ C, but involves integration over  ′ ∈ R. The fact that the complex-valued Hilbert transform only requires knowledge of  ( ′ ) for  ′ ∈ R along the real axis will motivate our reduced model equation, in which Ĥ is subdominant for later components of the asymptotic expansion of .
The analytically continued equations are then given by (2.6b) for the solutions (  ) and  (  ), with  ∈ C. We note that the energy condition (2.3) is enforced along the real-valued free surface, upon which Im[  ] = 0.There are two ways to interpret system (2.6).
(i) First, we may treat (2.6) as an initial value problem for which the solution along the real-axis, (),  (), , and , 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  and  specified to find (  ) and  (  ) for  along paths in C.This method was implemented numerically with  = 0 by Crew & Trinh (2016) to find branch points in the analytic continuation of gravity waves, and with  ≠ 0 by Shelton & Trinh (2022) to observe the parasitic capillary ripple amplitude increasing in the complex domain.(ii) Second, we may solve (2.6) in conjunction with energy condition (2.3), for the solutions , , and .Numerically, this requires for the implementation of a coupled scheme, in which | =1 and | =1 are solutions of (2.6) with  = 1, and | =−1 and | =−1 are solutions with  = −1 (these are coupled as they both involve the same eigenvalue ).The energy condition (2.3) is then evaluated along the real-axis with  = (| =1 + | =−1 )/2 and  = (| =1 + | =−1 )/2.This is the method we use in §4 to numerically validate solutions of our model equation.

Development of the model equation
In the small surface tension limit of  → 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 .To detect this behaviour of the solution of equations (2.6a) and (2.6b), we now consider truncated asymptotic expansions of the form In the above, the expansions are truncated at  (  −1 ), and the remainders q, θ, and F will be of  (  ).
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  =  (1/), for which  → ∞ as  → 0, the remainders in (3.1) are exponentially-small.Equations for the remainder functions q, θ, and 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 q by eliminating θ, an expression for which is found by rearranging (2.6b) and substituting this into (2.6a).This yields which is a linear second-order differential equation for q with eigenvalue F. The forcing term on the right-hand side of (3.2a) depends on the regular asymptotic series of , , and  introduced in (3.1), which we have denoted by In defining Bernoulli's equation (2.6a) and the boundary integral equation (2.6b) evaluated on these regular expansions as When solving equation (3.2a), we regard each order of  reg ,  reg , and  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  bern and  int in (3.2b) is identically zero up to and including  (  −1 ).Thus, the forcing term R is of  (  ).The exponentiallysmall components of the asymptotic solutions were derived by Shelton & Trinh (2022), who solved (3.2a) with  =  (1/), for which the forcing term R was shown to be exponentially-small in  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  = 2. Shelton & Trinh (2022) Using the methodology of exponential asymptotics, Shelton & Trinh (2022) solved equation (3.2a) under the limit of  → 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: (i)

The asymptotic theory of
Solve for the divergence of   and   in the limit of  → ∞.This is performed by approximating the  (  ) terms in equations (2.6a) and (2.6b) in the limit of  → ∞, for which a differential-difference equation emerges.The singular perturbative nature of Bernoulli's equation (2.6a) results in an  (  ) equation for which   is determined from  ′ −1 , and this leads to the solution diverging as  → ∞ in the factorial-over-power manner of Here, Λ is a constant of integration and  * is the location of the closest branch point of  0 to the real-axis.The divergence of the  expansion may then be found from Truncate optimally when  =  (1/).
The base asymptotic series reorders in the limit of  → ∞ when     ∼  +1  +1 , due to the divergence of   in (3.3a).This yields the optimal truncation point  =  + | |/, where 0 ≤  < 1 ensures that  is an integer, and the singulant  = i 2 0 ∫    *  0 d  ′ is the function in the denominator of (3.3a).The forcing term, R, from (3.2b) may then be simplified further by noting that: as  → 0, R is of  (  ); as  → ∞, only one component, − 2 0  ′′  −1 , of this order is dominant due to the solution divergence; and when  ∼ | |/ this term may be expanded to find R =  (e − |  |/ ).This yields those that cross the real-valued domain.These lie along the imaginary  -axis, and intersect at the wave crest,  = 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  = 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.

The model 𝑁 = 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  = 2. Second, in the coefficients of q′′ , q′ , q, and F, terms of orders  3 ,  2 , , and , respectively, are neglected.Third, only one term in the  ( 2 ) component of the forcing term is retained.This yields our model  = 2 equation with solution q (  ) and eigenvalue F. To emphasise the dependence of q on the direction of analytic continuation  = ±1, we have introduced above the notation q .We note that while (3.4a) is linear, the coefficients contain  0 ,  0 ,  0 ,  1 ,  1 , and  1 , which are solutions to the nonlinear equations found at  (1) and  () 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 F to be an eigenvalue.This is determined by enforcing an amplitude condition for q, the remainder energy equation, which is given by (3.4b) 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  -plane, the amplitude condition (3.4b) is evaluated along the real -axis where Im[  ] = 0. Thus, the unknown functions in (3.4b), q1 and q−1 , which denote solutions of (3.4a) with  = 1 or  = −1, are evaluated upon Im[  ] = 0.

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  = 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  → 0 are determined correctly.Our model  = 2 equation also produces solutions which violate our initial assumption of q ≪ 1; these are multiple-scales solutions which do not agree with those of the fully nonlinear equations.

Numerical methodology
We employ a spectral method to evaluate equations (3.4a) and (3.4b) on the real-valued periodic domain −1/2 ≤  ≤ 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,  0 ,  0 ,  0 ,  1 ,  1 , and  1 must be found for a specified energy, ℰ.These are solutions to nonlinear integro-differential equations, which may be evaluated in Fourier space by using the Fourier multipliers for differentiation, 2i, and the Hilbert transform, i • sgn(), where  is the wavenumber.The solutions found in this section have an energy of ℰ = 0.4.
We then find q1 , q−1 , and 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 .The domain is discretised into  collocation points given by   = −1/2 + (  − 1)/ for  = 1, . . ., , upon which we evaluate the spectral representation of the solution to find values for q , q′  , and q′′  .For specified initial guesses of q1 , q−1 , and F (each either zero, or a previous numerical solution found with a different value of ), we calculate the energy (3.4b) and evaluate equation (3.4a) at   , both for  = 1 and  = −1.This yields 2 + 1 unknowns for 2 + 1 equations, specified by 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 (  ) to (), towards smaller .Our exponential asymptotic results in §5 will be valid across the center of this branch, where the ripple amplitude is small.

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, ℰ = 0.4, we calculate the forcing terms  0 ,  0 ,  0 ,  1 ,  1 , and  1 with the code from figure C.3.Then, we solve for q1 , q−1 , and F with a specified Bond number, .The code in figure C.2 does this for a range of values of .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  = 2 system (3.4) found by numerical continuation are shown in figure 2. These are shown in the (, F)-plane in figure 2(), and the (,  0 +  1 + F)-plane in figure 2().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  = 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 q and F are small), for which the assumptions made for the  = 2 model are valid.
For the solutions presented in this section, we first calculate ( q1 + q−1 )/2, which produces a realvalued contribution along the free surface, and then plot  =  0 +  1 + ( q1 + 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 () to () in the middle of this branch are those whose ripple amplitude is exponentially small as  → 0, and are therefore correctly predicted by the model equation.In the derivation of the  = 2 equation, we neglected terms of  ( q2 ).Thus, solutions () and (  ) 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, 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 () and ().This profile, which has  = 0.001330, is shown in figure 4() alongside analogous solutions from other solution branches.As the surface tension decreases ( → 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 [ q(1/2)] vs 1/ plot of figure 4().In this plot, a straight line with gradient Δ would correspond to the exponential scaling of e Δ/ as  → 0. This value may be estimated from the slope in figure 4() to be Δ ≈ −0.0066, which one may compare to the later exponential-asymptotic prediction of Δ ≈ −0.007718 from equation (5.9).Furthermore, we see in figure 4() 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.

Comparing solutions of the 𝑁 = 2 model to the full problem
We now compare numerical solutions of the  = 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 () and (), and convert these into the surface speed, (), by the relation Comparison between the solution profiles, (), of the full problem and the  = 2 equation occurs in figure 5(, ).In figure 5, we also plot the capillary ripple amplitude, q(1/2), for numerical and asymptotic  = 2 solutions, and numerical solutions of the full problem.To determine q(1/2) for the numerical solutions of the full problem, such as the grey profiles in figure 5(, ), we have evaluated 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  = 2 reduced problem.Regarding the numerical comparisons shown in figure 5(, ), we see that while the functional behaviour of the parasitic capillary ripples looks similar in both systems, the  = 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) ≈ −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  = 2 solutions (black) and the full solutions (grey).

Exponential asymptotics
In this section, we perform a beyond-all-orders study of the  = 2 model equation (3.4a) in the limit of  → 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, q(  ).This study of divergent eigenvalues is more tractable in the model  = 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.

An initial asymptotic expansion
We now calculate the exponentially-small component of solutions to the  = 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, q and F, as  → 0, approximate later orders of these expansions, and then resolve the associated Stokes phenomenon.
In substituting expansions of the form into the  = 2 equation (3.4a), we find at  (  ) for  = 2 and  ≥ 4, the equations (5.2b) In (5.2a) above, we determine the solution  2 and eigenvalue  2 by also enforcing the  ( 2 ) component of the  = 2 energy condition (3.4b).The  ( 3 ) equation has not been given above as only knowledge of   as  → ∞ and the singular behaviour of  2 will be required.Rather than solve (5.2b) with a fixed value of , for which  −1 and  −2 are assumed known and   unknown, we solve under the limit of  → ∞ for which a differential-difference equation emerges.
5.1.1.Late-order divergence Due to singular behaviour of  2 at  =   * , and the fact that later orders of the asymptotic solution are determined from repeated differentiation of previous orders, the power of the singularity of   will increase as  → ∞.Furthermore, factorial divergence will emerge due to this differentiation of increasingly singularity behaviour.We thus posit a factorial-over-power ansatz of the form (5.3) for this divergence.Substitution of (5.3) into the  (  ) equation (5.2b) then yields at leading order as  → ∞ an equation for the singulant (  ), and at the next order an equation for the amplitude function (  ).The constants of integration from these two equations, and , will then be determined by matching this outer solution for   to an inner solution at the singular points,  =   * .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 where the constant of integration for  has been specified by the matching condition   (  * ) = 0.In (5.4) above, we have introduced the notation   and   to emphasise the dependence of these quantities on the direction of analytic continuation,  = ±1.
The constants  and Λ  are found by matching with an inner solution at  =   * , which is performed in Appendix B. This yields Here,   is a constant from the singular behaviour  0 ∼   (  −   * ) 1/4 , and P (  * ), which is defined in equation (A7) of Shelton & Trinh (2022), arises from taking the inner limit of   (  ) in (5.4b).Additionally, q is determined from recurrence relation (B.4).Evaluation of each component in (5.5b) is difficult, so we focus here on the difference between the value of Λ  in the  = 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 (5.6) In the above, we have taken the ratio of Λ (  =2)  from (5.5b) to Λ (full)  from equation (6.13) of Shelton & Trinh (2022).The superscript notation | (  =2) and | (full) has been introduced to distinguish similar variables between these works.The quantity q(=2) −2 is determined by evaluating recurrence relation (B.4) numerically, and q(full)  is found from recurrence relation (B14) of Shelton & Trinh (2022).Since lim →∞ q(=2) −2 /Γ( + ) ≈ 0.7393, and lim →∞ q(full)  /Γ( + ) ≈ −1.616, we have demonstrated in (5.6) that our model  = 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  = 2 model that was observed in figure 5.
Essentially, this incorrect amplitude factor is the result of a second-order truncation at  = 2 in the asymptotic expansions.In order to obtain the correct beyond-all-orders behaviour one must truncate optimally with  ∼  (1/).One would expect this disparity to improve as more terms are retained in the model, but this is inconvenient from a numerical perspective.

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  → 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 e − 2 /2 d  1 ()e −  1 ( )/ + .. (5.7) The last component in (5.7) above arises from the lower-half plane contribution with  = −1, which is the complex-conjugate of the first component when  is specified to take real-values along the free surface.
The constant  is given by (5.8) Analytically, the Stokes phenomenon is captured by the integral in (5.7).The behaviour of this depends on the sign on arg (  1 ) in the upper limit of integration; if this is negative, the integral yields zero, and when positive, the integral gives 2i.There is a boundary layer with width of  ( 1/2 ) at arg (  1 ) = 0 [the Stokes line], across which this change smoothly occurs.The constants − cos ()/sin () − i in (5.7) were found by enforcing periodicity on this solution.We note that for certain values of , the constant  in (5.8) is zero.This results in  ripples growing without bound as these values of  are approached, due to the factor of 1/sin () 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 (, )-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  = 2 equation.Since it is assumed that one already knows  0 ,  0 ,  0 ,  1 , and  1 along the real-axis as part of this model (the code in figure C.3 calculates these), the constant , and the function  1 (), may be evaluated straightforwardly.The difficulty is in calculating  1 (), as the expression in (5.4a) involves integration of  0 through the analytically continued domain.To evaluate this term, we split the path of integration into a component along the imaginary axis between  * and 0, and another along the real axis between 0 and  by writing (5.9) The first component of (5.9), a constant found by integrating  0 through the analytically continued domain, corresponds to the exponentially-small scaling of the parasitic ripples.When ℰ = 0.4, we may evaluate the integral in (5.9) with  * ≈ 0.07411i, which yields Re[  1 ] ≈ 0.007718.This is the asymptotic prediction for the gradient shown in the log (ripple amplitude) vs 1/ plots in figure 1, figure 4(), and figure 5. We have used the numerical method of Crew & Trinh (2016) to determine values of the analytically continued Stokes wave,  0 .Comparison between the current exponential asymptotic theory and numerical solutions to the  = 2 problem is shown in figure 5, where we have used equation (5.7) to find  ripples (1/2).

The divergent eigenvalue expansion
In our formulation, we have imposed an amplitude condition (3.4b), for which the Froude number, F, is obtained as an eigenvalue of the problem.When considering an asymptotic expansion of q for small Bond number, , 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 F = ∞ =2     in equation (5.1b).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,  =   * , 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 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 , there is a contribution of  exp =  ( −  e −Δ/ ) that should be considered.We note that the constants  and Δ can take complex values; analogous to (5.10) there will be multiple exponentially-small contributions that sum so that the final (optimally truncated)  is real-valued.Our focus here will be on deriving the components of the eigenfunction, q, produced by the inclusion of the divergent eigenvalue at  (  ).In the late-term equation (5.2b),   appears as a forcing term, and the particular solution associated with this may be found by substituting the divergent solutions for   and   from (5.3) and (5.10).At leading-order as  → ∞, one finds  = Δ.At the next order in , the equation is obtained for the amplitude function, (  ).Solving this equation yields the contribution to the late-order solution as (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  ( −  e −Δ/ ).Since Δ here is constant, this particular solution will not be oscillatory across the real-valued domain, ; therefore it does not form part of the spatially-varying parasitic capillary ripple solution.

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  = 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  = 2 model equation.
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.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(𝑏)]; time-evolution of viscous surface waves by Hung & Tsai (2009) [fig. 6(𝑐)]; time-dependant standing gravity-capillary waves by Shelton et al. (2023b) [fig. 6(𝑑)]; and solitary waves on finite depth by Champneys et al. (2002) [fig. 6(𝑔)].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(𝑒)]; and large-amplitude gravity-capillary waves with trapped bubbles by Debiane & Kharif (1996) [fig. 6( 𝑓 )].
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.
The inner equation is then found by substituting definitions (B.1) into the  = 2 equation (3.4a), and then retaining only the leading-order components as  → 0. This yields We can then match the  (  ) component of (B.5) with the inner limit of the outer divergent solution (5.3) to find the unknown constants  and Λ  .An explicit expression for Λ  is given in equation (5.5b), which is primarily used to find the difference between this constant (which corresponds to a scaling factor of the parasitic capillary ripples) in the  = 2 model equation, and the corresponding constant from the full problem determined by Shelton & Trinh (2022).This is performed in §5.1.1.

C. Numerical codes
In this section, the MATLAB code used to find numerical solutions to the  = 2 model equation (3.4a) is provided.The main code that implements this is given in figure C.1.First, the forcing terms must be found.These are solutions to the  (1) and  () components of the gravity-capillary wave equations (2.1), and are found numerically by Newton iteration with the code in figure C.3.For the amplitude used here, ℰ = 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 ℰ.The main code then determines the solution q for different values of the surface tension, , using the code in figure C.2.In this example, we began with an initial guess of q () = 0, but to produce the branches in figure 2 a continuation routine was used instead.

Fig. 1 .
Fig. 1.Numerical solutions of the full gravity-capillary wave problem are shown at fixed energy ℰ = 0.4 for different values of the surface-tension (or Bond) parameter, .The main plot shows the logarithm of the capillary-ripple amplitude as a function of 1/.The approximate straight-line fit is indicative that this amplitude is of  (e −const./ ) as  → 0. The insets show the physical wave profile,  =  (  ), at two selected points of the main diagram.These solutions are computed using the algorithm discussed inShelton et al. (2021).A comparison between these amplitudes, and those predicted both numerically and asymptotically by our reduced model occurs in figure5.
Fig.3.A single solution branch from figure 2() with ℰ = 0.4 is shown in the ( , F )-plane.Insets ()-(  ) 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 (  ) to (), towards smaller .Our exponential asymptotic results in §5 will be valid across the center of this branch, where the ripple amplitude is small.

Fig. 4 .
Fig.4.Subfigure () shows the capillary ripple amplitude found for solutions across the branches shown in figure2.Fifteen numerical solutions have been selected across each branch.Each of these branches has a solution for which the ripple amplitude, q (1/2) is minimal, and these solutions are shown in subfigure ().

Fig. 5 .
Fig. 5. Comparison is shown between numerical solutions of the full problem (2.1), and numerical and analytical solutions to the  = 2 equation (3.4a).The main figure plots the capillary ripple amplitude, q (1/2), for different values of the surface tension, .The black dots show numerical values for the  = 2 equation from figure4.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  = 2 equation.Numerical solutions to the model (black) and full system (grey) are shown in subfigures () and ().The  = 2 equation is seen to underpredict the amplitude of the parasitic capillary ripples.

Fig. 6 .
Fig.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 (), is an inviscid, irrotational, and incompressible water wave on infinite depth which includes the effects of gravity and surface tension.
consider solutions under the outer limit of  → ∞ by specifying qinner () = of (B.3) into the inner equation (B.2) gives the solutions q0 = 1, q1 = 171/55, and for  ≥ 2 determine the outer limit of the inner solution by writing (B.3) in outer variables, which yields q ∼ − 9  10 (  −   * ) Capture the Stokes phenomenon by localising in a boundary layer about the Stokes lines, where Im[ ] = 0 and Re[ ] ≥ 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