2. Governing equations

Free-surface irrotational water waves with a flat bottom, depth $h$ , satisfy the following equations in the absence of surface tension.

1. (i) Ideal flow:

   (2.1) \begin{equation} \Delta \phi = 0, \quad -h\lt z\lt \eta (\boldsymbol {r},t), \end{equation}

2. (ii) no flow through the bottom:

   (2.2) \begin{equation} \frac {\partial \phi }{\partial z} = 0, \quad \text{ on } \quad z=-h, \end{equation}

3. (iii) Bernoulli (pressure) equation:

   (2.3) \begin{equation} \frac {\partial \phi }{\partial t} + \frac {1}{2}|{\boldsymbol{\nabla }}\phi |^2 + \frac {1}{2}\left ( \frac {\partial \phi }{\partial z} \right )^2 + g\eta = 0 \quad \text{ on } \quad z=\eta (\boldsymbol {r},t), \end{equation}

4. (iv) kinematic boundary condition:

   (2.4) \begin{equation} \frac {\partial \phi }{\partial z} = \frac {\partial \eta }{\partial t} + \boldsymbol{\nabla }\phi \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }\eta , \quad \text{ on } \quad z=\eta (\boldsymbol {r},t), \end{equation}

where $\boldsymbol{\nabla }= \partial _x \mathrm{i} + \partial _y \mathrm{j}$ is the transverse gradient, $\varDelta = \partial _x^2 + \partial _y^2 +\partial _z^2$ is the three dimensional Laplacian and $g$ is the acceleration due to gravity. These four equations constitute the classical equations for irrotational water waves. Here, the unknowns are: $\phi ({\mathrm{\boldsymbol{r}}},z,t)$ the velocity potential and $\eta (\boldsymbol{r},t)$ the surface wave elevation; $\boldsymbol {r}=(x,y)$ is the set of horizontal coordinates, $z$ is the vertical coordinate and $t$ is time. This is a free-boundary problem for the unknowns $\phi (\boldsymbol {r},z,t)$ and $\eta (\boldsymbol{r},t)$ .

In Ablowitz, Fokas & Musslimani (Reference Ablowitz, Fokas and Musslimani2006), the water-wave problem was reformulated as a non-local differential–integral system for two surface unknowns, $\eta (\boldsymbol{r},t)$ and the horizontal velocity potential $q=q(\boldsymbol{r},t)= \phi (\boldsymbol {r},\eta (\boldsymbol {r},t),t)$ . These equations, which we will call the AFM equations, are given by

(2.5a) \begin{align} \iint _{\mathbb{R}^2} {\rm d} \boldsymbol {r} {\rm e}^{ - i \boldsymbol{k} \,\boldsymbol{\cdot }\,\boldsymbol {r} } \left ( i \eta _t \cosh [k (\eta +h)] - \frac {\boldsymbol {k} \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }q}{k} \sinh \left [ k (\eta +h)\right ]\right ) = 0 , \end{align}

(2.5b) \begin{align}\,\,\,\quad\quad\quad\quad\quad\quad\quad q_t + \frac {1}{2} | \boldsymbol{\nabla }q|^2 + g\eta - \frac {(\eta _t + \boldsymbol{\nabla }q \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }\eta )^2} {2(1+|\boldsymbol{\nabla }\eta |^2)}=0, \end{align}

where $k = | \boldsymbol {k}| = \sqrt {k_x^2 + k_y^2} \geqslant 0$ . It is assumed that $\eta , \boldsymbol{\nabla }q, q_t$ decay rapidly to zero at infinity. Equation (2.5b ) is Bernoulli’s equation on the free surface. Equations (2.5a )–(2.5b ) provide an explicit formulation of the surface variables and the non-local formulation lends itself well to asymptotic calculations, such as those considered in this paper. In the infinite depth limit, $h \to \infty$ , (2.5a ) reduces to

(2.6) \begin{equation} \iint _{\mathbb{R}^2} {\rm d} \boldsymbol {r} {\rm e}^{ - i \boldsymbol {k} \,\boldsymbol{\cdot }\,\boldsymbol {r} } {\rm e}^{k\eta } \left ( i \eta _t - \frac { \boldsymbol {k}}{k} \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }q \right ) = 0 . \end{equation}

We consider the case of weakly nonlinear waves, for which it is convenient to let $\eta \to \varepsilon \eta ,\ q \to \varepsilon q$ and assume $|\varepsilon | \ll 1$ . Doing so, and formally expanding the hyperbolic functions in (2.5a ) to order $\varepsilon$ , we find: (i) for finite depth,

(2.7) \begin{equation} \iint _{\mathbb{R}^2} {\rm d} \boldsymbol {r} {\rm e}^{-i \boldsymbol {k} \,\boldsymbol{\cdot }\,\boldsymbol {r} } \left ( i \eta _t (1 + \varepsilon k \tanh\, (k h) \eta ) - \frac { \boldsymbol {k}}{k} \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }q (\tanh\, (k h) + \varepsilon k \eta ) + \cdots \right ) =0 , \end{equation}

and (ii) for infinite depth, or $h \to \infty$ ,

(2.8) \begin{equation} \iint _{\mathbb{R}^2} {\rm d} \boldsymbol {r} {\rm e}^{-i \boldsymbol {k} \,\boldsymbol{\cdot }\,\boldsymbol {r} } \left ( \left(i \eta _t -\frac { \boldsymbol {k}}{k} \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }q\right) (1+ \varepsilon k \eta )+ \cdots \right ) =0 . \end{equation}

Note that these two equations differ only by simple factors. The free-surface Bernoulli equation (2.5b ) is unchanged regardless of finite or infinite depth; this equation to order $\varepsilon$ reads

(2.9) \begin{equation} q_t=-g\eta +\frac {\varepsilon }{2} \big(\eta _t^2-|\boldsymbol{\nabla }q|^2\big)+ \cdots . \end{equation}

The weakly nonlinear equations are considered in more detail in § 4.

We define the Fourier and inverse Fourier transforms as

