2 Surface waves on a flow

We consider surface waves propagating on a flow of constant depth $h_{0}$ . The flow is assumed incompressible ( $\unicode[STIX]{x1D70C}=\text{const}.$ ) and irrotational, so that the velocity flow is derived from a potential $\boldsymbol{v}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}$ . On top of a background (mean) flow $\boldsymbol{v}_{0}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}_{0}$ , small perturbations $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D6F7}-\unicode[STIX]{x1D6F7}_{0}$ propagate according to the wave equation (Milewski & Keller Reference Milewski and Keller1996)

(2.1) $$\begin{eqnarray}\displaystyle {\mathcal{D}}_{t}^{2}\unicode[STIX]{x1D719}-\text{i}(g\unicode[STIX]{x1D735}-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D735}^{3})\tanh (-\text{i}h_{0}\unicode[STIX]{x1D735})\unicode[STIX]{x1D719}+2\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}{\mathcal{D}}_{t}\unicode[STIX]{x1D719}=0, & & \displaystyle\end{eqnarray}$$

where ${\mathcal{D}}_{t}=\unicode[STIX]{x2202}_{t}+\boldsymbol{v}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ is the material derivative, $g$ is the gravitational acceleration, $\unicode[STIX]{x1D6FE}$ is the ratio of surface tension and the fluid density, $\unicode[STIX]{x1D708}$ is the viscosity of the fluid, and $\text{i}$ is the imaginary number such that $\text{i}^{2}=-1$ . In this equation, we included the effect of bulk dissipation, governed by the last term, but other sources of dissipation are generally present (Hunt Reference Hunt1964; Przadka et al. Reference Przadka, Cabane, Pagneux, Maurel and Petitjeans2012; Robertson & Rousseaux Reference Robertson and Rousseaux2017). To keep the discussion general, we introduce a pair of functions, $F$ and $\unicode[STIX]{x1D6E4}$ , encoding the dispersion relation and dissipation, in order to write the wave equation in the following form:

(2.2) $$\begin{eqnarray}\displaystyle {\mathcal{D}}_{t}^{2}\unicode[STIX]{x1D719}+F(-\text{i}\unicode[STIX]{x1D735})\unicode[STIX]{x1D719}-2\unicode[STIX]{x1D708}\unicode[STIX]{x1D6E4}(-\text{i}\unicode[STIX]{x1D735}){\mathcal{D}}_{t}\unicode[STIX]{x1D719}=0. & & \displaystyle\end{eqnarray}$$

In the absence of background flow, solutions of (2.1) are given by plane waves $\unicode[STIX]{x1D719}=\text{Re}[A\exp (-\text{i}\unicode[STIX]{x1D714}t+\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x})]$ of amplitude $A$ , angular frequency $\unicode[STIX]{x1D714}$ and wavevector $\boldsymbol{k}$ (of norm $k=\Vert \boldsymbol{k}\Vert$ ). Then $\unicode[STIX]{x1D714}$ and $\boldsymbol{k}$ obey the dispersion relation

(2.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}^{2}=F(\boldsymbol{k})-2\text{i}\unicode[STIX]{x1D708}\unicode[STIX]{x1D714}\unicode[STIX]{x1D6E4}(\boldsymbol{k}), & & \displaystyle\end{eqnarray}$$

with the functions

(2.4a,b ) $$\begin{eqnarray}\displaystyle F(\boldsymbol{k})=(gk+\unicode[STIX]{x1D6FE}k^{3})\tanh (h_{0}k)\quad \text{and}\quad \unicode[STIX]{x1D6E4}(k)=k^{2}. & & \displaystyle\end{eqnarray}$$

Note that the method we shall present stays identical for general functions $F$ and $\unicode[STIX]{x1D6E4}$ , the only necessary assumptions being that they are analytic and even. In the absence of background flow, a wave packet propagates at a group velocity $\boldsymbol{v}_{g}$ given by

(2.5) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{g}\doteq \unicode[STIX]{x1D735}_{k}\unicode[STIX]{x1D714}=\unicode[STIX]{x1D735}_{k}\sqrt{F}. & & \displaystyle\end{eqnarray}$$

In the presence of a constant flow, one must distinguish between two notions of (angular) frequency: the intrinsic (fluid-frame) frequency $\unicode[STIX]{x1D6FA}$ and the laboratory frequency $\unicode[STIX]{x1D714}$ . The intrinsic frequency is the one satisfying the dispersion relation, while the laboratory one is conserved due to stationarity. The two are related by a Doppler shift: $\unicode[STIX]{x1D6FA}=\unicode[STIX]{x1D714}-\boldsymbol{v}_{0}\boldsymbol{\cdot }\boldsymbol{k}$ . As a result, the group velocity in the laboratory frame is defined as

(2.6) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{g}^{lab}=\unicode[STIX]{x1D735}_{k}\unicode[STIX]{x1D714}, & & \displaystyle\end{eqnarray}$$

and obtained from the group velocity in the fluid frame using

(2.7) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{g}^{f}\doteq \unicode[STIX]{x1D735}_{k}\sqrt{F}=\boldsymbol{v}_{g}^{lab}-\boldsymbol{v}_{0}. & & \displaystyle\end{eqnarray}$$

Now when the background flow varies spatially, solutions of the wave equation become much more involved. This is why we shall now solve (2.1) using a ray-tracing method, to characterize the propagation of waves around vortices.

3 Ray-tracing equations

