2.1. Quasi-periodic functions and the Hilbert transform

A function $u(\alpha )$ is quasi-periodic if there exists a continuous, periodic function $\tilde u(\boldsymbol \alpha )$ defined on the $d$-dimensional torus $\mathbb{T}^d$ such that

(2.1a,b)\begin{equation} u(\alpha) = \tilde u(\boldsymbol{k} \alpha), \quad \tilde u(\boldsymbol\alpha) = \sum_{\boldsymbol{j}\in\mathbb{Z}^d}\hat{u}_{\boldsymbol{j}} \,\textrm{e}^{\textrm{i}\langle\,\boldsymbol{j},\,\boldsymbol{\alpha}\rangle}, \quad \alpha\in\mathbb{R}, \ \boldsymbol\alpha,\boldsymbol k \in \mathbb{R}^d. \end{equation}

We generally assume $\tilde u(\boldsymbol \alpha )$ is real analytic, which means the Fourier modes satisfy the symmetry condition $\hat u_{-\boldsymbol j}=\overline {\hat u_{\boldsymbol j}}$ and decay exponentially as $|\,\boldsymbol j|\rightarrow \infty$, i.e. $|\hat u_{\boldsymbol j}|\le M\, \textrm {e}^{-\sigma |\,\boldsymbol j|}$ for some $M,\sigma >0$. Entries of the vector $\boldsymbol {k}$ are required to be linearly independent over $\mathbb {Z}$. Fixing this vector $\boldsymbol k$, we define two versions of the Hilbert transform, one acting on $u$ (the quasi-periodic version) and the other on $\tilde u$ (the torus version)

(2.2a,b)\begin{equation} H[u](\alpha) = \frac{1}{\rm \pi}\operatorname{PV} \int_{-\infty}^\infty\frac{u(\xi)}{\alpha-\xi}\,\textrm{d} \xi, \quad H[\tilde u](\boldsymbol\alpha) = \sum_{\boldsymbol{j}\in\mathbb{Z}^d} ({-}i)\operatorname{sgn}(\langle\,\boldsymbol{j},\boldsymbol{k}\rangle) \hat{u}_{\boldsymbol{j}}\, \textrm{e}^{\textrm{i} \langle\,\boldsymbol{j},\,\boldsymbol{\alpha}\rangle}. \end{equation}

Here, $\operatorname {sgn}(q)\in \{1,0,-1\}$ 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(\mathbb{T}^d)$ that depends on $\boldsymbol k$. It is shown in Wilkening & Zhao (Reference Wilkening and Zhao2020) that

(2.3)\begin{equation} H[u](\alpha)=H[\tilde u](\boldsymbol k\alpha), \end{equation}

and 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

(2.4)\begin{equation} f(w) = \hat u_{\boldsymbol 0} + i\hat v_{\boldsymbol 0} + \sum_{\langle\,\boldsymbol j,\boldsymbol k\rangle<0} 2 \hat u_{\boldsymbol j}\,\textrm{e}^{\textrm{i}\langle\,\boldsymbol j,\boldsymbol k\rangle w}, \quad (w=\alpha+i\beta, \beta\le0), \end{equation}

where $\hat v_{\boldsymbol 0}$ is an arbitrary constant and the sum is over all $\boldsymbol j\in \mathbb{Z}^d$ satisfying $\langle \,\boldsymbol j,\boldsymbol k\rangle <0$. The imaginary part of $f$ on the real axis is then given by $v=\hat v_{\boldsymbol 0}-H[u]$. Similarly, given $v$, the most general bounded analytic function $f(w)$ in the lower half-plane whose imaginary part agrees with $v$ on the real axis has the form (2.4) with $u=\hat u_{\boldsymbol 0} + H[v]$, where $\hat u_{\boldsymbol 0}$ is an arbitrary constant. This analytic extension is quasi-periodic on slices of constant depth, i.e.

(2.5)\begin{equation} f(w) = \tilde f(\boldsymbol k\alpha,\beta), \quad (w=\alpha+i\beta,\beta\le0), \end{equation}

where $\tilde f(\boldsymbol \alpha ,\beta ) = \hat u_{\boldsymbol 0} + i\hat v_{\boldsymbol 0} + \sum _{\langle \,\boldsymbol j,\boldsymbol k\rangle <0} 2[\hat u_{\boldsymbol j}\, \textrm {e}^{-\langle \,\boldsymbol j,\boldsymbol k\rangle \beta }]\, \textrm {e}^{\textrm {i}\langle \,\boldsymbol j,\boldsymbol \alpha \rangle }$ is periodic in $\boldsymbol \alpha$ for fixed $\beta \le 0$. The torus version of the bounded analytic extension corresponding to $\tilde u(\boldsymbol \alpha +\boldsymbol \theta )$ is simply $\tilde f(\boldsymbol \alpha +\boldsymbol \theta ,\beta )$, which has imaginary part $\tilde v(\boldsymbol \alpha +\boldsymbol \theta )$ on the real axis. As a result, the Hilbert transform commutes with the shift operator,

(2.6)\begin{equation} H[\tilde u(\boldsymbol{\cdot}+\boldsymbol\theta)](\boldsymbol\alpha) = H[\tilde u](\boldsymbol\alpha+\boldsymbol\theta), \end{equation}

which can also be checked directly from (2.2a,b). We also define quasi-periodic and torus versions of two projection operators,

(2.7)\begin{equation} P = \operatorname{id} - P_0, \quad P_0 [u] = P_0 [\tilde u] = \hat{u}_{\boldsymbol{0}} = \frac{1}{(2{\rm \pi})^d} \int_{\mathbb{T}^d} \tilde u(\boldsymbol{\alpha})\, \textrm{d} \alpha_1\ldots \, \textrm{d} \alpha_d, \end{equation}

where $P_0[u]$ is a constant function on $\mathbb{R}$, $P_0[\tilde u]$ is a constant function on $\mathbb{T}^d$ and $P[u]$ has zero mean on $\mathbb{R}$ in the sense that its torus representation, $P[\tilde u]$, which satisfies $P[u](\alpha )=P[\tilde u](\boldsymbol k\alpha )$, has zero mean on $\mathbb{T}^d$.

2.2. A quasi-periodic conformal mapping

For the general initial value problem (Wilkening & Zhao Reference Wilkening and Zhao2020), we consider a time-dependent conformal map $z(w,t)$ that maps the lower half-plane

(2.8)\begin{equation} \mathbb{C}^{-} = \{w=\alpha+i\beta: \alpha\in\mathbb{R},\beta<0\} \end{equation}

to the fluid domain $\varOmega _f(t)$ that lies below the free surface in physical space. At each time $t$, we assume $z(w,t)$ extends continuously to $\overline {\mathbb{C}^-}$, and in fact is analytic on a slightly larger half-plane $\mathbb{C}^-_\varepsilon =\{w: \operatorname {Im} w<\varepsilon \}$, where $\varepsilon >0$ could depend on $t$. The free surface $\varGamma (t)$ is parameterized by

(2.9)\begin{equation} \zeta(\alpha,t)=\xi(\alpha,t)+i\eta(\alpha,t), \quad (\alpha\in\mathbb{R}, t \,{fixed}), \quad \zeta = z\vert_{\beta=0}. \end{equation}

We assume $\alpha \mapsto \zeta (\alpha ,t)$ is injective but do not assume $\varGamma (t)$ is the graph of a single-valued function of $x$ in the derivation. An example of a time-dependent spatially quasi-periodic overturning water wave is computed by Wilkening & Zhao (Reference Wilkening and Zhao2020). In future work we will study travelling quasi-periodic perturbations of the overhanging periodic travelling water waves computed by Akers et al. (Reference Akers, Ambrose and Wright2014).

The conformal map is required to remain a bounded distance from the identity map in the lower half-plane. Specifically, we require that

(2.10)\begin{equation} |z(w,t)-w|\le M(t) \quad (w=\alpha+i\beta, \beta\le0), \end{equation}

where $M(t)$ is a uniform bound that could vary in time. The Cauchy integral formula implies that $|z_w-1|\le M(t)/|\beta |$, so at any fixed time,

(2.11)\begin{equation} z_w\rightarrow1 \text{ as } \beta\to-\infty. \end{equation}

In this paper and its companion (Wilkening & Zhao Reference Wilkening and Zhao2020), we assume $\eta$ has two spatial quasi-periods, i.e. at any time it has the form (2.1a,b) with $d=2$ and $\boldsymbol {k}=[k_1,k_2]^T$. This is a major departure from previous work (Meiron, Orszag & Israeli Reference Meiron, Orszag and Israeli1981; Dyachenko et al. Reference Dyachenko, Kuznetsov, Spector and Zakharov1996a; Zakharov et al. Reference Zakharov, Dyachenko and Vasilyev2002; Dyachenko, Lushnikov & Korotkevich Reference Dyachenko, Lushnikov and Korotkevich2016), where $\eta$ is assumed to be periodic. Through non-dimensionalization, we may set $k_1=1$ and $k_2 = k$, where $k$ is irrational

(2.12a,b)\begin{equation} \eta(\alpha,t) = \tilde\eta(\alpha,k\alpha,t), \quad \tilde\eta(\alpha_1,\alpha_2,t) = \sum_{j_1, j_2\in\mathbb{Z}} \hat{\eta}_{j_1, j_2}(t)\,\textrm{e}^{\textrm{i}(\,j_1\alpha_1+j_2\alpha_2)}. \end{equation}

Here, $\hat {\eta }_{-j_1, -j_2}(t) = \overline {\hat {\eta }_{j_1,j_2}(t)}$ since $\tilde \eta (\alpha _1,\alpha _2,t)$ is real valued. Since $w\mapsto [z(w,t)-w]$ is bounded and analytic on $\mathbb{C}^-$ and its imaginary part agrees with $\eta$ on the real axis, there is a real number $x_0$ (possibly depending on time) such that