(2.10a) \begin{equation} \hat {f}(\boldsymbol {k})= \frac {1}{2\pi } \iint _{\mathbb{R}^2} {\rm d} \boldsymbol {r}\kern-1pt {f}\kern-1pt(\boldsymbol {r} ){\rm e}^{-i \boldsymbol {k} \,\boldsymbol{\cdot }\,\boldsymbol {r} } , \\ \end{equation}

(2.10b) \begin{equation} f(\boldsymbol {r})= \frac {1}{2\pi } \iint _{\mathbb{R}^2} {\rm d} \boldsymbol {k}\kern-1pt\hat{f}\kern-1pt(\boldsymbol {k} ){\rm e}^{i \boldsymbol {k} \,\boldsymbol{\cdot }\,\boldsymbol {r} } . \end{equation}

Then, from (2.7) and (2.9), the linearised equations can be written as

(2.11a) \begin{equation} i\hat {\eta }_t -\frac { \boldsymbol {k}}{k} \,\boldsymbol{\cdot }\,\widehat { \boldsymbol{\nabla }q} \tanh\, (k h)=0 \end{equation}

and

(2.11b) \begin{equation} \widehat {\boldsymbol{\nabla }q}_t =-g\widehat { \boldsymbol{\nabla }\eta } , \end{equation}

where $\widehat { \boldsymbol{\nabla }\!f} = i \boldsymbol {k} \hat {f}$ . Differentiating (2.11a ) with respect to time and combining these equations leads to

(2.12) \begin{equation} \hat {\eta }_{tt}+ \omega ^2(k) \hat {\eta }=0, \end{equation}

where $\omega (k)$ is the linear dispersion relation for two-dimensional water waves in finite depth given by

(2.13) \begin{equation} \omega ^2(k)= gk\tanh\, (k h). \end{equation}

In the infinite depth limit, $h \to \infty$ , the dispersion relation is given by

(2.14) \begin{equation} \omega ^2(k)= gk. \end{equation}

The linear surface wave $\eta (x,y,t)$ can be obtained by solving (2.12) and then taking its inverse Fourier transform. When we mention the exact solution below, this is to what we are referring.

3. Linear spiral waves

In this section, linear spiral waves are observed for certain initial conditions and analysed via the method of stationary phase. Before we do this, as an illustrative example, using Fourier transforms we solve for the wave elevation $\eta (\boldsymbol {r},t)$ from (2.12) with the following initial conditions:

(3.1) \begin{equation} \eta (x,y,t=0)= x {\rm e}^{- \big( \frac{r}{L_0} \big)^2}, \quad \eta _t(x,y,t=0)= \frac {y}{T_0} {\rm e}^{- \big( \frac{r}{L_0} \big)^2}, \quad r^2=x^2+y^2, \end{equation}

where $L_0$ is a characteristic length scale and $T_0$ is a characteristic time scale. Physically, this corresponds to spatial and velocity profiles which are odd about the origin and $90^\circ$ out-of-phase with each other. This configuration provides the initial ‘twist’ sufficient to generate spiral waves. We also find that spirals evolve from a much broader set of initial conditions of the form

(3.2) \begin{equation} \eta (x,y,t=0)= P_1(x,y) f(r), \quad \eta _t(x,y,t=0)= P_2(x,y) f(r), \end{equation}

where $P_{\!j}(x,y)$ are suitable polynomials and $f(r)$ is a rapidly decaying radially symmetric function (this is discussed further in Appendix A). We highlight the initial conditions in (3.1) due to their rather simple form that generates spiral waves.

Typical evolutions of the linear water-wave equations in dimensional units appear in figure 1 for different fluid depths. In each case, we observe two spiral arms: a (positive) arm of elevation and a (negative) arm of depression. Examining the differences, the shallower depths have more identifiable spirals, but have lost more amplitude. We do not show it here, but the finite-depth evolutions for $h \gt 7.5 \,\mathrm{m}$ converge asymptotically to the infinite-depth limit dynamics (2.14). These results clearly indicate the presence of spiral waves in linear water-wave systems. We note that these are induced by the initial conditions and not by external forcing.

Figure 1.

Evolution of the linearised water-wave equation (2.12) with dispersion relation (2.13) for initial conditions (3.1). Here $t=5\, \mathrm{s}$ , $g = 9.8 \,\mathrm{m}\,\mathrm{s}^{-2}$ , $L_0 = 2.45$ m and $T_0 = 0.5$ s.

Figure 1. Long description

The heat map displays the evolution of the linearised water-wave equation with a dispersion relation for initial conditions. The map consists of three subplots, each representing different water depths: 7.5 meters, 2.5 meters, and 1.25 meters. Each subplot is a square grid with x and y axes ranging from -20 to 20 meters. The color scale on the right of each subplot ranges from -0.2 to 0.2, indicating the intensity of the wave. The colors transition from blue (negative values) through green (neutral values) to yellow (positive values). The central regions of each subplot show spiral patterns, with the intensity decreasing outward. The patterns indicate the distribution and propagation of wave frequencies over time for the given initial conditions and water depths.

Having observed these spirals numerically, we seek to describe them analytically. Here, it is convenient to transform to non-dimensional variables. This is accomplished by scaling by a characteristic length $L_0$ , $\boldsymbol{r} \rightarrow L_0\mathrm{\boldsymbol{r}}$ , $h \rightarrow L_0 h$ , chosen such that $\mathrm{\boldsymbol{r}} \sim O(1)$ , where the initial conditions are localised and characteristic time scale $T_0$ , $t\rightarrow T_0 t$ , chosen such that the non-dimensional parameter related to gravity,

(3.3) \begin{equation} \tilde {g} = \frac {T_0^2}{L_0} g, \end{equation}

is of order $1$ . In the examples below we take $\tilde {g} = 1$ and use the dimensionless initial conditions

(3.4) \begin{equation} \eta (x,y,t=0)= x {\rm e}^{-r^2}, \quad \eta _t(x,y,t=0)= y {\rm e}^{-r^2}, \quad r^2=x^2+y^2. \end{equation}

From here on, all variables presented are non-dimensional.

Motivated by previous work used to describe spiral motion in a Klein–Gordon equation (Ablowitz et al. Reference Ablowitz, Cole and Nixon2025), we implement the method of stationary phase to describe the structure of these spirals. The second-order linear water-wave equation (2.12) has the exact solution for decaying data