We shall now solve (2.1) using a gradient expansion – a common method in wave physics (see e.g. Bühler Reference Bühler2014). This means that we assume that the background flow changes over a scale significantly larger than the wavelength. To materialize this, a convenient procedure is to rescale the derivatives in (2.1) as $\unicode[STIX]{x2202}\rightarrow \unicode[STIX]{x1D716}\unicode[STIX]{x2202}$ , and solve the equation perturbatively in $\unicode[STIX]{x1D716}$ . We then look for wave solutions with a rapidly varying phase and a slowly varying amplitude, that is, we seek solutions of the form

(3.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}=A(\boldsymbol{x})\exp \!\left(\text{i}\frac{S(\boldsymbol{x})}{\unicode[STIX]{x1D716}}\right), & & \displaystyle\end{eqnarray}$$

where $A$ is the local amplitude and $S$ is the local phase, both real-valued functions (of course, ultimately, one should take the real part of (3.1) to obtain the velocity potential fluctuation field). We now expand $S$ and $A$ in powers of $\unicode[STIX]{x1D716}$ , and obtain an equation for the phase and one for the amplitude (in appendix A, we detail how to obtain them). At leading order in $\unicode[STIX]{x1D716}$ , the phase $S_{0}$ obeys the Hamilton–Jacobi equation,

(3.2) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x2202}_{t}S_{0}+\boldsymbol{v}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}S_{0})^{2}-F(\unicode[STIX]{x1D735}S_{0})=0. & & \displaystyle\end{eqnarray}$$

This first-order partial differential equation can be solved by the method of (bi)characteristics. The characteristics are the rays of the waves. Substituting the $\unicode[STIX]{x1D735}S_{0}$ by the wavevector $\boldsymbol{k}$ and $\unicode[STIX]{x2202}_{t}S_{0}$ by $-\unicode[STIX]{x1D714}$ , we obtain the Hamiltonian of our system,

(3.3) $$\begin{eqnarray}\displaystyle {\mathcal{H}}=-{\textstyle \frac{1}{2}}(\unicode[STIX]{x1D714}-\boldsymbol{v}_{0}\boldsymbol{\cdot }\boldsymbol{k})^{2}+{\textstyle \frac{1}{2}}F(\boldsymbol{k}). & & \displaystyle\end{eqnarray}$$

The characteristics are now obtained through Hamilton’s equations:

(3.4a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\dot{t}}=-\frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}},\quad \dot{\unicode[STIX]{x1D714}}=\frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}t}, & \displaystyle\end{eqnarray}$$

(3.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{x_{j}}=\frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}k_{j}}\quad \text{and}\quad \dot{k_{j}}=-\frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}x_{j}}, & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D706}$

$j$

$(x,y)$

$t$

$-\unicode[STIX]{x1D714}$

$\unicode[STIX]{x1D714}$

(3.5) $$\begin{eqnarray}\displaystyle {\mathcal{H}}=0, & & \displaystyle\end{eqnarray}$$

is satisfied for all solutions. A satisfying feature of this constrained Hamiltonian formalism is that it is independent of the choice of parameter. One can multiply the Hamiltonian ${\mathcal{H}}$ by any non-vanishing function $N(\unicode[STIX]{x1D706})$ , and one will find the same set of rays, instead parametrized by $\tilde{\unicode[STIX]{x1D706}}$ , where

(3.6) $$\begin{eqnarray}\displaystyle \text{d}\tilde{\unicode[STIX]{x1D706}}=N(\unicode[STIX]{x1D706})\,\text{d}\unicode[STIX]{x1D706}. & & \displaystyle\end{eqnarray}$$

The next-to-leading-order expansion of (2.1) in $\unicode[STIX]{x1D716}$ gives us the equation for the amplitude. This takes the form of a transport equation,

(3.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}(\unicode[STIX]{x1D6FA}_{0}A_{0}^{2})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FA}_{0}A_{0}^{2}\boldsymbol{v}_{g}^{lab})=-2\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FA}_{0}A_{0}^{2}\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D735}S_{0}). & & \displaystyle\end{eqnarray}$$

Note that this equation can be solved once one has a solution $S_{0}$ of the Hamilton–Jacobi equation (3.2). In particular, the group velocity in (3.7) is a function of only the position: $\boldsymbol{v}_{g}(\boldsymbol{x},p)\rightarrow \boldsymbol{v}_{g}(\boldsymbol{x},\unicode[STIX]{x1D735}S_{0}(\boldsymbol{x}))$ . The quantity $A_{0}^{2}\unicode[STIX]{x1D6FA}_{0}$ is in fact the wave action (Bühler Reference Bühler2014), up to a factor of $\unicode[STIX]{x1D70C}/g$ . This equation implies that the wave action is moving with the velocity $\boldsymbol{v}_{g}^{lab}$ . The dissipative term simply encodes the loss of wave action carried along the propagation. Since dissipation does not affect the equations (3.4) that determine the paths of the rays, we neglect it in the following.

Notice that we have so far considered a constant surface height. However, it is formally simple to add a varying height $h_{0}(x,y)$ at the level of Hamilton’s equations (3.4). Such replacement will be valid if the variations of the surface height are slow enough. By this we mean, first, that the relative variations of the height are small $|\unicode[STIX]{x1D735}h_{0}|\ll 1$ , so that its gradient can be neglected in the derivation of the wave equation (2.1); and, second, that the wavelength is small compared to the variation scale $|\unicode[STIX]{x1D735}h_{0}/(kh_{0})|$ , so that the eikonal approximation stays valid. In the experiment of § 4.2.1, the assumption of constant height is well satisfied outside a small radius near the centre, and hence we shall assume it in the rest of the paper.

4 Rays around a vortex

We shall now apply the ray-tracing equations to study the propagation of waves around a vortex. For this we consider the simplest model of an irrotational vortex with a drain. This model assumes that the surface is flat, and that the flow is homogeneous in the vertical direction and symmetric around the vertical axis. Under these assumptions, the irrotational and incompressibility conditions give us the profile

