Spatially Quasi-Periodic Water Waves of Infinite Depth

We formulate the two-dimensional gravity-capillary water wave equations in a spatially quasi-periodic setting and present a numerical study of solutions of the initial value problem. We propose a Fourier pseudo-spectral discretization of the equations of motion in which one-dimensional quasi-periodic functions are represented by two-dimensional periodic functions on a torus. We adopt a conformal mapping formulation and employ a quasi-periodic version of the Hilbert transform to determine the normal velocity of the free surface. Two methods of time-stepping the initial value problem are proposed, an explicit Runge–Kutta (ERK) method and an exponential time-differencing (ETD) scheme. The ETD approach makes use of the small-scale decomposition to eliminate stiffness due to surface tension. We perform a convergence study to compare the accuracy and efficiency of the methods on a traveling wave test problem. We also present an example of a periodic wave profile containing vertical tangent lines that is set in motion with a quasi-periodic velocity potential. As time evolves, each wave peak evolves differently, and only some of them overturn. Beyond water waves, we argue that spatial quasi-periodicity is a natural setting to study the dynamics of linear and nonlinear waves, offering a third option to the usual modeling assumption that solutions either evolve on a periodic domain or decay at infinity.


Introduction
Linear and nonlinear wave equations are generally studied under the assumption that the solution is spatially periodic or decays to zero at infinity (Lax, 1976). Beginning with Berenger (1994), a great deal of effort has been devoted to developing perfectly matched layer (PML) techniques for imposing absorbing boundary conditions over a finite computational domain to simulate wave propagation problems on unbounded domains. However, in many situations, assuming the waves decay to zero at infinity is not a realistic model. For example, a large body of water such as the ocean is often covered in surface waves in every direction over vast distances. But assuming spatial periodicity may limit one's ability to observe interesting dynamics. In this paper, we formulate the initial value problem of the surface water wave equations in a spatially quasi-periodic setting, design numerical algorithms to compute such waves, and study their properties.
Since the pioneering work of Benjamin and Feir (1967) and Zakharov (1968), it has been recognized that water waves exhibit interesting nonlinear interactions between component waves of different wavelength. For example, in oceanography, modulational instabilities of periodic wavetrains introduce perturbations that lead to spatially quasi-periodic dynamics and are believed to be one of the mechanisms responsible for the formation of rogue waves (Osborne et al. 2000;Onorato et al. 2006;Ablowitz and Horikis 2015). These instabilities have been studied extensively using a variety of techniques, summarized below, including linearization using Bloch stability theory, evolving the nonlinear equations on a larger periodic domain, developing coupled weakly nonlinear models, and solving weakly nonlinear models via the inverse scattering transform. However, it has not been known how to formulate or compute fully nonlinear water waves in a spatially quasi-periodic setting. We show that a conformal mapping formulation of the water wave equations, introduced by Dyachenko et al. (1996a) and further developed by many authors (Dyachenko et al. 1996b;Choi and Camassa 1999;Dyachenko 2001;Zakharov et al. 2002;Li et al. 2004;Milewski et al. 2010;Wang et al. 2013;Viotti et al. 2014;Dyachenko and Newell 2016;Turner and Bridges 2016), extends nicely to this setting via a quasi-periodic generalization of the Hilbert transform. Currently our method is limited to two-dimensional fluids, but we formulate the equations of motion and discuss computational challenges of 3D quasi-periodic water waves in Appendix D.
The Bloch stability approach can be carried out by linearizing the full water wave equations about a traveling Stokes wave (Longuet-Higgins 1978;McLean 1982;MacKay and Saffman 1986;Deconinck and Oliveras 2011;Trichtchenko et al. 2016;Wilkening 2021) or within a weakly nonlinear model such as the nonlinear Schrödinger (NLS) equation (Benney and Newell 1967;Zakharov 1968). A major drawback is that unstable modes grow exponentially forever and eventually leave the realm of validity of the linearization. Osborne et al. (2000) have observed that if nonlinear effects are taken into account in this scenario, the perturbation often exhibits Fermi-Pasta-Ulam recurrence (Berman and Izrailev 2005). They solve a 2+1-dimensional NLS equation on a domain that is 10 times larger than the wavelength of the carrier wave and look for rogue-wave formation and recurrence over long simulation times. Similarly, Bryant and Stiassnie (1994) study recurrence using both a weakly nonlinear model (Zakharov's equation) and the full water wave equations in the context of standing water waves when the wavelength of the subharmonic perturbation is 9 times that of the unperturbed standing wave. The main drawback of this approach is that the larger periodic computational domain must be an integer multiple of both the base wave and the perturbation, which requires that the ratio of their wavelengths be a rational number with a small numerator and denominator. We propose a method below that allows for more general perturbations and plan to investigate the long-time nonlinear dynamics of unstable subharmonic perturbations of traveling and standing waves in future work.
An alternative approach that does not require rationally related wave numbers is to model the interaction of two periodic wavetrains as a coupled weakly nonlinear system. This is particularly useful for studying the interaction between oblique waves on the surface of a three-dimensional fluid. For example, Bridges and Laine-Pearson (2001) studied coupled NLS equations and extended the theory to analyze the stability of short-crested waves (Bridges and Laine-Pearson 2005). More recently, Ablowitz and Horikis (2015) showed that some propagation angles enhance the number and amplitude of rogue wave events in a coupled NLS water wave model. Some weakly nonlinear models are completely integrable and can be studied using the inverse scattering transform (IST) (Ablowitz and Segur 1981). Osborne et al. (2000) discuss IST results for the 1+1-dimensional NLS equation in the context of rogue waves. Other equations such as the Korteweg-deVries and Benjamin-Ono equations are meant to model wave dynamics in shallow water (Ablowitz and Segur 1981) and internal waves in a stratified fluid (Ono 1975), respectively. Solving these equations using the inverse scattering transform (Ablowitz and Segur 1981) leads to infinite hierarchies of exact spatio-temporal quasi-periodic solutions (Flaschka et al. 1980;Dobrokhotov and Krichever 1991).
In these and other examples involving weakly nonlinear theory, it is natural to seek analogous quasi-periodic solutions of the Euler equations in regimes where the model equations are intended to be accurate. It is also of interest to search for new regimes and behavior not predicted by model water wave equations. Even within weakly nonlinear theory, except for the exact quasi-periodic solutions obtained via the IST, spatial quasi-periodicity has only been approximated by embedding in a larger periodic domain or by introducing coupling terms between two or more single-mode NLS equations. The framework we propose below for water waves, which involves representing quasi-periodic functions as periodic functions on a higher dimensional torus and using a spectral method to solve a torus version of the equations of motion, could also be used to find true quasi-periodic solutions of weakly nonlinear equations without introducing systems of coupled equations.
Only recently have quasi-periodic dynamics of water waves been studied mathematically. Berti and Montalto (2016) and Baldi et al. (2018) used Nash-Moser theory to prove the existence of small-amplitude temporally quasi-periodic gravity-capillary standing waves. With different assumptions on the form of solutions, Berti et al. (2020) have proved the existence of time quasi-periodic gravity-capillary waves with constant vorticity while Feola and Giuliani (2020) have proved the existence of time quasi-periodic irrotational gravity waves. New families of relative-periodic (Wilkening 2015) and traveling-standing (Wilkening 2021) water wave solutions have been computed by Wilkening. As with Berti and Montalto (2016); Baldi et al. (2018); Berti et al. (2020); Feola and Giuliani (2020), these solutions are quasi-periodic in time rather than space.
Another motivation for studying spatially quasi-periodic water waves is the work of Wilton (1915), who observed that a resonance can occur that causes Stokes' regular perturbation expansion for traveling water waves to break down (Akers and Gao 2012;Akers and Nicholls 2020;Vanden-Broeck 2010;Trichtchenko et al. 2016;Akers et al. 2019). Suppose k 1 and k 2 are both roots of the dispersion relation for linearized water waves of infinite depth, c 2 = gk −1 + τ k. (1.1) Here the gravitational acceleration g, wave speed c, and surface tension τ are held constant when solving for the wave numbers k 1 and k 2 . If k 2 = K k 1 with K ≥ 2 an integer, the K th harmonic will enter a modified Stokes expansion for traveling waves of wavelength 2π/k 1 at order ε max(K-2,1) instead of ε K , where ε is the expansion coefficient of the fundamental mode. This resonance occurs because the two waves travel at the same speed under the linearized water wave equations. Bridges and Dias (1996) consider a generalization in which the wave numbers k 1 and k 2 are irrationally related. They use 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. This inspired us to develop a conformal mapping framework for computing spatially quasi-periodic, fully nonlinear traveling gravitycapillary waves, which is the topic of the companion paper, namely Wilkening and Zhao (2021). We show in (Wilkening and Zhao 2021) that these spatially quasi-periodic traveling waves come in two-parameter families in which the amplitudes of the base modes with wave numbers k 1 and k 2 serve as bifurcation parameters. The wave speed and surface tension depend nonlinearly on these parameters as well. Akers et al. (2019) have proved existence of similar two-parameter families of traveling waves for the case of a two-fluid hydro-elastic interface. They develop an integral equation formulation of the equations governing traveling hydro-elastic waves such that the linearization about any state is a compact perturbation of the identity and use global bifurcation theory to establish existence and uniqueness results. They show that the nullspace of the linearized operator about the flat rest state has dimension one or two, and is two if and only if the non-dimensionalized wave numbers k 1 < k 2 that travel with a given speed are integers. When they are integers, Akers et al. distinguish resonant and non-resonant cases depending on whether k 2 /k 1 is an integer. The non-resonant case leads to a smooth two-parameter family of traveling waves with wave speed and surface tension depending nonlinearly on the amplitude parameters. This is the case most analogous to the spatially quasi-periodic traveling waves that we compute in (Wilkening and Zhao 2021).
The present paper focuses on the more general spatially quasi-periodic initial value problem, which we use to validate the traveling wave computations of Wilkening and Zhao (2021) and explore new dynamic phenomena. In recent years, conformal mapping methods have proved useful for studying two-dimensional traveling water waves (Choi and Camassa 1999;Milewski et al. 2010;Wang et al. 2013;) and time-dependent water waves (Dyachenko et al. 1996a;Dyachenko 2001;Dyachenko et al. 1996b;Dyachenko and Newell 2016;Li et al. 2004;Milewski et al. 2010;Turner and Bridges 2016;Zakharov et al. 2002) with periodic boundary conditions. We introduce a Hilbert transform for quasi-periodic functions to compute the normal velocity and maintain a conformal parametrization of the free surface. This leads to a numerical method to compute the time evolution of solutions of the Euler equations from arbitrary quasi-periodic initial data. Following the definitions in Moser (1966); Dynnikov and Novikov (2005), we represent a general quasi-periodic function u(α) in one dimension by a periodic functionũ(α) on a d-dimensional torus, i.e. u(α) =ũ(kα) for α ∈ R, where k = (k 1 , . . . , k d ). The k i are assumed to be linearly independent over the integers. We take these basic wave numbers k i and the initial conditions on the torus as given, focusing on the d = 2 case. This leaves open the important question of how best to measure a one-dimensional wave profile or velocity potential and identify its quasi-periods and corresponding torus function.
We present two variants of the numerical method, one in a high-order explicit Runge-Kutta framework and one in an exponential time-differencing (ETD) framework. The former is suitable for the case of zero or small surface tension while the latter makes use of the small-scale decomposition (Hou et al. 1994(Hou et al. , 2001 to eliminate stiffness due to surface tension. The conformal mapping method has not been implemented in an ETD framework before, even for periodic boundary conditions. We present a convergence study of the methods as well as a large-scale computation of a quasi-periodic wave in which some of the wave crests overturn when evolved forward in time while others do not. Due to the torus representation of solutions, there are infinitely many wave crests and no two of them evolve in exactly the same way. The computation involves over 33 million degrees of freedom evolved over 5400 time steps to maintain double-precision accuracy. We include four appendices that cover various technical aspects of this work. In Appendix A, we prove a theorem establishing sufficient conditions for an analytic function z(w) to map the lower half-plane topologically onto a semi-infinite region bounded above by a parametrized curve and for 1/|z w | to be uniformly bounded. In Appendix 1, we study families of quasi-periodic solutions obtained by introducing phases in the reconstruction formula for extracting 1D quasi-periodic functions from periodic functions on a torus. This enables us to prove that if all the solutions in the family are single-valued and have no vertical tangent lines, the solutions are also quasi-periodic in the original graph-based formulation of the Euler equations. We also present a simple procedure for computing the change of variables from the conformal representation to the graph representation. This appears to be a new result even for periodic boundary conditions. In Appendix C we provide details on how to implement the equations of motion in an exponential time-differencing framework to avoid stepsize limitations due to stiffness caused by surface tension. And in Appendix D, we discuss the equations of motion for spatially quasi-periodic water waves in three dimensions and outline possible alternatives to the conformal mapping approach.

Mathematical Formulation
In this section, we review the governing equations for gravity-capillary waves in both physical space and conformal space. We then extend the conformal mapping framework to allow for spatially quasi-periodic solutions. For simplicity, we initially assume the wave profile η(x, t) remains single-valued. This assumption is relaxed when discussing the conformal formulation, and an example of a wave in which some of the peaks overturn as time advances is presented in Sect. 4.2.

Governing Equations in Physical Space
Gravity-capillary waves of infinite depth are governed by the two-dimensional freesurface Euler equations (Zakharov 1968;Craig and Sulem 1993) where x is the horizontal coordinate, y is the vertical coordinate, t is the time, (x, y, t) is the velocity potential in the fluid, η(x, t) is the free surface elevation, is the boundary value of the velocity potential on the free surface, g is the vertical acceleration due to gravity and τ is the coefficient of surface tension. Following Zakharov (1968); Craig and Sulem (1993), only the surface variables η and ϕ are evolved in time; the velocity potential in the bulk fluid is reconstructed from η and ϕ by solving (2.2), which causes the problem to be nonlocal. The function C(t) in the Bernoulli condition (2.4) is an arbitrary integration constant that is allowed to depend on time but not space. When the domain is periodic or quasi-periodic, one can choose C(t) so that the mean value of ϕ(x, t) remains constant in time, where the mean is defined as lim a→∞

The Quasi-Periodic Hilbert Transform
We find that a conformal mapping representation of the free surface greatly simplifies the solution of the Laplace equation for the velocity potential in the quasi-periodic setting. In this section, we establish the properties of the Hilbert transform that will be needed to study quasi-periodic water waves in a conformal mapping framework.
As defined in Moser (1966); Dynnikov and Novikov (2005), a quasi-periodic, real analytic function u(α) is a function of the form where ·, · denotes the standard inner product in R d andũ is a periodic, real analytic function defined on the d-dimensional torus Entries of the vector k are called the basic wave numbers (or basic frequencies) of u and are required to be linearly independent over Z. If k is given, one can reconstruct the Fourier coefficientsû j from u viâ A similar averaging formula holds for functions in the more general class of almost periodic functions (Moser 1966;Bohr 2018;Fink 1974;Amerio and Prouse 1971;Hino et al. 2002), which is the closure with respect to uniform convergence on R of the set of trigonometric polynomials p(x) = N n=1 c n e iκ n x . Before taking limits to obtain the closure, this set includes polynomials of any degree and there is no restriction on the real numbers κ n . Within the framework of almost periodic functions, one obtains quasi-periodic functions if one assumes the κ n in the approximating polynomials are integer linear combinations of a fixed, finite set of basic wave numbers k 1 , . . . , k d .
We have not attempted to formulate the water wave problem in the full generality of almost periodic functions, and instead assume the basic wave numbers are given and the torus representation (2.6) is available. Thus, the average over R on the right-hand side of (2.8) can be replaced by the simpler Fourier coefficient formulâ Our assumption thatũ(α) is real analytic is equivalent to the conditions thatû − j = u j for j ∈ Z d and there exist positive numbers M and σ such that |û j | ≤ Me −σ j , i.e. the Fourier modesû j decay exponentially as j → ∞. This is proved, e.g. in Lemma 5.6 of Broer and Takens (2011).
Next we define the projection operators P and P 0 that act on u andũ via (2.10) Note that P projects onto the space of zero-mean functions while P 0 returns the mean value, viewed as a constant function on R or T d . There are two versions of P and P 0 , one acting on quasi-periodic functions defined on R and one acting on torus functions defined on T d . Given u(α) as in (2.6), the most general bounded analytic function f (w) in the lower half-plane whose real part agrees with u on the real axis has the form wherev 0 ∈ R and the sum is over all j ∈ Z d satisfying j , k < 0. The imaginary part of f (z) on the real axis is given by (2.12) where sgn(q) ∈ {1, 0, −1} depending on whether q > 0, q = 0 or q < 0, respectively. Similarly, given v(α) and requiring (Im f )| R = v yields (2.11) withû j replaced by iv j . We introduce a quasi-periodic Hilbert transform to compute is a free parameter when computing v or u, respectively. H returns the "zero-mean" solution, i.e. P 0 H [u] = 0. Uniqueness of the bounded extension from u or v to f up to the additive constant iv 0 orû 0 follows from two well-known results: the only bounded solution of the Laplace equation on a half-space satisfying homogeneous Dirichlet boundary conditions is identically zero (Axler et al. 1992), and the harmonic conjugate of the zero function on a connected domain is constant.

The Conformal Mapping
We consider a time dependent conformal mapping that maps the conformal domain to the fluid domain This conformal mapping, denoted by z(w, t), is assumed to extend continuously to C − and maps the real line β = 0 to the free surface (2.20) We express z(w, t) as We also introduce the notation ζ = z| β=0 , ξ = x| β=0 and η = y| β=0 so that the free surface is parametrized by This allows us to denote a generic field point in the physical fluid by z = x + iy while simultaneously discussing points ζ = ξ + iη on the free surface. To avoid ambiguity, we will henceforth denote the free surface elevation function from the previous section by η phys (x, t). Thus, (2.23) The parametrization (2.22) is more general than (2.23) in that it allows for overturning waves. In deriving the equations of motion for ζ(α, t) and ϕ(α, t) in Sect. 2.5 below, we will indicate the modifications necessary to handle the case of overturning waves.
In particular, as discussed in Appendix A, (t) is defined in this case as the image of ζ(·, t), which is assumed to be injective on R, and (t) can be obtained from (t) using the Jordan curve theorem. The conformal map is required to remain a bounded distance from the identity map in the lower half-plane. Specifically, we require that where M(t) is a uniform bound that could vary in time. The Cauchy integral formula implies that |z w − 1| ≤ M(t)/|β|, so at any fixed time, Our goal is to investigate the case when the free surface is quasi-periodic in α. This differs from conformal mappings discussed in Meiron et al. (1981); Dyachenko et al. (1996a); Dyachenko (2001); Zakharov et al. (2002); Li et al. (2004) and Milewski et al. (2010), where it is assumed to be periodic. In the present work, η is assumed to have two spatial quasi-periods, i.e. at any time it has the form (2.6) with d = 2 and k = (k 1 , k 2 ). Since k 1 and k 2 are irrationally related, we assume without loss of generality that k 1 = 1 and k 2 = k, where k is irrational: (2.26) 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.27) Using (2.25) and Im[z w | β=0 ] = η α or differentiating (2.27) gives (2.28) We use a tilde to denote the periodic functions on the torus that correspond to the quasi-periodic parts of ξ , ζ and z, (2.30) While the mean surface height remains constant in physical space,η 0,0 (t) generally varies in time. Since the modesη j 1 , j 2 are assumed to decay exponentially, there is a uniform bound M(t) such that |z(α 1 , α 2 , β, t)| ≤ M(t) for (α 1 , α 2 ) ∈ T 2 and β ≤ 0. In Appendix A we show that as long as the free surface ζ(α, t) does not self-intersect at a given time t, the mapping w → z(w, t) is an analytic isomorphism of the lower half-plane onto the fluid region.

The Complex Velocity Potential
Let phys (x, y, t) denote the velocity potential in physical space from Sect. 2.1 above and let W phys (x + iy, t) = phys (x, y, t) + i phys (x, y, t) be the complex velocity potential, where phys is the stream function. Using the conformal mapping (2.21), we pull back these functions to the lower half-plane and define We also define ϕ = | β=0 and ψ = | β=0 and use (2.2) and (2.23) to obtain We assume ϕ is quasi-periodic with the same quasi-periods as η, The fluid velocity ∇ phys (x, y, t) is assumed to decay to zero as y → −∞ (since we work in the lab frame). From (2.25) and the chain rule (see (2.33) (2.32) Here we have set the integration constant to zero and assumed P 0 [ϕ] =φ 0,0 (t) = 0 and P 0 [ψ] =ψ 0,0 (t) = 0, which is allowed since and can be modified by additive constants (or functions of time only) without affecting the fluid motion.

Governing Equations in Conformal Space
Following Dyachenko et al. (1996a); Choi and Camassa (1999); Zakharov et al. (2002); Li et al. (2004); Viotti et al. (2014) and Turner and Bridges (2016), we present a derivation of the equations of motion for surface water waves in a conformal mapping formulation, modified as needed to handle quasi-periodic solutions. We also justify the assumption that z t /z α remains bounded in the lower half-plane, which we have not seen discussed previously in the literature.
From the chain rule, (2.33) Evaluating (2.33) on the free surface gives (2.34) Using (2.23) and (2.34) in (2.3) and multiplying by ξ α , we obtain This states that the normal velocity of the free surface is equal to the normal velocity of the fluid,n · (ξ t , η t ) =n · ∇ phys , wheren = (−η α , ξ α )/ √ J . This can also be obtained by tracking a fluid particle x p (t) + iy p (t) = ζ(α p (t), t) on the free surface. We haveẋ p = ξ αα p + ξ t = phys x andẏ p = η αα p + η t = phys y , which leads to (2.35) after eliminatingα p . This argument does not assume the free surface is a graph, i.e. (2.35) is also valid for overturning waves.
Next we define a new function, Since q is quasi-periodic in α and extends analytically to the lower half-plane via z t /z α , the real and imaginary part of q can be related by the Hilbert transform. Here we have assumed that z t /z α is bounded, which will be justified below. Thus, where C 1 is an arbitrary integration constant that may depend on time but not space. Letŝ = (ξ α , η α )/ √ J denote the unit tangent vector to the curve. Equation (2.37) prescribes the tangential velocityŝ · (ξ t , η t ) of points on the curve in terms of the normal velocity in order to maintain a conformal parametrization. Note that the tangent velocity of the curve differs from that of the underlying fluid particles. This is similar in spirit to a method of Hou et al. (1994Hou et al. ( , 2001, who proposed a tangential velocity that maintains a uniform parametrization of the curve (rather than a conformal one); see also (Choi and Camassa 1999;Turner and Bridges 2016;Ambrose et al. 2021;Wilkening and Zhao 2021). Combining (2.35) and (2.37), we obtain the kinematic boundary conditions in conformal space, (2.38) The right-hand side can be interpreted as complex multiplication of z α with z t /z α . Since both functions are analytic in the lower half-plane, their product is, too. Thus, ξ t is related to η t via the Hilbert transform (up to a constant). The constant is determined by comparing (2.27) with (2.38), which gives (2.39) The three most natural choices of C 1 are (2.40) In options (b) and (c), the evolution equation ensures that dx 0 /dt = 0 and ξ t (0, t) = 0, respectively; we have assumed the initial conditions satisfy x 0 (0) = 0 or ξ(0, 0) = 0, respectively. Option (c) amounts to setting x 0 (t) = −H [η](0, t) in (2.27). This arguably leads to the most natural parametrization, but would have a problem if the vertical part of an overturning wave crosses α = 0. Indeed, such a crossing would lead to ξ α (0, t) = 0 at some time t in the denominator of (2.40 c). We recommend option (b) in this scenario.

Remark 2.3
In infinite depth as considered here, iff andg are torus functions on T d and f (α) =f (kα), g(α) =g(kα), the identity is easily proved, where the sum is over j ∈ Z d satisfying j , k = 0. In the periodic case with d = 1 and k = (1) or the quasi-periodic case with k = (k 1 , . . . , k d ) and the k i linearly independent over the integers, the right-hand side of (2.41) isf 0ĝ0 . This simplifies (2.40 b) to C 1 = 0 and (2.39) to dx 0 /dt = C 1 , i.e. cases (a) and (b) in (2.40) coincide. However, this only works in infinite depth as the finite-depth Hilbert transform with symbol −i tanh j , k h does not satisfy (2.41). Here h is the fluid depth in conformal space, which evolves in time to maintain constant fluid depth in physical space (Li et al. 2004;Turner and Bridges 2016). Finite-depth quasi-periodic water waves will be investigated in future work.
Next we evaluate the Bernoulli equation phys t + 1 2 ∇ phys 2 + p ρ + gy = C 2 at the free surface to obtain an evolution equation for ϕ(α, t). Here C 2 is an arbitrary integration constant that may depend on time but not space. The pressure at the free surface is determined by the Laplace-Young condition, p = p 0 − ρτ κ, where κ is the curvature, ρτ is the surface tension, and p 0 is a constant that can be absorbed into C 2 and set to zero. From (2.33) or (2.34), we know ∇ phys 2 = (ϕ 2 α + ψ 2 α )/J on the free surface. Finally, differentiating ϕ(α, t) = phys (ξ(α, t), η(α, t), t) and using (2.34) and (2.38), we obtain (2.42) We choose C 2 so that P 0 [ϕ t ] = 0. In conclusion, we obtain the following governing equations for spatially quasi-periodic gravity-capillary waves in conformal space (2.43) Note that these equations govern the evolution of x 0 , η and ϕ, which determine the state of the system. The functions ξ , ψ, J , χ and κ are determined at any moment by x 0 , η and ϕ through the auxiliary equations in (2.43). We emphasize that C 1 can be chosen arbitrarily as long as dx 0 /dt satisfies (2.39). The special cases (2.40 b) and (2.40 c) lead to nice formulas for x 0 (t) without having to evolve (2.39) numerically. An alternative approach was proposed by Li et al. (2004), who set C 1 = 0 (by not introducing it) and avoid writing down a differential equation for x 0 (t) by instead solving both the ξ t and η t equations in (2.38).
In deriving (2.43) from (2.35) and (2.42), we had to assume z t /z α remains bounded in the lower half-plane. Conditions that ensure the boundedness of 1/z α are given in Appendix A. We note that z t /z α is automatically bounded in the converse direction, where (2.35) and (2.42) are derived from (2.43). In more detail, when solving (2.43), z t /z α is constructed first, before z t , as the bounded extension of the quasiperiodic function with imaginary part (−ψ α /J ) to the lower half-plane. Equation (2.38) then defines z t as the product of this function by z α , which is also bounded since ξ α = 1 + H [η α ]. Thus, the first component of each side of (2.38) is related to the corresponding second component by the Hilbert transform, up to a constant. Since the second components are equal (i.e. the η t equation holds), the ξ t equation also holds -the constants are accounted for by (2.39). Left-multiplying (2.38) by the row vector (−η α , ξ α ) gives the kinematic condition (2.35), as required.
Equations (2.43) break down if J becomes zero somewhere on the curve. Such a singularity would arise, for example, if the wave profile were to form a corner in finite time. To our knowledge, it remains an open question whether the free-surface Euler equations can form such a corner.

Remark 2.4 Equations
which is the one place this notation becomes awkward. Using (2.27) and (2.29),ζ is completely determined by x 0 (t) andη, so only these have to be evolved -the formula forξ t in (2.38) is redundant as long as (2.39) is satisfied. If both components ofζ (α 1 , α 2 , t) are given, we say that (ζ ,φ) satisfy the torus version of (2.43) if there is a continuously differentiable function and if x 0 ,η andφ satisfy the torus version of (2.43) with We show in Appendix B that solving the torus version of (2.43) yields a threeparameter family of one-dimensional solutions of the form We also show that if all the waves in this family are single-valued and have no vertical tangent lines, there is a corresponding family of solutions of the Euler equations in the original graph-based formulation of (2.1-2.4) that are quasi-periodic in physical space. A precise statement is given in Theorem B.2 and the discussion that follows.

Numerical Method
In this section, we describe a pseudo-spectral time-stepping strategy for evolving water waves with spatially quasi-periodic initial conditions. The evolution equations (2.43) for η and ϕ are nonlinear and involve computing derivatives, antiderivatives and Hilbert transforms of quasi-periodic functions. Let f denote one of these functions (e.g. η, ϕ or χ ) and letf denote the corresponding periodic function on the torus, We propose a pseudo-spectral method in which each such f that arises in the formulas (2.43) is represented by the values off at M 1 × M 2 equidistant gridpoints on the torus T 2 , We visualize a 90 • rotation between the matrixF holding the entriesf m 1 ,m 2 and the collocation points in the torus. The columns of the matrix correspond to horizontal slices of gridpoints while the rows of the matrix correspond to vertical slices indexed from bottom to top. The nonlinear operations in (2.43) consist of products, powers and division; they are carried out pointwise on the grid. Derivatives and the Hilbert transform are computed in Fourier space via (3.2). To plot the solution, we also need to compute an antiderivative to get ξ from f = ξ α . This involves dividingf j 1 , j 2 by i( j 1 + j 2 k) when ( j 1 , j 2 ) = (0, 0) and adjusting the (0, 0) mode to obtain ξ(0, t) = 0.
Since the functions f that arise in the computation are real-valued, we use the realto-complex ('r2c') version of the two-dimensional discrete Fourier transform. The 'r2c' transform of a one-dimensional array of length M (assumed even) is given by In practice, theĝ j are computed simultaneously in O(M log M) time rather than by this formula. The fully complex ('c2c') transform of this (real) data would give additional valuesĝ j with M/2 + 1 ≤ j ≤ M − 1. These extra entries are actually aliased values of negative-index modes; they are redundant due toĝ j =ĝ j−M =ĝ M− j . Since the imaginary components ofĝ 0 andĝ M/2 are zero, the number of real degrees of freedom on both sides of (3.4) is M. The Nyquist modeĝ M/2 requires special attention. Settinĝ g M/2 = 1 and the other modes to zero yields g m = cos(π m) = (−1) m . The derivative and Hilbert transform of this mode are taken to be zero since they would involve evaluating sin(Mα/2) at the gridpoints α m = 2π m/M. The two-dimensional 'r2c' transform can be computed by applying one-dimensional 'r2c' transforms in the x-direction (i.e. to the columns ofF) followed by onedimensional 'c2c' transforms in the y-direction (i.e. to the rows ofF): (3.5) The 'r2c' routine in the FFTW library actually returns the index range 0 ≤ j 2 < M 2 , but we usef j 1 , j 2 −M 2 =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 byf This imposes additional constrains on the computed Fourier modes, namely Im{f 0,0 } = 0, Im{f M 1 /2,0 } = 0, Im{f 0,M 2 /2 } = 0, Im{f M 1 /2,M 2 /2 } = 0, 7) where we also usedf −M 1 /2, j 2 =f M 1 /2, j 2 . This reduces the number of real degrees of freedom in the complex (M 1 /2 + 1) × M 2 array of Fourier modes to M 1 M 2 . When computing f α and H [ f ] via (3.2), the Nyquist modes with j 1 = M 1 /2 or j 2 = M 2 /2 are set to zero. Otherwise the formulas (3.2) respect the constraints (3.7) and the 'c2r' transform reconstructs real-valued functions f α and H [ f ] from their Fourier modes.
The evolution Eq. (2.43) are not stiff when the surface tension parameter is small or vanishes, but become moderately stiff for larger values of τ . We find that the 5th and 8th order explicit Runge-Kutta methods of Dormand and Prince (Hairer et al., 2000) work well for smaller values of τ , and exponential time-differencing (ETD) methods (Cox and Matthews, 2002;Kassam and Trefethen, 2005;Berland et al., 2005;Whalen et al., 2015;Chen and Wilkening, 2021) work well generally. This will be demonstrated in Sects. 4.1 and 4.2 below. In the ETD framework, we follow the basic idea of the small-scale decomposition for removing stiffness from interfacial flows (Hou et al. 1994(Hou et al. , 2001 and write the evolution Eq. (2.43) in the form where P is the projection in (2.10), H is the Hilbert transform in (2.14), and Note that N is obtained by subtracting the terms included in L from (2.43).
In particular, ψ α in (3.9) is −H ∂ α ϕ from (3.8). The eigenvalues of L are ±i | j 1 + j 2 k| g + τ ( j 1 + j 2 k) 2 , so the leading source of stiffness is dispersive. This 3/2 power growth rate of the eigenvalues of the leading dispersive term with respect to wave number is typical of interfacial fluid flows with surface tension (Hou et al. 1994(Hou et al. , 2001. For stiffer problems such as the Benjamin-Ono and KdV equations, the growth rate is faster (quadratic and cubic, respectively) and it becomes essential to use a semi-implicit or exponential time-differencing scheme to avoid severe timestepping restrictions. Here it is less critical, but still useful. Further details on how to implement (3.8) and (3.9) in the ETD framework are given in Appendix C.
In both the explicit Runge-Kutta and ETD methods, as explained above, the functions evolved in time areη(α 1 , α 2 , t) andφ(α 1 , α 2 , t), sampled on the uniform M 1 ×M 2 grid covering T 2 . At the end of each time step, we apply a 36th order filter (Hou and Li 2007;Hou et al. 1994) with Fourier multiplier 36 otherwise.
(3.10) In all the computations reported below, we used the same number of gridpoints in the α 1 and α 2 -directions, M 1 = M 2 = M. It is easy to check a-posteriori that the Fourier modes decay sufficiently (e.g. to machine precision) by the time the filter deviates appreciably from 1. If they do not, the calculation can be repeated with a larger value of M. This will be demonstrated in Sect. 4.2 below.

Numerical Results
In this section, we compute spatially quasi-periodic solutions of the initial value problem (2.43) with k = 1/ √ 2 and g normalized to 1. First we validate the traveling wave computations of Wilkening and Zhao (2021) and compare the accuracy and efficiency of the ERK and ETD schemes in a convergence study. We then consider more complex dynamics in which some of the wave peaks overturn. Here we plotted every 30th step of the 5th order, 6 stage explicit Runge-Kutta method of Dormand and Prince (Hairer et al. 2000), so the Runge-Kutta stepsize was t = 0.1/30. The initial condition is plotted with a thick blue line, and the wave travels right at constant speed c = 1.552197 in physical space. The solution is plotted over the representative interval 0 ≤ x ≤ 12π , though it extends in both directions to ±∞ without exactly repeating. Panel (b) shows the error in time-stepping this traveling wave solution from t = 0 to t = 3 using the 5th and 8th order explicit Runge-Kutta methods of Dormand and Prince (Hairer et al. 2000), the 4th order ETD scheme of Cox and Matthews (2002), and the 5th order ETD scheme of Whalen et al. (2015). These errors compare the numerical solution from time-stepping the initial condition to the exact formula of how a quasi-periodic traveling wave should evolve under (2.43) and (2.44), which is worked out in Wilkening and Zhao (2021). If the initial wave profile has the torus representationη 0 (α 1 , α 2 ), we defineξ 0 = H [η 0 ] andφ 0 = cξ 0 . By construction,η 0 is an even function of α = (α 1 , α 2 ), soξ 0 andφ 0 are odd. The exact traveling solution is thenη t) is a periodic function on T 2 defined implicitly by (B.7) below. To derive (4.1), one makes use of the change of variables formula (B.8) from physical space, where the wave speed is constant, to conformal space; see Wilkening and Zhao (2021). Note that the waves in (4.1) do not change shape as they move through the torus in the direction k, but the traveling speed α 0 (t) in conformal space varies in time in order to maintainξ(0, 0, t) = 0 via (2.44).

Traveling Waves
The errorssfig plotted in panel (b) is the discrete norm at the final time computed, T = 3: The surface tension in this example (τ = 1.410902) is high enough that once the stepsize is sufficiently small for the Runge-Kutta methods to be stable, roundoff error dominates truncation error. So the errors suddenly drop from very large values (10 20 or more) to machine precision. By contrast, the error in the ETD methods decreases steadily as the stepsize is reduced, indicating that the small-scale decomposition introduced in (3.8) is successful in removing stiffness from the equations of motion (Hou et al. 1994(Hou et al. , 2001. A contour plot ofη(α 1 , α 2 , T ) at t = T = 3 is shown in panel (c) of Fig. 1. The dashed line shows the trajectory from t = 0 to t = 3 of the wave crest that begins at (0, 0) and continues along the path α 1 = α 0 (t), α 2 = kα 0 (t). We use Newton's method to solve the implicit equation (B.7) for A(x, 0) at each point of the pseudo-spectral grid. We then use FFTW to compute the 2D Fourier representation of A(x, 0), which can then be used to quickly evaluate the function at any point. We find that the Fourier modes of A(x, 0) decay to machine precision on the M × M grid with M = 60, corroborating the assertion in Theorem B.2 below that A(x 1 , x 2 , t) is real analytic in x 1 and x 2 .

Overturning Waves
Next we present a spatially quasi-periodic water wave computation in which some of the wave peaks overturn as they evolve while others do not. Conformal mapping methods have been used previously to compute overturning waves. For example, Dyachenko and Newell (2016) use this approach to study whitecapping in the ocean and Wang et al. (2013) use it to compute solitary and periodic overturning traveling flexural-gravity waves. The novelty of our work is the computation of a spatially quasi-periodic water wave in which every wave peak evolves differently, and only some of them overturn.
Since torus functions are involved, the number of degrees of freedom is squared, leading to a large-scale computation. For simplicity, we set the surface tension parameter, τ , to zero.
We first seek spatially periodic dynamics in which the initial wave profile has a vertical tangent line that overturns when evolved forward in time and flattens out when evolved backward in time. Through trial and error, we selected the following parametric curves for the initial wave profile and velocity potential of this auxiliary periodic problem: Note that ξ 1 (σ ) = 0 when σ ∈ π + 2π Z, and otherwise ξ 1 (σ ) > 0. Thus, vertical tangent lines occur where ξ 1 (σ ) ∈ π + 2π Z and η 1 (σ ) = −0.5 cos(1.4π) = 0.154508; see Fig. 2.
As shown in Figure 2, the initial conditions η 2 (α) and have the desired property that the wave overturns when evolved forward in time and flattens out when evolved backward in time. In other words, the wave becomes less steep in the neighborhood of the initial vertical tangent line when time is reversed. However, it does not evolve backward to a flat state. Instead, a secondary wave crest forms to the right of the initial wave crest and grows in amplitude as t decreases. This secondary wave crest resembles the early stages of the fluid jets that were observed by Aurther et al. (2019) to form in the wave troughs when the initial condition η 0 (x) = 1 3 sin x + 1 6 sin 2x + 1 3 sin 3x) is evolved from rest in the graph-based formulation (2.1).
The blue markers in panels (a) and (b) of Fig. 2 show the location of the vertical tangent line in physical space at t = 0. The blue marker in panel (c) shows the corresponding point in conformal space. When the wave overturns for t > 0 in physical space, it is because α → ξ(α, t) no longer increases monotonically. Indeed, we see in panel (c) that η(α, t) remains single-valued as a function of the conformal variable α The initial condition (t = 0) is shown in blue in each plot but becomes very steep. This causes the Fourier mode amplitudes in panel (d) to decay more slowly as t increases. We used different mesh sizes and timesteps in the regions 0 ≤ t ≤ 0.3, 0.3 ≤ t ≤ 0.45 and 0 ≥ t ≥ −0.45 to maintain spectral accuracy; details are given below when discussing the quasi-periodic calculation. The drop-off in |η j | from 10 −14 to 10 −18 as j approaches the Nyquist frequency M/2 is due to the 1D version of the filter (3.10), which is applied after each timestep. Floating point errors of size 10 −14 occur in the discretization of the equations of motion while errors of size 10 −18 are due to having computed the inverse FFT of the filtered data to get back to real space before taking the FFT again to plot the Fourier data.
We turn the solution of this auxiliary periodic problem into a spatially quasi-periodic solution by defining initial conditions on the torus of the form where q is a free parameter that we choose heuristically to be q = 0.6kπ = 1.3329 in order to make the first wave crest to the right of the origin behave similarly to the periodic 1D solution of Fig. 2. (This will be explained below). The results of the quasi-periodic calculation are summarized in Figs. 3 and 4. Panel (a) of Fig. 3 shows snapshots of the solution at t = ( /6)T for 0 ≤ ≤ 6 over the range 0 ≤ ξ(α) ≤ 16π , where T = 0.225. The initial wave profile, ζ 0 (α) = ξ 0 (α) + iη 0 (α) with η 0 (α) =η 0 (α, kα), is plotted with a thick blue line. The wave profile is plotted with a thick black line at t = T and with thin grey lines at intermediate times. Panel Snapshots in time of a spatially quasi-periodic water wave with a periodic initial wave profile with vertical tangent lines at ξ = π + 2π n, n ∈ Z. A quasi-periodic initial velocity potential causes some of the peaks to overturn for t > 0 while others do not. Panels (a) and (f) show η(α, t) and ϕ(α, t) versus ξ(α, t) over 0 ≤ x ≤ 16π and 0 ≤ t ≤ T = 0.225. Panels (b-e) show the results of panel (a) in more detail. The blue arrows show the direction of travel of the wave at various locations (b) zooms in on the first wave in panel (a), which overturns as the wave crest moves up and right while the wave trough moves down and left, as indicated by the blue arrows. This is very similar (by design) to the forward evolution of the auxiliary periodic wave of Fig. 2, with initial conditions η 2 (α), ϕ 2 (α). Panels (c) and (d) zoom in on two other wave crests from panel (a) that flatten out (rather than overturn) as t advances from 0 to T . Panel (e) shows another type of behavior in which the wave overturns due to the wave trough moving down and left faster than the wave crest moves down and left. Panel (f) shows the evolution of the velocity potential ϕ(α, t) over 0 ≤ t ≤ T . Unlike η 0 (α), the initial velocity potential ϕ 0 (α) =φ 0 (α, kα) is not 2π -periodic due to the factor of cos(α 2 − q) in (4.5).
Panels (a)  The rapid dropoff inη(α 1 , α 2 , t) over the window 0.6π ≤ α 1 ≤ 0.667π persists from the initial state in whichη 0 (α 1 , α 2 ) does not depend on α 2 . Panel (e) shows the exponential decay of Fourier modes with respect to the shell index s at different times in panels (b) and (c). Initially,η(α 1 , α 2 , 0) depends only on α 1 ; however, by t = T , the dependence on α 2 is clearly visible. Although the waves overturn in some places when η(α, t) =η(α, kα, t) is plotted parametrically versus ξ(α, t) with t > 0 held fixed, bothη andφ are single-valued functions of α 1 and α 2 at all times. Nevertheless, throughout the evolution,η(α 1 , α 2 , t) has a steep dropoff over a narrow range of values of α 1 . Initially,η 0 (α 1 , α 2 ) = η 2 (α 1 ) = η 1 (α 1 + B 2 (α 1 )) and the rapid dropoff occurs for α 1 near the solution of α 1 + B 2 (α 1 ) = π (since the vertical tangent line occurs at ξ 1 (σ ) + iη 1 (σ ) with σ = π ). Using Newton's method, we find that this occurs at α 1 = 0.634185π . The blue curve in panel (c) of Fig. 2 gives η 2 (α). If one zooms in on this plot, one finds that η 2 (α) decreases rapidly by more than half its crest-to-trough height over the narrow range 0.6π ≤ α ≤ 0.667π . At later times, η(α 1 , α 2 , t) continues to drop off rapidly when α 1 traverses this narrow range in spite of the dependence on α 2 . This can be seen in panel (b) of Fig. 4, where there is a high clustering of nearly vertical contour lines separating the yellow-orange region from the blue region. Over this narrow window,φ(α 1 , α 2 , t) also varies rapidly with respect to α 1 . Many gridpoints are needed to resolve these rapid variations with spectral accuracy. Although ξ 1 (σ ), η 1 (σ ) and ϕ 1 (σ ) involve only a few nonzero Fourier modes, conformal reparametrization via (4.3) vastly increases the Fourier content of the initial condition. We used M = 6144 gridpoints to evolve the periodic auxiliary problem of Fig. 2 from t = 0 to t = 0.3 using the 8th order Runge-Kutta method of Dormand and Prince (Hairer et al. 2000) with stepsize t = 2.08333 × 10 −5 . We then switched to M = 12288 gridpoints to evolve from t = 0.3 to t = 0.45 with t = 7.5 × 10 −6 . In the reverse direction, we used M = 4096 gridpoints to evolve from t = 0 to t = −0.45 with t = −4.6875 × 10 −5 . Studying the Fourier modes in panel (d) of Fig. 2, it appears that 4096 gridpoints (2048 modes) are sufficient to maintain double-precision accuracy forward or backward in time to t = ±0.225. Using this as a guideline for the quasi-periodic calculation, we evolved (2.43) on a 4096 × 4096 spatial grid using the 8th order explicit Runge-Kutta method described in Sect. 3. The calculation involved 5400 time steps from t = 0 to t = T = 0.225, which took 2.5 days on 12 threads running on a server with two 3.0 GHz Intel Xeon Gold 6136 processors. Additional threads had little effect on the running time as the FFT calculations require a lot of data movement relative to the number of floating point operations involved.
Panel (e) of Fig. 4 shows the 2 average of the Fourier mode amplitudes |η j 1 , j 2 | in each shell of indices satisfying max(| j 1 |, | j 2 |) = s for 1 ≤ s ≤ 2048. Sincê η − j 1 ,− j 2 =η j 1 , j 2 , we can discard half the modes and sweep through the lattice along straight lines from (0, s) to (s, s) to (s, −s) to (1, −s), which sweeps out 4s index pairs. (The same ordering is used to enumerate the unknowns in the nonlinear least squares method proposed in Wilkening and Zhao (2021) to compute quasi-periodic traveling water waves.) We see in panel (e) that as time increases, the modes continue to decay at an exponential rate with respect to s, but the decay rate is slower at later times. The rapid dropoff in the mode amplitudes for s ≥ 1536 is due to the Fourier filter. At the final time t = T = 0.225, the modes still decay by 12 orders of magnitude from s = 1 to s = 1536, so we believe the solution is correct to 10-12 digits. A finer grid would be required to maintain this accuracy over longer times. As in Fig. 2, the additional drop-off in the amplitude of Fourier modes from s = 1536 to s = 2048 in panel (e) is due applying the filter (3.10) to the solution after every timestep.
Beyond monitoring the decay of Fourier modes, as an additional check of accuracy, we compute the average energy, mass and momentum of the solution of Fig. 4 (4.6) which may be shown to be conserved quantites under the free-surface Euler equations following the derivations in Zakharov et al. (2002); Dyachenko et al. (2019). The only changes required for the quasi-periodic case are that derivatives and the Hilbert transform are replaced by their torus versions, and the integrals over R or T are replaced by integrals over T 2 . We also divide by (2π) 2 to obtain average values over the torus. This scaling has the advantage that E, M and P x do not suddenly jump by a factor of 2π when periodic functions are viewed as quasi-periodic functions that depend on α 1 only. The equations of motion are Hamiltonian (Zakharov 1968;Dyachenko et al. 2019) whether or not the energy is scaled by 1/(2π) 2 -one just has to multiply the symplectic 2-form (Abraham et al., 1988) by the same factor.
The numerical results in Table 1 show that energy is conserved to a relative error of 1.5 × 10 −13 /0.0975 = 1.5 × 10 −12 ; mass is conserved to a relative error of The values of P x and the digits highlighted in bold fluctuate due to discretization and floating-point errors 3.3 × 10 −13 /0.0464 = 7.1 × 10 −12 ; and momentum is conserved to an absolute error of 4.6 × 10 −15 over the course of the numerical computation of Fig. 4. This gives further evidence thatη andφ are accurate to 10-12 digits. The mass being negative is an artifact of the choice of ξ 1 (σ ) and η 1 (σ ) in (4.2). While η 1 (σ ) has zero mean as a function of σ , its average value with respect to x is (2π) −1 2π 0 η 1 (σ )ξ 1 (σ ) dσ = −(3π/40)( √ 5 − 1) = −0.04635254915624. If we had added a constant toη to make M zero initially, it would have remained zero up numerical errors, similar to P x , which is initially zero sinceη andφ in (4.5) are independent of α 2 except for the factor of cos(α 2 − q). The constant value of the energy would also change if a constant were added toη.
The rationale for setting q = 0.6kπ in (4.5) is that cos(α 2 − q) ≈ 1 where the characteristic line (α, kα) crosses the dropoff in the torus near α 1 = 0.6π for the first time. Locally, ϕ 0 (α) =φ 0 (α, kα) is close to ϕ 2 (α), the initial condition of the auxiliary periodic problem, so we expect the quasi-periodic wave to evolve similarly to the periodic wave near x = π for a short time. (Here z = x + iy describes physical space). This is indeed what happens, which may be seen by comparing panel (b) of Fig. 2 to panel (b) of Fig. 3, keeping in mind that t ∈ [−0.45, 0.45] in the former plot and t ∈ [0, 0.225] in the latter plot. Advancing α from 0.6π to 10.6π causes the characteristic line (α, kα) to cross a periodic image of the dropoff at α 2 = 10.6kπ , where cos(α 2 − q) = −0.9752 ≈ −1. Locally, ϕ 0 (α) is close to −ϕ 2 (α), the initial condition of the time-reversed auxiliary periodic problem. Thus, we expect the quasiperiodic wave to evolve similarly to the time-reversed periodic wave near x = 11π . (Recall that ξ 0 (0.634185π) = π , so ξ 0 (10.634185π) = 11π ). Comparing panel (b) of Fig. 2 to panel (d) of Fig. 3 confirms that this does indeed happen. At most wave peaks, the velocity potential of the quasi-periodic solution is not closely related to that of the periodic auxiliary problem since the cosine factor is not near a relative maximum or minimum, where it is flat. As a result, the wave peaks of the quasi-periodic solution evolve in many different ways as α varies over the real line.

Conclusion
In this work, we have formulated the two-dimensional, infinite depth gravity-capillary water wave problem in a spatially quasi-periodic, conformal mapping framework. We developed two time-stepping strategies for solving the quasi-periodic initial value problem, an explicit Runge-Kutta method and an exponential time differencing scheme. We numerically verified a result in Wilkening and Zhao (2021) that quasiperiodic traveling waves evolve in time on the torus T 2 in the direction (1, k) without changing form, though their speed is non-uniform in conformal space if the conditionξ(0, 0, t) = 0 is imposed via (2.44). We then performed a convergence study to demonstrate the effectiveness of the small-scale decomposition at removing stiffness from the evolution equations when the surface tension is large. Finally, we presented the results of a large-scale computation of a spatially quasi-periodic overturning water wave for which the wave peaks exhibit a wide array of dynamic behavior.
In the appendices, we establish minimal conditions to ensure that a quasi-periodic analytic function z(w) maps the lower half-plane topologically onto a region bounded above by a curve, and that 1/|z w | is bounded. We also show that if all the solutions in the family are single-valued and have no vertical tangent lines, the corresponding solutions of the original graph-based formulation (2.1-2.4) of the Euler equations are quasi-periodic in physical space. This analysis includes a change of variables formula from conformal space to physical space, which to our knowledge is also new in the periodic setting. We then provide details on implementing the exponential timedifferencing scheme and discuss generalizations to quasi-periodic water waves at the surface of a 3D fluid, where conformal mapping methods are no longer applicable.
We believe that spatial quasi-periodicity is a natural setting to study the dynamics of linear and nonlinear waves, and has largely been overlooked as a possible third option to the usual modeling assumption that the solution either evolves on a periodic domain or decays at infinity. In the future, we plan to develop numerical methods to compute temporally quasi-periodic water waves with n = 3 quasi-periods and to study subharmonic instabilities (Longuet-Higgins 1978;MacKay and Saffman 1986;Mercer and Roberts 1992;Deconinck and Oliveras 2011;Trichtchenko et al. 2016) of periodic traveling waves and standing waves as well as the long-time dynamics of spatially quasi-periodic perturbations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. We decompose n( a , z 0 ) = n 2 (z 0 , a) − n 1 (z 0 , a), where n 1 (z 0 , a). We will show that Thus, ζ(α)− z 0 does not cross the principal branch cut of the logarithm. As a result, Since |ζ(α) − α| ≤ M < A, as a → ∞ we obtain A similar argument shows that n 1 (i A) = 1/2. Let [z 0 , z 1 ] denote the line segment in C connecting two points z 0 and z 1 , and suppose this line segment does not intersect the curve ζ(α). We claim that if the limit defining n 1 (z 0 ) exists, the limit defining n 1 (z 1 ) also exists, and n 1 (z 1 ) = n 1 (z 0 ). First note that Since ζ(α) does not cross [z 0 , z 1 ], ζ(α)−z 1 ζ(α)−z 0 does not cross the negative real axis or 0. Thus, Taking the limit as a → ∞ gives n 1 (z 1 ) = n 1 (z 0 ) + 0, as claimed. Every point of is connected to −i A by a polygonal path that remains inside , and every point in v is connected to i A by a polygonal path that remains inside v . The result (A.4) follows.
Next consider n 2 (z 0 , a) in (A.3). Since |z(w) − w| ≤ M, the Cauchy integral formula gives |z thus the modulus of the integrand in the formula for n 2 (z 0 , a) in (A.3) is bounded uniformly by 2(1 + M/ε) for large a. For fixed θ ∈ (−π, 0), z (ae iθ ) → 1 and [z(ae iθ ) − z 0 ]/a → e iθ , so the integrand approaches i pointwise on the interior of the integration interval as a → ∞. By the dominated convergence theorem, lim a→∞ n 2 (z 0 , a) = 1/2. Combining these results, we find that lim a→∞ n( a , z 0 ) = 1/2 − n 1 (z 0 ). But since n( a , z 0 ) is constant for a > |z 0 | + M, we conclude that This shows that when solving the equation z(w) = z 0 , if z 0 ∈ , there is precisely one solution w 0 in the lower half-plane, and if z 0 ∈ v , there are no solutions w 0 in the lower half-plane. Since z(w) is an open mapping, it cannot map a point in the lower half-plane to the boundary , since a nearby point would then have to be mapped to v . It follows that z(w) is a 1-1 mapping of {Im w < 0} onto . It is then a standard result that z (w) has no zeros in the lower half-plane and the inverse function w(z) exists and is analytic on .

Appendix B. Quasi-Periodic Families of Solutions
In this appendix we explore the effect of introducing phases in the reconstruction of one-dimensional quasi-periodic solutions of (2.43) from solutions of the torus version of these equations. This ultimately makes it possible to show that if all the solutions in the family are single-valued and have no vertical tangent lines, the corresponding solutions of the original graph-based formulation (2.1-2.4) of the Euler equations are quasi-periodic in physical space.
Here we note that a phase shift does not affect the mean of a periodic function on the torus, i.e. P 0 [S θũ ] = P 0 [ũ] where S θ [ũ](α) =ũ(α + θ ). We have assumed that in the 1D calculation, C 1 is chosen to agree with that of the torus calculation.
Since C 1 only affects the tangential velocity of the interface parametrization, it can be specified arbitrarily. Left-multiplying (2.38) by (−η α , ξ α ) eliminates C 1 and yields the kinematic condition (2.35). Since the Bernoulli Eq. (2.42) holds on the torus, it also holds in the 1D calculation, as claimed.
For each solution in the family (B.1), there are many others that represent identical dynamics up to a spatial phase shift or α-reparametrization. Changing δ merely shifts the solution in physical space. In fact, δ does not appear in the equations of motion (2.38) -it is only used to reconstruct the curve via (B.2). The relations show that shifting α by α 0 leads to another solution already in the family. This shift reparametrizes the curve but has no effect on its evolution in physical space. If we identify two solutions that differ only by a spatial phase shift or α-reparametrization, the parameters (θ 1 , θ 2 , δ) become identified with (0, θ 2 − kθ 1 , 0). Every solution is therefore equivalent to one of the form Within this smaller family, two values of θ lead to equivalent solutions if they differ by 2π(n 1 k + n 2 ) for some integers n 1 and n 2 . This equivalence is due to solutions "wrapping around" the torus with a spatial shift, (B.5) Here n 2 is chosen so that 0 ≤ θ + 2π(n 1 k + n 2 ) < 2π and we used periodicity of ζ(α, t ; θ 1 , θ 2 , δ) with respect to θ 1 and θ 2 . It usually suffices to restrict attention to α ∈ [0, 2π) by making use of (B.5). One exception is determining whether the curve self-intersects. In that case it is more natural to tile the plane with periodic copies of the torus and consider the straight line parametrization of (B.4). Indeed, it is conceivable that ζ(α, t; 0, θ, 0) = ζ(β, t; 0, θ, 0) (B.6) with |α − β| as large as 2M, where M is a bound on |ζ | over T 2 , and the condition (B.6) becomes hard to understand if (B.5) is used to map α and β back to [0, 2π) with different choices of n 1 or n 2 . We now show that η phys (x, t ; θ 1 , θ 2 , δ) and ϕ phys (x, t ; θ 1 , θ 2 , δ) can be defined and computed easily from ζ(α, t ; θ 1 , θ 2 , δ) and ϕ(α, t ; θ 1 , θ 2 ) if all of the waves in the family (B.4) are single-valued and have no vertical tangent lines, and that η phys and ϕ phys are quasi-periodic functions of x. To simplify notation, let α = (α 1 , α 2 ), x = (x 1 , x 2 ) and k = (1, k).
defines a unique function A(x 1 , x 2 , t) on T 2 that is periodic and real analytic in x 1 and x 2 . The inverse of the change of variables x = α + kξ(α, t) on T 2 is given by Proof First we check that if A satisfies (B.7), then (B.8) is the inverse of the change of variables x = α + kξ(α, t). Given x ∈ T 2 , define α by (B.8). Then (B.9) as required. Next we show existence and uniqueness of a solution A of (B.7) under the assumed hypotheses. Given α = (α 1 , α 2 ) ∈ T 2 , the definition (B.1) gives where the left-hand side means (d/dα) α=α 1 ξ(α, t ; 0, α 2 − kα 1 , 0). We know the right-hand side is periodic and continuous on T 2 while the left-hand is positive on the primitive cell {(α 1 , α 2 ) : 0 ≤ α 1 < 2π , kα 1 ≤ α 2 < kα 1 + 2π }. Therefore, both sides of (B.10) are bounded below by some ε(t) > 0 that does not depend on α ∈ T 2 . Let M(t) be a bound on |ξ(α, t)| over T 2 . Then for fixed x ∈ R 2 (with t also fixed), the function g(α) = g(α; x, t) = α +ξ(x + kα, t) is strictly monotonically increasing on R (as g (α) ≥ ε(t)) and satisfies g(−M(t)) ≤ 0 and g(M(t)) ≥ 0. Thus, we can define A(x, t) as the unique solution of g(α) = 0. It follows that |A(x, t)| ≤ M(t). If n 1 and n 2 are integers, replacing x in (B.7) by y = (y 1 , y 2 ) = (x 1 +2π n 1 , x 2 +2π n 2 ) and using periodicity ofξ(α, t) gives Since the solution of this equation is unique, A( y, t) = A(x, t). This shows that A(x, t) is periodic in x, and hence well-defined on T 2 . It is also real analytic, which follows from the implicit function theorem, noting that g(α; x 1 , x 2 , t) is real analytic in α, x 1 and x 2 for fixed t and ∂ g/∂α is never zero. For the same reason, A(x, t) will depend as smoothly on t asξ(α, t) does.

Appendix C. Details on Implementing the Exponential Time Differencing Schemes
In this section we summarize how to solve the evolution equations (3.8-3.9) using the 4-stage 4th order ETD scheme of Cox and Matthews (2002) or the 6-stage 5th order ETD scheme of Whalen et al. (2015). When using an s-stage ETD scheme to solve the ODE u t = Lu + N (t, u), (C.1) the numerical solution u n is advanced from t n to (t n + h) via U r = e c r hL u n + h s j=1 a r j (hL)N j , The Butcher array for the Cox-Matthews scheme may be written 0 0 (C.4) The Butcher array for the Whalen-Brio-Moloney scheme is given in Table 2 of Whalen et al. (2015). Both of these schemes are explicit methods, i.e. A(z) is strictly lowertriangular. Moreover, c r = s j=1 A r j (0). It follows that c 1 = 0, U 1 = u n , and for r ≥ 2, the upper limit of the sum in the formula for U r can be replaced by r − 1 so that only previously computed values of N j are needed to evaluate U r . The only implicit ETD schemes we are aware of are the high-order collocation methods of Chen and Wilkening (2021); Aurther et al. (2019).
Let V denote the space of real-valued functions defined on a uniform M 1 × M 2 grid overlaid on the torus T 2 via (3.3). When solving (3.8-3.9), the state space for (C.1) consists of state vectors u = (η,φ) ∈ V 2 . The nonlinear function N (u) in (3.9) does not explicitly depend on time and is computed using the pseudo-spectral approach described in Sect. 3, i.e. derivatives and the Hilbert transform are computed in Fourier space while the quadratic nonlinearities in (3.9) are evaluated pointwise on the grid. Evaluation of e c r hL u n , a r j (hL)N j , e hL u n and b r (hL)N r in (C.2) are all of the form ψ(hL)v where ψ(z) is an entire function involving the φ k (z) functions and v is a state vector. Since the two-dimensional FFT diagonalizes L (see below), the numerical instabilities discussed by Kassam and Trefethen (2005) are avoided by using the series expansion of (C.4) for |z| ≤ 1, which is ∞ l=k (k!/l!)z l−k . Replacing the upper limit by k + 19 is sufficient to achieve double-precision accuracy. There is no catastrophic cancellation of digits since the leading terms of e z are eliminated analytically before evaluating the series numerically. It is therefore not necessary to use the contour integral approach advocated in Kassam and Trefethen (2005). For |z| > 1, one can just evaluate (k!/z k )(e z − k−1 l=0 z l /l!) as written, or use the scaling and modified squaring algorithm of Skaflestad and Wright (2009); see also (Sison et al. 2021;Chen and Wilkening 2021).
Next we explain how to compute ψ(hL)v. We temporarily allow functions in V to be complex-valued and regard it as a complex vector space. In the end, only realvalued functions in V will actually arise in the calculation. Let F denote the 2D FFT, and consider the spaceV = l 2 (J ) of functionsf taking valuesf j 1 , j 2 on the discrete Fourier lattice (C.5) The operator F = diag(F, F) is an isomorphism of the state space V 2 ontoV 2 , and L is block-diagonalized by F. In more detail L = F −1 SF, where S leaves invariant the two-dimensional Fourier subspaces (V 2 ) j = span{e + j , e − j } ⊂V 2 , where j ∈ J and Here 0 ∈V is the zero function (0 l = 0 for l ∈ J ), and δ j ∈V is a lattice version of the Kronecker delta, i.e. (δ j ) l = 1 if l = j and 0 if l ∈ J \ { j }. The restriction of S to (V 2 ) j has the following matrix representation with respect to the basis {e + j , e − j }: j , otherwise. (C.7) When j = (0, 0), S j is the zero matrix, so it is diagonalized by Let and E be the operators onV 2 that leave the subspaces (V 2 ) j invariant and have matrix representations j and E j with respect to {e + j , e − j }. Then Q = F −1 E diagonalizes L and hence ψ(hL): (C.9) The functions ψ(z) that arise in (C.2) all have the property that ψ(z) = ψ(z). Let v ∈ V 2 be a real-valued state vector and denote the intermediate steps in computing ψ(hL)v as F (C.10) We denote the components ofv in the e ± j basis byv ± j , with similar notation forŵ,x, andŷ. Since v is real-valued and a − j = a j , b − j = b j in (C.7), we havê v − j =v j , (C.11) Inspecting the formulas in (C.8), we see that conjugating the inputs of E −1 j gives the conjugates of the outputs in the opposite order. (When computing b = Ax, we refer to the components of x as inputs to A and those of b as outputs.) It then follows from (C.11) that Conjugating and reversing the order of the inputs to E j gives the conjugates of the outputs in the original order. Thus,ŷ − j =ŷ j (C.13) and ψ(hL)v is real. This also justifies using the 'r2c' version of the two-dimensional FFT to computev = Fv, which only returns values of v j with j = ( j 1 , j 2 ) ∈ J and j 1 ≥ 0. We then only computeŵ,x andŷ for these values of j . The missing entries are assured to satisfy (C.13), which is the assumption needed to apply the 'c2r' version of the inverse FFT to obtain ψ(hL)v = F −1ŷ . The most expensive steps of evaluating ψ(hL)v are the FFTs in F and F −1 . To implement (C.2), one has to apply F to u n and N 1 , . . . , N s and apply F −1 to obtain U 2 , . . . , U s and u n+1 . Since F and F −1 each involve 2 FFT's, the functional calculus steps of (C.2) involve a total of 4s + 2 FFT's. Meanwhile, evaluating N 1 , . . . , N s using the pseudo-spectral method to compute derivatives and Hilbert transforms in (3.9) involves 10s FFTs. These would have to be computed using a Runge-Kutta method anyway, so the cost of an s-stage ETD method is approximately 40% higher than that of an s-stage RK method with the same stepsize h. But as seen in Fig. 1, significantly larger steps can often be taken with the ETD method for a given accuracy goal, making the ETD method more efficient in spite of the additional cost per step.
We remark that we use the formulas for a j and b j in (C.7) even for the Nyquist modes j = ( j 1 , j 2 ) with j 1 = M 1 /2 or j 2 = M 2 /2. We avoid setting a j = 0 and b j = g for these modes, which would have been consistent with our treatment of H , ∂ 2 α and P in the definition of L in (3.8), as it would lead to a Jordan block in the diagonalization of L. It makes little difference since the Nyquist modes are intended to remain close to roundoff-level values throughout the computation, and in fact are set to zero at the end of each timestep by the filter (3.10). In general, if a Jordan block arises in the diagonalization of L and the corresponding eigenspace contains "lowfrequency" modes that have to be resolved to get an accurate result, one can transfer a term from L to N in the decomposition (C.1) so that the modified L is diagonalizable. This technique is demonstrated in Aurther et al. (2019).
where C(t) is an arbitrary function of time, ∇ x = ∂ x 1 , ∂ x 2 T , and is the mean curvature. We have also introduced the Dirichlet-Neumann operator (Craig and Sulem 1993), where (x 1 , x 2 , y, t) is the solution of (D.6) Note that t could be dropped from the notation in ϕ, η and when defining G as time is frozen when solving the auxiliary problem of reconstructing the velocity potential in the fluid from its boundary value ϕ on the free surface η and computing the scaled normal derivative (D.5). Using x i = ϕ x i − η x i y in (D.5), we see that (D.3) is equivalent to which can also be derived easily from the 3D analog of (2.4). Following the definitions in Dynnikov and Novikov (2005), we consider quasiperiodic functions on R 2 with d quasi-periods, which have the form u(x) =ũ(K x) = j ∈Z dû j e i j T K x , x = (x 1 ; x 2 ) ∈ R 2 , K ∈ Mat d×2 (R).
Hereũ : R d → R is real analytic and periodic (so well-defined on T d ),û j are its Fourier modes, and a semicolon separates entries of a column vector. The rows of K may be assumed to be linearly independent over Z since otherwise a new functionṽ : T d−1 → R and matrix L ∈ Mat (d−1)×2 (R) can be constructed so that u(x) =ṽ(Lx). Indeed, if the rows of K are linearly dependent over Z, there is a unimodular integer matrix J such that J K = (L; 0), i.e. the last row of J K contains only zeros and the first d − 1 rows define L. One can then definev l = l d ∈Zû j(l,l d ) and confirm that u(x) = l∈Z d−1vl e i l T Lx , where j (l, l d ) = J T (l 1 ; . . . ; l d−1 ; l d ) for l ∈ Z d−1 and l d ∈ Z. One also finds thatṽ(y) =ũ J −1 (y; 0) for y ∈ T d−1 . We do not know a reference for these calculations, but they are straightforward. The procedure can be repeated until the minimal d is found. For the wave to be genuinely two-dimensional and quasi-periodic, we require K to have full rank and d ≥ 3. For example, one choice of K when d = 3 is Just as in Remark 2.4, substitution of the quasi-periodic functions η(x, t) =η(K x, t), ϕ(x, t) =φ(K x, t) (D.8) into (D.2-D.4) gives evolution equations forη(x, t) andφ(x, t) on the torus T d . One just has to replace ∇ x by (k 1 · ∇x; k 2 · ∇x), where k 1 , k 2 ∈ R d are the columns of K andx = (x 1 ; . . . ;x d ) ∈ T d . This amounts to using the chain rule, D x η = Dxη · K , where D x = (∂ x 1 , ∂ x 2 ) = (∇ x ) T . It is also necessary to define a quasi-periodic Dirichlet-Neumann operator viaG Givenη andφ in (D.8) with t fixed, one would construct the solution˜ (x, y, t) of (k 1 · ∇x) 2˜ + (k 2 · ∇x) 2˜ + (∂ y ) 2˜ = 0, −∞ < y <η(x, t),x ∈ R d , (x,η(x, t), t) =φ(x, t), y =η(x, t), ∂ y˜ → 0, y → −∞ (D.10) and then evaluatẽ G(η)φ = ∂ y˜ − (k 1 · ∇xη)(k 1 · ∇x˜ ) − (k 2 · ∇xη)(k 2 · ∇x˜ ) y=η .
Whereas the conformal mapping approach of Sect. 2.5 employs a quasi-periodic Hilbert transform to efficiently compute the Dirichlet-Neumann operator for a 2D fluid, it remains an open problem to devise and implement an efficient method for solving (D.10) for a 3D fluid with quasi-periodic boundary conditions. In a finite-depth variant of the problem, the finite element approach of Rycroft and Wilkening (2013) and the Transformed Field Expansion method of Reitich (2001, 2006); Qadeer and Wilkening (2019) are candidate approaches that have been used successfully for periodic boundary conditions. Both approaches could work in principle for quasiperiodic boundary conditions but will suffer from the curse of dimensionality as the (d + 1)-dimensional region η,h,t = (x; y) : −h < y <η(x, t) ,x ∈ T d has to be discretized, where h is the fluid depth. We do not know what to expect for the condition number of a linear system that discretizes (D.10), which is elliptic with respect to (x; y) on quasi-periodic slices through η,h,t but not with respect to (x; y) on η,h,t . The infinite depth case can often be dealt with by introducing a transparent boundary condition along a fictitious interface at some depth y = −h below the free surface, discretizing the fluid region η,h,t above the interface, and using series expansions in the unbounded region below the interface (Nicholls and Reitich 2006). We have not worked out the details in a quasi-periodic setting.
A final option is to avoid the Dirichlet-Neumann operator altogether by using a weakly nonlinear model water wave equation. One could introduce torus versions of the equations of motion to obtain genuine quasi-periodic solutions rather than solving a system of coupled single-mode nonlinear Schrödinger equations (Bridges and Laine-Pearson 2001;Ablowitz and Horikis 2015). The curse of dimensionality will also be a challenge for weakly nonlinear models with this torus approach, especially if larger values of d are considered.