(3.5) \begin{equation} \eta (x,y,t) = I_+(\mathrm{\boldsymbol{r}},t) + I_-(\mathrm{\boldsymbol{r}},t), \end{equation}

where

(3.6) \begin{align} \begin{aligned} I_{\pm }(\mathrm{\boldsymbol{r}}, t) &= \frac {1}{2\pi } \iint _{\mathbb{R}^2} \widehat {A}_{\pm }(\mathrm{\boldsymbol{k}})\, \mathrm{e}^{{i}(\mathrm{\boldsymbol{k}} \,\boldsymbol{\cdot }\,\mathrm{\boldsymbol{r}} \pm \omega (k) \,t)}\mathrm{d} k_x \, \mathrm{d} k_y \\ & =\frac {1}{2\pi } \iint _{\mathbb{R}^2} \widehat {A}_{\pm }(\rho ,\varphi )\, \mathrm{e}^{{i}(\rho \cos\, \varphi \, x + \rho \sin\, \varphi \, y \pm \omega (\rho ) \, t)} \rho \, \mathrm{d} \varphi \, \mathrm{d} \rho . \end{aligned} \end{align}

In the polar version of the inverse Fourier transform, we take $k_x = \rho \cos\, \varphi$ and $k_y = \rho \sin\, \varphi$ . We intentionally take $(\rho , \varphi )$ to denote the spectral polar coordinates to distinguish them from the standard labels $(r, \theta ), \text{ where } x=r \cos\, \theta ,\ y=r \sin\, \theta$ . Notice that the dispersion relation only depends on the modulus of the Fourier wavenumbers. The coefficients $\widehat {A}_{\pm }$ are related to the initial conditions by

(3.7) \begin{equation} \widehat {A}_{\pm }(\rho , \varphi ) = \frac {1}{2}\left ( \widehat {\eta }(\rho , \varphi ,0) \pm \frac {1}{\mathrm{i} \omega (\rho )}\,\frac {\partial \widehat {\eta }}{\partial t}(\rho ,\varphi ,0) \right )\!, \end{equation}

where $\widehat {A}_+=\widehat {A}^*_-$ are complex conjugates for real initial conditions. As indicated above, the long-time, $t\gg 1$ , asymptotic approximation of the exact solutions is obtained through a stationary-phase analysis. This is done by first approximating in $\varphi$ , then $\rho$ . Since the dispersion relation only depends on $\rho$ , this approach simplifies the calculation. Redefining the spatial variables, $x$ , $y$ , $r$ , we create $O(1)$ variables

(3.8) \begin{equation} \overline {x} = \frac {x}{t}, \quad \overline {y} = \frac {y}{t}, \quad \overline {r} = \frac {r}{t}. \end{equation}

First, isolate the integral in $\varphi$ :

(3.9) \begin{equation} J_{\pm }(\rho , \overline {r},\theta , t) = \int _{-\pi }^{\pi } \widehat {A}_{\pm }(\rho ,\varphi )\, \mathrm{e}^{{i} (\rho \cos\, \varphi \, \overline {x} + \rho \sin\, \varphi \, \overline {y})\, t} \mathrm{d} \varphi , \end{equation}

which has two stationary points of the Fourier phase, $\varphi _{\!j}$ $j = 1,2$ , in the interval $[-\pi ,\pi ]$ that satisfy the equation

(3.10) \begin{equation} \varphi _{1} = \theta ,\quad \varphi _{2} = \theta + \pi . \end{equation}

Here it appears that there are a total of four stationary points: two for the positive case and two for the negative. However, as we show below, only two of these values correspond to a stationary point of the solution (3.6). For the other two, the monotonicity of $\omega (\rho )$ means the phase is stationary only with respect to the angular variable, $\varphi$ , but not with respect to the radial variable, $\rho$ .

The second derivative of the phase is given by

(3.11a) \begin{align} \frac {{\rm d} ^2}{\mathrm{d} \varphi ^2}\big (\rho \cos\, \varphi \, \overline {x} + \rho \sin\, \varphi \, \overline {y}\, \big )\big |_{\varphi _{\!j}} & = \rho \big (\!- \cos\, \varphi _{\!j}\, \overline {x} -\sin\, \varphi _{\!j} \, \overline {y}\, \big ) \\[-12pt]\nonumber \end{align}

(3.11b) \begin{align} & =(-1 )^{j} \rho \, \overline {r}. \end{align}

Using this, we find the stationary-phase approximation:

(3.12) \begin{equation} J_{\pm }(\rho , \overline {r}, \theta , t) \approx \sqrt {\frac {2 \pi }{\rho \overline {r} \, t}} \widehat {A}_{\pm }(\rho , \varphi _1) \mathrm{e}^{{i}\, \rho \overline {r} \, t - {i} {\pi }/{4}} + \sqrt {\frac {2 \pi }{\rho \overline {r} \, t}} \widehat {A}_{\pm }(\rho , \varphi _2) \mathrm{e}^{-\, ({i}\,\rho \overline {r} \, t - {i} {\pi }/{4} )} . \end{equation}

The integral in $\rho$ can now be expressed as

(3.13) \begin{align} I_{\pm }(\overline {r},\theta , t) & \approx \frac {1}{2\pi } \int _0^{\infty } J_{\pm }(\rho , \overline {r}, t) \mathrm{e}^{\pm {i} \omega (\rho ) t} \rho \, \mathrm{d} \rho \nonumber\\ & = \mathrm{e}^{-{i} \frac{\pi}{4}}\frac {1}{\sqrt {2 \pi \overline {r} \, t}} \int _0^{\infty } \, \widehat {A}_{\pm }(\rho , \varphi _1) \mathrm{e}^{{i}\, ( \rho \,\overline {r} \pm \omega (\rho ) ) \, t} \sqrt {\rho }\,\mathrm{d} \rho \nonumber\\ & \quad + \mathrm{e}^{{i} \frac{\pi}{4}}\frac {1}{\sqrt {2 \pi \overline {r} \, t}} \int _0^{\infty } \, \widehat {A}_{\pm }(\rho , \varphi _2) \mathrm{e}^{{i}\, (- \, \rho \,\overline {r} \pm \omega (\rho ) ) \, t} \sqrt {\rho }\, \mathrm{d} \rho . \notag \end{align}