(4.1) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{0}=-\frac{D}{r}\boldsymbol{e}_{r}+\frac{C}{r}\boldsymbol{e}_{\unicode[STIX]{x1D703}}, & & \displaystyle\end{eqnarray}$$

written in polar coordinates $(r,\unicode[STIX]{x1D703})$ (see figure 1), and where $C$ and $D$ are the circulation and draining constants (positive), respectively. The schematic of the flow and the coordinate systems are presented in figure 1. This model is not only simple to analyse, it also provides a realistic description of a wide variety of situations. The reason is that, far from the centre (where the water depth drops), and far from the edges of the fluid tank (where circular symmetry breaks down), the free surface tends to be flat and the flow to resemble that of (4.1). As we will see, the effect we are after is contained in this region (see also figure 4 in Torres et al. (Reference Torres, Patrick, Coutant, Richartz, Tedford and Weinfurtner2017)). This model also corresponds (in the shallow-water regime) to a water-wave analogue of a rotating black hole (Visser Reference Visser1998; Schützhold & Unruh Reference Schützhold and Unruh2002).

Figure 1. Schematic of the flow and associated coordinate systems, from the perspective of an observer looking down on the wavetank. The vortex considered here is rotating anticlockwise when $C>0$ (as shown by the bold blue arrows). The light grey lines represent streamlines of the flow. A point M on the surface is located using either Cartesian coordinates $(x,y)$ or polar coordinates $(r,\unicode[STIX]{x1D703})$ .

4.1 Circular orbits

Although we shall have in mind the flow profile of (4.1), our results stay conceptually identical for a large class of flows. To see this, the following discussion is done using a general symmetric flow of the form

(4.2) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{0}=v_{r}(r)\boldsymbol{e}_{r}+v_{\unicode[STIX]{x1D703}}(r)\boldsymbol{e}_{\unicode[STIX]{x1D703}}. & & \displaystyle\end{eqnarray}$$

Using the symmetry of the vortex flow, we first switch to polar coordinates. We build the conjugate variables as before, that is, $k_{r}=\unicode[STIX]{x2202}_{r}S$ and $m=\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}S$ . As we shall see, depending on the context, $m$ will be a real number (interpreted as an angular momentum) or an integer (interpreted as the azimuthal number of the wave). We also note that $\unicode[STIX]{x1D714}$ is conserved since the flow is stationary. For simplicity, we assume $\unicode[STIX]{x1D714}>0$ (but negative values can be obtained via the symmetry $\unicode[STIX]{x1D714}\rightarrow -\unicode[STIX]{x1D714}$ and $m\rightarrow -m$ ).

A circular orbit is an equilibrium point in the radial direction. This means that it is a critical point of the Hamiltonian for $(r,k_{r})$ . It thus satisfies the conditions