(2.13a,b)\begin{equation} \xi(\alpha, t) = \alpha + x_0(t) + H[\eta](\alpha, t), \quad \xi_\alpha(\alpha, t) = 1 + H[\eta_\alpha](\alpha, t). \end{equation}

We use a tilde to denote the periodic functions on the torus that correspond to the quasi-periodic parts of $\xi$, $\zeta$ and $z$,

(2.14)\begin{equation} \left. \begin{gathered} \xi(\alpha,t) = \alpha + \tilde\xi(\alpha,k\alpha,t), \quad \zeta(\alpha,t) = \alpha + \tilde\zeta(\alpha,k\alpha,t), \\ z(\alpha+i\beta,t)=(\alpha+i\beta )+ \tilde z(\alpha,k\alpha,\beta,t), \quad (\beta\le 0). \end{gathered} \right\} \end{equation}

Specifically, $\tilde \xi = x_0(t) + H[\tilde \eta ]$, $\tilde \zeta = \tilde \xi + i\tilde \eta$, and

(2.15)\begin{equation} \tilde z(\alpha_1,\alpha_2,\beta,t) = x_0(t)+i\hat\eta_{0,0}(t)+\sum_{j_1+j_2k<0} (2i\hat\eta_{j_1,j_2}(t)\, \textrm{e}^{-(\,j_1+j_2k)\beta}) \, \textrm{e}^{\textrm{i}(\,j_1\alpha_1+j_2\alpha_2)}. \end{equation}

Since the modes $\hat \eta _{j_1,j_2}$ are assumed to decay exponentially, there is a uniform bound $M(t)$ such that $|\tilde z(\alpha _1,\alpha _2,\beta ,t)|\le M(t)$ for $(\alpha _1,\alpha _2)\in \mathbb{T}^2$ and $\beta \le 0$. Wilkening & Zhao (Reference Wilkening and Zhao2020) show that as long as the free surface $\zeta (\alpha ,t)$ does not self-intersect at a given time $t$, the mapping $w\mapsto z(w,t)$ 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 Wilkening & Zhao (Reference Wilkening and Zhao2020), let $\varPhi ^{phys}(x,y,t)$ denote the velocity potential in physical space and let $W^{phys}(x+iy,t)= \varPhi ^{phys}(x,y,t)+i\varPsi ^{phys}(x,y,t)$ denote the complex velocity potential, where $\varPsi ^{phys}$ is the streamfunction. Using the conformal mapping $z(w,t)$, we pull back these functions to the lower half-plane and define

(2.16)\begin{equation} W(w,t) = \varPhi(\alpha,\beta,t)+i\varPsi(\alpha,\beta,t) = W^{phys}(z(w,t),t), \quad (w=\alpha+i\beta). \end{equation}

We also define

(2.17a,b)\begin{equation} \varphi=\varPhi\vert_{\beta=0}, \quad \psi=\varPsi\vert_{\beta=0}. \end{equation}

We assume $\varphi$ is quasi-periodic with the same quasi-periods as $\eta$,

(2.18a,b)\begin{equation} \varphi(\alpha, t) = \tilde\varphi(\alpha,k\alpha,t), \quad \tilde\varphi(\alpha_1,\alpha_2,t) = \sum_{j_1, j_2\in\mathbb{Z}} \hat{\varphi}_{j_1, j_2}(t) \, \textrm{e}^{\textrm{i}(\,j_1\alpha_1+j_2\alpha_2)}. \end{equation}

The fluid velocity $\boldsymbol {\nabla }\varPhi ^{phys}(x,y,t)$ is assumed to decay to zero as $y\rightarrow -\infty$ (since we work in the laboratory frame). Since $\textrm {d} W/ \textrm {d} w=(\textrm {d} W^{phys}/\textrm {d} z)(\textrm {d} z/ \textrm {d} w)$, it follows from (2.11) that $\textrm {d} W/ \textrm {d} w\rightarrow 0$ as $\beta \to -\infty$. Thus,

(2.19a,b)\begin{equation} \psi_\alpha ={-}H[\varphi_\alpha], \quad \psi(\alpha, t) ={-}H[\varphi](\alpha, t). \end{equation}

Here, we have assumed $P_0[\varphi ] = \hat {\varphi }_{0,0}(t) = 0$ and $P_0[\psi ] = \hat {\psi }_{0,0}(t) = 0$, which is allowed since $\varPhi$ and $\varPsi$ can be modified by additive constants (or functions of time only) without affecting the fluid motion.

Let $U$ and $V$ denote the normal and tangential velocities of the curve parameterization, respectively; let $s_\alpha =|\zeta _\alpha |=(\xi _\alpha ^2+\eta _\alpha ^2)^{1/2}$ denote the rate at which arclength increases as the curve $\alpha \mapsto \zeta (\alpha ,t)$ is traversed; and let $\theta$ denote the tangent angle of the curve relative to the horizontal. Then

(2.20a,b)\begin{equation} \zeta_\alpha = s_\alpha\, \textrm{e}^{\textrm{i}\theta}, \quad \zeta_t = (V+iU)\,\textrm{e}^{\textrm{i}\theta}. \end{equation}

Tracking a fluid particle $x_p(t)+iy_p(t)= \zeta (\alpha _p(t),t)$ on the free surface, we find that

(2.21a,b)\begin{equation} \dot x_p=\xi_\alpha\dot\alpha_p+\xi_t=\varPhi^{phys}_x, \quad \dot y_p=\eta_\alpha\dot\alpha_p+\eta_t=\varPhi^{phys}_y. \end{equation}

Eliminating $\dot \alpha _p$ gives the kinematic condition

(2.22)\begin{equation} U = \zeta_t\boldsymbol{\cdot}\hat{\boldsymbol n} = \boldsymbol{\nabla}\varPhi^{phys}\boldsymbol{\cdot}\hat{\boldsymbol n}, \end{equation}

where $\hat {\boldsymbol n}=(-\eta _\alpha ,\xi _\alpha )/s_\alpha$ is the outward unit normal to $\varGamma$ and we have identified $\zeta _t$ with the vector $(\xi _t,\eta _t)$ in $\mathbb{R}^2$. The general philosophy proposed by Hou, Lowengrub & Shelley (HLS) (Hou et al. Reference Hou, Lowengrub and Shelley1994, Reference Hou, Lowengrub and Shelley1997) is that while (2.22) 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_\alpha (t)$ independent of $\alpha$, we interpret the work of Dyachenko et al. (Reference Dyachenko, Kuznetsov, Spector and Zakharov1996a) and subsequent authors (Choi & Camassa Reference Choi and Camassa1999; Dyachenko Reference Dyachenko2001, Reference Dyachenko2019; Zakharov et al. Reference Zakharov, Dyachenko and Vasilyev2002) as choosing $V$ to maintain a conformal representation. Briefly, since $\varPhi ^{phys}$ and $\varPsi ^{phys}$ satisfy the Cauchy–Riemann equations, we have

(2.23)\begin{equation} -\frac{\psi_\alpha}{s_\alpha} ={-}\frac{\varPsi^{phys}_x\xi_\alpha + \varPsi^{phys}_y\eta_\alpha}{s_\alpha} = \frac{\varPhi^{phys}_y\xi_\alpha - \varPhi^{phys}_x\eta_\alpha}{s_\alpha} = \boldsymbol{\nabla}\varPhi^{phys}\boldsymbol{\cdot} \hat{\boldsymbol n} = U. \end{equation}

Assuming $z_t/z_\alpha$ is bounded and analytic in the lower half-plane (justified below),

(2.24)\begin{equation} \left. \frac{z_t}{z_\alpha}\right\vert_{\beta=0} = \frac{\zeta_t}{\zeta_\alpha} = \frac{V+iU}{s_\alpha} \Rightarrow \frac{V}{s_\alpha} = H\left(\frac{U}{s_\alpha}\right)+C_1 ={-}H\left(\frac{\psi_\alpha}{s_\alpha^2}\right)+C_1, \end{equation}

where $C_1$ is an arbitrary constant (in space) that we are free to choose. For any differentiable function $\alpha _0(t)$, replacing $C_1(t)$ by $C_1(t)-\alpha _0'(t)$ will cause a reparameterization of the solution with $\alpha$ replaced by $\alpha -\alpha _0(t)$; see Appendix A. The tangential and normal velocities can be rotated back to obtain $\xi _t$ and $\eta _t$ via

(2.25)\begin{equation} \begin{pmatrix} \xi_t \\ \eta_t \end{pmatrix} = \begin{pmatrix} \xi_\alpha & -\eta_\alpha \\ \eta_\alpha & \xi_\alpha \end{pmatrix} \begin{pmatrix} V/s_\alpha \\ U/s_\alpha \end{pmatrix}, \end{equation}

which can be interpreted as the real and imaginary parts of the complex multiplication $\zeta _t = (\zeta _\alpha )(\zeta _t/\zeta _\alpha )$. As explained in Wilkening & Zhao (Reference Wilkening and Zhao2020), the first equation of (2.25) is automatically satisfied if the second equation holds and $\xi$ is reconstructed from $\eta$ via (2.13a,b), provided $x_0(t)$ satisfies

(2.26)\begin{equation} \frac{\textrm{d}\kern0.7pt x_0}{\textrm{d} t} = P_0\left[ \xi_\alpha \frac V{s_\alpha} - \eta_\alpha \frac U{s_\alpha} \right]. \end{equation}

The equations of motion for water waves in the conformal framework may now be written