Defining the phase term

(3.14) \begin{equation} Q_{\!j}^{\pm }(\rho , \overline {r}) = (-1)^{j+1} \rho \,\overline {r} \pm \omega \big (\rho \big ), \quad j = 1,2, \end{equation}

the stationary points occur for values of $\overline {r}$ that satisfy

(3.15) \begin{equation} \frac {\mathrm{d} Q_{\!j}^{\pm }}{\mathrm{d} \rho } = (-1)^{j+1} \overline {r} \pm \frac {\partial \omega }{\partial \rho } = 0 . \end{equation}

For water waves, $\omega '(\rho )\gt 0$ for $\rho \gt 0$ ; thus solutions to (3.15) only occur when $j=1$ for the $-$ case or $j=2$ for the $+$ case. In either case, the radial component of the stationary point is given by

(3.16) \begin{equation} \overline {r} = \frac {\partial \omega }{\partial \rho }. \end{equation}

We define the solution to (3.16) as $\rho _0 (\overline {r} )$ , which represents the only stationary point for both integrals in (3.13). Thus, stationary points of the phase in solution (3.6) occur at $ (\rho _0 (\overline {r} ),\theta )$ and $ (\rho _0 (\overline {r} ),\theta +\pi )$ . Checking the second derivative of the phase term, we have

(3.17) \begin{equation} \mathrm{sgn}\left [\frac {\partial ^2 Q_{\!j}^{\pm }}{\partial \rho ^2}\left (\rho _0 \right )\right ] = \mathrm{sgn}\left [\pm \frac {\partial ^2\omega }{\partial \rho ^2}\left (\rho _0 \right )\right ] = \mp \,1, \end{equation}

since $\omega ''(\rho )\lt 0$ . Thus, for $t \gg 1$ , we have the final asymptotic approximation

(3.18) \begin{align} \eta (\overline {r},\theta , t) & = I_-(\overline {r},\theta ,t)+I_+(\overline {r},\theta ,t) \nonumber \\[3pt] & \approx \frac {\sqrt {\rho _0}}{ \sqrt {\overline {r} \,|\omega ''(\rho _0)|} \, t} \big( \, \widehat {A}_{-}(\rho _0, \varphi _1) \mathrm{e}^{{i}\, ( \rho _0 \,\overline {r} - \omega (\rho _0) ) \, t} + \text{c.c.}\big) , \notag \\[5pt] & \approx \frac {\rho _0^{3/2}\mathrm{e}^{- \frac{\rho _0^2}{4}}}{4 \sqrt {\overline {r} \,|\omega ''(\rho _0)|} \, t} \Bigg [ \frac { \sin\, \theta }{\omega (\rho _0)} \cos\, \big ( \rho _0 \,\overline {r}t - \omega (\rho _0)t\big ) - \cos\, \theta \sin\, \big ( \rho _0 \,\overline {r}t - \omega (\rho _0)t\big ) \Bigg ], \end{align}

where $\rho _0$ is the function of $\overline {r}$ that solves (3.16) and c.c. is the complex conjugate. We note that this unsteady solution decays like $t^{-1}$ for $\overline {r} = {O}(1)$ . This leading-order unsteady solution has three distinctive parts: (i) linear decay in time; (ii) envelope; and (iii) spiral structure inside the brackets. In general, we cannot write down an explicit formula for $\rho _0$ that satisfies (3.16) and instead use a numerical approximation.

Figure 2.

(a) Approximation given by the stationary-phase formulae (3.18) of a finite-depth spiral at $t=20$ . (b) Relative norm squared error (3.19) between the exact solution and the approximation is observed to converge like $O(t^{-2})$ as $t \rightarrow \infty$ . For both figures $h =1$ and $\tilde {g}=1$ .

Figure 2. Long description

The image contains two graphs. The first graph on the left shows an approximation of a finite-depth spiral wave using stationary-phase formulae. The x-axis ranges from negative 20 to 20, and the y-axis ranges from negative 20 to 20. The spiral pattern is color-coded, with a color bar indicating values from negative 0.03 to 0.03. The second graph on the right depicts the relative norm squared error between the exact solution and the approximation. The x-axis represents time (t) ranging from 0 to 40, and the y-axis represents relative error ranging from 0 to 0.30. The error decreases as time increases, showing convergence. All values are approximated.

Figure 3.

Approximation given by the stationary-phase formulae (3.18) of a finite-depth spiral; here $t=20$ , $h =0.5$ and $\tilde {g}=1$ . (a) Contour plot. (b) Comparison of the solution obtained via Fourier transform, (2.12), and the approximation at cross-section $y=0$ .

Figure 3. Long description

The image contains two graphs. The left graph, labeled (a), is a contour plot showing the approximation of a finite-depth spiral wave with a step size of 0.5. The x and y axes range from -15 to 15, and the color scale indicates values from -0.04 to 0.04. The right graph, labeled (b), compares the exact solution with the approximate solution at a cross-section where y equals 0 and t equals 20. The x-axis ranges from -15 to 15, and the y-axis represents the variable eta, ranging from -0.08 to 0.06. The exact solution is represented by a solid black line, while the approximate solution is shown by a dashed cyan line. The graph illustrates how closely the approximation matches the exact solution across the x-axis.

A comparison between the exact solution obtained via Fourier transform from (2.12) and the stationary-phase approximation in (3.18) is shown in figures 2 and 3 (here and below in dimensionless units) for water depths of $h=1$ and $h=0.5$ , respectively. To quantify how accurate the stationary-phase approximation is, we compute the relative norm squared error:

(3.19) \begin{equation} \mathrm{Relative \ error} = \frac {\iint \left |\eta _{\textit{Exact}} - \eta _{\textit{Approx}} \right |^2 \mathrm{d} A} {\iint \left | \eta _{\textit{Exact}} \right |^2 \mathrm{d} A}, \end{equation}

and observe that this converges like $O(t^{-2})$ as $t$ increases. Note that the modulus of the stationary point, $\rho _0$ , is found by numerically solving (3.16). The circumferential part, for a given $\overline {x},\overline {y}$ , is obtained from (3.10). We note that the approximation is valid in the domain $x^2 + y^2 \lt$ $\tilde {g}$ $h t^2$ , which can be seen by noting that