(4.3a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{r}{\mathcal{H}}=0, & \displaystyle\end{eqnarray}$$

(4.3b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{k_{r}}{\mathcal{H}}=0. & \displaystyle\end{eqnarray}$$

In the bathtub profile of (4.1), equations (4.3) give two relations between four unknowns. For a fixed pair $(\unicode[STIX]{x1D714},m)$ , equations (4.3) determine a unique pair $(r_{\star },k_{r\star })$ . Furthermore, imposing the Hamiltonian constraint ${\mathcal{H}}=0$ gives a relation between $\unicode[STIX]{x1D714}$ and $m$ on the circular orbit, analogous to a dispersion relation, which we write as

(4.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{\star }(m). & & \displaystyle\end{eqnarray}$$

This relation can be read in either of two ways: at fixed angular frequency $\unicode[STIX]{x1D714}$ , it gives the value of $m$ such that one has a circular orbit; at fixed azimuthal number $m$ , it yields the particular wave frequency corresponding to the circular orbit.

4.1.1 Shallow-water regime

In the limit of shallow water, where all modes propagate at the same speed, the radius of the circular orbit is independent of the angular frequency $\unicode[STIX]{x1D714}$ , and furthermore the relationship between $\unicode[STIX]{x1D714}$ and $m$ is linear (see e.g. Dolan, Oliveira & Crispino Reference Dolan, Oliveira and Crispino2012; Dempsey Reference Dempsey2017):

(4.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{\star }=\hat{\unicode[STIX]{x1D6FA}}^{\pm }m. & & \displaystyle\end{eqnarray}$$

Here $\hat{\unicode[STIX]{x1D6FA}}^{\pm }\equiv \dot{\unicode[STIX]{x1D703}}/{\dot{t}}$ is the orbital frequency of the circular orbit. In general, there exist two circular orbits, one co-rotating ( $\dot{\unicode[STIX]{x1D703}}>0$ , $\hat{\unicode[STIX]{x1D6FA}}^{+}>0$ ) and one counter-rotating ( $\dot{\unicode[STIX]{x1D703}}<0$ , $\hat{\unicode[STIX]{x1D6FA}}^{-}<0$ ) relative to the circulating flow. The frequencies are given by

(4.6) $$\begin{eqnarray}\displaystyle \hat{\unicode[STIX]{x1D6FA}}^{\pm }=\pm \frac{c^{2}\sqrt{C^{2}+D^{2}}}{B_{\pm }^{2}}. & & \displaystyle\end{eqnarray}$$

The radii of these circular orbits (which we shall now refer to as ‘orbital radii’) are then given by

(4.7) $$\begin{eqnarray}\displaystyle r_{\star }^{\pm }=\frac{B_{\pm }}{c}, & & \displaystyle\end{eqnarray}$$

where $c=\sqrt{gh}$ is the propagation speed of shallow-water waves, and a $+$ ( $-$ ) sign denotes the co-rotating (counter-rotating) case. The parameter $B_{\pm }$ is defined as

(4.8) $$\begin{eqnarray}\displaystyle B_{\pm }=(2(C^{2}+D^{2})\mp 2C\sqrt{C^{2}+D^{2}}\,)^{1/2}. & & \displaystyle\end{eqnarray}$$

As we shall see, $B_{\pm }$ also governs the orbital radius for deep-water waves. An immediate conclusion of (4.7) is that the co-rotating circular orbit is closer to the centre of the vortex, while the counter-rotating orbit is further out. Hence, the counter-rotating orbit is in general more visible (this is further confirmed in § 5 when studying damped resonances).

4.1.2 Deep-water regime

For further insight into circular orbits in the presence of dispersion, it is instructive to look at another simple case where analytic formulae can be obtained. Here we consider deep water without capillarity, and thus a group velocity given by

(4.9) $$\begin{eqnarray}\displaystyle v_{g}(k)=\frac{1}{2}\sqrt{\frac{g}{k}}. & & \displaystyle\end{eqnarray}$$

In this regime it is possible to derive a simple formula for the radius of the orbits as a function of the frequency. We can directly express the norm of the wavevector in terms of the radius as

(4.10) $$\begin{eqnarray}\displaystyle |\boldsymbol{k}|=\frac{gr_{\star }^{2}}{4B_{\pm }^{2}} & & \displaystyle\end{eqnarray}$$

(with $\pm$ indices omitted on $k$ and $r_{\star }$ for clarity). With the other conditions for the circular orbits, equations (4.3), we can express both $k_{r\star }$ and $m$ as functions of $r_{\star }$ , viz.,

(4.11a,b ) $$\begin{eqnarray}\displaystyle k_{r\star }=\frac{Dgr_{\star }^{2}}{4B_{\pm }^{3}}\quad \text{and}\quad m=\frac{gr_{\star }^{3}}{4B_{\pm }^{2}}\sqrt{1-\frac{D^{2}}{B_{\pm }^{2}}}. & & \displaystyle\end{eqnarray}$$

Combining these expressions and using the Hamiltonian constraint, we obtain the radii of the unstable orbits as a function of the frequency and the effective dispersion relation (see appendix B for detail):

(4.12) $$\begin{eqnarray}\displaystyle r_{\star }=\frac{8B_{\pm }}{3g}\unicode[STIX]{x1D714}_{\star } & & \displaystyle\end{eqnarray}$$

and

(4.13) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{\star }(m)=\frac{3}{8}\left(\frac{4g^{2}}{\sqrt{B_{\pm }^{2}-D^{2}}}\right)^{1/3}m^{1/3}. & & \displaystyle\end{eqnarray}$$

Note that the wave angular frequency $\unicode[STIX]{x1D714}_{\star }$ scales linearly with $m$ in shallow water, but with $m^{1/3}$ in deep water. This is a direct consequence of the dispersive nature of the waves.

4.1.3 General case

For any nonlinear dispersion relation, equations (4.3) imply the relation

(4.14) $$\begin{eqnarray}\displaystyle r_{\star }=\frac{B_{\pm }}{v_{g}^{f}}. & & \displaystyle\end{eqnarray}$$

This equation shows that it is the change in the group velocity $v_{g}^{f}$ that shifts the location of the circular orbit. Here $v_{g}^{f}$ is evaluated with the local wavevector on the orbit $k^{2}=k_{r\star }^{2}+m^{2}/r_{\star }^{2}$ , and hence the equation above is implicit. To find $r_{\star }$ and $\unicode[STIX]{x1D714}_{\star }(m)$ explicitly, one must insert (4.14) into (4.3) and the Hamiltonian constraint.

Figure 2 shows numerically obtained solutions for the orbital radius as a function of wave frequency, for a vortex with $C/hc=D/hc=0.9$ . The two branches correspond to the co-rotating ( $+$ ) and counter-rotating ( $-$ ) orbits. At low frequency (long wavelength), the orbital radius is anticipated by the shallow-water result, equation (4.5). At higher frequencies, the circular orbits migrate away from the vortex core, as anticipated in (4.12). Even though the ray structure will qualitatively be similar for various values of the parameters $C$ and $D$ , one can distinguish two interesting limits. First, when $C\gg D$ , we can see that the radius of the co-rotating orbit shrinks to zero as $B_{+}\rightarrow 0$ . In contrast, when $C\ll D$ , the two circular orbits asymptotically approach one another as $B_{+}\rightarrow B_{-}$ .

At this level, it is instructive to discuss the case of a vortex without a drain. By taking the limit $D\rightarrow 0$ in (4.7) and (4.12), we see that the co-rotating orbit degenerates to $r=0$ , but the counter-rotating one still exists at a non-zero radius (e.g.  $C/2c$ in the shallow-water regime). For sound waves, it was shown by Nazarenko (Reference Nazarenko1994) and Nazarenko et al. (Reference Nazarenko, Zabusky and Scheidegger1995) that rays approaching $r=0$ of an ideal non-draining vortex (equation (4.1) with $D=0$ ) see their wavelength decrease to zero, and ultimately reach the dissipative scale. As such rays do not escape, it suggests that the scattering will be similar to a vortex with small drain (that is $D\neq 0$ but $D\ll C$ ). In particular, we expect to have a counter-rotating orbit, as well as the associated resonance effects (see § 5). (We thank an anonymous referee for having stimulated this discussion, and for pointing out references Nazarenko (Reference Nazarenko1994) and Nazarenko et al. (Reference Nazarenko, Zabusky and Scheidegger1995).)

Figure 2. Radius of the unstable orbit as a function of the angular frequency of the co- and counter-rotating rays (in purple and orange, respectively) and comparison with the shallow- and deep-water regimes (in solid and dashed black, respectively). We have chosen the flow parameters such that $C/hc=D/hc=0.9$ (chosen to make the transition from shallow to deep water more visible), where $c=\sqrt{gh}$ is the celerity of shallow-water waves. On the vertical axis, $R_{rel}^{\pm }$ is the value of $r_{\star }^{\pm }$ in the shallow-water (or relativistic) regime. (We recall that, when the frequency varies, the corresponding orbital $m$ also does, so that (4.4) is satisfied.)

4.2 Ray scattering

We now study a family of characteristics (rays) encroaching upon the vortex from the right ( $x\rightarrow +\infty$ ) with angular frequency $\unicode[STIX]{x1D714}$ and various impact parameters $b$ , which is the ray equivalent of sending a plane wave. Here, $m$ is a continuous parameter that is given in terms of the impact parameter via the relation

(4.15) $$\begin{eqnarray}\displaystyle b=\frac{m}{k_{in}}+\frac{C}{v_{g}^{in}}. & & \displaystyle\end{eqnarray}$$

Here $k_{in}$ is the incoming wavenumber satisfying the dispersion relation (or the Hamiltonian constraint) at infinity for a given $\unicode[STIX]{x1D714}$ and $v_{g}^{in}$ is the group velocity at infinity (as the flow is negligible at infinity, the group velocities in the fluid frame and in the laboratory frame are identical). We numerically solve Hamilton’s equations (3.4) to find the characteristics. A standard fourth-order Runge–Kutta (RK4) scheme is used. To ensure the validity of our numerical simulation, the step size is chosen such that along the rays the adimensional quantity $|{\mathcal{H}}/\unicode[STIX]{x1D714}^{2}|$ is smaller than $10^{-7}$ . Our solution therefore satisfies the Hamiltonian constraint. To numerically compute the amplitude of the wave, we use (3.7) with $\unicode[STIX]{x1D708}=0$ . This shows that the flux of wave action is conserved along a tube of rays. Numerically, we use this conservation to obtain the change in amplitude between neighbouring points along two neighbouring rays, by ensuring that the product of the distance between the rays and the wave action is constant. This means that, if the rays converge (diverge), the amplitude will increase (decrease).

Figure 3(a) depicts a congruence of rays at fixed angular frequency $\unicode[STIX]{x1D714}=19.8~\text{rad}~\text{s}^{-1}$ impinging on a draining bathtub vortex with parameters $C=1.6\times 10^{-2}~\text{m}^{2}~\text{s}^{-1}$ , $D=1\times 10^{-3}~\text{m}^{2}~\text{s}^{-1}$ and $h=0.06~\text{m}$ , corresponding to the laboratory study of Torres et al. (Reference Torres, Patrick, Coutant, Richartz, Tedford and Weinfurtner2017). We can clearly distinguish two types of rays. The first type are rays that are able to escape to infinity (co-rotating in red and counter-rotating in dashed blue). The impact parameter for those rays are beyond some critical values $b_{c}^{\pm }(\unicode[STIX]{x1D714})$ . We note that $|b_{c}^{-}(\unicode[STIX]{x1D714})|>|b_{c}^{+}(\unicode[STIX]{x1D714})|$ , implying that co-rotating rays are able to go closer to the vortex than counter-rotating ones. The second type of rays are those that fall into the vortex core (in dotted brown). We note from (4.15) that the impact parameter depends on the group velocity and therefore on the frequency of the wave. This behaviour differs from the shallow-water regime where the group velocity matches the phase velocity and is constant for all frequency.

In the process of obtaining the rays by solving Hamilton’s equations, we also compute the phase of the wave along each trajectory. It is then possible to reconstruct the eikonal wavefronts by finding the constant phase points along each ray (see Dempsey & Dolan (Reference Dempsey and Dolan2016) for a shallow-water study). The wavefronts are presented in figure 3(b). The coloured dots represent the location of the constant phase points and the colour scale gives the amplitude of the wave at these points. Far from the vortex, where the flow becomes negligible, we can see that the wavefronts are orthogonal to the rays (black lines). On the contrary, closer to the vortex, we clearly observe the non-orthogonality of the wavefronts and the rays where the flow becomes more rapid.

Figure 3. (a) A congruence of rays incoming from far right on a draining vortex flow. The flow parameters are $C=1.6\times 10^{-2}~\text{m}^{2}~\text{s}^{-1}$ , $D=1\times 10^{-3}~\text{m}^{2}~\text{s}^{-1}$ and $h=0.06~\text{m}$ . They correspond to the experimental flow of Torres et al. (Reference Torres, Patrick, Coutant, Richartz, Tedford and Weinfurtner2017). The angular frequency of the wave is $\unicode[STIX]{x1D714}=19.8~\text{rad}~\text{s}^{-1}$ . For these parameters, we have $hk_{in}\simeq 2.4$ and hence are in the deep-water regime. The red (blue) curves represent the co-rotating (counter-rotating) rays escaping to infinity. The dotted brown lines are the rays falling into the hole. The two black circles represent the orbital radius for co- and counter-rotating rays, respectively, at $r=1$  cm and $r=15.3~\text{cm}$ . (b) The eikonal wavefront reconstructed from the phase along the rays. The colour scale represents the amplitude of the wave normalized to one initially. The black lines are some rays from (3 a).

When varying the angular frequency $\unicode[STIX]{x1D714}$ , the rays display a similar pattern as in figure 3(a), but the orbital radii $r_{\star }^{\pm }$ vary. When increasing the angular frequency, they interpolate between their shallow-water value to the deep-water behaviour of (4.12). To illustrate this, we plotted the dependence of the orbital radii with the frequency in figure 2.

4.2.1 Experimental validation

In figure 4 we compare our ray-tracing results with experimental data taken from the wavetank experiment of Torres et al. (Reference Torres, Patrick, Coutant, Richartz, Tedford and Weinfurtner2017) with a circulation-to-draining ratio of approximately $C/D\simeq 16$ . One can see that, broadly, the eikonal wavefronts agree well with the experiment. We observe small deviations in two regions: close to the centre of the vortex, and after the wave has propagated through the vortex (on the left of the image). The eikonal approximation we used is expected to degrade close to the centre of the vortex, where the flow varies more rapidly. Moreover, cumulative errors might explain the shift between our simulation and the data after the wave has propagated away from the vortex (on the far left side of figure 4).

Figure 4. Comparison between eikonal wavefront computed numerically and experimental data. The bright yellow dots represent the eikonal wavefront reconstructed from the phase along the rays. The two black circles are the unstable orbits, also present in figure 3(a). The background image is a measurement of the free surface of the water. The flow parameters and wave frequency are the same as in figure 3. The frequency of the eikonal wave ( $3.15$  Hz) is chosen so as to obtain the best fit with the $3.27$  Hz wave of figure 1 in Torres et al. (Reference Torres, Patrick, Coutant, Richartz, Tedford and Weinfurtner2017). The two frequencies agree within error bars (approximately $4\,\%$ ). The colour bar represents the amplitude of the wave in metres.

The good agreement between our numerical solution and the experimental data suggests that circular rays are present in the system (see § 4.1). While the inner one (co-rotating with the vortex) is very close to the centre, in a region where the vortex model of (4.1) becomes questionable, the outer one (counter-rotating) lies in a region where our method provides an accurate description. To estimate the order of the error, we compare the gradient scale of the background to the frequency of the wave: $|\unicode[STIX]{x2202}_{r}v_{\unicode[STIX]{x1D703}}|/\unicode[STIX]{x1D714}$ . Near the outer light ring $r_{\star }^{+}$ , this ratio reduces to the inverse of the azimuthal number $m_{\star }$ and is approximately $0.1$ . As a last remark, we point out that, within our approximation, the locations of the rays are usually more accurate than the phase itself. This is because the phase is rapidly varying, but the local momentum is slowly varying (see the discussion in § 4.3.4 of Bühler (Reference Bühler2014)). Hence, this fact and the good agreement in figure 4 strengthen our conclusion regarding the presence of circular orbits around the vortex of figure 4.

5 Characteristic damped resonances of a draining vortex: the quasi-normal mode spectrum

We anticipate that the presence of ‘rings’ (i.e. circular rays/characteristics) in the vortex system will influence the time-dependent response to a small wave excitation, through certain characteristic damped oscillations. In the black hole case, the time-dependent response of a black hole to a perturbation typically bears the imprint of one or more quasi-normal modes (QNMs) (Berti et al. Reference Berti, Cardoso and Starinets2009). The QNMs are a discrete set of modes with complex frequencies, with the real part determining the oscillation frequency, and the (negative) imaginary part determining the damping rate. Formally, QNMs are defined via a pair of boundary conditions to be those modes which are ingoing at the horizon and outgoing far from the black hole. In the eikonal limit, one finds that the spectrum is governed chiefly by the properties of the circular orbits (see e.g. Cardoso et al. Reference Cardoso, Miranda, Berti, Witek and Zanchin2009). More precisely, the eikonal complex frequencies are given by

(5.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{QNM}(m)=\unicode[STIX]{x1D714}_{\star }(m)-\text{i}\unicode[STIX]{x1D6EC}(m)(n+{\textstyle \frac{1}{2}}), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{\star }(m)$ is the circular orbit wave frequency of § 4.1, $n$ is an integer called the overtone index, and $\unicode[STIX]{x1D6EC}$ is the decay rate of these (unstable) orbits (Lyapunov exponent), to be defined below.

The notion of QNMs extends to many different contexts. QNMs typically arise in open systems, or systems with leakage (in our case, the drain plays that role). As we shall see, around a draining vortex, surface wave decay is also driven by a discrete set of complex frequencies, which can be obtained approximately using the structure of circular orbits.

Owing to the presence of circular orbits, we expect that, if one perturbs the system initially, the perturbations will spread out rapidly, except where the perturbation lingers around the circular orbit. Indeed, a wave packet passing close to such a ‘ring’ can hover around it for a longer time before dispersing. To understand this behaviour more precisely, we consider a wave packet in the neighbourhood of a circular orbit, that is, around the family of rays $(r_{\star }+\unicode[STIX]{x1D6FF}r,~k_{r\star }+\unicode[STIX]{x1D6FF}k,~\unicode[STIX]{x1D714}_{\star }+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D714})$ at a fixed azimuthal number $m$ . To obtain the effective dynamics of these rays, we expand the Hamiltonian in the vicinity of the critical point $(r_{\star },k_{r\star })$ . This gives

(5.2) $$\begin{eqnarray}\displaystyle {\mathcal{H}}=\unicode[STIX]{x2202}_{\unicode[STIX]{x1D714}}{\mathcal{H}}\,\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D714}+{\textstyle \frac{1}{2}}\boldsymbol{X}^{\text{T}}\boldsymbol{\cdot }[\text{d}^{2}{\mathcal{H}}]\boldsymbol{\cdot }\boldsymbol{X}, & & \displaystyle\end{eqnarray}$$

where we have defined the Hessian matrix

(5.3a,b ) $$\begin{eqnarray}\displaystyle [\text{d}^{2}{\mathcal{H}}]=\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x2202}_{k}^{2}{\mathcal{H}} & \unicode[STIX]{x2202}_{k}\unicode[STIX]{x2202}_{r}{\mathcal{H}}\\ \unicode[STIX]{x2202}_{k}\unicode[STIX]{x2202}_{r}{\mathcal{H}} & \unicode[STIX]{x2202}_{r}^{2}{\mathcal{H}}\end{array}\right)\quad \text{and}\quad \boldsymbol{X}=\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D6FF}k\\ \unicode[STIX]{x1D6FF}r\end{array}\right). & & \displaystyle\end{eqnarray}$$