(2.27) \begin{equation} \left. \begin{gathered} \xi_\alpha = 1 + H[\eta_\alpha], \quad \psi ={-}H[\varphi], \quad J = \xi_\alpha^2 + \eta_\alpha^2, \quad \chi = \dfrac{\psi_\alpha}{J}, \\ \text{choose }C_1 \text{ (see below), compute } \dfrac{\textrm{d}\kern0.7pt x_0}{\textrm{d} t} \text{ in (2.26)} \text{ if necessary},\\ \eta_t ={-}\eta_\alpha H[\chi] - \xi_\alpha\chi + C_1\eta_\alpha, \quad \kappa = \dfrac{\xi_\alpha\eta_{\alpha\alpha} - \eta_\alpha\xi_{\alpha\alpha}}{J^{3/2}}, \\ \varphi_t = P\left[\dfrac{\psi_\alpha^2 - \varphi_\alpha^2}{2J} - \varphi_\alpha H[\chi] + C_1\varphi_\alpha - g\eta + \tau\kappa\right], \end{gathered} \right\} \end{equation}

where the last equation comes from the unsteady Bernoulli equation and the Laplace–Young condition for the pressure. These equations govern the evolution of $x_0$, $\eta$ and $\varphi$. The full curve $\zeta =\xi +i\eta$ and its analytic extension $z$ to the lower half-plane can be reconstructed from $\eta$ using (2.13a,b) and (2.14). Doing so ensures that $z$ is injective and that $z_t/z_\alpha$ remains bounded in the lower half-plane provided that the free surface does not self-intersect and $J$ remains non-zero on the surface; see Wilkening & Zhao (Reference Wilkening and Zhao2020) for details.

As noted by Wilkening & Zhao (Reference Wilkening and Zhao2020), (2.27) can be interpreted as an evolution equation for the functions $\tilde \zeta (\alpha _1,\alpha _2,t)$ and $\tilde \varphi (\alpha _1,\alpha _2,t)$ on the torus $\mathbb{T}^2$. The $\alpha$-derivatives are replaced by the directional derivatives $[\partial _{\alpha _1}+k\partial _{\alpha _2}]$ and the quasi-periodic Hilbert transform is replaced by its torus version, i.e. $H[\tilde u]$ in (2.2a,b) above. The pseudo-spectral method proposed in Wilkening & Zhao (Reference Wilkening and Zhao2020) is based on this representation. A convenient choice of $C_1$ is

(2.28)\begin{equation} C_1 = \left[H\left(\frac{\tilde\psi_\alpha}{\tilde J}\right) - \frac{\tilde\eta_\alpha\tilde\psi_\alpha}{ (1+\tilde\xi_\alpha)\tilde J} \right]_{(\alpha_1,\alpha_2)=(0,0)}, \end{equation}

which causes $\tilde \xi (0,0,t)$ to remain constant in time, alleviating the need to evolve $x_0(t)$ explicitly. Here $\tilde J=(1+\tilde \xi _\alpha )^2 + \tilde \eta _\alpha ^2$, and all instances of $\xi _\alpha$ in (2.27) must be replaced by

(2.29)\begin{equation} \widetilde{\xi_\alpha} = 1 + \tilde\xi_\alpha, \end{equation}

since the secular growth term $\alpha$ is not part of $\tilde \xi$ in (2.14). Using (2.13a,b) and (2.14), $\tilde \zeta$ is completely determined by $x_0(t)$ and $\tilde \eta$, so only these have to be evolved – the formula for $\tilde \xi _t$ in (2.25) is redundant as long as (2.26) is satisfied. Other choices of $C_1$ are considered in Appendix A.

Wilkening & Zhao (Reference Wilkening and Zhao2020) show that solving the torus version of (2.27) yields a three-parameter family of one-dimensional solutions of the form

(2.30)\begin{equation} \begin{aligned} \zeta(\alpha,t; \theta_1,\theta_2,\delta) & = \alpha + \delta + \tilde\zeta( \theta_1+\alpha,\theta_2+k\alpha,t), \\ \varphi(\alpha,t; \theta_1,\theta_2) & = \tilde\varphi( \theta_1+\alpha,\theta_2+k\alpha,t), \end{aligned} \quad \left( \begin{aligned} \alpha\in\mathbb{R}, \ t\ge0 \\ \theta_1,\theta_2,\ \delta\in\mathbb{R} \end{aligned}\right). \end{equation}

The parameters $(\theta _1,\theta _2,\delta )$ lead to the same solution as $(0,\theta _2-k\theta _1,0)$ up to a spatial phase shift and $\alpha$-reparameterization. Thus, every solution is equivalent to one of the form

(2.31)\begin{equation} \begin{aligned} \zeta(\alpha,t; 0,\theta,0) & = \alpha + \tilde\zeta(\alpha,\theta+k\alpha,t), \\ \varphi(\alpha,t; 0,\theta) & = \tilde\varphi(\alpha,\theta+k\alpha,t) \end{aligned} \quad \alpha\in\mathbb{R}, \ t\ge 0,\ \theta\in[0,2{\rm \pi}). \end{equation}

Within this smaller family, two values of $\theta$ lead to equivalent solutions if they differ by $2{\rm \pi} (n_1k+n_2)$ for some integers $n_1$ and $n_2$. This equivalence is due to solutions ‘wrapping around’ the torus with a spatial shift,

(2.32)\begin{align} \zeta(\alpha+2{\rm \pi} n_1,t; 0,\theta,0) = \zeta(\alpha,t; 0,\theta+2{\rm \pi}(n_1k+n_2),2{\rm \pi} n_1), \quad (\alpha\in[0,2{\rm \pi}),\ n_1\in\mathbb{Z}). \end{align}

Here, $n_2$ is chosen so that $0\le [\theta +2{\rm \pi} (n_1k+n_2)]<2{\rm \pi}$ and we used periodicity of $\zeta (\alpha ,t; \theta _1,\theta _2,\delta )$ with respect to $\theta _1$ and $\theta _2$. Wilkening & Zhao (Reference Wilkening and Zhao2020) also show that if all the waves in the family (2.31) 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 (Zakharov Reference Zakharov1968; Craig & Sulem Reference Craig and Sulem1993; Johnson Reference Johnson1997) that are quasi-periodic in physical space.

2.4. Quasi-periodic travelling water waves

We now specialize to the case of quasi-periodic travelling waves and derive the equations of motion in a conformal mapping framework. One approach (see e.g. Milewski et al. (Reference Milewski, Vanden-Broeck and Wang2010) for the periodic case) is to write down the equations of motion in a graph-based representation of the surface variables $\eta ^{phys}(x,t)$ and $\varphi ^{phys}(x,t)=\varPhi ^{phys}(x,\eta (x,t),t)$ and substitute $\eta ^{phys}_t = -c\eta ^{phys}_x$, $\varphi ^{phys}_t = -c\varphi ^{phys}_x$ to solve for the initial condition of a solution of the form

(2.33a,b)\begin{equation} \eta^{phys}(x,t)=\eta^{phys}_0(x-ct), \quad \varphi^{phys}(x,t)=\varphi^{phys}_0(x-ct). \end{equation}

We present below an alternative derivation of the equations in Milewski et al. (Reference Milewski, Vanden-Broeck and Wang2010) that is more direct and does not assume the wave profile is single valued. Other systems of equations have also been derived to describe travelling water waves, e.g. by Nekrasov (Reference Nekrasov1921), Milne-Thomson (Reference Milne-Thomson1968) and Dyachenko et al. (Reference Dyachenko, Lushnikov and Korotkevich2016).

Recall the kinematic condition (2.23) that the normal velocity of the curve is given by $\zeta _t\boldsymbol {\cdot }\hat {\boldsymbol n} = U = -\psi _\alpha /s_\alpha$. Since the wave travels at constant speed $c$ in physical space, there is a reparameterization $\beta (\alpha ,t)$ such that $\zeta (\alpha ,t)=\zeta (\beta (\alpha ,t),0)+ct$. Since $\zeta _\alpha$ is tangent to the curve, the normal velocity is simply $\zeta _t\boldsymbol {\cdot }\hat {\boldsymbol n}=(c,0)\boldsymbol {\cdot }\hat {\boldsymbol n}=-c\eta _\alpha /s_\alpha$, where we used $\hat {\boldsymbol n}=(-\eta _\alpha ,\xi _\alpha )/s_\alpha$. We conclude that

(2.34a,b)\begin{equation} \psi_\alpha = c\eta_\alpha, \quad \varphi_\alpha = H[\psi_\alpha] = cH[\eta_\alpha] = c(\xi_\alpha -1). \end{equation}

This expresses $\psi$ and $\varphi$ (up to additive constants) in terms of $\eta$ and $\xi =\alpha +x_0+H[\eta ]$, leaving only $\eta$ to be determined. As in the graph-based approach of (2.33a,b) above, it suffices to compute the initial wave profile at $t=0$ to know the full evolution of the travelling wave under (2.27); however, the wave generally travels at a non-uniform speed in conformal space in order to travel at constant speed in physical space. This is demonstrated in § 4.2 and proved in Appendix A.

The two-dimensional velocity potential $\varPhi ^{phys}(x,y,t)$ may be assumed to exist even if the travelling wave possesses overhanging regions that cause the graph-based representation via $\eta ^{phys}(x,t)$ and $\varphi ^{phys}(x,t)$ to break down. In a moving frame travelling at constant speed $c$ along with the wave, the free surface will be a streamline. Let $\breve{z}=z-ct$ denote position in the moving frame and note that the complex velocity potential picks up a background flow term, $\breve{W}^{phys}(\breve{z},t)=W^{phys}(\breve{z}+ct,t)-c\breve{z}$, and becomes time independent. We drop $t$ in the notation and define $\breve{W}(w)=\breve{W}^{phys}(\breve{z}(w))$, where $\breve{z}(w)=z(w,0)$ conformally maps the lower half-plane onto the fluid region of this stationary problem. We assume $W^{phys}(\breve{z}(\alpha ),0)$ is quasi-periodic with exponentially decaying mode amplitudes, so