(3.20) \begin{equation} \lim _{\rho \rightarrow 0} \frac {\mathrm{d} \omega }{\mathrm{d} \rho } = \sqrt {\tilde {g} h} \end{equation}

is the supremum of $\omega '$ for finite depth. Thus, no solution to (3.15) exists for $\overline {r}\geqslant \mathrm{sup}\ \omega '(\rho )$ . Within this region, there is good agreement with the exact solution. Other integral approximations, like steepest descent, are required outside of this region. Finally, we remark that these spirals are quite robust to perturbations and persist in the nonlinear regime (as seen in § 4).

3.1. Deep water

Thus far, we have found an asymptotic approximation for the spiral structures, which shows that the phenomena can be expressed as a slowly varying, radially symmetric, envelope over an expression for the underlying spiral structure. In the deep-water limit, we further investigate this spiral structure to obtain a more explicit formula for the shape of the spiral expressed as $r$ as a function of $\theta$ , which traces the path the spiral takes as it expands out from the origin. Here, as $h\rightarrow \infty$ , $\omega$ is given by (2.14), or in polar form

(3.21) \begin{equation} \omega (\rho ) = \tilde {g}^{1/2} \rho ^{1/2}. \end{equation}

The equation for the stationary points, (3.16), now becomes

(3.22) \begin{equation} \overline {r} = \frac {1}{2} \tilde {g}^{1/2} \rho _0^{-1/2} \end{equation}

with solutions

(3.23) \begin{equation} \rho _{0} = \frac {\tilde {g}}{4 \overline {r}^2}. \end{equation}

Substituting this into (3.18), we find

(3.24) \begin{equation} \eta (\overline {r}, \theta , t) \approx - \frac {\tilde {g}}{t} \, \frac {1}{ 2 \sqrt {2}\,\overline {r}^{3}} \big ( \widehat {A}_{-}(\rho _0, \varphi _1) \mathrm{e}^{-{i}( {\tilde {g}}/{4 \overline {r}) } \, t} + \text{c.c.} \big ) \end{equation}

with initial conditions (3.4) from § 2. After the Fourier coefficients in (3.7) are evaluated, we have

(3.25) \begin{equation} \eta (\overline {r}, \theta ,t) \approx \frac {\tilde {g}^2 \mathrm{e}^{- \frac{\tilde {g}^2}{64 \overline {r}^4}}}{32 \sqrt {2}\, \overline {r}^5 \, t} \left ( \frac {2 \overline {r}}{\tilde {g}} \sin\, (\theta ) \cos\, \left (\frac {\tilde {g}t}{4 \overline {r}} \right ) - \cos\, (\theta ) \sin\, \left (\frac {\tilde {g} t}{4 \overline {r}} \right )\right ) \!. \end{equation}

With regard to the envelope structure, we observe a dominant exponential decay to zero when $\overline {r} \rightarrow 0$ , but an algebraic (quintic) decay as $\overline {r} \rightarrow \infty$ . The spiral structure is further analysed in the next section.

Figure 4.

Comparison of the deep-water (a) exact solution (2.12) computed via Fourier transforms and the (b) approximation given by the stationary-phase formulae (3.25). (c) Cross-section at $y=0$ comparing the exact solution and approximation in (a,b). (d) Relative norm squared error (3.19) between the exact solution and the approximation; it converges like $O(t^{-2})$ as $t \rightarrow \infty$ . Here $\tilde {g}=1$ .

Figure 4. Long description

The image contains four graphs comparing the exact solution and approximation of deep-water waves. The first graph (a) shows the exact solution computed via Fourier transforms, displaying a spiral pattern. The second graph (b) shows the approximation given by the stationary-phase formulae, with a similar spiral pattern but slight differences. The third graph (c) is a cross-section at t equals 15, y equals 0, comparing the exact solution and approximation, showing two lines that closely follow each other with minor deviations. The fourth graph (d) depicts the relative norm squared error between the exact solution and the approximation, showing a decreasing trend over time. The error converges as t increases.

Figure 5.

Comparison of the equation for the spiral shape (3.29) found in the deep-water limit (dashed line) and the exact solution obtained via Fourier transform. The black dashed curve tracks the zeros of the spiral.

Figure 5. Long description

The heat map displays three spiral patterns at different time points, specifically at t equals 10, 15, and 20. Each subplot shows a square grid with x and y axes ranging from negative 10 to 10. The color scale on the right of each subplot indicates the value intensity, with colors ranging from blue (negative values) to yellow (positive values). The dashed lines in each subplot track the zeros of the spiral. The central region of each spiral shows higher intensity values, gradually decreasing outward. The patterns exhibit a concentric structure, with the intensity values becoming less pronounced as they move away from the center.

A comparison between the exact solution obtained from (2.12) with the deep-water dispersion and the stationary-phase approximation in (3.25) is shown in figure 4. We observe that both qualitatively and quantitatively, the approximation describes the spiral wave phenomena.

3.1.1. Deep-water spiral structure

We now isolate the underlying spiral structure from the slowly varying envelope in order to describe the fundamental shape of the deep-water spirals. This spiral shape can be derived from (3.25). Since the surface spiral consists of regions of elevation (positive) and depression (negative), we can trace the shape of the spiral along the zeros of the solution, which occur at

(3.26) \begin{equation} \frac {2 \overline {r}}{\tilde {g}} \sin\, (\theta ) \cos\, \left (\frac {\tilde {g} t}{4 \overline {r}} \right ) - \cos\, (\theta ) \sin\, \left (\frac {\tilde {g} t}{4 \overline {r}} \right ) = 0 \end{equation}

or, in terms of $r = t \overline {r}$ (see (3.8)),

(3.27) \begin{equation} r \tan (\theta ) = \frac {\tilde {g} t}{2} \tan \left (\frac {\tilde {g} t^2}{4 r} \right )\!. \end{equation}

To better visualise the shape of the spiral, we seek an approximate solution for $r$ in terms of $\theta$ . If we look for the zeros ( $\theta = 0$ ) that lie along the positive $x$ axis, we find that $n{\mathrm{th}}$ zero has a radius of