Since $[\text{d}^{2}{\mathcal{H}}]$ is a symmetric and real matrix, it can be diagonalized. We call $(\unicode[STIX]{x1D707}_{1},\unicode[STIX]{x1D707}_{2})$ the eigenvalues, and $(R,K)$ the components in the eigenbasis. (In addition, the eigenbasis is orthogonal, which in two dimensions is enough to ensure a canonical transformation $(\unicode[STIX]{x1D6FF}r,\unicode[STIX]{x1D6FF}k)\rightarrow (R,K)$ .) In this representation, the Hamiltonian becomes

(5.4) $$\begin{eqnarray}\displaystyle {\mathcal{H}}=\unicode[STIX]{x2202}_{\unicode[STIX]{x1D714}}{\mathcal{H}}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D714}+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D707}_{1}K^{2}+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D707}_{2}R^{2}. & & \displaystyle\end{eqnarray}$$

This is the Hamiltonian of a harmonic oscillator. The relative sign of the coefficients $\unicode[STIX]{x1D707}_{1}$ and $\unicode[STIX]{x1D707}_{2}$ will determine the stability of the orbits. In our case, we always have $\unicode[STIX]{x1D707}_{1}\unicode[STIX]{x1D707}_{2}<0$ and hence the orbit is unstable. In its vicinity, wave dynamics is described by an inverted harmonic oscillator. To obtain the Lyapunov exponent $\unicode[STIX]{x1D6EC}$ , we solve Hamilton’s equations with this reduced Hamiltonian. This gives