(2.35) \begin{equation} |\breve{W}(w)+cw|\le|W^{phys}(\breve{z}(w),0)|+c|\breve{z}(w)-w|, \end{equation}

is bounded in the lower half-plane. Since the streamfunction $\operatorname {Im}\{\breve{W}^{phys}(\breve{z})\}$ is constant on the free surface, we may assume $\operatorname {Im}\{\breve{W}(\alpha )\}=0$ for $\alpha \in \mathbb{R}$. The function $\operatorname {Im}\{\breve{W}(w)+cw\}$ is then bounded and harmonic in the lower half-plane and satisfies homogeneous Dirichlet boundary conditions on the real line, so it is zero (Axler, Bourdon & Ramey Reference Axler, Bourdon and Ramey1992). Up to an additive real constant,

(2.36)\begin{equation} \breve{W}(w)={-}cw. \end{equation}

Thus, $|\breve{\boldsymbol {\nabla }}\breve{\varPhi} ^{phys}|^2=|\breve{W}'(w)/\breve{z}'(w)|^2=c^2/J$. Since the free surface is a streamline in the moving frame, the steady Bernoulli equation $(1/2)|\breve{\boldsymbol {\nabla }}\breve{\varPhi}^{phys}|^2+g\eta +p/\rho =C$ together with the Laplace–Young condition $p=p_0-\rho \tau \kappa$ on the pressure gives

(2.37)\begin{equation} \left. \begin{gathered} \xi_\alpha = 1 + H[\eta_\alpha], \quad J = \xi_\alpha^2 + \eta_\alpha^2, \\ \kappa = \dfrac{\xi_\alpha\eta_{\alpha\alpha} - \eta_\alpha\xi_{\alpha\alpha}}{J^{3/2}}, \quad P\left[\dfrac{c^2}{2J} + g\eta -\tau\kappa\right] = 0, \end{gathered}\right\} \end{equation}

which is the desired system of equations for $\eta$.

In the quasi-periodic travelling wave problem, we seek a solution of (2.37) of the form (2.12a,b), except that $\tilde \eta$ and its Fourier modes will not depend on time. Like the initial value problem, (2.37) can be interpreted as a nonlinear system of equations for $\tilde \eta (\alpha _1,\alpha _2)$ defined on $\mathbb{T}^2$, where the $\alpha$-derivatives are replaced by $[\partial _{\alpha _1} + k\partial _{\alpha _2}]$ and the Hilbert transform is replaced by its torus version in (2.2a,b). Without loss of generality, we assume

(2.38)\begin{equation} \hat{\eta}_{0,0} = 0. \end{equation}

We also assume that $\tilde \eta$ is an even, real function of $(\alpha _1,\alpha _2)$ on $\mathbb{T}^2$. Hence, in our set-up, the Fourier modes of $\tilde \eta$ satisfy

(2.39a,b)\begin{equation} \hat{\eta}_{{-}j_1, -j_2} = \bar{\hat{\eta}}_{j_1, j_2}, \quad \hat{\eta}_{{-}j_1, -j_2} = \hat{\eta}_{j_1, j_2}, \quad (\,j_1, j_2)\in \mathbb{Z}^2. \end{equation}

This implies that all the Fourier modes $\hat {\eta }_{j_1, j_2}$ are real, and causes $\eta (\alpha )=\tilde \eta (\alpha ,k\alpha )$ to be even as well, which is compatible with the symmetry of (2.37). However, as in (2.30), there is a larger family of quasi-periodic travelling solutions embedded in this solution, namely

(2.40)\begin{equation} \eta(\alpha;\theta) = \tilde\eta(\alpha,\theta+k\alpha). \end{equation}

As in (2.32), two values of $\theta$ lead to equivalent solutions (up to $\alpha$-reparameterization and a spatial phase shift) if they differ by $2{\rm \pi} (n_1k+n_2)$ for some integers $n_1$ and $n_2$. In general, $\eta (\alpha -\alpha _0;\theta )$ will not be an even function of $\alpha$ for any choice of $\alpha _0$ unless $\theta =2{\rm \pi} (n_1k+n_2)$ for some integers $n_1$ and $n_2$. In the periodic case, symmetry breaking travelling water waves have been found by Zufiria (Reference Zufiria1987), though most of the literature is devoted to periodic travelling waves with even symmetry.

2.5. Linear theory of quasi-periodic travelling waves

Linearizing (2.37) around the trivial solution $\eta (\alpha ) = 0$, we obtain,

(2.41)\begin{equation} c^2H[\delta{\eta}_\alpha] - g\delta{\eta} + \tau\delta{\eta}_{\alpha\alpha} = 0, \end{equation}

where $\delta \eta$ denotes the variation of $\eta$. Substituting (1.2a,b) into (2.41), we obtain a resonance relation for the Fourier modes of $\delta \eta$

(2.42)\begin{equation} (c^2|\,j_1k_1+j_2k_2| - g - \tau(\,j_1k_1+j_2k_2)^2) \widehat{\delta\eta}_{j_1,j_2} = 0, \quad (\,j_1, j_2) \in \mathbb{Z}^2. \end{equation}

Note that, $j_1k_1+j_2k_2$, which appears in the exponent of the Fourier plane wave representation (1.2a,b), plays the role of $k$ in the dispersion relation (1.1). Many families of quasi-periodic travelling wave solutions bifurcate from the trivial solution even after specifying $k_1$ and $k_2$. Selecting a branch amounts to choosing two of the modes $\widehat {\delta \eta }_{j_1,j_2}$ to bring in at linear order and setting the others to zero in (2.42). In this paper, we focus on the case in which $\hat \eta _{1,0}$ and $\hat \eta _{0,1}$ enter linearly. This gives the first-order resonance conditions

(2.43a,b)\begin{equation} c^2k_1 - g - \tau k_1^2 = 0, \quad c^2k_2 - g - \tau k_2^2 = 0, \end{equation}

where $k_1=1$ and $k_2=k$ in our non-dimensionalized setting. For right-moving waves, we then have $c=\sqrt {g/k_1 + g/k_2}$ and $\tau =g/(k_1k_2)$. Any superposition of waves with dimensionless wave numbers $k_1=1$ and $k_2=k$ travelling with speed $c=c_{lin}$ will solve the linearized problem (2.41) for $\tau =\tau _{lin}$. Here, we have introduced the notation $c_{lin}=\sqrt {g+g/k}$ and $\tau _{lin}=g/k$ to facilitate the discussion of nonlinear effects below.

4.1. Spatially quasi-periodic travelling waves

We now present a detailed numerical study of solutions of (2.37) with $k=1/\sqrt 2$ and $g=1$ on three continuation paths corresponding to $\gamma \in \{5,1,0.2\}$, where $\gamma =\hat \eta _{1,0}/\hat \eta _{0,1}$ is the amplitude ratio of the prescribed base modes. In each case, we vary the larger of $\hat \eta _{1,0}$ and $\hat \eta _{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.13), which by (3.10) corresponds to

(4.1a–c)\begin{equation} p^{(0)}_1=\tau^{(0)}=\sqrt2, \quad p^{(0)}_3=b^{(0)}=1+\sqrt2, \quad p^{(0)}_j=0, \quad j\not\in\{1,3\}. \end{equation}