(3.28) \begin{equation} \frac {\tilde {g} t^2}{4 r_n} = n \pi , \quad n \in \mathbb{Z} . \end{equation}

Interpolating these zeros on the positive $x$ axis for general $\theta$ , it suggests that the shape of the spiral is given by

(3.29) \begin{equation} r = \frac {\tilde {g} t^2}{4 \theta }. \end{equation}

For a fixed time $t$ , this is a hyperbolic spiral. A comparison of the curve given in (3.29) superimposed on top of the exact solution obtained from (2.12) in the deep-water limit is shown in figure 5. The curve fits the zeros of the spiral with good agreement between the two. While it is not shown here, these curves also correspond to the extrema of the surface velocity profile, $\eta _t(x,y,t)$ . It is interesting to note that hyperbolic spirals have been used to describe certain galaxies where the radial spiral arm grows as it moves away from the centre (Kennicutt Reference Kennicutt1981).

4. Weakly nonlinear spiral waves

While the spiral waves discussed above are a linear phenomenon, these waves are found to be robust structures that remain intact in the presence of weak nonlinearity. Next we derive a numerical scheme for evolving spiral waves in the weakly nonlinear limit by taking the Taylor expansion of (2.5) to find suitable approximations that can be solved numerically using spectral methods to approximate the spatial derivatives and standard fourth-order Runge–Kutta methods for the time evolution. Transforming the non-local system (2.5) to non-dimensional variables and considering the amplitude limit, $\eta \rightarrow \varepsilon \eta$ and $q \rightarrow \varepsilon q,\ |\varepsilon | \ll 1$ , yields

(4.1a) \begin{equation} \iint _{\mathbb{R}^2} {\rm d} \mathrm{\boldsymbol{r}} {\rm e}^{ - i \boldsymbol {k} \,\boldsymbol{\cdot }\,\boldsymbol {r} } \left ( i \eta _t \cosh \big [k h +\varepsilon k \eta \big ] - \frac {\boldsymbol {k} \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }q}{k} \sinh \big [k h +\varepsilon k \eta \big ]\right ) = 0 , \end{equation}

(4.1b) \begin{equation} q_t + \tilde {g}\eta + \frac {\varepsilon }{2} | \boldsymbol{\nabla }q|^2 - \varepsilon \frac {(\eta _t + \varepsilon \boldsymbol{\nabla }q \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }\eta )^2} {2(1+ \varepsilon ^2 |\boldsymbol{\nabla }\eta |^2)}=0. \end{equation}

Then, expanding (4.1a ) in powers of $\varepsilon$ , we find

(4.2) \begin{align} 0 = &\,\widehat {\eta }_t + \tanh\, (k h ) \frac {\mathrm{i\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\boldsymbol{\nabla }q}}{k} + \varepsilon \big (k \tanh\, (k h ) \widehat {\eta \eta _t} + \mathrm{i\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\eta \boldsymbol{\nabla }q} \big ) \notag \\ &\quad \quad + \frac {\varepsilon ^2}{2} \big (k^2 \widehat {\eta ^2 \eta _t} + k \big ( \mathrm{i\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\eta ^2 \boldsymbol{\nabla }q }\big ) \tanh\, (k h ) \big ) + \cdots . \end{align}

In order to evolve the system, we need to solve for $\eta _t$ . However, in this form, $\eta _t$ is nested within the nonlinear terms. To derive a tractable weakly nonlinear system, we expand the governing equations in powers of $\varepsilon$ and solve perturbatively for $\eta _t$ :

(4.3) \begin{equation} \eta _t = H_0 + \varepsilon H_1 + \cdots . \end{equation}

Substituting this into (4.2) and collecting powers of $\varepsilon$ gives

(4.4a) \begin{align} H_0 & = - \frac {1}{2\pi } \iint _{\mathbb{R}^2} {\rm d} \mathrm{\boldsymbol{k}} {\rm e}^{ i \boldsymbol {k} \,\boldsymbol{\cdot }\,\boldsymbol {r} } \,\tanh\, (k h ) \frac {\mathrm{i\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\boldsymbol{\nabla }q}}{k} , \\[-12pt]\nonumber \end{align}

(4.4b) \begin{align} H_1 & = - \frac {1}{2\pi } \iint _{\mathbb{R}^2} {\rm d} \mathrm{\boldsymbol{k}} {\rm e}^{ i \boldsymbol {k} \,\boldsymbol{\cdot }\,\boldsymbol {r} } \, \left [ k \tanh\, (k h ) \widehat {\eta H_0} + \mathrm{i\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\eta \boldsymbol{\nabla }q } \right ]\! . \end{align}

Equation (4.2) may now be written in explicit terms as

(4.5) \begin{align} \widehat {\eta }_t \approx &\, - \tanh\, (k h ) \frac {\mathrm{i\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\boldsymbol{\nabla }q}}{k} - \varepsilon \big ( k \tanh\, (k h ) \widehat {\eta H_0} + \mathrm{i\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\eta \boldsymbol{\nabla }q }\big ) \notag\\ &\quad - \varepsilon ^2 \Bigg (k \tanh\, (k h ) \widehat {\eta H_1} + \frac {k^2}{2} \widehat {\eta ^2 H_0} + \frac {k}{2} \big ( \mathrm{i\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\eta ^2 \boldsymbol{\nabla }q }\big ) \tanh\, (k h ) \Bigg ), \end{align}

where we truncate the expansion at $O(\varepsilon ^2)$ .

Applying expansion (4.3) in (4.1b ) now gives

(4.6a) \begin{align} 0 & = q_t + \tilde {g}\,\eta + \frac {\varepsilon }{2} \big |\boldsymbol{\nabla }q\big |^2 - \frac {\varepsilon }{2} \frac {\big (\eta _t + \varepsilon \boldsymbol{\nabla }\eta \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }q\big )^2}{1+\varepsilon ^2\big | \boldsymbol{\nabla }\eta \big |^2}, \\[-12pt]\nonumber \end{align}