(5.5a-c ) $$\begin{eqnarray}\displaystyle {\dot{R}}=\unicode[STIX]{x1D707}_{1}K,\quad \dot{K}=-\unicode[STIX]{x1D707}_{2}R,\quad {\dot{t}}=-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D714}}{\mathcal{H}}. & & \displaystyle\end{eqnarray}$$

The general solution is a linear superposition of decaying and growing exponentials. At late times, the growing behaviour dominates, and we have

(5.6a,b ) $$\begin{eqnarray}\displaystyle R\sim A\text{e}^{(\sqrt{-\unicode[STIX]{x1D707}_{1}\unicode[STIX]{x1D707}_{2}}/\unicode[STIX]{x2202}_{\unicode[STIX]{x1D714}}{\mathcal{H}})t}\quad \text{and}\quad K\sim B\text{e}^{(\sqrt{-\unicode[STIX]{x1D707}_{1}\unicode[STIX]{x1D707}_{2}}/\unicode[STIX]{x2202}_{\unicode[STIX]{x1D714}}{\mathcal{H}})t}, & & \displaystyle\end{eqnarray}$$

where $A$ and $B$ are constant coefficients such that (5.5) is satisfied. The Lyapunov exponent is the rate of exponential growth away from the orbit, and is given by

(5.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6EC}=\frac{\sqrt{-\unicode[STIX]{x1D707}_{1}\unicode[STIX]{x1D707}_{2}}}{|\unicode[STIX]{x2202}_{\unicode[STIX]{x1D714}}{\mathcal{H}}|}=\frac{\sqrt{-\text{det}[\text{d}^{2}{\mathcal{H}}]}}{|\unicode[STIX]{x2202}_{\unicode[STIX]{x1D714}}{\mathcal{H}}|}. & & \displaystyle\end{eqnarray}$$