As noted already, the amplitudes $\hat \eta _{1,0}$ and $\hat \eta _{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=24$ for the Fourier cutoff in each dimension, where we drop the subscripts when $M_1=M_2$ and $N_1=N_2$. The nonlinear least-squares problem involves $M^2=3600$ equations in $N_{tot}=1200$ unknowns.

Figure 1 shows the initial conditions $\eta$ and $\varphi$ for the last solution on each continuation path (with $max\{\hat \eta _{1,0}, \hat \eta _{0,1}\}=0.01$). Panels (a), (b) and (c) correspond to $\gamma =5, 1$ and $0.2$, respectively. The solution in all three cases is quasi-periodic, i.e. $\eta$ and $\varphi$ never exactly repeat themselves; we plot the solution from $x=0$ to $x=36{\rm \pi}$ as a representative snapshot. For these three solutions, the objective function $f$ in (3.11), which is a squared error, was minimized to $6.05\times 10^{-28}$, $9.28\times 10^{-28}$ and $4.25\times 10^{-28}$, respectively, with similar or smaller values for lower-amplitude solutions on each path. For each of the 30 solutions computed on these paths, only one Jacobian evaluation and 3–5 $f$ evaluations were needed to achieve roundoff-error accuracy. In our computations, $\eta$ and $\varphi$ are represented by $\tilde \eta (\alpha _1, \alpha _2)$ and $\tilde \varphi (\alpha _1, \alpha _2)$, which are defined on the torus $\mathbb {T}^2$. In figure 2, we show contour plots of $\tilde \eta (\alpha _1, \alpha _2)$ and $\tilde \varphi (\alpha _1, \alpha _2)$ corresponding to the final solution on each path. Following the dashed lines through $\mathbb{T}^2$ in figure 2 leads to the plots in figure 1. By construction in (2.39a,b), $\tilde \eta (-\boldsymbol \alpha )=\tilde \eta (\boldsymbol \alpha )$ while $\tilde \varphi (-\boldsymbol \alpha )=-\tilde \varphi (\boldsymbol \alpha )$.

Figure 1. Spatially quasi-periodic travelling solutions in the laboratory frame at $t = 0$. The wave height $\eta (\alpha )$ (solid red line) and velocity potential $\varphi (\alpha )$ (dashed blue line) are plotted parametrically against $\xi (\alpha )$ to show the wave in physical space.

Figure 2. Contour plots of $\tilde \eta$ and $\tilde \varphi$ on $\mathbb {T}^2$. The dashed lines show $(\alpha , k\alpha )$ and its periodic images with $0\le \alpha \le 10{\rm \pi}$ and $k=1/\sqrt {2}$. Evaluating $\tilde \eta$ and $\tilde \varphi$ at these points gives $\eta$ and $\varphi$ in (2.12a) and (2.18a), which were plotted in figure 1.

The amplitude ratio, $\gamma := \hat {\eta }_{1,0} / \hat {\eta }_{0,1}$, determines the bulk shape of the solution. If $\gamma \gg 1$, the component wave with wavenumber 1 will be dominant; if $\gamma \ll 1$, the component wave with wavenumber $k=1/\sqrt 2$ will be dominant; and if $\gamma$ 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 $\gamma =5$, $1$ and $0.2$ in panels (a), (b) and (c) of figure 1. Panels (a,c) show a clear dominant mode with visible variations in the amplitude. The oscillations are faster in panel (a) than in (c) since $1>k\approx 0.707$. By contrast, in panel (b), there is no single dominant wavelength.

This can also be understood from the contour plots of figure 2. In case (a), $\gamma \gg 1$ and the contour lines of $\tilde \eta$ and $\tilde \varphi$ are perturbations of sinusoidal waves depending only on $\alpha _1$. The unperturbed waves would have vertical contour lines. The $\alpha _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 $\alpha _1$ direction. Case (c) is similar, but the contour lines are elongated horizontally since the dominant variation is in the $\alpha _2$ direction. Following the dashed lines in figure 2, a cycle of $\alpha _1$ is completed before a cycle of $\alpha _2$ (since $k<1$). In case (a), a cycle of $\alpha _1$ traverses the dominant variation of $\tilde \eta$ and $\tilde \varphi$ on the torus, whereas in case (c), this is true of $\alpha _2$. So the waves in figure 1 appear to oscillate faster in case (a) than case (c). In the intermediate case (b) with $\gamma =1$, the contour lines of the crests and troughs are nearly circular, but not perfectly round. The amplitude of the waves in figure 1 are largest when the dashed lines in figure 2 pass near the extrema of $\tilde \eta$ and $\tilde \varphi$, and are smallest when the dashed lines pass near the zero level sets of $\tilde \eta$ and $\tilde \varphi$.

Next we examine the behaviour of the Fourier modes that make up these solutions. Figure 3 shows two-dimensional plots of the Fourier modes $\hat \eta _{j_1,j_2}$ for the 3 cases above, with $\gamma \in \{5,1,0.2\}$ and $max\{\hat \eta _{1,0},\hat \eta _{0,1}\}= 0.01$. Only the prescribed modes and the modes that were optimized by the solver (see (3.10)) are plotted, which have indices in the range $0\le j_1\le N$ and $-N\le j_2\le N$, excluding $j_2\le 0$ when $j_1=0$. The other modes are determined by the symmetry of (2.39a,b) and by zero padding $\hat \eta _{j_1,j_2}=0$ if $N < j_1\le M/2$ or $N<|\,j_2|\le M/2$. We used $N=24$ and $M=60$ in all 3 calculations. One can see that the fixed Fourier modes $\hat {\eta }_{1,0}$ and $\hat {\eta }_{0,1}$ are the two highest-amplitude modes in all three cases. In this sense, our solutions of the nonlinear problem (2.37) are small-amplitude perturbations of the solutions (3.13) of the linearized problem. However, in the plots of figure 3, there are many active Fourier modes other than the four modes $\textrm {e}^{\pm \textrm {i}\alpha _1}$, $\textrm {e}^{\pm \textrm {i} \alpha _2}$ from linear theory. In this sense, these solutions have left the linear regime. Carrying out a weakly nonlinear Stokes expansion to high enough order to accurately predict all these modes would be difficult due to the two-dimensional array of unknown Fourier modes, which would complicate the analysis of the periodic Wilton ripple problem (Vanden-Broeck Reference Vanden-Broeck2010; Trichtchenko et al. Reference Trichtchenko, Deconinck and Wilkening2016; Akers et al. Reference Akers, Ambrose and Sulon2019; Akers & Nicholls Reference Akers and Nicholls2020). Steeper waves that are well outside of the linear regime will be computed in § 4.3.

Figure 3. Two-dimensional Fourier modes of $\tilde \eta$ for the $k=1/\sqrt 2$ solutions plotted in figures 1 and 2: (a) $\gamma =5$; (b,d) $\gamma =1$; (c) $\gamma =0.2$. In all three cases, the modes decay visibly slower along the line $j_1+j_2k=0$, indicating the presence of resonant mode interactions.

In panels (a,b,c) of figure 3, 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 $\gamma =1$. (The other two cases are similar.) The direction along which the modes decay less rapidly appears to coincide with the line $\{(\,j_1,j_2):j_1+j_2k=0\}$, which is plotted in red. A partial explanation is that when $j_1+j_2k$ is close to zero, the corresponding modes $\textrm {e}^{\textrm {i}(\,j_1+j_2k)\alpha }$ in the expansion of $\eta (\alpha )$ in (2.12a,b) have very long wavelengths. 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.37). We believe this would lead to a small divisor problem that would complicate a rigorous proof of existence of quasi-periodic travelling water waves. Similar small divisor problems arise in proving the existence of standing water waves (Plotnikov & Toland Reference Plotnikov and Toland2001; Iooss et al. Reference Iooss, Plotnikov and Toland2005), three-dimensional travelling gravity waves (Iooss & Plotnikov Reference Iooss and Plotnikov2009) and two-dimensional time quasi-periodic gravity–capillary waves (Berti & Montalto Reference Berti and Montalto2016; Baldi et al. Reference Baldi, Berti, Haus and Montalto2018; Berti et al. Reference Berti, Franzoi and Maspero2020), where small divisors are tackled using a Nash–Moser iterative scheme.

Next we show that $\tau$ and $c$ depend nonlinearly on the amplitude of the Fourier modes $\hat {\eta }_{1,0}$ and $\hat {\eta }_{0,1}$. Panels (a) and (b) of figure 4 show plots of $\tau$ and $c$ versus $\hat \eta _{max}:=max(\hat \eta _{1,0},\hat \eta _{0,1})$ for 9 values of $\gamma =\hat \eta _{1,0}/\hat \eta _{0,1}$, namely $\gamma =0.1, 0.2, 0.5, 0.8, 1, 1.25, 2, 5, 10$. On each curve, $\hat \eta _{max}$ varies from 0 to $0.01$ in increments of $0.001$. At small amplitude, linear theory predicts $\tau =g/k=1.41421$ and $c=\sqrt {g(1+1/k)}=1.55377$. This is represented by the black marker at $\hat \eta _{max}=0$ in each plot. For each value of $\gamma$, the curves $\tau$ and $c$ are seen to have zero slope at $\hat \eta _{max}=0$, and can be concave up or concave down depending on $\gamma$. This can be understood from the contour plots of panels (e,f). Both $\tau$ and $c$ appear to be even functions of $\hat \eta _{1,0}$ and $\hat \eta _{0,1}$ when the other is held constant. Both plots have a saddle point at the origin, are concave down in the $\hat \eta _{1,0}$ direction holding $\hat \eta _{0,1}$ fixed, and are concave up in the $\hat \eta _{0,1}$ direction holding $\hat \eta _{1,0}$ fixed. The solid lines in the first quadrant of these plots are the slices corresponding to the values of $\gamma$ plotted in panels (a,b). The concavity of the one-dimensional plots depends on how these lines intersect the saddle in the two-dimensional plots.

Figure 4. Surface tension, wave speed, energy and momentum of small-amplitude quasi-periodic water waves with $k=1/\sqrt 2$. (a–d) Plots of $\tau$, $c$, $E$ and $P_x$ versus $\hat {\eta }_{max}=max\{\hat \eta _{1,0},\hat \eta _{0,1}\}$ holding $\gamma =\hat \eta _{1,0}/\hat \eta _{0,1}$ fixed. The black arrow in each plot shows how the curves change as $\gamma$ increases from $0.1$ to $10$. (e, f) Contour plots of $\tau$ and $c$ and the rays of constant $\gamma$ corresponding to (a,b). (g) Mode amplitudes of a two-dimensional Chebyshev expansion of $c(\hat \eta _{1,0},\hat \eta _{0,1})$ over the rectangle $-0.01\le \hat \eta _{1,0},\hat \eta _{0,1}\le 0.01$.

The contour plots of panels (e,f) of figure 4 were made by solving (2.37) with $(\hat \eta _{1,0},\hat \eta _{0,1})$ ranging over a uniform $26\times 26$ grid on the square $[-0.01,0.01]\times [-0.01,0.01]$. Using an even number of grid points avoids the degenerate case where $\hat \eta _{1,0}$ or $\hat \eta _{0,1}$ is zero. At those values, the two-dimensional family of quasi-periodic solutions meets a sheet of periodic solutions where $\tau$ or $c$ becomes a free parameter. Alternative techniques would be needed in these degenerate cases to determine the value of $\tau$ or $c$ from which a periodic travelling wave in the nonlinear regime bifurcates to a quasi-periodic wave. In panel (g), we plot the magnitude of the Chebyshev coefficients in the expansion

(4.2)\begin{equation} c(\hat\eta_{1,0},\hat\eta_{0,1}) = \sum_{m=0}^{15}\sum_{n=0}^{15} \hat c_{mn}T_m(100\hat\eta_{1,0})T_n(100\hat\eta_{0,1}), \quad -0.01\le \hat\eta_{1,0},\hat\eta_{0,1}\le 0.01. \end{equation}