(4.6b) \begin{align} 0 & = q_t + \tilde {g}\,\eta + \frac {\varepsilon }{2} \big |\boldsymbol{\nabla }q\big |^2 \notag \\ &\quad - \frac {\varepsilon }{2} \big ( \eta _t^2 + 2\varepsilon \eta _t \boldsymbol{\nabla }\eta \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }q + \cdots \big )\big (1 - \varepsilon ^2 \big | \boldsymbol{\nabla }\eta \big |^2 + \cdots \big ), \\[-12pt]\nonumber \end{align}

(4.6c) \begin{align} 0 & = q_t + \tilde {g}\,\eta + \frac {\varepsilon }{2} \big |\boldsymbol{\nabla }q\big |^2 - \frac {\varepsilon }{2} \eta _t^2 - \varepsilon ^2 \eta _t \boldsymbol{\nabla }\eta \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }q + \cdots , \\[-12pt]\nonumber \end{align}

(4.6d) \begin{align} 0 & = q_t + \tilde {g}\,\eta + \frac {\varepsilon }{2} \big |\boldsymbol{\nabla }q\big |^2 - \frac {\varepsilon }{2} (H_0^2 + 2 \varepsilon H_0 H_1) - \varepsilon ^2 H_0\boldsymbol{\nabla }\eta \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }q - \cdots ,\end{align}

where $H_0$ and $H_1$ are taken from the $\eta _t$ expansion (4.4). Hence,

(4.7) \begin{equation} q_t \approx - \tilde {g}\,\eta + \frac {\varepsilon }{2} \big (H_0^2 - \big |\boldsymbol{\nabla }q\big |^2 \big ) + \varepsilon ^2 \big ( H_0 H_1 + H_0 \boldsymbol{\nabla }\eta \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }q \big ) . \end{equation}

The final system describing the evolution of weakly nonlinear waves is given by

(4.8a) \begin{align} \widehat {\eta }_t &\approx - \tanh\, (k h ) \frac {\mathrm{i\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\boldsymbol{\nabla }q}}{k} - \varepsilon \big ( k \tanh\, (k h ) \widehat {\eta H_0} + \mathrm{i\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\eta \boldsymbol{\nabla }q }\big ) \notag \\ &\quad - \varepsilon ^2 \Bigg (k \tanh\, (k h ) \widehat {\eta H_1} + \frac {k^2}{2} \widehat {\eta ^2 H_0} + \frac {k}{2} \big ( \mathrm{i\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\eta ^2 \boldsymbol{\nabla }q }\big ) \tanh\, (k h ) \Bigg ), \\[-12pt]\nonumber \end{align}

(4.8b) \begin{align} \big(\widehat {\boldsymbol{\nabla }q}\big)_t &\approx \mathrm{i\boldsymbol{k}} \Bigg (\!- \tilde {g}\,\widehat {\eta } + \frac {\varepsilon }{2} \big (\widehat {H_0^2} -\widehat {\big |\boldsymbol{\nabla }q\big |^2} \big ) + \varepsilon ^2 \big ( \widehat {H_0 H_1} + \widehat {H_0\boldsymbol{\nabla }\eta \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }q} \big ) \Bigg ), \end{align}

where $H_0,H_1$ are given by (4.4). Or, in terms of the dispersion relation, $\omega (k)$ , we can write these equations as

(4.9a) \begin{align} H_0 & = - \frac {1}{2\pi \tilde {g}} \iint _{\mathbb{R}^2} {\rm d} \mathrm{\boldsymbol{k}} {\rm e}^{ i \boldsymbol {k} \,\boldsymbol{\cdot }\,\boldsymbol {r} } \,\omega ^2(k) \frac {\mathrm{i\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\boldsymbol{\nabla }q}}{k^2} , \\[-12pt]\nonumber \end{align}

(4.9b) \begin{align} H_1 & = - \frac {1}{2\pi \tilde {g}} \iint _{\mathbb{R}^2} {\rm d} \mathrm{\boldsymbol{k}} {\rm e}^{ i \boldsymbol {k} \,\boldsymbol{\cdot }\,\boldsymbol {r} } \, \big[\omega ^2(k) \widehat {\eta H_0} + \mathrm{i} \tilde {g} \mathrm{\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\eta \boldsymbol{\nabla }q } \big], \\[-12pt]\nonumber \end{align}

(4.9c) \begin{align} \widehat {\eta }_t &\approx - \omega ^2(k) \frac {\mathrm{i\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\boldsymbol{\nabla }q}}{\tilde {g} k^2} - \varepsilon \Bigg ( \frac {1}{ \tilde {g}}\omega ^2(k) \widehat {\eta H_0} + \mathrm{i\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\eta \boldsymbol{\nabla }q }\Bigg ) \notag \\ &\quad - \varepsilon ^2 \Bigg (\frac {1}{\tilde {g}}\omega ^2(k) \widehat {\eta H_1} + \frac {k^2}{2} \widehat {\eta ^2 H_0} + \frac {1}{2 \tilde {g}}\omega ^2(k) \big ( \mathrm{i\boldsymbol{k}}\,\boldsymbol{\cdot }\,\widehat {\eta ^2 \boldsymbol{\nabla }q }\big ) \Bigg ), \\[-12pt]\nonumber \end{align}

(4.9d) \begin{align} \big(\widehat {\boldsymbol{\nabla }q}\big)_t &\approx \mathrm{i\boldsymbol{k}} \Bigg (\!- \tilde {g}\,\widehat {\eta } + \frac {\varepsilon }{2} \big (\widehat {H_0^2} -\widehat {\big |\boldsymbol{\nabla }q\big |^2} \big ) + \varepsilon ^2 \big ( \widehat {H_0 H_1} + \widehat {H_0\boldsymbol{\nabla }\eta \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }q} \big ) \Bigg ). \end{align}

We note that, to $O(\varepsilon ^2)$ , these equations are equivalent to those derived in Craig & Sulem (Reference Craig and Sulem1993) where the expansion is found using the Dirichlet to Neumann operator (see e.g. Wang & Milewski Reference Wang and Milewski2012).