The Lyapunov exponent governs the decay of waves close to the circular orbit. One can see this by inspecting the amplitude governed by the transport equation (3.7), namely,

(5.8) $$\begin{eqnarray}\displaystyle A_{0}(t)\propto \text{e}^{-\unicode[STIX]{x1D6EC}t/2}. & & \displaystyle\end{eqnarray}$$

However, this argument misses the overtone frequencies (i.e.  $n\neq 0$ in (5.1)).

A more precise argument is to lift the reduced Hamiltonian (5.4) at the level of the wave equation. To see this, we consider a wave packet of fixed $m$ , localized around the circular orbit, written as

(5.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}(t,r_{\star }+\unicode[STIX]{x1D6FF}r)\sim \unicode[STIX]{x1D713}(t,\unicode[STIX]{x1D6FF}r)\text{e}^{-\text{i}\unicode[STIX]{x1D714}_{\star }t+\text{i}k_{r\star }r}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D713}$ is a slowly varying envelope. Formally, the wave equation (2.1) can be written as

(5.10) $$\begin{eqnarray}\displaystyle \hat{{\mathcal{H}}}(\unicode[STIX]{x1D714}=\text{i}\unicode[STIX]{x2202}_{t};r,k_{r}=-\text{i}\unicode[STIX]{x2202}_{r})\unicode[STIX]{x1D719}=0. & & \displaystyle\end{eqnarray}$$