This was done by evaluating $c$ on a cartesian product of two 16-point Chebyshev–Lobatto grids over $[-0.01,0.01]$ 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\ge 10$ or so, and only even modes $m$ and $n$ are active. The plot for $|\hat \tau _{mn}|$ is very similar, so we omit it. These plots confirm the visual observation from the contour plots that $\tau$ and $c$ are even functions of $\hat \eta _{1,0}$ and $\hat \eta _{0,1}$ when the other is held constant. These properties of $\tau$ and $c$ make them unsuitable as continuation parameters near the trivial solution, as discussed in § 3.

In panels (c,d) of figure 4, we show the energy $E$ and momentum $P_x$ of waves in the above two-parameter family of quasi-periodic solutions,

(4.3)\begin{equation} \left. \begin{array}{c@{}} \begin{aligned} E & = \int_{\mathbb{T}^2} \tfrac{1}{2} \tilde{\psi}(\partial_{\alpha_1} + k\partial_{\alpha_2})\tilde{\varphi} + \tfrac{1}{2} g\tilde{\eta}^2(1+ (\partial_{\alpha_1} + k\partial_{\alpha_2})\tilde{\xi}) \\ & \quad + \tau\left(\sqrt{(1+ (\partial_{\alpha_1} + k\partial_{\alpha_2})\tilde{\xi})^2 + ((\partial_{\alpha_1} + k\partial_{\alpha_2})\tilde{\eta})^2} - 1\right)\, \textrm{d} \alpha_1\,\textrm{d} \alpha_2, \end{aligned}\\ P_x ={-}\displaystyle\int_{\mathbb{T}^2} \tilde{\varphi}(\partial_{\alpha_1} + k\partial_{\alpha_2})\tilde{\eta}\, \textrm{d} \alpha_1 \,\textrm{d} \alpha_2. \end{array}\right\} \end{equation}

These formulas are derived by Zakharov et al. (Reference Zakharov, Dyachenko and Vasilyev2002) and Dyachenko, Lushnikov & Zakharov (Reference Dyachenko, Lushnikov and Zakharov2019) in the conformal mapping framework for a water wave of infinite depth. The only modification needed for spatially quasi-periodic waves with $d$ quasi-periods is that integrals over $\mathbb{R}$ or $\mathbb{T}$ are replaced by integrals over $\mathbb{T}^d$. Wilkening & Zhao (Reference Wilkening and Zhao2020) confirm that $E$ and $P_x$ in (4.3) are conserved quantities under the evolution equations (2.27). We see in figure 4 that the energy and momentum of the quasi-periodic waves are positively correlated. In particular, the quasi-periodic wave family with $\gamma = 1$ possesses the largest energy and momentum when $\hat {\eta }_{max}$ is fixed, even though it does not have the highest wave speed. Energy and momentum can both be regarded as measures of the amplitude of the wave. Unlike the wave speed, they are both zero at the flat rest state. We note that $\gamma =1$ corresponds to maximizing both $|\eta _{1,0}|$ and $|\eta _{0,1}|$ to have the value $\hat {\eta }_{max}$, and also leads to the largest-amplitude oscillations in figure 1. The Hamiltonian structure of the equations of motion could be useful e.g. in generalizing the time quasi-periodic results of Berti et al. (Reference Berti, Franzoi and Maspero2020) to the spatially quasi-periodic setting.

4.2. Time evolution of spatially quasi-periodic travelling waves

In this section, we confirm that the quasi-periodic solutions we obtain by minimizing the objective function (3.11) are indeed travelling waves under the evolution equations (2.27). This allows us to measure the accuracy of our independent codes for solving these two problems by comparing the numerical results. An interesting feature of the conformal mapping formulation arises in this comparison, namely that for most choices of $C_1$ in (2.27), travelling waves move at a non-uniform speed through conformal space in order to travel at constant speed in physical space. This is discussed in this section and proved in Appendix A.

In figure 5, we plot the time evolution of $\zeta (\alpha ,t)$ in the laboratory frame from $t=0$ to $t=3$. The initial conditions, plotted with thick blue lines, are those of the travelling waves computed in figures 1 and 2 above by minimizing the objective function (3.11). The grey curves give snapshots of the solution at uniformly sampled times with ${\rm \Delta} t=0.1$. They were computed using the fifth-order explicit Runge–Kutta method described by Wilkening & Zhao (Reference Wilkening and Zhao2020) with a step size of 1/300, so there are 30 Runge–Kutta steps between snapshots in the figure. The solutions are plotted over the representative interval $0\le x\le 12{\rm \pi}$, though they extend in both directions to $\pm \infty$ without exactly repeating. The initial condition and time evolution were computed on the torus and then sampled along the $(1,k)$ direction to extract the data for these one-dimensional plots.

Figure 5. Time evolution of the travelling wave profiles, $\zeta (\alpha ,t)$, from $t=0$ to $t=3$ in the laboratory frame. The thick blue lines correspond to the initial conditions.

For quantitative comparison, let $\tilde \eta _0(\boldsymbol \alpha )$ denote the initial condition on the torus, which is computed numerically by minimizing (3.11). We then compute $\tilde \xi _0=H[\tilde \eta _0]$ and $\tilde \varphi _0=c\tilde \xi _0$, which are odd functions of $\boldsymbol \alpha =(\alpha _1,\alpha _2)\in \mathbb{T}^2$ since $\tilde \eta$ is even. From Corollary A.5 of Appendix A, we define the ‘exact solution’ of the time evolution of the travelling wave under (2.27) and (2.28) with these initial conditions as

(4.4)\begin{equation} \left. \begin{gathered} \tilde\eta_{exact}(\boldsymbol\alpha,t) = \tilde\eta_0(\boldsymbol\alpha - \boldsymbol k \alpha_0(t)), \\ \tilde\varphi_{exact}(\boldsymbol\alpha,t) = \tilde\varphi_0(\boldsymbol\alpha - \boldsymbol k\alpha_0(t)), \end{gathered}\right\} \end{equation}

where $\boldsymbol k=(1,k)$, $\alpha _0(t)=ct - \mathcal{A}_0(-\boldsymbol kct)$ and $\mathcal{A}_0(x_1,x_2)$ is a periodic function on $\mathbb{T}^2$ defined implicitly by (A12) below. We see in (4.4) that the waves do not change shape as they move through the torus along the characteristic direction $\boldsymbol k$, but the travelling speed $\alpha _0'(t)$ in conformal space varies in time in order to maintain $\tilde \xi (0,0,t)=0$ via (2.28). By Corollary A.5, the exact reconstruction of $\tilde \xi _{exact}$ from $\tilde \eta _{exact}$ is

(4.5)\begin{equation} \tilde\xi_{exact}(\boldsymbol\alpha,t) = \tilde\xi_0(\boldsymbol\alpha - \boldsymbol k\alpha_0(t)) + \delta_0(t), \end{equation}

where $\delta _0(t) = ct-\alpha _0(t) = \mathcal{A}_0(-\boldsymbol kct)$ measures the deviation in position from travelling at the constant speed $ct$ in conformal space. The defining property (A12) of $\mathcal{A}_0(x_1,x_2)$ ensures that $\tilde \xi _{exact}(0,0,t)=0$.

The significance of $\mathcal{A}_0$ is that the inverse of the mapping $\boldsymbol x = \boldsymbol \alpha + \boldsymbol k\tilde \xi _0(\boldsymbol \alpha )$ on $\mathbb{T}^2$, assuming it is single valued, is

(4.6)\begin{equation} \boldsymbol \alpha = \boldsymbol x + \boldsymbol k\mathcal{A}_0(\boldsymbol x). \end{equation}

This result of Wilkening & Zhao (Reference Wilkening and Zhao2020) allows us to express quasi-periodic solutions of the initial value problem in conformal space as quasi-periodic functions in physical space. In the travelling case considered here, the exact solutions on the torus in physical space are $\tilde \eta _0^{phys}(\boldsymbol x-\boldsymbol kct)$ and $\tilde \varphi _0^{phys}(\boldsymbol x-\boldsymbol kct)$, where e.g. $\tilde \eta _0^{phys}(\boldsymbol x) = \tilde \eta _0(\boldsymbol x +\boldsymbol k \mathcal{A}_0(\boldsymbol x))$. We know this already on physical grounds, but it also follows from (4.4) and (4.5) using

(4.7a,b)\begin{equation} \tilde\eta^{phys}_{exact}(\boldsymbol x,t) = \tilde\eta_{exact}(\boldsymbol x + \boldsymbol k\mathcal{A}(\boldsymbol x,t), t), \quad \tilde\varphi^{phys}_{exact}(\boldsymbol x,t) = \tilde\varphi_{exact}(\boldsymbol x + \boldsymbol k\mathcal{A}(\boldsymbol x,t), t), \end{equation}

where $\mathcal{A}(\boldsymbol x,t) = \mathcal{A}_0(\boldsymbol x-\boldsymbol kct) - \mathcal{A}_0(-\boldsymbol kct)$ satisfies the time-dependent analogue of (A12).

Figure 6 shows contour plots of the torus version of the $\gamma =5$ and $\gamma =0.2$ solutions shown in panels (a,c) of figure 5 at the final time computed, $T=3$. A similar plot of the $\gamma =1$ solution is given in Wilkening & Zhao (Reference Wilkening and Zhao2020). The dashed lines show the trajectory from $t=0$ to $t=T$ of the wave crest that begins at $(0,0)$ and continues along the path $\alpha _1=\alpha _0(t),$ $\alpha _2=k\alpha _0(t)$ through the torus in (4.4). The following table gives the phase speed, $c$, surface tension, $\tau$, translational shift in conformal space at the final time computed, $\alpha _0(T)$, and deviation from steady motion in conformal space, $\delta _0(T)$, for these three finite-amplitude solutions (recall that $max\{\hat \eta _{1,0},\hat \eta _{0,1}\}=0.01$ and $\hat \eta _{1,0}/\hat \eta _{0,1}=\gamma$) as well as for the zero-amplitude limit