In these equations we evolve $\widehat {\eta }$ , $\widehat {q_x}$ and $\widehat {q_y}$ ; note that the last line has $\boldsymbol{\nabla }q$ , not $q$ . This has two benefits, one physical and one numerical. First, $\boldsymbol{\nabla }q$ represents the surface fluid velocity, which is an important physical quantity. Second, an important source of numerical instability arises from the high-frequency ( $k\gg 1$ ) modes, and in this form the scaling of the nonlinear terms with respect to $k$ is reduced in order, i.e. terms of order $\varepsilon$ scale like $O(k)$ instead of $O(k^2)$ and terms of order $\varepsilon ^2$ scale like $O(k^2)$ instead of $O(k^3)$ . The numerical scheme used to solve (4.8) is a fourth-order Runge–Kutta method with integrating factor. Details are given in Appendix B.

In figure 6, we plot $\eta (\mathrm{\boldsymbol{r}},t)$ overlaid with the vector field defined by $\boldsymbol{\nabla }q(\mathrm{\boldsymbol{r}},t)$ ; this illustrates how the initial spiral structure develops from the initial conditions (3.4). In the initial conditions, the velocity field is characterised by a sink and source that are located along the axis separating the initial regions of depression and elevation (in this example, the y axis). As the system is evolved, this sink and source begin rotating around each other. For long-time evolutions, this develops into a radially oriented velocity field that is inwardly directed for regions of depression and outwardly oriented for regions of elevation.

Figure 6.

Evolution of the finite-depth nonlinear water-wave equations (4.8) with initial conditions (3.4). Surface velocity $\boldsymbol{\nabla }q(\mathrm{\boldsymbol{r}},t)$ is displayed as a vector field overlaying the plot of $\eta (\mathrm{\boldsymbol{r}},t)$ . Here $h=1,\ \varepsilon = 0.25$ and $\tilde {g}=1$ .

Figure 6. Long description

A heat map displays the evolution of finite-depth nonlinear water-wave equations with initial conditions. The surface velocity is represented as a vector field overlaying the plot of eta. The heat map is divided into three panels, each representing different time points: t equals 0, t equals 3, and t equals 6. Each panel shows a grid with x and y axes ranging from negative 5 to 5. The color scale on the right of each panel indicates the magnitude of eta, with colors ranging from blue (negative values) to yellow (positive values). The maximum absolute value of eta decreases over time, from 0.429 at t equals 0 to 0.232 at t equals 3, and further to 0.149 at t equals 6. The vector field shows the direction and magnitude of surface velocity, with vectors pointing outward in a spiral pattern from the center. The overall trend indicates a spreading and decreasing intensity of the water wave over time.

Figure 7.

Evolution of the finite-depth nonlinear water-wave equations (4.8) with initial conditions (3.4). Here $\varepsilon = 0.25$ , $h=0.7$ and $\tilde {g}=1$ . This can be compared with the linear results shown in figure 3.

Figure 7. Long description

The heat map displays the evolution of the finite-depth nonlinear water-wave equations over time. The map is divided into three panels, each representing a different time point: t equals 0, t equals 10, and t equals 20. Each panel shows a grid with x and y axes ranging from negative 20 to positive 20. The color scale on the right of each panel indicates the magnitude of the variable eta, with colors ranging from blue (negative values) to yellow (positive values). At t equals 0, the maximum absolute value of eta is 0.429, with a concentrated spot in the center. At t equals 10, the maximum absolute value of eta decreases to 0.11, showing a circular pattern. At t equals 20, the maximum absolute value of eta further decreases to 0.06, with a more pronounced spiral pattern. The overall trend shows a decrease in the maximum absolute value of eta over time, along with the development of more complex patterns.

Using the same initial conditions (3.4) prescribed in the linear system described § 3, we find that linear spiral waves constitute a robust family of solutions and adding weak nonlinearity does not have significant impact on the formation and persistence of the spiral waves observed in the linear problem. Snapshots of the nonlinear evolution for finite-depth waves are shown in figure 7. These results can be compared against the linear evolution shown in figure 3. The difference between the two is small (see also figure 12) with the nonlinear version showing a few extra spirals near the centre and a slightly larger peak amplitude. This makes intuitive sense, since the spirals are composed of non-stationary waves and the weak nonlinearity does not grow in strength due to nonlinear feedback.

Deep-water spirals in the presence of weak nonlinearity are highlighted in figure 8. These results can be compared with the linearised version in figure 4. Comparing the two at $t = 10$ , we see that the spiral phenomenon persists and the hyperbolic spiral shape, as given in § 3.1.1, remains unchanged. Similar to the finite-depth case, the nonlinear evolution has a slightly larger amplitude than the linear version. For the problems considered here, the weak nonlinearity has modest impact on the evolution of the linear spiral waves. We also find numerically that adding a small amount of noise, of the order of 1 %–5 %, to the initial conditions has little effect; i.e. the spiral waves appear to be robust.

Figure 8.

Evolution of the nonlinear water-wave equations (4.8) in infinite depth with initial conditions (3.4). Here $\varepsilon = 0.25$ and $\tilde {g}=1$ . This can be compared with the linear results shown in figure 4.

Figure 8. Long description

The heat map displays the evolution of nonlinear water-wave equations in infinite depth with specific initial conditions. The map is divided into three panels representing different time points: t equals 0, t equals 10, and t equals 20. Each panel shows a grid layout with x and y axes ranging from negative 10 to positive 10. The color scale on the right of each panel indicates the magnitude of the variable eta, with colors ranging from blue (negative values) to yellow (positive values). At t equals 0, the maximum absolute value of eta is 0.428, with a concentrated spot in the center. At t equals 10, the maximum absolute value of eta decreases to 0.137, forming a spiral pattern. At t equals 20, the maximum absolute value of eta further decreases to 0.072, showing a more diffused spiral pattern. The overall trend indicates a decrease in the maximum absolute value of eta over time, with the initial concentrated spot evolving into a spiral pattern that becomes more diffused.

We also note the similarity between the spiral problem analysed here and the well-known problem of a small object, e.g. a pebble, thrown into still water. We refer to this as the ‘pebble problem’. The only difference is that the pebble problem arises from radially symmetric initial conditions whereas spiral waves evolve from initial conditions of the form given in (3.2), e.g. those given by (3.4). This is also discussed further in Appendix A. Based on this analogy we expect the spiral waves discussed here to be observable in various settings.