Quasi-periodic traveling gravity-capillary waves

We present a numerical study of spatially quasi-periodic traveling waves on the surface of an ideal fluid of infinite depth. This is a generalization of the classic Wilton ripple problem to the case when the ratio of wave numbers satisfying the dispersion relation is irrational. We develop a conformal mapping formulation of the water wave equations that employs a quasi-periodic variant of the Hilbert transform to compute the normal velocity of the fluid from its velocity potential on the free surface. We develop a Fourier pseudo-spectral discretization of the traveling water wave equations in which one-dimensional quasi-periodic functions are represented by two-dimensional periodic functions on the torus. This leads to an overdetermined nonlinear least squares problem that we solve using a variant of the Levenberg-Marquardt method. We investigate various properties of quasi-periodic traveling waves, including Fourier resonances and the dependence of wave speed and surface tension on the amplitude parameters that describe a two-parameter family of waves.


I
Traveling water waves have long played a central role in the field of fluid mechanics. Spatially periodic traveling waves, dating back to Stokes [11,33], have been studied extensively [6,[21][22][23][24]29,30,34]. However, little research has been done on spatially quasi-periodic water waves in spite of their abundance in integrable model water wave equations such as the Korteweg-de Vries equation and the nonlinear Schrödinger equation. On the theoretical side, Bridges and Dias [8] used a spatial Hamiltonian structure to construct weakly nonlinear approximations of spatially quasi-periodic traveling gravity-capillary waves for two special cases: deep water and shallow water. The existence of such waves in the fully nonlinear setting is still an open problem. In this paper, we formulate the quasi-periodic traveling wave problem in a conformal mapping framework, demonstrate their existence numerically, and explore their properties.
To motivate our work, recall the dispersion relation for linearized traveling gravitycapillary waves in deep water: Here c is the phase speed, k is the wave number, g is the acceleration due to gravity and τ is the coefficient of surface tension. Notice that c " a pg{kq`τk has a positive minimum, denoted by c crit . For any fixed phase speed c ą c crit , there are two distinct positive wave numbers satisfying the dispersion relation (1.1), denoted k 1 and k 2 . Any traveling solution of the linearized problem with this speed can be expressed as a superposition of waves with these wave numbers. If k 1 and k 2 are rationally related, the motion is spatially periodic This work was supported in part by the National Science Foundation under award number DMS-1716560 and by the Department of Energy, Office of Science, Applied Scientific Computing Research, under award number DE-AC02-05CH11231. and corresponds to the well-known Wilton ripples [1,35,39]. However, if k 1 and k 2 are irrationally related, the motion will be spatially quasi-periodic.
Recently, Berti et al [5,7] have proved the existence of small-amplitude temporally quasiperiodic gravity-capillary waves using Nash-Moser theory. They show that solutions of the linearized standing water wave problem can be combined and perturbed to obtain temporally quasi-periodic solutions of the nonlinear problem. Following the same philosophy, we look for spatially quasi-periodic solutions of the traveling water wave equations that are perturbations of solutions of the linearized problem. The velocity potential can be eliminated from the Euler equations when looking for traveling solutions, so our goal is to study traveling waves with height functions of the form (1.2) ηpαq "ηpk 1 α, k 2 αq,ηpα 1 , α 2 q " ÿ pj 1 , j 2 qPZ 2η j 1 , j 2 e ip j 1 α 1`j2 α 2 q .
Hereη is real-valued and defined on the torus T 2 " R 2 {2πZ 2 , and α parametrizes the free surface in such a way that the fluid domain is the image of the lower half-plane tw " α`iβ : β ă 0u under a conformal map zpwq whose imaginary part on the upper boundary is Imtz| β"0 u " η. The leading term here is η lin pαq " 2 Retη 1,0 e ik 1 α`η 0,1 e ik 2 α u, which will be a solution of the linearized problem.
Unlike [8], as noted above, we use a a conformal mapping formulation [9,14,25] of the gravity-capillary water wave problem. This makes it possible to compute the normal velocity of the fluid from the velocity potential on the free surface via a quasi-periodic variant of the Hilbert transform. As in the periodic case, the Hilbert transform is a Fourier multiplier operator, but now acts on functions defined on a higher-dimensional torus. In a companion paper [37], we use this idea to develop a numerical method to compute the time evolution of solutions of the Euler equations from arbitrary quasi-periodic initial data. The present paper focuses on traveling waves in this framework.
We formulate the traveling wave computation as a nonlinear least-squares problem and use the Levenberg-Marquardt method to search for solutions. This approach builds on the overdetermined shooting methods developed by Wilkening et al [2,3,17,32,38] to compute standing waves and other time-periodic solutions. Specifically, we fix the ratio k 2 {k 1 , denoted by k, and solve simultaneously for the phase speed c, the coefficient of surface tension τ, and the unknown Fourier modesη j 1 , j 2 in (1.2) subject to the constraint thatη 1,0 andη 0,1 have prescribed values. In Section 3, we discuss the merits of these bifurcation parameters over, say, prescribing τ andη 1,0 and solving forη 0,1 along with c and the other unknown Fourier modes. While the numerical method is general and can be used to search for solutions for any irrational k, for brevity we present results only for k " 1{ ?
2. In future work we plan to extend our results to the case of finite-depth water waves and analyze the stability of solutions [12,26,35].
In Section 2, we define a quasi-periodic Hilbert transform, derive the equations of motion governing quasi-periodic traveling water waves, and summarize the main results and notation introduced in [37] on the more general spatially quasi-periodic initial value problem. In Section 3, we design a Fourier pseudo-spectral method to numerically solve the torus version of the quasi-periodic traveling wave equations. The discretization leads to an overdetermined nonlinear least-squares problem that we solve using a variant of the Levenberg-Marquardt method [31,38]. In Section 4, we present a detailed numerical study of a two-parameter family of quasi-periodic traveling waves with k " 1{ ?
2 and g " 1. In Section 5, we summarize the results and discuss the effects of floating point arithmetic and whether solutions might exist for rational values of k. Finally, in Appendix A, we recall a theorem proved in [37] establishing sufficient conditions for an analytic function zpwq to map the lower half-plane topologically onto a semi-infinite region bounded above by a parametrized curve. We also discuss conditions that ensure 1{|z w | is uniformly bounded in the lower half-plane. We then study the dynamics of traveling waves in conformal space for various choices of a free parameter in the equations of motion that controls the tangential velocity of the surface parametrization. We show that the waves maintain a permanent form but generally travel at a non-uniform speed in conformal space as they evolve.

P
The primary goal of this paper is to study quasi-periodic traveling water waves using a conformal mapping framework. In this section, we establish notation; review the properties of the quasi-periodic Hilbert transform; discuss quasi-periodic conformal maps and complex velocity potentials; propose a synthesis of viewpoints between the Hou, Lowengrub and Shelley formalism for evolving interfaces [18,19] and the conformal mapping method of Dyachenko and Zakharov [9,13,41]; summarize the one-dimensional (1d) and torus versions of the equations of motion for the spatially quasi-periodic initial value problem [37]; discuss families of 1d quasi-periodic solutions corresponding to a single solution of the torus version of the problem; derive the equations governing traveling waves; and review the linear theory of quasi-periodic traveling waves.

Quasi-periodic functions and the Hilbert transform.
A function upαq is quasi-periodic if there exists a continuous, periodic functionũpαq defined on the d-dimensional torus T d such that (2.1) upαq "ũpkαq,ũpαq " ÿ jPZ dû j e ixj, αy , α P R, α, k P R d .
We generally assumeũpαq is real analytic, which means the Fourier modes satisfy the symmetry conditionû´j "û j and decay exponentially as |j| Ñ 8, i.e. |û j | ď Me´σ |j| for some M, σ ą 0. Entries of the vector k are required to be linearly independent over Z.
Fixing this vector k, we define two versions of the Hilbert transform, one acting on u (the quasi-periodic version) and the other onũ (the torus version): Here sgnpqq P t1, 0,´1u depending on whether q ą 0, q " 0 or q ă 0, respectively. Note that the torus version of H is a Fourier multiplier on L 2 pT d q that depends on k. It is shown in [37] that Hruspαq " Hrũspkαq, and the most general bounded analytic function f pwq in the lower half-plane whose real part agrees with u on the real axis has the form 2û j e ix j,kyw , pw " α`iβ , β ď 0q, wherev 0 is an arbitrary constant and the sum is over all j P Z d satisfying xj, ky ă 0. The imaginary part of f on the real axis is then given by v "v 0´H rus. Similarly, given v, the most general bounded analytic function f pwq in the lower half-plane whose imaginary part agrees with v on the real axis has the form (2.4) with u "û 0`H rvs, whereû 0 is an arbitrary constant. This analytic extension is quasi-periodic on slices of constant depth, i.e.
The torus version of the bounded analytic extension corresponding toũpα`θq is simplỹ f pα`θ, βq, which has imaginary partṽpα`θq on the real axis. As a result, the Hilbert transform commutes with the shift operator, which can also be checked directly from (2.2). We also define quasi-periodic and torus versions of two projection operators, where P 0 rus is a constant function on R, P 0 rũs is a constant function on T d , and Prus has zero-mean on R in the sense that its torus representation, Prũs, which satisfies Pruspαq " Prũspkαq, has zero mean on T d .

A quasi-periodic conformal mapping.
For the general initial value problem [37], we consider a time-dependent conformal map zpw, tq that maps the lower half-plane to the fluid domain Ω f ptq that lies below the free surface in physical space. At each time t, we assume zpw, tq extends continuously to C´, and in fact is analytic on a slightly larger half-plane Cέ " tw : Im w ă εu, where ε ą 0 could depend on t. The free surface Γptq is parametrized by (2.9) ζpα, tq " ξpα, tq`iηpα, tq, pα P R , t fixedq, ζ " z| β"0 .
We assume α Þ Ñ ζpα, tq is injective but do not assume Γptq is the graph of a single-valued function of x; an example of a spatially quasi-periodic overturning wave is computed in [37]. The conformal map is required to remain a bounded distance from the identity map in the lower half-plane. Specifically, we require that where Mptq is a uniform bound that could vary in time. The Cauchy integral formula implies that |z w´1 | ď Mptq{|β|, so at any fixed time, In this paper and its companion [37], we assume η has two spatial quasi-periods, i.e. at any time it has the form (2.1) with d " 2 and k " rk 1 , k 2 s T . The is a major departure from previous work [13,16,27,41], where η is periodic. Through non-dimensionalization, we may assume k 1 " 1 and k 2 " k, where k is irrational: (2.12) ηpα, tq "ηpα, kα, tq,ηpα 1 , α 2 , tq " ÿ j 1 , j 2 PZη Hereη´j 1 ,´j 2 ptq "η j 1 , j 2 ptq sinceηpα 1 , α 2 , tq is real-valued. Since w Þ Ñ rzpw, tq´ws is bounded and analytic on C´and its imaginary part agrees with η on the real axis, there is a real number x 0 (possibly depending on time) such that (2.13) ξpα, tq " α`x 0 ptq`Hrηspα, tq, ξ α pα, tq " 1`Hrη α spα, tq.
Specifically,ξ " x 0 ptq`Hrηs,ζ "ξ`iη, and Since the modesη j 1 , j 2 are assumed to decay exponentially, there is a uniform bound Mptq such that |zpα 1 , α 2 , β, tq| ď Mptq for pα 1 , α 2 q P T 2 and β ď 0. In [37], we show that as long as the free surface ζpα, tq does not self-intersect at a given time t, the mapping w Þ Ñ zpw, tq is an analytic isomorphism of the lower half-plane onto the fluid region.
2.3. The complex velocity potential and equations of motion for the initial value problem. Adopting the notation of [37], let Φ phys px, y, tq denote the velocity potential in physical space and let W phys px`i y, tq " Φ phys px, y, tq`iΨ phys px, y, tq denote the complex velocity potential, where Ψ phys is the stream function. Using the conformal mapping zpw, tq, we pull back these functions to the lower half-plane and define Wpw, tq " Φpα, β, tq`iΨpα, β, tq " W phys pzpw, tq, tq, pw " α`iβq.
Eliminating 9 α p gives the kinematic condition wheren " p´η α , ξ α q{s α is the outward unit normal to Γ and we have identified ζ t with the vector pξ t , η t q in R 2 . The general philosophy proposed by Hou, Lowengrub and Shelley (HLS) [18,19] is that while (2.20) constrains the normal velocity U of the curve to match that of the fluid, the tangential velocity V can be chosen arbitrarily to improve the mathematical properties of the representation or the accuracy and stability of the numerical scheme. Whereas HLS propose choosing V to keep s α ptq independent of α, we interpret the work of Zakharov and Dyachenko [9,13,41] as choosing V to maintain a conformal representation. Briefly, since Φ phys and Ψ phys satisfy the Cauchy-Riemann equations, Assuming z t {z α is bounded and analytic in the lower half-plane (see Appendix ??), where C 1 is an arbitrary constant that we are free to choose. The tangential and normal velocities can be rotated back to obtain ξ t and η t via which can be interpreted as the real and imaginary parts of the complex multiplication ζ t " pζ α qpζ t {ζ α q. As explained in [37], the first equation of (2.23) is automatically satisfied if the second equation holds and ξ is reconstructed from η via (2.13), provided x 0 ptq satisfies The equations of motion for water waves in the conformal framework may now be written where the last equation comes from the unsteady Bernoulli equation and the Laplace-Young condition for the pressure; see [37] for details. As noted in [37], equations (2.25) can be interpreted as an evolution equation for the functionsζpα 1 , α 2 , tq andφpα 1 , α 2 , tq on the torus T 2 . The α-derivatives are replaced by the directional derivatives rB α 1`k B α 2 s and the quasi-periodic Hilbert transform is replaced by its torus version, i.e. Hrũs in (2.2) above. The pseudo-spectral method proposed in [37] is based on this representation. A convenient choice of C 1 is , which causesξp0, 0, tq to remain constant in time, alleviating the need to evolve x 0 ptq explicitly. HereJ " p1`ξ α q 2`η2 α . Note that ξ α in (2.25) is replaced by since the secular growth term α is not part ofξ in (2.14). Using (2.13) and (2.14),ζ is completely determined by x 0 ptq andη, so only these have to be evolved -the formula for ξ t in (2.23) is redundant as long as (2.24) is satisfied.
It is shown in [37] that solving the torus version of (2.25) yields a three-parameter family of one-dimensional solutions of the form The parameters pθ 1 , θ 2 , δq lead to the same solution as p0, θ 2´k θ 1 , 0q up to a spatial phase shift and α-reparametrization. Thus, every solution is equivalent to one of the form Within this smaller family, two values of θ lead to equivalent solutions if they differ by 2πpn 1 k`n 2 q for some integers n 1 and n 2 . This equivalence is due to solutions "wrapping around" the torus with a spatial shift, (2.30) ζpα`2πn 1 , t ; 0, θ, 0q " ζpα, t ; 0, θ`2πpn 1 k`n 2 q, 2πn 1 q,`α P r0, 2πq, n 1 P Z˘.
It is shown in [37] that if all the waves in the family (2.29) are single-valued and have no vertical tangent lines, there is a corresponding family of solutions of the Euler equations in a standard graph-based formulation [10,21,40] that are quasi-periodic in physical space.

Quasi-periodic traveling water waves.
We now specialize to the case of quasi-periodic traveling waves and derive the equations of motion in a conformal mapping framework. One approach (see e.g. [28] for the periodic case) is to write down the equations of motion in a graph-based representation of the surface variables η phys px, tq and ϕ phys px, tq " Φ phys px, ηpx, tq, tq and substitute η We present below an alternative derivation of the equations in [28] that is more direct and does not assume the wave profile is single-valued. Other systems of equations have also been derived to describe traveling water waves, e.g. by Nekrasov [29,30] and Dyachenko et. al. [15].
This expresses ψ and ϕ (up to additive constants) in terms of η and ξ " α`x 0`H rηs, leaving only η to be determined. As in the graph-based approach of (2.31) above, it suffices to compute the initial wave profile at t " 0 to know the full evolution of the traveling wave under (2.25); however, the wave generally travels at a non-uniform speed in conformal space in order to travel at constant speed in physical space; see Appendix A.
The two-dimensional velocity potential Φ phys px, y, tq may be assumed to exist even if the traveling wave possesses overhanging regions that cause the graph-based representation via η phys px, tq and ϕ phys px, tq to break down. In a moving frame traveling at constant speed c along with the wave, the free surface will be a streamline. Letz " z´ct denote position in the moving frame and note that the complex velocity potential picks up a background flow term,W phys pz, tq " W phys pz`ct, tq´cz, and becomes time-independent. We drop t in the notation and defineWpwq "W phys pzpwqq, wherezpwq " zpw, 0q conformally maps the lower half-plane onto the fluid region of this stationary problem. We assume W phys pzpαq, 0q is quasi-periodic with exponentially decaying mode amplitudes, so |Wpwq`cw| ď |W phys pzpwq, 0q|`c|zpwq´w| is bounded in the lower half-plane. Since the stream function ImtW phys pzqu is constant on the free surface, we may assume ImtWpαqu " 0 for α P R. The function ImtWpwq`cwu is then bounded and harmonic in the lower half-plane and satisfies homogeneous Dirichlet boundary conditions on the real line, so it is zero [4]. Up to an additive real constant, Thus, |∇Φ phys | 2 " |W 1 pwq{z 1 pwq| 2 " c 2 {J. Since the free surface is a streamline in the moving frame, the steady Bernoulli equation p1{2q|∇Φ phys | 2`g η`p{ρ " C together with the Laplace-Young condition p " p 0´ρ τκ on the pressure gives which is the desired system of equations for η.
In the quasi-periodic traveling wave problem, we seek a solution of (2.34) of the form (2.12), except thatη and its Fourier modes will not depend on time. Like the initial value problem, (2.34) can be interpreted as a nonlinear system of equations forηpα 1 , α 2 q defined on T 2 , where the α-derivatives are replaced by rB α 1`k B α 2 s and the Hilbert transform is replaced by its torus version in (2.2). Without loss of generality, we assume We also assume thatη is an even, real function of pα 1 , α 2 q on T 2 . Hence, in our setup, the Fourier modes ofη satisfy This implies that all the Fourier modesη j 1 , j 2 are real, and causes ηpαq "ηpα, kαq to be even as well, which is compatible with the symmetry of (2.34). However, as in (2.28), there is a larger family of quasi-periodic traveling solutions embedded in this solution, namely (2.37) ηpα; θq "ηpα, θ`kαq.
As in (2.30), two values of θ lead to equivalent solutions (up to α-reparametrization and a spatial phase shift) if they differ by 2πpn 1 k`n 2 q for some integers n 1 and n 2 . In general, ηpα´α 0 ; θq will not be an even function of α for any choice of α 0 unless θ " 2πpn 1 k`n 2 q for some integers n 1 and n 2 . In the periodic case, symmetry breaking traveling water waves have been found by Zufiria [42], though most of the literature is devoted to periodic traveling waves with even symmetry. Equations (2.34) were derived from the requirement that their solutions travel at a constant speed in physical space. In Appendix A, we consider their evolution in conformal space under (2.25) for various choices of C 1 . The 1d waves maintain a permanent form as functions of α that travel at a generally non-uniform speed, and the torus version of the waves maintain a permanent two-dimensional form that travels through the torus in the p1, kq direction at a speed that generally varies in time. A particular choice of C 1 causesη andφ to remain stationary in time, though it is not the choice (2.26) in whichξp0, 0, tq " 0.
2.5. Linear theory of quasi-periodic traveling waves. Linearizing (2.34) around the trivial solution ηpαq " 0, we obtain, where δη denotes the variation of η. Substituting (2.12) into (2.38), we obtain a resonance relation for the Fourier modes of δη: Note that j 1`j2 k, which appears in the exponent of (2.12), plays the role of k in the dispersion relation (1.1). In the numerical scheme, we assume that both of the base modesη 1,0 ,η 0,1 are nonzero. (If either is zero, there is another family of periodic solutions bifurcating from the quasi-periodic family of interest here.) Setting pj 1 , j 2 q to p1, 0q and p0, 1q, respectively, gives the first-order resonance conditions These are dimensionless equations, where the wave number k 1 of the first wave has been set to 1, and k 2 " k 1 k. For right-moving waves, we then have c " ? g`τ and k " g{τ.
Any superposition of waves with wave numbers 1 and k traveling with speed c will solve the linearized problem (2.38). We introduce the notation c lin " a g`g{k and τ lin " g{k to facilitate the discussion of nonlinear effects below.

N M
Equations (2.34) involve computing derivatives and Hilbert transforms of quasi-periodic functions that arise in intermediate computations. Let f pαq denote one of these functions, and letf denote the corresponding periodic function on the torus, Eachf that arises is represented by its values on a uniform M 1ˆM2 grid on the torus T 2 , Products, powers and quotients in (2.34) are evaluated pointwise on the grid while derivatives and the Hilbert transform are computed in Fourier space via We use the 'r2c' version of the 2d FFTW library to rapidly compute the forward and inverse transform given by (3.4) The FFTW library actually returns the index range 0 ď j 2 ă M 2 , but we usef j 1 , j 2´M2 "f j 1 , j 2 to de-alias the Fourier modes and map the indices j 2 ą M 2 {2 to their correct negative values. The missing entries with´M 1 {2 ă j 1 ă 0 are determined implicitly by When computing f α and Hr f s via (3.3), the Nyquist modes with are set to zero, which ensures that the 'c2r' transform reconstructs real-valued functions Ă f α and Ć Hr f s from their Fourier modes. Further details on this pseudo-spectral representation are given in [37] in the context of timestepping the dynamic equations (2.25).
In [38], an overdetermined shooting algorithm based on the Levenberg-Marquardt method [31] was proposed for computing standing water waves accurately and efficiently.
Here we adapt this method to compute quasi-periodic traveling waves instead of standing waves. We first formulate the problem in a nonlinear least-squares framework. We consider τ, c 2 (which we denote as b) and η as unknowns in (2.34) and define the residual function Here,η represents the Fourier modes of η, which are assumed real via (2.36); J and κ depend on η through the auxiliary equations of (2.34); and a tilde indicates that the function is represented on the torus, T 2 , as in (3.1). We also define the objective function Note that solving (2.34) is equivalent to finding a zero of the objective function Frτ, b,ηs.
The parameter k in (3.1) is taken to be a fixed, irrational number when searching for zeros of F.
In the numerical computation, we truncate the problem to finite dimensions by varying only the leading Fourier modesη j 1 , j 2 with 0 ď |j 1 |, | j 2 | ď N{2. We evaluate the residual R (and compute the Fourier transforms) on an MˆM grid, where M ą N. The resulting nonlinear least squares problem is overdetermined because we zero-pad the Fourier modes η j 1 , j 2 when |j 1 | or |j 2 | is larger than N{2. Assuming theη j 1 , j 2 are real (i.e. that η is even) also reduces the number of unknowns relative to the number of equations, which are enumerated by the M 2 gridpoints without exploiting symmetry. According to the linear theory of Section 2.5, we fix the two base Fourier modesη 1,0 andη 0,1 at nonzero amplitudes; these amplitudes are chosen independently. It might seem more natural to prescribe τ and η 1,0 and solve forη 0,1 along with b and the other unknown Fourier modes of η. However, linearization about the flat state leads to BR{Bτ " 0 (since η " 0 ñ κ " 0). This prevents the use of the implicit function theorem to solve the system in terms of τ andη 1,0 and would also cause difficulties for the numerical solver. Note that by (2.40) above, to linear order we have τ " g{k and c " a g`g{k. Variations in τ and c enter at higher order when the two amplitude parametersη 1,0 andη 0,1 are perturbed from 0, as shown below. The Levenberg-Marquardt solver requires a linear ordering of the unknowns. We enumerate theη j 1 , j 2 so that lower-frequency modes appear first. As the "shell index" s ranges from 1 to N{2, we enumerate all the index pairs pj 1 , j 2 q with maxp|j 1 |, |j 2 |q " s before increasing s. Within shell s, we proceed clockwise, along straight lines through the lattice, from p0, sq to ps, sq to ps,´sq to p1,´sq. The other Fourier modes are known from (2.35) and (2.36). Shell s contains 4s index pairs, so the total number of independent modesη j 1 , j 2 with maxp| j 1 |, | j 2 |q ď N{2 is ř N{2 s"1 4s " NpN{2`1q. We replaceη 1,0 by τ andη 0,1 by b in the list of unknowns to avoid additional shuffling of the variables when the prescribed base modes are removed from the list. Eventually there are NpN{2`1q parameters to compute: The objective function F is evaluated numerically by the trapezoidal rule approximation over T 2 , which is spectrally accurate: The parameters p j are chosen to minimize f ppq using the Levenberg-Marquardt method [31,38]. The method requires a Jacobian matrix pBr m {Bp j q m j , which we compute by solving the following variational equations: (3.10) δξ α " Hrδη α s, δ J " 2 pξ α δξ α`ηα δη α q , δκ "´3 2 κ δ J J`1 J 3{2´δ ξ α η αα`ξα δη αα´δ η α ξ αα´ηα δξ αα¯, We then have Br m Bp j " δRrτ, b,ηspα m 1 , α m 2 q, where m " 1`m 1`M m 2 and the jth column of the Jacobian corresponds to setting the perturbation δτ, δb or δη j 1 , j 2 corresponding to p j in (3.8) to 1 and the others to 0.
Like Newton's method, the Levenberg-Marquardt method generates a sequence of approximate solutions p p0q , p p1q , etc., which terminate when the residual drops below the desired tolerance or fails to decrease sufficiently. If no other solutions have been computed, we use the solution of the linearized problem as an initial guess: (3.11) η p0q pαq "η 1,0 pe iα`e´iα q`η 0,1 pe ikα`e´ikα q, τ p0q " τ lin " g{k, b p0q " c 2 lin " g`g{k. After computing two small-amplitude solutions, we use numerical continuation to increase the amplitude beyond the applicability of linear theory. In the present work, we hold the ratio γ "η 1,0 {η 0,1 constant to explore one-dimensional slices (or paths) through the twodimensional family of quasi-periodic traveling waves. We find that linear extrapolation from the previous two solutions on a path works well as the starting guess for the next Levenberg-Marquardt solve. Details of our Levenberg-Marquardt implementation, including stopping criteria and a strategy for delaying the re-computation of the Jacobian, are given in [38].

N R
We now present a detailed numerical study of solutions of (2.34) with k " 1{ ? 2 and g " 1 on three continuation paths corresponding to γ P t5, 1, 0.2u, where γ "η 1,0 {η 0,1 is the amplitude ratio of the prescribed base modes. In each case, we vary the larger of η 1,0 andη 0,1 from 0.001 to 0.01 in increments of 0.001. The initial guess for the first two solutions on each path are obtained using the linear approximation (3.11), which by (3.8) corresponds to As noted already, the amplitudesη 1,0 andη 0,1 are prescribed -they are not included among the unknowns. The initial guess for the remaining 8 solutions on each continuation path are obtained from linear extrapolation from the previous two computed solutions. In all cases, we use M " 60 for the grid size and N " 48 for the Fourier cutoff in each dimension. Figure 1 shows the initial conditions η and ϕ for the last solution on each continuation path (with maxtη 1,0 ,η 0,1 u " 0.01). Panels (a), (b) and (c) correspond to γ " 5, 1, and 0.2, respectively. The solution in all three cases is quasi-periodic, i.e. η and ϕ never exactly repeat themselves; we plot the solution from x " 0 to x " 36π as a representative snapshot. For these three solutions, the objective function f in (3.9) was minimized to 6.05ˆ10´2 8 , 9.28ˆ10´2 8 and 4.25ˆ10´2 8 , respectively, with similar or smaller values for lower-amplitude solutions on each path. The number of Jacobian evaluations in the Levenberg-Marquardt method for each of the 30 solutions computed on these paths never exceeded 5, and is typically 3 or 4. In our computations, η and ϕ are represented bỹ ηpα 1 , α 2 q andφpα 1 , α 2 q, which are defined on the torus T 2 . In Figure 2, we show contour plots ofηpα 1 , α 2 q andφpα 1 , α 2 q corresponding to the final solution on each path. Following the dashed lines through T 2 in Figure 2 leads to the plots in Figure 1. By construction in (2.36),ηp´αq "ηpαq whileφp´αq "´φpαq.
The amplitude ratio, γ :"η 1,0 {η 0,1 , determines the bulk shape of the solution. If γ " 1, the wave with wave number 1 will be dominant; if γ ! 1, the wave with wave number k " 2´1 {2 will be dominant; and if γ is close to 1, both waves together will be dominant over higher-frequency Fourier modes (at least in the regime we study here). This is demonstrated with γ " 5, 1 and 0. This can also be understood from the plots in Figure 2. In case (a), γ " 1 and the contour lines ofη andφ are perturbations of sinusoidal waves depending only on α 1 . The unperturbed waves would have vertical contour lines. The α 2 -dependence of the perturbation causes local extrema to form at the crest and trough. As a result, the contour lines join to form closed curves that are elongated vertically since the dominant variation is in the α 1 direction. Case (c) is similar, but the contour lines are elongated horizontally since the dominant variation is in the α 2 direction. Following the dashed lines in Figure 2, a cycle of α 1 is completed before a cycle of α 2 (since k ă 1). In case (a), a cycle of α 1 traverses the dominant variation ofη andφ on the torus, whereas in case (c), this is true of α 2 . So the waves in Figure 1 appear to oscillate faster in case (a) than case (c). In the intermediate case (b) with γ " 1, the contour lines of the crests and troughs are nearly circular, but not perfectly round. The amplitude of the waves in Figure 1  in Figure 2 pass near the extrema ofη andφ, and are smallest when the dashed lines pass near the zero level sets ofη andφ. If the slope of the dashed lines were closer to 1 and the functionsη andφ were to remain qualitatively similar to the results of panel (b) of Figure 2, the waves would have a beating pattern with many cycles with larger amplitude followed by many cycles with smaller amplitude. The former would occur when the dashed lines pass near the diagonal from p0, 0q to p2π, 2πq, which passes over the peaks and troughs of η andφ, while the latter would occur when the dashed lines pass near the lines connecting pπ, 0q to p2π, πq and p0, πq to pπ, 2πq, whereη andφ are close to zero. The dashed lines would linger in each regime over many cycles if k were close to 1.
In Figure 3, we plot the time evolution of ζpα, tq in the lab frame from t " 0 to t " 3 using the timestepping algorithm described in [37]. The initial conditions, plotted with thick blue lines, are those of the traveling waves computed in Figures 1 and 2 above. The grey curves F 3. Time evolution of the traveling wave profiles, ζpα, tq, from t " 0 to t " 3 in the lab frame. The thick blue lines correspond to the initial conditions.
give snapshots of the solution at uniformly sampled times with ∆t " 0.1. The solutions are plotted over the representative interval 0 ď x ď 12π, though they extend in both directions to˘8 without exactly repeating. Note that the solutions appear to propagate to the right at constant speed without changing shape. Our next goal is to verify this quantitatively to confirm that the quasi-periodic solutions we obtained by minimizing the objective function (3.9) are indeed traveling waves under the evolution equations (2.25).
Quantitative comparison requires an "exact" solution, which we take to be the numerically computed traveling wave, spatially shifted according to the exact time evolution derived in Corollary A.5 of Appendix A. In more detail, minimizing the objective function (3.9) gives the torus version of the traveling wave profileη 0 pα 1 , α 2 q, the surface tension τ, and the wave speed c such that p p η 0 q 1,0 and p p η 0 q 0,1 have prescribed values at t " 0. We then computeξ 0 " Hrη 0 s andφ 0 " cξ 0 , which are odd functions of α " pα 1 , α 2 q sinceη is even. From Corollary A.5, the time evolution of the traveling wave with these initial conditions under the torus version of (2.25) and (2.26) is given by where α 0 ptq " ct´Ap´ct,´kctq and Apx 1 , x 2 q is a periodic function on T 2 defined implicitly by (A.12) below. We see in (4.2) that the waves do not change shape as they move through the torus along the characteristic direction p1, kq, but the traveling speed α 1 0 ptq in conformal space varies in time in order to maintainξp0, 0, tq " 0 via (2.26). By Corollary A.5, the exact reconstruction ofξ exact fromη exact is (4.3)ξ exact pα 1 , α 2 , tq "ξ 0`α1´α0 ptq, α 2´k α 0 ptq˘`δ 0 ptq, where δ 0 ptq " ct´α 0 ptq " Ap´ct,´kctq measures the deviation in position from traveling at the constant speed ct in conformal space. The defining property (A.12) of Apx 1 , x 2 q ensures thatξ exact p0, 0, tq " 0. Figure 4 shows contour plots of the torus version of the γ " 5 and γ " 0.2 solutions shown in panels (a) and (c) of Figure 3 at the final time computed, T " 3. A similar plot of the γ " 1 solution is given in [37]. The dashed lines show the trajectory from t " 0 to t " T of the wave crest that begins at p0, 0q and continues along the path α 1 " α 0 ptq, α 2 " kα 0 ptq through the torus in (4.2). The following table gives the phase speed, c, surface tension, τ, translational shift in conformal space at the final time computed, α 0 pTq, and deviation from steady motion in conformal space, δ 0 pTq, for these three finite-amplitude solutions (recall that maxtη 1,0 ,η 0,1 u " 0.01 andη 1,0 {η 0,1 " γ) as well as for the zero-amplitude limit: In Figure 5, we plot δ 0 ptq for 0 ď t ď T (solid lines) along with pc´c lin qt (dashed and dotted lines) for the three finite-amplitude solutions in this table. Writing α 0 ptq " c lin t`rpc´c lin qt´δ 0 ptqs, we see that the deviation of α 0 ptq from linear theory over this time interval is due mostly to fluctuations in δ 0 ptq rather than the steady drift pc´c lin qt due to the change in phase speed c of the finite-amplitude wave.
Computing the exact solution (4.2) requires evaluating δ 0 ptq " Ap´ct,´kctq. We use Newton's method to solve the implicit equation (A.12) for Apx 1 , x 2 q at each point of a uniform MˆM grid, with M as in Section 3. We then use FFTW to compute the 2d Fourier representation of Apx 1 , x 2 q, which is used to quickly evaluate the function at any point. It would also have been easy to compute Ap´ct,´kctq directly by Newton's method, but the Fourier approach is also very fast and gives more information about the function Apx 1 , x 2 q. In particular, the modes decay to machine roundoff on the grid, corroborating the assertion in [37] that A is real analytic. We use the exact solution to compute the error in timestepping (2.25) and (2.26) from t " 0 to t " T, A detailed convergence study is given in [37] to compare the accuracy and efficiency of the Runge-Kutta and exponential time differencing schemes proposed in that paper using the γ " 1 traveling solution above as a test case. Here we report the errors for all three waves plotted in Figure 3 γ " 5 γ " 1 γ " 0.2 err 1.04ˆ10´1 6 1.16ˆ10´1 6 7.38ˆ10´1 7 using the simplest timestepping method proposed in [37] to solve (2.25), namely a 5th order explicit Runge-Kutta method using 900 uniform steps from t " 0 to t " 3. These errors appear to mostly be due to roundoff error in floating-point arithmetic, validating the accuracy of both the timestepping algorithm of [37] and the traveling wave solver of Section 3, which was taken as the exact solution. Evolving the solutions to compute these errors took less than a second on a laptop (with M 2 " 3600 gridpoints and 900 timesteps), while computing the traveling waves via the Levenberg-Marquardt method took 30-40 seconds on a laptop and only 3 seconds on a server (Intel Xeon Gold 6136, 3GHz) running on 12 threads (with M 2 " 3600 gridpoints and NpN{2`1q " 1200 unknowns). Next we examine the behavior of the Fourier modes that make up these solutions. Figure 6 shows two-dimensional plots of the Fourier modesη j 1 , j 2 for the 3 cases above, with γ P t5, 1, 0.2u and maxtη 1,0 ,η 0,1 u " 0.01. Only the prescribed modes and the modes that were optimized by the solver (see (3.8)) are plotted, which have indices in the range 0 ď j 1 ď N{2 and´N{2 ď j 2 ď N{2, excluding j 2 ď 0 when j 1 " 0. The other modes are determined by the symmetry of (2.36) and by zero-paddingη j 1 , j 2 " 0 if N{2 ă j 1 ď M or N{2 ă |j 2 | ď M. We used N " 48 and M " 60 in all 3 calculations. One can see that the fixed Fourier modesη 1,0 andη 0,1 are the two highest-amplitude modes in all three cases. In this sense, our solutions of the nonlinear problem (2.34) are small-amplitude perturbations of (3.11). However, there are many active Fourier modes, so these solutions are well outside of the linear regime. Carrying out a weakly nonlinear Stokes expansion to high enough order to accurately predict all these modes would be difficult, especially considering the degeneracies that arise already in the periodic Wilton ripple problem [35,36].
In panels (a), (b) and (c) of Figure 6, the modes appear to decay more slowly in one direction than in other directions. This is seen more clearly when viewed from above, as shown in panel (d) for the case of γ " 1. (The other two cases are similar). The direction along which the modes decay less rapidly appears to coincide with the line  tp j 1 , j 2 q : j 1`j2 k " 0u, which is plotted in red. A partial explanation is that when j 1`j2 k is close to zero, the corresponding modes e ip j 1`j2 kqα in the expansion of ηpαq in (2.12) have very long wavelength. Slowly varying perturbations lead to small changes in the residual of the water wave equations, so these modes are not strongly controlled by the governing equations (2.34). We believe this would lead to a small divisor problem [20] that would complicate a rigorous proof of existence of quasi-periodic traveling water waves.
Next we show that τ and c depend nonlinearly on the amplitude of the Fourier modesη 1,0 andη 0,1 . Panels (a) and (b) of Figure 7 show plots of τ and c versusη max :" maxpη 1,0 ,η 0,1 q for 9 values of γ "η 1,0 {η 0,1 , namely γ " 0.1, 0.2, 0.5, 0.8, 1, 1.25, 2, 5, 10. On each curve, η max varies from 0 to 0.01 in increments of 0.001. At small amplitude, linear theory predicts τ " g{k " 1.41421 and c " a gp1`1{kq " 1.55377. This is represented by the black marker atη max " 0 in each plot. For each value of γ, the curves τ and c are seen to have zero slope atη max " 0, and can be concave up or concave down depending on γ. This can be understood from the contour plots of panels (c) and (d). Both τ and c appear to be even functions ofη 1,0 andη 0,1 when the other is held constant. Both plots have a saddle point at the origin, are concave down in theη 1,0 direction holdingη 0,1 fixed, and are concave up in theη 0,1 direction holdingη 1,0 fixed. The solid lines in the first quadrant of these plots are the slices corresponding to the values of γ plotted in panels (a) and (b). The concavity of the 1d plots depends on how these lines intersect the saddle in the 2d plots.
The contour plots of panels (c) and (d) of Figure 7 were made by solving (2.34) with pη 1,0 ,η 0,1 q ranging over a uniform 26ˆ26 grid on the square r´0.01, 0.01sˆr´0.01, 0.01s. Using an even number of gridpoints avoids the degenerate case whereη 1,0 orη 0,1 is zero. At those values, the two-dimensional family of quasi-periodic solutions meets a sheet of periodic solutions where τ or c becomes a free parameter. Alternative techniques would be needed in these degenerate cases to determine the value of τ or c from which a periodic traveling wave in the nonlinear regime bifurcates to a quasi-periodic wave. In panel (e), we plot the magnitude of the Chebyshev coefficients in the expansion This was done by evaluating c on a cartesian product of two 16-point Chebyshev-Lobatto grids over r´0.01, 0.01s and using the one-dimensional Fast Fourier Transform in each direction to compute the Chebyshev modes. We see that the modes decay to machine precision by the time m`n ě 10 or so, and only even modes m and n are active. The plot for |τ mn | is very similar, so we omit it. These plots confirm the visual observation from the contour plots that τ and c are even functions ofη 1,0 andη 0,1 when the other is held constant. In summary, over the range´0.01 ďη 1,0 ,η 0,1 ď 0.01 considered here, τ and c show interesting nonlinear effects that would be difficult to model using weakly nonlinear theory since polynomials of degree 10 are needed to represent τ and c accurately to machine precision. Also, as seen in Figures 5 and 6 above, other aspects of the solution such as the deviation δ 0 ptq from traveling at a constant speed in conformal space and higher-frequency Fourier modesη j 1 , j 2 show greater sensitivity to nonlinear effects than c and τ do.

C
In this work, we have formulated the two-dimensional, infinite depth gravity-capillary traveling wave problem in a spatially quasi-periodic, conformal mapping framework. We have numerically demonstrated the existence of traveling solutions, which are a quasiperiodic generalizations of Wilton's ripples. To compute them, we adapted an overdetermined nonlinear least squares technique introduced in [38] for a different problem. For each solution computed, the value of k and the amplitudes of two base Fourier modesη 1,0 andη 0,1 are fixed while τ, c and the other Fourier modesη j 1 , j 2 are varied to search for solutions of (2.34). Before minimizing (3.9), the initial guess for each solution is computed using either the linear approximation (3.11) or numerical continuation. We validate the accuracy of the traveling solutions using the timestepping algorithm of [37]. To evolve at constant speed in physical space, we have shown that the 2d representation of the quasiperiodic waves travel at a nonuniform speed through the torus. We observed resonance effects in the Fourier modesη j 1 , j 2 along the line j 1`j2 k " 0 and computed the nonlinear dependence of phase speed and surface tension for the two-dimensional family with amplitude parameters in the range maxtˇˇη 1,0ˇ,ˇη0,1ˇu ď 0.01.
The question of what happens in our framework if k is rational is interesting. We believe the initial value problem (2.25) could still be solved, though in that case solving the torus version of the equations is equivalent to simultaneously computing a family of 1d solutions on a periodic domain. Families of 1d waves corresponding to a single solution of the torus problem are discussed in detail in [37], and take the form (2.28) above. If k " q{p with p and q relatively prime integers, the waves in this family all have period 2πp. The traveling wave problem becomes degenerate if k is rational -solutions of the torus version of (2.34) may still exist (we do not know), but if so, they are not unique. Indeed, if k " q{p as above andη 1 solves the torus version of (2.34), then for any 2π-periodic, real analytic function α 0 prq, ill also be a solution of (2.34) since the corresponding 1d functions passing through the torus along characteristic lines are related by a simple reparametrization, Another degeneracy is that the modesη j 1 , j 2 of a solution of (2.34) with j 1`k j 2 " 0 and pj 1 , j 2 q ‰ p0, 0q can be modified arbitrarily (maintainingη´j 1 ,´j 2 "η j 1 , j 2 ) to obtain additional solutions of (2.34). These modes are plane waves that only affect the 1d functions passing through the torus along characteristic lines by an additive constant. The resonance phenomenon observed in the Fourier modes in Figure 6 is presumably a small-divisor phenomenon [20] in the irrational case related to this degeneracy. If solutions for rational k exist, a natural open question is whether they can be selected to fit together continuously with solutions for nearby irrational wave numbers. In floating point arithmetic, irrational wave numbers are approximated by rational ones. We did not encounter difficulties with this, presumably because the above degeneracies are not visible with the grid resolution used. More work is needed to understand this rigorously.
Our results show that the amplitude ratio γ "η 1,0 {η 0,1 plays an important role in determining the shapes of solutions. As seen in Figures 1 and 3, the quasi-periodic features of the solutions are most evident when γ « 1. In the future, we plan to study the behavior of different perturbation families, e.g. fixing the amplitudes of different base Fourier modes in (2.39) such asη 1,0 andη 1,1 . We also aim to use this methodology to compute spatially quasiperiodic traveling gravity-capillary waves of finite depth, to compute the time evolution of solutions of the finite depth quasi-periodic initial value problem, and to study the stability of spatially quasi-periodic water waves along the lines of what has been done for periodic traveling waves [12,26] and Wilton ripples [35].
In this section we study the dynamics of the traveling waves computed in Section 2.4 under the evolution equations (2.25) for various choices of C 1 . We show that the waves maintain a permanent form but generally travel at a non-uniform speed in conformal space. We start by showing that there is a choice of C 1 for which η and ϕ remain stationary in time. We then show how C 1 changes when the waves are phase shifted by α 0 ptq, and how to determine α 0 ptq so that C 1 takes the value in (2.26). The evolution of the torus version of (2.34) under (2.25) is also worked out. We will need the following theorem and corollary, proved in [37]: Theorem A.1. Suppose ε ą 0 and zpwq is analytic on the half-plane Cέ " tw : Im w ă εu. Suppose there is a constant M ą 0 such that |zpwq´w| ď M for w P Cέ , and that the restriction ζ " z| R is injective. Then the curve ζpαq separates the complex plane into two regions, and zpwq is an analytic isomorphism of the lower half-plane onto the region below the curve ζpαq.
Corollary A.2. Suppose k ą 0 is irrational,ηpα 1 , α 2 q " ř p j 1 , j 2 qPZ 2η j 1 , j 2 e ipj 1 α 1`j2 α 2 q , and there exist constants C and ε ą 0 such that where K " maxpk, 1q. Let x 0 be real and defineξ " x 0`H rηs,ζ "ξ`iη and where the sum is over all integer pairs pj 1 , j 2 q satisfying the inequality. Suppose also that for each fixed θ P r0, 2πq, the function α Þ Ñ ζpα; θq " α`ζpα, θ`kαq is injective from R to C and ζ α pα; θq ‰ 0 for α P R. Then for each θ P R, the curve ζpα; θq separates the complex plane into two regions and is an analytic isomorphism of the lower half-plane onto the region below ζpα; θq. Moreover, there is a constant δ ą 0 such that |z w pw; θq| ě δ for Im w ď 0 and θ P R.