(4.8) \begin{align} \begin{array}{c|c|c|c|c} & \gamma=5 & \gamma=1 & \gamma = 0.2 & \text{linear theory} \\\hline c & \phantom{-} 1.552175 & \phantom{-} 1.552197 & \phantom{-} 1.553743 & c_{lin} = 1.553774 \\ \tau & \phantom{-} 1.409665 & \phantom{-} 1.410902 & \phantom{-} 1.415342 & \tau_{lin} = 1.414214 \\ \alpha_0(T) & \phantom{-} 4.677416 & \phantom{-} 4.681174 & \phantom{-} 4.668757 & c_{lin}T = 4.661322 \\ \delta_0(T) & -0.020890 & -0.024583 & -0.007527 & 0\end{array} \quad \begin{array}{c} \\ \\ \\ (T=3). \end{array} \end{align}

In figure 7, we plot $\delta _0(t)$ for $0\le t\le T$ (solid lines) along with $(c-c_{lin})t$ (dashed and dotted lines) for the three finite-amplitude solutions in this table. Writing $\alpha _0(t) = c_{lin}t + [(c-c_{lin})t - \delta _0(t)]$, we see that the deviation of $\alpha _0(t)$ from linear theory over this time interval is due mostly to fluctuations in $\delta _0(t)$ rather than the steady drift $(c-c_{lin})t$ due to the change in phase speed $c$ of the finite-amplitude wave.

Figure 6. Contour plots of the numerical solution $\tilde \eta (\alpha _1,\alpha _2,T)$ on the torus corresponding to the quasi-periodic solutions $\eta (\alpha ,t)$ of panels (a,c) of figure 5 at the final time shown, $t=T=3$. The dashed lines show the trajectory of the wave crest from $t=0$ to $t=T$.

Figure 7. Plots of $\delta _0(t)=ct-\alpha _0(t)$ in (4.4) and $(c-c_{lin})t$ for the solutions of figure 5.

Computing the exact solution (4.4) requires evaluating $\delta _0(t) = \mathcal{A}_0(-ct,-kct)$. We use Newton's method to solve the implicit equation (A12) for $\mathcal{A}_0(x_1,x_2)$ at each point of a uniform $M\times M$ grid, with $M_1=M_2=M$ in the notation of § 3. We then use FFTW to compute the two-dimensional Fourier representation of $\mathcal{A}_0(x_1,x_2)$, which is used to quickly evaluate the function at any point. It would also have been easy to compute $\mathcal{A}_0(-ct,-kct)$ directly by Newton's method, but the Fourier approach is also very fast and gives more information about the function $\mathcal{A}_0(x_1,x_2)$. In particular, the modes decay to machine roundoff on the grid, corroborating the assertion of Wilkening & Zhao (Reference Wilkening and Zhao2020) that $\mathcal{A}_0$ is real analytic. We use the exact solution to compute the error in time stepping (2.27) and (2.28) from $t=0$ to $t=T$,

(4.9a,b)\begin{align} \text{err} = \sqrt{\|\tilde\eta-\tilde\eta_{{exact}}\|^2 + \|\tilde\varphi-\tilde\varphi_{{exact}}\|^2}, \quad \|\tilde\eta\|^2 = \frac1{M_1M_2}\sum_{m_1,m_2} \tilde\eta\left(\frac{2{\rm \pi} m_1}{M_1},\frac{2{\rm \pi} m_2}{M_2},T\right)^2. \end{align}

A detailed convergence study is given in Wilkening & Zhao (Reference Wilkening and Zhao2020) to compare the accuracy and efficiency of the Runge–Kutta and exponential time differencing schemes proposed in that paper using the $\gamma =1$ travelling solution above as a test case. Here we report the errors for all three waves plotted in figure 5

(4.10)\begin{equation} \begin{array}{c|c|c|c} & \gamma=5 & \gamma=1 & \gamma=0.2 \\\hline \text{err} & 1.04\times10^{{-}16} & 1.16\times10^{{-}16} & 7.38\times10^{{-}17} \end{array}, \end{equation}

using the simplest time stepping method of Wilkening & Zhao (Reference Wilkening and Zhao2020) to solve (2.27), namely a fifth-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 time stepping algorithm of Wilkening & Zhao (Reference Wilkening and Zhao2020) and the travelling wave solver of § 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$ grid points and 900 time steps), while computing the travelling waves via the Levenberg-Marquardt method took 7 seconds on a laptop and only 0.9 seconds on a server (Intel Xeon Gold 6136, 3GHz) running on 12 threads (with $M^2=3600$ grid points and $N_{tot}=1200$ unknowns).

4.3. Larger-amplitude gravity–capillary waves

In the previous sections we studied the full two-parameter family of quasi-periodic travelling waves with $k=1/\sqrt 2$, varying both $\hat \eta _{1,0}$ and $\hat \eta _{0,1}$ over the range $[-0.01,0.01]$. Here we search for larger-amplitude waves along the path $\gamma =1$, where $\hat \eta _{1,0}=\hat \eta _{0,1}=\hat \eta _{max}$ serves as an amplitude parameter. The calculations are done on an $M\times M$ grid with Fourier cutoff $N$. As the amplitude increases with $M$ and $N$ fixed, the Fourier modes outside of the cutoff region $max(|\,j_1|,|\,j_2|)\le N$ eventually grow in magnitude to exceed $\varepsilon \hat \eta _{max}$, where $\varepsilon$ is machine precision. Because we formulate the problem as an overdetermined least-squares problem, it ceases to be possible to satisfy all the equations with the limited number of Fourier degrees of freedom, and the minimum value of the objective function begins to grow rapidly with amplitude.

This is demonstrated in figure 8 using five grids ranging from $(M,N)\,{=}\,(48,11)$ to $(M,N)=(240{,}100)$ and $\hat \eta _{max}$ ranging from $0.001$ to $0.029$. Because the objective function $f$ is a squared error, if the solution has 14 digits of accuracy the objective function will be around $10^{-28}$. The coarsest grid becomes under-resolved for $\hat \eta _{max}>0.013$ while the finest grid becomes under-resolved for $\hat \eta _{max}\ge 0.0281$. One can see from the two-dimensional Fourier plots in figure 3 that $N=24$ was overkill at the amplitude $\hat \eta _{max}=0.01$ since the modes have decayed below $\varepsilon \hat \eta _{max}=1.11\times 10^{-18}$ by the time $max(|\,j_1|,|\,j_2|)\ge 11$. But we see in figure 8 that once $\hat \eta _{max}$ reaches 0.021, it becomes necessary to increase $M$ and $N$ to maintain accuracy. At this amplitude, a two-dimensional Fourier plot (not shown) contains larger-amplitude modes extending all the way to the boundary of $max(|\,j_1|,|\,j_2|)\le 24$.

Figure 8. Minimum value of the objective function $f$ for different values of $M$, $N$ and amplitude, $\hat \eta _{max}$. Each curve is labelled by two numbers, $M$ and $N$, with $N$ the smaller one. The objective function grows rapidly with $\hat \eta _{max}$ once there are not enough Fourier modes to represent the solution to machine precision.

The running time grows rapidly with grid size, with each calculation on the grids in figure 8 requiring an average of

(4.11)\begin{equation} \begin{array}{r|c|c|c|c|c} (M,N) & (48,11) & (72,24) & (96,40) & (192,75) & (240,100) \\\hline \text{running time} & 0.1\ \text{s} & 0.9\ \text{s} & 10.2\ \text{s} & 8.5\ \text{min} & 55.3\ \text{min} \end{array}, \end{equation}

on a 3 GHz server with 24 cores. The memory requirements also grow rapidly as several matrices of size $M^2\times N_{tot}$ are computed in the Levenberg–Marquardt algorithm, namely the Jacobian and its (reduced) singular value decomposition. In the $(M,N)=(240,100)$ case, each of these matrices requires 9.3 GB of storage, and we are not able to increase the problem size further due to hardware limitations. As a possible future research direction, one can try to improve the performance of the Levenberg–Marquardt method for high-dimensional problems using Krylov subspace approximations (Lin, O'Malley & Vesselinov Reference Lin, O'Malley and Vesselinov2016) without computing the entire Jacobian or its singular value decomposition (SVD).

In figure 9(a), we plot the largest-amplitude ‘fully resolved’ solution in figure 8 with $(M,N)=(240,100)$ and $(\hat {\eta }_{1, 0}, \hat {\eta }_{0, 1}) = (0.028, 0.028)$. The solid black curve is the nonlinear travelling wave, which has a maximum slope of $0.107$ over the representative interval $[0,10{\rm \pi} ]$ shown in the plot, while the dashed red curve is the linear prediction $\eta (\alpha ) = 0.056\cos (\alpha ) + 0.056\cos (\alpha /\sqrt {2})$. At this amplitude, there is a visible difference between the nonlinear and linear quasi-periodic waves, especially near the peaks and troughs of the waves. However, the difference is not large since the two base modes are still the dominant Fourier modes: the amplitudes of the other Fourier modes are less than 0.0035, which is $1/8$ of the base modes.

Figure 9. Plots of higher-amplitude quasi-periodic travelling waves. Panels (a,c) show the initial conditions $\eta$ over $[0,10{\rm \pi} ]$. Panels (b,d) show the amplitudes of Fourier modes along different directions versus the magnitude of the two-dimensional mode index $(\,j_1,j_2)$.