(This writing is, of course, only formal, since there is a priori no unique way to promote the Hamiltonian function ${\mathcal{H}}(\unicode[STIX]{x1D714},r,k_{r})$ into an operator $\hat{{\mathcal{H}}}(\text{i}\unicode[STIX]{x2202}_{t};r,-\text{i}\unicode[STIX]{x2202}_{r})$ . The ambiguities essentially come from the various possible orderings between functions of $r$ and $\unicode[STIX]{x2202}_{r}$ . However, since we work here at the level of the eikonal approximation, different choices of ordering lead to the same result.) Assuming $\unicode[STIX]{x1D713}$ in (5.9) is slowly varying, one can again rescale $\unicode[STIX]{x2202}\rightarrow \unicode[STIX]{x1D716}\unicode[STIX]{x2202}$ and expand in $\unicode[STIX]{x1D716}$ . Doing so, the effective wave equation in the vicinity of the circular orbit is given by

(5.11) $$\begin{eqnarray}\displaystyle -\text{i}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D714}}{\mathcal{H}}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D713}=-{\textstyle \frac{1}{2}}\unicode[STIX]{x2202}_{k}^{2}{\mathcal{H}}\unicode[STIX]{x2202}_{r}^{2}\unicode[STIX]{x1D713}-{\textstyle \frac{1}{2}}\text{i}\unicode[STIX]{x2202}_{k}\unicode[STIX]{x2202}_{r}{\mathcal{H}}(r\unicode[STIX]{x2202}_{r}+\unicode[STIX]{x2202}_{r}r)\unicode[STIX]{x1D713}+{\textstyle \frac{1}{2}}\unicode[STIX]{x2202}_{r}^{2}{\mathcal{H}}r^{2}\unicode[STIX]{x1D713}. & & \displaystyle\end{eqnarray}$$

(Note that the cross-term is chosen to be symmetric.) This equation is a Schrödinger equation in an upside-down harmonic potential. To see this, we use the fact that there is a canonical transformation to diagonalize the Hamiltonian (5.2). At the level of the Schrödinger equation (5.11), this means that there is a unitary transformation $\unicode[STIX]{x1D713}\rightarrow \tilde{\unicode[STIX]{x1D713}}$ such that

(5.12) $$\begin{eqnarray}\displaystyle -\text{i}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D714}}{\mathcal{H}}\unicode[STIX]{x2202}_{t}\tilde{\unicode[STIX]{x1D713}}=-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D707}_{1}\unicode[STIX]{x2202}_{R}^{2}\tilde{\unicode[STIX]{x1D713}}+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D707}_{2}R^{2}\tilde{\unicode[STIX]{x1D713}}. & & \displaystyle\end{eqnarray}$$

In the case of (5.12) the time-dependent response is well described with the help of quasi-normal modes, defined as the solutions of $\text{i}\unicode[STIX]{x2202}_{t}\tilde{\unicode[STIX]{x1D713}}=\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D714}\tilde{\unicode[STIX]{x1D713}}$ with purely outgoing boundary conditions (see e.g. Kokkotas & Schmidt Reference Kokkotas and Schmidt1999). Such solutions are analytically known, and can be expressed using parametric cylinder functions. This selects a discrete set of complex values for $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D714}$ , which are the quasi-normal frequencies for $\tilde{\unicode[STIX]{x1D713}}$ and hence for $\unicode[STIX]{x1D713}$ . Finally, using (5.9), we see that the quasi-normal frequencies for $\unicode[STIX]{x1D719}$ are given by (5.1).

Alternative boundary conditions, such as those for exotic black hole mimickers, can alter the global structure of the quasi-normal mode spectrum, without significantly altering the local ‘ringing’ phenomena associated with the light ring (see e.g. Cardoso & Pani Reference Cardoso and Pani2017a ,Reference Cardoso and Pani b ).

We have numerically computed the fundamental ( $n=0$ ) eikonal quasi-normal frequencies in the full dispersive regime for a selection of azimuthal numbers and for several values of the ratio $C/D$ . The result and the comparison with the linear regime are presented in figure 5. We can see that the behaviour of the counter-rotating modes ( $m<0$ ) is very similar in the two different regimes. As the ratio increases, the radius of the unstable orbits of these modes will increase too. For fixed $m$ , the increase of the radius will increase the wavelength (decrease the frequency) and therefore will bring the mode closer to the linear regime. On the other side, the behaviour for co-rotating ( $m>0$ ) modes is very different when dispersion is included. Indeed, in the relativistic regime, the lifetime of the co-rotating modes decreases while it appears to increase after a short drop as an effect of dispersion.

In general, it is noticeable that the co-rotating modes are damped quicker than the counter-rotating ones (see figure 5). Hence, if the initial excitation has a similar amplitude on both radii, the dominant signal at large time will be that of the counter-rotating orbit. This is the case, for instance, if the system is excited by sending a plane wave, which contains as much $m>0$ as $m<0$ (as in figure 4). In addition, we expect the counter-rotating mode to have a larger spatial extension than the co-rotating one since the radius of the counter-rotating orbit is larger that the co-rotating one, and hence will be more visible.

Figure 5. Fundamental quasi-normal mode frequencies $\unicode[STIX]{x1D714}_{QNM}$ of different $m$ modes, $m=\pm 1$ (a) and $m=\pm 2$ (b), for various ratios $C/D$ . The dashed line is the shallow-water result, and the solid one the (general) dispersive regime. In order to make the transition more transparent, the depth is chosen so that $gh^{3}/D^{2}=0.04$ . Each point represents a different ratio for the flow parameters. The negative $m$ modes have a similar behaviour in the two regimes, namely their lifetime increases as the ratio $C/D$ increases. On the other hand, the lifetime of the positive $m$ modes increases too with the ratio $C/D$ in the dispersive regime while it monotonically decrease in the absence of dispersion.