One of the main obstacles to computing high-amplitude solutions numerically is the slow decay of Fourier modes along certain resonant directions. To demonstrate this, we plot in figure 9(b) the amplitudes of the Fourier modes of $\tilde {\eta }$ along 7 directions: $j_1+akj_2\approx 0$ with $a\in \{1,1.1,0.9,1.5,0.5,0\}$ and $j_2 = 0$. Since $k$ is irrational, in direction $j_1 + akj_2 \approx 0$ we choose $j_1$ to be $min\{{floor}(-akj_2), N\}$ with $j_2\in \{-1,\ldots ,-N\}$. As shown in the figure, the Fourier modes decay more slowly when the ratio $-j_1/j_2$ is close to $k$. Even though $j_1 + k j_2 = 0$ is the resonant condition for linear quasi-periodic waves, the effects of this resonance persist into the nonlinear regime. Along the direction $j_1+kj_2\approx 0$, the mode with $(\,j_1,j_2)$ farthest from the origin is $\hat {\eta }_{70, -100}$ and its amplitude is $8.3\times 10^{-12}$, which is the point where floating-point error and finite $N$ truncation effects are roughly equal in this large-scale optimization problem.

In non-resonant directions, the modes decay faster, often remaining smaller than $\varepsilon \hat \eta _{max}$. For example, curves (4)–(7) in figure 9(b) drop below $10^{-20}$ for $(\,j_1^2+j_2^2)\ge 50$, whereas $\varepsilon \hat \eta _{max}=(2^{-53})(0.028)=3.1\times 10^{-18}$. This may also be observed in the two-dimensional Fourier plots of figure 3. Presumably the columns $\partial \mathcal{R}/\partial p_j$ of the Jacobian corresponding to these modes remain nearly orthogonal to the residual $\mathcal{R}$ throughout the computation, so the Levenberg–Marquardt algorithm brings them into the calculation with very small coefficients. Since increasing the amplitude beyond $\hat \eta _{max}=0.028$ leads to loss of spectral accuracy, this is the largest-amplitude wave of this type that we can compute.

Next we look for steeper waves by modifying the surface tension parameter $\tau$ and wavenumber ratio $k$. So far we have only shown calculations with $k = 1/\sqrt {2}$, which was an arbitrary choice. In ocean waves, the characteristic wavelength of gravity waves is larger than that of capillary waves by several orders of magnitude. Here we increase $k$ modestly to $\sqrt {151} \approx 12.29$, which is still much smaller than occurs in the ocean but could be relevant to a laboratory experiment. The case of pure gravity waves, which is more relevant to the ocean, will be undertaken in future work. Some comments on this were given in the introduction.

Rather than explore the two-parameter family of quasi-periodic water waves with $k=\sqrt {151}$ near the trivial solution or follow a path holding $\gamma =\hat \eta _{1,0}/\hat \eta _{0,1}$ constant, we attempt to compute steep quasi-periodic travelling waves as small quasi-periodic perturbations of large-amplitude periodic waves, which are comparatively inexpensive to compute (Dyachenko et al. Reference Dyachenko, Lushnikov and Korotkevich2016; Trichtchenko et al. Reference Trichtchenko, Deconinck and Wilkening2016). In panel (e) of figure 4 above, the contour plot of $\tau (\hat \eta _{1,0},\hat \eta _{0,1})$ represents a surface in three-dimensional $(\hat \eta _{1,0},\hat \eta _{0,1},\tau )$ space. The coordinate planes $\hat \eta _{0,1}=0$ and $\hat \eta _{1,0}=0$ in this space are two additional surfaces representing travelling waves, the first of periodic waves of wavelength $2{\rm \pi}$ and the second of periodic waves of wavelength $2{\rm \pi} /k$. The two parameters on the $\hat \eta _{0,1}=0$ surface are $\tau$ and $\hat \eta _{1,0}$. This surface intersects the $\tau (\hat \eta _{1,0},\hat \eta _{0,1})$ surface along a curve $\tau (\hat \eta _{1,0},0)$ where it is possible to bifurcate from periodic travelling waves to quasi-periodic travelling waves. As explained after (4.2) above, $\tau (\hat \eta _{1,0},0)$ is an even function of $\hat \eta _{1,0}$, so its deviation from $\tau _{lin}$ is a second-order correction.

In the present case of $k=\sqrt {151}$, we hold the surface tension fixed at $\tau =\tau _{lin}=1/\sqrt {151} \approx 0.0814$ and use the Levenberg–Marquardt method to compute the resulting one-parameter family of $2{\rm \pi}$-periodic travelling waves, denoted as $\eta _{per}$, over the range $0\le \hat \eta _{1,0} \le 0.1$. At the amplitude $\hat \eta _{1,0}=0.1$, the Fourier modes $\hat \eta _{j_1}$ of $\eta _{per}(\alpha )$ decay to machine precision around $j_1=200$. The maximum slope of this wave in physical space is 0.332, which is about 3 times steeper than the quasi-periodic wave computed above with $k=1/\sqrt 2$ and $\hat \eta _{1,0}=\hat \eta _{0,1}=0.028$. Rather than search within the family of periodic waves for the bifurcation point $\tau (0.1,0)$, we attempt to jump directly onto the family of quasi-periodic waves from the periodic wave with $\tau =\tau _{lin}$. As an initial guess for the Levenberg–Marquardt method, we set

(4.12a–c)\begin{equation} \tilde\eta^{(0)}(\alpha_1,\alpha_2) = \eta_{per}(\alpha_1) + \hat\eta_{0,1}( \textrm{e}^{\textrm{i}\alpha_2} + \textrm{e}^{-\textrm{i}\alpha_2}), \quad \tau^{(0)}=\tau_{lin}, \quad b^{(0)} = c_{per}^2, \end{equation}

in (4.1a–c), where $c_{per}$ is the wave speed of $\eta _{per}$. We succeeded in minimizing the objective function to $f=8.4\times 10^{-29}$ holding $(\hat \eta _{1,0},\hat \eta _{0,1})$ fixed at $(0.1,0.00003)$ and using $(N_1,N_2)=(216,8)$ for the Fourier cutoffs on an $M_1\times M_2=576\times 24$ grid. Using smaller values $N_2 < N_1$ and $M_2 < M_1$ is possible since the unperturbed wave is independent of $\alpha _2$, and is required to make the problem computationally tractable. We then use numerical continuation to increase $\hat \eta _{0,1}$ to $0.0003$ in increments of $0.00003$, holding $\hat \eta _{1,0}=0.1$ fixed. Polynomial interpolation of $\tau (0.1,\hat \eta _{0,1})$ from the points $\hat \eta _{0,1}\in \{\pm 0.00003m: 1\le m\le 5\}$ gives the value $\tau (0.1,0)=0.0807311$ for the surface tension of the periodic travelling wave where the bifurcation to quasi-periodicity occurs. This is only $0.8\,\%$ smaller than $\tau _{lin}$, which explains why it was possible to find nearby quasi-periodic waves to the $\tau =\tau _{lin}$ periodic waves even though this is not the precise location of the bifurcation.

The last solution on this path, with ($\hat \eta _{1,0},\hat \eta _{0,1})=(0.1,0.0003)$, is shown in panels (c) and (d) of figure 9. We had to increase the Fourier cutoffs $(N_1,N_2)$ to $(350,30)$ and the grid to $M_1\times M_2 = 720\times 64$ to achieve spectral accuracy. The objective function for this quasi-periodic solution has been minimized to $f=1.8\times 10^{-25}$ and the maximum slope over the representative interval $[0,10{\rm \pi} ]$ is $0.448$, so this wave is 35 % steeper than $\eta _{per}$ and 4.2 times steeper than the $k=1/\sqrt 2$ wave of panels (a,b) of the figure. Hardware limitations prevented increasing $\hat \eta _{0,1}$ further since there are already $M_1M_2=46\,080$ nonlinear equations in $N_{tot}=21{\kern1pt}380$ unknowns. The wave speed and surface tension of this quasi-periodic wave are $\tau =0.0809677$ and $c=1.072419$, which are close to the values $\tau (0.1,0)=0.0807311$ and $c(0.1,0)=1.071972$ of the periodic wave at the bifurcation.

Panel (c) shows the nonlinear periodic and quasi-periodic travelling waves as well as the linear quasi-periodic travelling wave $\eta = 0.2 \cos (\alpha ) + 0.0006 \cos (k\alpha )$ over the representative interval $[0,10{\rm \pi} ]$. Both nonlinear waves deviate from the linear wave by more than 50 % of the amplitude of the linear wave, which shows that these solutions are well outside of the linear regime. The difference between the periodic wave and the quasi-periodic wave is also visible, with the wave peak at $\xi =0$ perturbed upward and the others perturbed upward or downward and left or right, asymmetrically, in a non-repeating pattern. The small oscillations in the trough also change aperiodically from one trough to the next, which shows that some of the modes $\hat \eta _{j_1,j_2}$ with $j_2\ne 0$ are comparable in size to the modes of the periodic wave responsible for the capillary ripples in the troughs.

Panel (d) shows the Fourier mode amplitudes $\hat \eta _{j_1,j_2}$ along various directions in the $(\,j_1,j_2)$ lattice. Along the direction $j_1 + akj_2 \approx 0$ we choose $j_1$ to be $min\{\,{floor}(-akj_2), N_1\}$ with $j_2\in \{-1,\ldots ,-N_2\}$. One can see that the Fourier modes decay more slowly along directions $j_1 + akj_2 \approx 0$ when $a\in \{1,1.1,0.9,1.5\}$ than when $a\in \{0.5,0\}$ or when ${j_2\approx 0}$. Thus, the linear resonance condition $j_1 + k j_2 = 0$ continues to have a large effect on the Fourier modes in the nonlinear regime. As noted in § 4.1 above, we believe this is because the corresponding modes $\textrm {e}^{\textrm {i}(\,j_1+j_2k)\alpha }$ in the expansion of $\eta (\alpha )$ in (2.12a,b) have long wavelengths and are not as strongly controlled by the governing equations (2.37) as other modes, which leads to greater sensitivity to nonlinear interactions among the Fourier modes.