1. Introduction
Rotating flows are rather ubiquitous in nature. Examples range from domestic flows such as the bathtub vortex, to geophysical flows and Jupiter's Great Red Spot, to astrophysical flows such as accretion disks. Rotating flows support the propagation of many different type of waves. These are often viewed as small perturbations to an equilibrium state, also known as a base flow, and the behaviour of such waves depends both on the particular velocity field of the base flow and on the geometry of the problem under consideration. One important class of such waves are interfacial surface waves. These are waves forming at the interface between two fluids with different physical properties, and they often are strongly localized close to the interface.
Our interest here is in free-surface waves on a vortex flow. In fact, our initial motivation was out of curiosity concerning vortices and waves interacting in a swimming pool. Changes in the surface height can be easily visualized by light and dark patterns on the swimming pool floor (for which the authors recommend a warm sunny climate and an outdoor swimming pool). An ‘experiment’ (Skipp Reference Skipp2020) created pairs of vortices by drawing a dinner plate through the water. The initial waves generated by this disturbance disperse rapidly, leaving remarkably long-lived vortices with surface waves that appear trapped in the vortex, but which propagate around the vortex in the opposite direction to the vortex flow. A photograph of the phenomenon, together with a schematic representation, is shown in figure 1; interested readers are encouraged to follow the reference for a video.
The physical scenario motivating this work. Surface waves are trapped within a rotating vortex flow, but propagate in the direction opposite to the direction of rotation. (a) Photograph of two vortices in a swimming pool (a still taken from the video of Skipp (Reference Skipp2020)). (b) Schematic of the flow and surface waves on the left-hand vortex. Black arrows show the base velocity field. Blue lines show the wave crests. Red arrows show the direction of propagation of the waves.
The study of surface waves scattered by a vortex flow in a horizontally unbounded domain in a two-dimensional shallow-water setting with a bathtub vortex as a background flow was considered by Patrick, Coutant & Richartz (Reference Patrick, Coutant and Richartz2018) and Patrick (Reference Patrick2019), in order to provide an analogy between fluid phenomena and black holes. Their study revealed the presence of modes that they called ‘quasibound states’: normal modes having a very slow decaying rate in time and a spatial structure which remains trapped within the vortex core region, similar to standing waves in a potential well. Such states are perhaps understandable for the bathtub vortex which includes an inflow towards the vortex. The term ‘quasibound’ references that temporally weakly damped modes are not strictly bounded, and in fact may have weak exponential growth at infinity in order to satisfy a causal radiation condition (Briggs Reference Briggs1964; Bers Reference Bers1983; Berti, Cardoso & Starinets Reference Berti, Cardoso and Starinets2009, pp. 17–18; Le Dizès & Billant Reference Le Dizès and Billant2009; Riedinger, Le Dizès & Meunier Reference Riedinger, Le Dizès and Meunier2010); despite this weak growth, such modes are highly localized to a finite region (as discussed in § 4.1), and here we simply refer to these as trapped modes. Any slow decay of trapped waves may be due to dissipation (such as through viscosity), or by the wave not being perfectly trapped and ‘leaking’ outgoing waves to infinity. Trapped waves could be excited either by an incoming propagating wave from infinity, or by any initial disturbance (such as a dinner plate in the video depicted in figure 1); here, we investigate the existence of trapped waves about a vortex, being agnostic as to the mechanism creating them in the first place. We do not consider the related wave scattering problem.
We will study free surface waves on rotating fluid whose velocity 
$\boldsymbol {U_0}$ is that of a Lamb–Oseen vortex
\begin{equation} \boldsymbol{U_0} = \boldsymbol{\hat{\theta}} \frac{\varGamma_0}{2{\rm \pi} r} (1-\exp\{{-r^2/a^2}\}), \end{equation}
where 
$r$ is the horizontal distance from the vortex centre, 
$\boldsymbol {\hat {\theta }}$ is a unit vector in the direction of rotation and 
$\varGamma _0$ and 
$a$ are constants (see, e.g. Drazin & Riley (Reference Drazin and Riley2006), for details of the derivation of this flow). The Lamb–Oseen vortex is in fact an exact solution to the viscous Navier–Stokes equations if we take 
$a^2 = 4\nu t$, where 
$\nu$ is the kinematic viscosity, although for short time scales compared with the viscous diffusion time scale, 
$a$ may be approximated as constant. Unlike the bathtub vortex, this base flow is not the result of vortex stretching by a downwards outflow. In that regard, it is more similar to a geophysical vortex, although here we do not consider a stratified fluid, nor do we assume the depth to be shallow. Unlike other vortices such as the Rankine vortex (made up of an inner solid body rotation and an outer potential vortex), which have been found to be unstable (e.g. Ford Reference Ford1994; Mougel, Fabre & Lacaze Reference Mougel, Fabre and Lacaze2014), the Lamb–Oseen vortex is usually considered to be a robust and stable coherent structure. Riedinger et al. (Reference Riedinger, Le Dizès and Meunier2010) showed that it can become unstable in a stratified environment, but with no stratification; Fabre, Sipp & Jacquin (Reference Fabre, Sipp and Jacquin2006) have shown that the Lamb–Oseen vortex is stable to axial-invariant perturbations. Their work highlighted the presence of critical layer waves, arising within a critical layer. However, because the critical layer provides a damping effect, they find the associated critical layer waves to be damped; we will not consider the critical layer further here.
Free surface rotating waves have previously been studied computationally by Mougel, Fabre & Lacaze (Reference Mougel, Fabre and Lacaze2015); Mougel et al. (Reference Mougel, Fabre, Lacaze and Bohr2017) in a confined geometry for background flows with particular forms, such as solid-body rotation or a potential vortex flow. However, since a potential vortex has an infinite velocity at its centre, and a solid-body rotation has an infinite velocity at infinity, neither solid-body rotation nor a potential vortex are appropriate for the case we consider here. Their studies reveal the presence of four main types of waves: surface gravity waves; inertial waves; Rossby waves; and centrifugal waves. The first three waves are observed for solid-body rotation, while centrifugal waves and gravity waves arise in a potential vortex flow. Depending on the problem, each of these waves can be easily recognized and distinguished by its spatial structure and the associated frequency of oscillation. Further examples of problems where surface gravity waves, inertial waves and Rossby waves appear can be found in Johnson (Reference Johnson1997), Greenspan (Reference Greenspan1969) and McWilliams (Reference McWilliams2006), respectively. An interesting feature of the study on the potential vortex flow is the possibility of an instability of the base solution – known as a polygon instability – due to an interaction between the centrifugal and gravity waves. This interaction mechanism was initially studied by Tophøj et al. (Reference Tophøj, Mougel, Bohr and Fabre2013) and subsequently by Mougel et al. (Reference Mougel, Fabre, Lacaze and Bohr2017). The stability of vortex flows in bounded domains is, however, markedly different to the stability of vortex flows in unbounded domains, such as considered here, as the instability in bounded domains found by Tophøj et al. (Reference Tophøj, Mougel, Bohr and Fabre2013) and Mougel et al. (Reference Mougel, Fabre, Lacaze and Bohr2017) involves the interaction of two types of wave modes, one of which is dependent on the outer boundary. Instabilities of vortices on horizontally unbounded domains have been found to occur within the shallow-water limit, as studied by Ford (Reference Ford1994). However, it is worth noting that the study of shallow-water rotating flows is rather different to the study of finite-depth rotating flows studied here, on two counts: firstly, shallow-water rotating flows can only rotate very slowly before the deformation in the free surface at the centre of the vortex touches the bottom boundary and a dry inner region is formed; and secondly, shallow-water waves are non-dispersive and have a fixed wave speed, while finite-depth water waves are dispersive such that any wave speed is available to the system, and for example a wave can exist whose speed matches the flow speed at a given location. Surface waves have been also studied by Hunt et al. (Reference Hunt, Parau, Vanden-Broeck and Papageorgiou2015), particularly considering the effects of an electric field on both linear and nonlinear inviscid, irrotational waves. Moreover, two-dimensional surface waves in electrohydrodynamics and magneohydrodynamics have been studied by Hunt & Vanden-Broeck (Reference Hunt and Vanden-Broeck2015), Hunt (Reference Hunt2019) and Hunt & Dutykh (Reference Hunt and Dutykh2021), although in the present work we do not consider any electric or magnetic field, but instead we consider the full three-dimensional linear response to a prescribed swirling flow.
The present work considers a laterally unbounded domain (i.e. a domain that is in principle of infinite radial extent), with a finite depth, and thus assumes neither a shallow water nor a deep-water limit. In order to numerically study wave propagation problems in an unbounded domain, a non-reflecting boundary condition (NRBC) needs to be imposed at infinity. This is relatively straightforward in a two-dimensional shallow-water setting, since the dispersion relation for surface waves is known, and hence a characteristics analysis (John Reference John1978) allows the imposition of an exact NRBC. Characteristic boundary conditions are, however, specific to classical hyperbolic wave equations (Engquist & Majda Reference Engquist and Majda1977; Bayliss & Turkel Reference Bayliss and Turkel1980; Keller & Givoli Reference Keller and Givoli1989; Grote & Keller Reference Grote and Keller1995; Grote Reference Grote2000). Such exact boundary conditions cannot be directly imposed on our three-dimensional problem, as our problem is not necessarily hyperbolic, and it is not possible to find real characteristics. Moreover, surface waves on non-shallow water are dispersive, with the propagation speed of each mode depending on its frequency, which significantly complicates a NRBC even for hyperbolic problems (e.g. Lindquist, Neta & Giraldo Reference Lindquist, Neta and Giraldo2012), and can limit the range of wavenumbers being correctly treated (e.g. Wellens & Borsboom Reference Wellens and Borsboom2020). One possible solution to this problem is provided by absorbing layer methods. The idea behind this class of techniques is to surround the physical region with computational layers through which waves are progressively damped such that they reach the edge of the computational domain with sufficiently small amplitude that any reflections back into the physical region will be negligible. The absorbing layers need to vary gradually, however, as otherwise waves can spuriously reflect from the edges of the absorbing layers. Absorbing layer methods have been applied initially to waves problems in electromagnetics (e.g. Berenger Reference Berenger1994) and then extended to other physical fields like elasticity (Chew & Liu Reference Chew and Liu1996; Hastings, Schneider & Broschat Reference Hastings, Schneider and Broschat1996) and fluid dynamics (Söderstrom, Karlsson & Museth Reference Söderstrom, Karlsson and Museth2010). A shortcoming of absorbing layer formulations lies in the number of new unknowns introduced in the problem; this may be reduced using the formulation of Sim (Reference Sim2010) in the case of second-order hyperbolic equations. Here, a related but different technique is developed for the absorbing layer formulation in the far field, without introducing further unknowns in the mathematical problem.
The paper is organized as follows. In § 2, a general description is given of the mathematical model used to study perturbations to a Lamb–Oseen swirling flow, including a mathematical description of NRBCs. The resulting eigenvalue problem is then solved using a numerical procedure described in § 3. The results of this numerical solution are described in § 4. Finally, resulting conclusions and possible future directions for continued research are highlighted in § 5.
2. Mathematical model
Assuming that viscosity is negligible over the time scales of interest here, the governing equations are the incompressible Euler equations,
\begin{equation} \frac{\partial\boldsymbol{U}}{\partial t} + \boldsymbol{U \cdot \nabla}\boldsymbol{U} + \frac{1}{\rho}\boldsymbol{\nabla} P + g\boldsymbol{\hat{z}} = 0,\quad \boldsymbol{\nabla \cdot U} = 0, \end{equation}
where 
$\boldsymbol {U}$ is the fluid velocity, 
$P$ is the fluid pressure, 
$\boldsymbol {\hat {z}}$ is a unit vector in the vertical direction and the constants 
$\rho$ and 
$g$ are the fluid density and the acceleration due to gravity, respectively. The fluid is contained between a bottom boundary at 
$z=0$ and an upper free surface at 
$z = H$. The fluid must satisfy no penetration through the bottom boundary, giving 
$\boldsymbol {U\cdot \hat {z}} = 0$ at 
$z=0$. Two boundary conditions, a kinematic and a dynamic boundary condition, must be satisfied along the free surface itself. Here, we assume the fluid above the free surface to be dynamically passive, and in particular, to have a constant pressure 
$\bar {P}$. Together, these give the boundary conditions
\begin{equation} \frac{\partial H}{\partial t} + \boldsymbol{U \cdot \nabla} H = \boldsymbol{U\cdot \hat{z}} \quad \text{and}\quad P = \bar{P}\quad \text{on}\ z = H. \end{equation} We split the overall velocity and pressure into a steady purely swirling base flow and a small (magnitude 
$\varepsilon$) time-dependent perturbation, where
\begin{equation} \boldsymbol{U} = U_0(r)\boldsymbol{\hat{\theta}} + \varepsilon(u_r\boldsymbol{\hat{r}} + u_\theta\boldsymbol{\hat{\theta}} + u_z\boldsymbol{\hat{z}}), \quad P = P_0(r,z) + \varepsilon p,\quad H = h_0(r) + \varepsilon h. \end{equation}2.1. The steady base flow solution
A purely swirling steady base flow has a velocity given by 
$\boldsymbol {U_{0}} = U_{0}(r)\boldsymbol {\hat {\theta }}$. The governing equations and boundary conditions for this steady flow are satisfied provided we take
$$\begin{gather} P_0(r, z) = \bar{P} + \rho g ( h_0(r) - z), \end{gather}$$
$$\begin{gather}h_{0}(r) = h_\infty - \frac{1}{g}\int_r^\infty \frac{U^{2}_{0}(r')}{r'}\,\mathrm{d} r', \end{gather}$$
where 
$h_\infty$ is the depth of the fluid at 
$r=\infty$. This holds for any velocity profile 
$U_0(r)$. For the specific case of the Lamb–Oseen vortex considered here, we have
\begin{equation} U_{0}(r) = \frac{\varGamma_{0}}{2{\rm \pi} r}(1-\exp({-}r^2/a^2)), \end{equation}
where 
$a$ sets the radial size of the core and 
$\varGamma _{0}$ sets the circulation of the vortex. Note that, for the Lamb–Oseen vortex, for small 
$r$, we have 
$U_0(r) \approx r\varGamma _0/2{\rm \pi} a^2$, so that (2.5) is a solid body rotation near the centre of the vortex, while for large 
$r$ we have 
$U_0(r)\approx \varGamma _0/2{\rm \pi} r$, so that (2.5) is a potential swirl far from the centre of the vortex. A typical steady base flow free surface for the Lamb–Oseen vortex is shown in figure 2.
A typical section through a Lamb–Oseen vortex, showing the steady base flow free surface.
2.2. Perturbation dynamics
Waves arise when a small perturbation is introduced to the steady base solution. By linearizing the governing equations (2.1) about the base solution 
$(U_{0}, P_{0}, h_{0})$ given above, the governing equations for the perturbation are
$$\begin{gather} D_{t}u_{r} - 2\varOmega_{0}(r)u_\theta + \frac{1}{\rho}\frac{\partial p}{\partial r} = 0, \end{gather}$$
$$\begin{gather}D_{t}u_{\theta} + \frac{1}{r}(rU_0(r))'u_{r} + \frac{1}{\rho r}\frac{\partial p}{\partial \theta} = 0, \end{gather}$$
$$\begin{gather}D_{t}u_{z} + \frac{1}{\rho}\frac{\partial p}{\partial z} = 0, \end{gather}$$
$$\begin{gather}\frac{1}{r}\frac{\partial}{\partial r}(ru_{r}) + \frac{1}{r}\frac{\partial u_{\theta}}{\partial \theta} + \frac{\partial u_{z}}{\partial z} = 0, \end{gather}$$
where 
$D_{t} = \partial _{t} + \varOmega _{0}(r)\partial _{\theta }$ is the convective derivative and 
$\varOmega _{0}(r) = U_{0}(r)/r$ is the steady base flow angular velocity. (Here and throughout, primes denote derivatives with respect to the argument for a function of only one variable.) Linearizing the boundary conditions (2.2) about the steady base flow leads to
$$\begin{gather} u_{z} = \frac{1}{\rho g}D_{t}p + h'_{0}u_{r} \quad\text{and}\quad h = \frac{p}{\rho g} \quad\text{on}\ z = h_{0}(r), \end{gather}$$
$$\begin{gather}u_{z} = 0 \quad \text{on} \ z = 0. \end{gather}$$ Since our focus is on trapped waves and their formation, in addition to these boundary conditions there is another implicit condition, which is that there are no waves entering the domain from 
$r=\infty$, and thus only outgoing waves are allowed at 
$r=\infty$. This rather subtle condition will become more concrete when we truncate the domain to finite 
$r$ in order to numerically solve the above equations.
Before we progress further, we now make two changes to the governing equations for the perturbation. Firstly, since the governing equations for the perturbation are linear, we may assume a modal solution of the form
$$\begin{gather} \boldsymbol{u} = [u(r, z)\boldsymbol{\hat{r}} + v(r, z)\boldsymbol{\hat{\theta}} + w(r, z)\boldsymbol{\hat{z}}]\exp(-\mathrm{i}\omega t + \mathrm{i} m\theta), \end{gather}$$
$$\begin{gather}p = \phi(r, z)\exp(-\mathrm{i}\omega t + \mathrm{i} m\theta), \end{gather}$$
where 
$m$ is restricted to integer values, since the solution must be 
$2{\rm \pi}$ periodic, while 
$\omega$ in general will be complex, with 
${\rm Re}(\omega )/2{\rm \pi}$ being the oscillation frequency and 
$-{\rm Im}(\omega )$ being the decay rate. The eigenvalue problem that will eventually result will have 
$\omega$ as the eigenvalue to be found.
Secondly, we rewrite the governing equations and boundary conditions in a non-dimensional form. To do so, we choose the reference length scale to be 
$a$, the scale of the vortex core ((1.1) and figure 2). The reference time scale is that given by this length scale and gravity, 
$\sqrt {a/g}$. The velocity scale is thus 
$\sqrt {ag}$. All dimensional variable may then be expressed in terms of a non-dimensional variable (denoted by a tilde), as
\begin{equation} \left.\begin{gathered} (r,\theta,z) = (a\tilde{r}, \tilde{\theta}, a\tilde{z}), \quad t = \sqrt{\frac{a}{g}}\tilde{t}, \quad U_{0}(r) = \frac{\varGamma_{0}}{2{\rm \pi} a}\widetilde{U}_{0}(\tilde{r}), \quad \varOmega_0 = \frac{\varGamma_{0}}{2{\rm \pi} a^2}\tilde{\varOmega_0},\\ h = a\tilde{h}, \qquad p = \rho ag\tilde{p}, \quad \omega = \sqrt\frac{g}{a}\tilde{\omega}, \quad \boldsymbol{u} = \sqrt{ag}\boldsymbol{\tilde{u}}. \end{gathered}\right\} \end{equation}There are two physical parameters that are not scaled to unity by this non-dimensionalization: these may be thought of as the strength of the vortex and the depth of the fluid at infinity, given, respectively, as
\begin{equation} F =\frac{\varGamma_0}{2{\rm \pi}\sqrt{ga^3}}, \quad \tilde{h}_\infty = \frac{h_\infty}{a}. \end{equation}
Here 
$F$ is the Froude number and sets the non-dimensionalized velocity of the vortex. Hence, 
$F \to 0$ corresponds to a slow vortex with negligible steady surface height variation, while 
$F \to \infty$ corresponds to a fast vortex with significant steady surface height variation, as can be seen from (2.13c) below. The dimensionless depth 
$\tilde {h}_\infty$ is exactly that, so that the limit 
$\tilde {h}_\infty \to 0$ corresponds to the shallow-water limit and 
$\tilde {h}_\infty \to \infty$ corresponds to the deep-water limit. Care is needed, however, in considering the shallow-water limit: in what follows, we will assume that the steady fluid height never reaches zero, so that the bottom stays wetted, and consequently the shallow water limit 
$\tilde {h}_\infty \to 0$ must be taken together with the slow wide vortex limit 
$F\to 0$ such that 
$\tilde {h}_\infty /F^2$ is bounded away from zero, as can also be seen from (2.13) below. The azimuthal wavenumber 
$m$ is also a dimensionless parameter, and represents the rotational symmetry of the solution being investigated.
Dropping the tildes, the complete non-dimensional eigenvalue problem is
$$\begin{gather} ({-\mathrm{i}}\omega + \mathrm{i} m F \varOmega_{0}(r))u - 2F\varOmega_{0}(r)v + \frac{\partial \phi}{\partial r} = 0, \end{gather}$$
$$\begin{gather}({-\mathrm{i}}\omega + \mathrm{i} m F \varOmega_{0}(r))v + \frac{F}{r}(rU_0(r))'u + \frac{\mathrm{i} m}{r}\phi = 0, \end{gather}$$
$$\begin{gather}({-\mathrm{i}}\omega + \mathrm{i} m F \varOmega_{0}(r))w + \frac{\partial \phi}{\partial z} = 0, \end{gather}$$
$$\begin{gather}\frac{1}{r}\frac{\partial}{\partial r}(ru) + \frac{\mathrm{i} m}{r}v + \frac{\partial w}{\partial z} = 0, \end{gather}$$
together with the boundary conditions of no incoming modes at 
$r=\infty$, and
$$\begin{gather} w = 0, \quad \text{on} \ z = 0, \end{gather}$$
$$\begin{gather}w = ({-\mathrm{i}}\omega + \mathrm{i} m F \varOmega_{0}(r))\phi + F^{2} r\varOmega^{2}_{0}(r)u, \quad \text{on} \ z = h_{0}(r). \end{gather}$$For the Lamb–Oseen vortex, in dimensionless terms
$$\begin{gather} U_0(r) = \frac{1 - \exp({-}r^2)}{r}, \quad \varOmega_0(r) = \frac{U_0(r)}{r}, \end{gather}$$
$$\begin{gather}h_{0}(r) = h_{\infty} - F^2 \int_r^\infty \frac{U^{2}_{0}(r')}{r'}\,\mathrm{d} r'. \end{gather}$$ One consequence of the harmonic assumption (2.8) is the creation of a so-called ‘critical layer’ at a radial location 
$r=r_c$, given by 
$D_t = -\mathrm {i}\omega + \mathrm {i} mF\varOmega _0(r_c) = 0$. It is not immediately obvious from (2.11) that anything particularly special occurs at the critical layer, as in no equation does it cause the highest derivative to vanish, and therefore our numerics described in § 3 have no difficulty in this case. However, in fact the critical layer is related to behaviour that is not of the harmonic form assumed in (2.8); investigating the effect of the critical layer requires a different mathematical and numerical technique (such as Frobenius series; see, by way of example, King et al. (Reference King, Brambley, Liupekevicius, Radia, Lafourcade and Shah2022)), which is beyond the scope of this paper. However, we comment in passing that Fabre et al. (Reference Fabre, Sipp and Jacquin2006) found that disturbances related to the critical layer are necessarily damped, suggesting that the undamped trapped wave modes we are interested in here are not related to the critical layer.
3. Numerical methods
3.1. Absorbing layer for three-dimensional incompressible Euler equations
In practice, we solve the eigenvalue problem in a finite computational domain and hence introduce an artificial boundary at finite radius 
$R$. The far-field boundary condition of no incoming waves becomes a NRBC at this boundary. However, exact NRBCs are generally derived using the method of characteristics, but this is not available to us since (2.11d) is not hyperbolic. A widely used alternative is based on damping layers or even perfectly matched layers. The idea is to introduce a ‘damping layer’ – also called an ‘absorbing layer’ – surrounding the physical outer boundary. In the damping layer, waves are progressively damped until their amplitude is sufficiently small that reflections from the outer boundary do not re-enter the physical domain. Outgoing waves are effectively absorbed in the layer. The boundary condition at the outer boundary of the computational domain is then unimportant, and is typically taken to be either a Dirichlet-type or Neumann-type boundary condition.
To add damping to our governing equations, we add a ‘damped compressibility’ term into the continuity equation (2.11d), which then becomes
\begin{equation} \xi(r)\phi + \frac{1}{r}\frac{\partial}{\partial r}(ru) + \frac{\mathrm{i} m}{r}v + \frac{\partial w}{\partial z} = 0, \end{equation}
where 
$\xi (r)$ is a sufficiently smooth function that is identically zero outside the damping region. This is motivated by the introduction of damping in the acoustic wave equation following Gao et al. (Reference Gao, Song, Zhang and Yao2017, pp. 81–82), and is explained in further detail in Appendix B. In what follows, we find good results using the simple form for 
$\xi (r)$ given by
\begin{equation} \xi(r) =\begin{cases} 0, & r < R_{c}, \\ \bar{\xi}\left(\dfrac{r - R_{c}}{R - R_{c}}\right)^{2}, & R_{c} \le r \le R. \end{cases} \end{equation}
This gives a physical region, 
$r < R_{c}$, with no damping, and computational region 
$r\in [R_c,R]$ with a damping strength governed by the constant 
$\bar {\xi }$. We then impose a Dirichlet boundary condition on the pressure at an artificial boundary 
$r = R \gg 1$. We show in § 4.2 that this leads to the correct causal behaviour. A representation of the damping function as well as the subdivision of the numerical domain is shown in figure 3.
(a) Profile of the quadratic damping function defined in (3.2) with parameters: 
$\bar {\xi } = 5$, 
$R_c = 5$ and 
$R = 10$. (b) Schematic of the computational domain in the 
$r$–
$z$ plane.
The governing equations (2.11a–c) together with the modified continuity equation (3.1) and boundary conditions (2.12) form an eigenvalue problem to solve for the allowable frequencies 
$\omega$ permitting a non-zero modal solution.
3.2. Numerical discretization
In order to solve (2.11a–c, 2.12, 3.1), we used a Galerkin spectral method and expand with a combination of Legendre polynomials as basis functions in both the radial and axial coordinate in order to satisfy the Dirichlet boundary conditions. We also remap the domain from the domain 
$D = [0, R]\times [0, h_{0}(r)]$ to the computational square domain 
$S = [-1, 1]\times [-1, 1]$ to account for the shape of the computational domain with a variable surface height 
$h_0(r)$. We then obtain the weak formulation of the problem. Full details are given in Appendix A. The discretized problem ends up being of the form
\begin{equation} \boldsymbol{\mathsf{A}}\boldsymbol{w} = \omega \boldsymbol{\mathsf{B}}\boldsymbol{w}, \end{equation}
with 
$\boldsymbol {w} = (u_{ij}, v_{ij}, w_{ij}, \phi _{ij})$ representing our array containing the spectral coefficients of each unknown, 
$\boldsymbol{\mathsf{A}}$ and 
$\boldsymbol{\mathsf{B}}$ being matrices of order 
$4N_{x}N_{y} \times 4N_{x}N_{y}$, where 
$N_x$ and 
$N_y$ are the number of Legendre polynomials used in the radial and vertical directions, respectively. This discretized problem may then be solved using any numerical eigenvalue solver; here, we use the eig solver in MATLAB. We sort eigenvalues by imaginary part and focus on the leading eigenvalues, i.e. those with the most positive imaginary part. These correspond to unstable modes, should there be any, or the least damped stable modes, if not.
Not all solutions to the discretized problem (3.3) correspond to solutions to the continuous problem being approximated, however. To remove under-resolved eigenmodes and spurious eigenvalues, the numerical solutions to (3.3) are filtered, as described in the Appendix A.2. The numerical scheme has been validated against three published cases: two with walls at finite radius (Mougel et al. Reference Mougel, Fabre and Lacaze2014, Reference Mougel, Fabre and Lacaze2015), and one unbounded but in shallow water (Ford Reference Ford1994). We obtain unstable modes where they are known to exist. Details are given in Appendix C.
3.3. Spurious reflected modes
The finite numerical domain and damping region introduces another source of spurious modes besides those coming from the numerical discretization, namely those modes which are well-resolved but which include a significant reflection from either the damping boundary at 
$r=R_c$ or the truncation boundary at 
$r=R$. Such spurious eigenmodes are affected by the values of 
$R_c$ and 
$R$, and by the amount of damping 
$\bar {\xi }$, whereas good approximations to the modes on the infinite domain should be insensitive to these values. Since 
$-\textrm {Im}(\omega )$ is the decay rate of the mode, variations in damping typically have a strong effect on 
$\textrm {Im}(\omega )$ for spurious reflected modes. We may therefore remove these spurious reflected modes by running our numerical code twice: the first time with a suitably chosen amount of damping 
$\bar {\xi }$, and the second time with twice that amount of damping 
$2\bar {\xi }$. Only those modes whose eigenvalues do not change significantly with the change in the damping coefficient are retained (as measured using the same metric (A11) used for the numerical resolution).
3.4. Numerical convergence study
The first convergence study involves varying the amount of damping in the damping layer and checking that the eigenvalues do not change, nor does the shape of the corresponding eigenfunctions in the physical domain 
$r < R_{c}$. Here we present the convergence results for 
$m = 7$ and 
$F = 0.3$, since for these parameters the system supports radiating (outwardly propagating) modes, and such modes provide the most stringent test of a NRBC. For varying magnitudes of damping 
$\bar {\xi }$, the radiating eigenfunction is plotted in figure 4, and the corresponding eigenvalues are given in table 1. Two domain sizes are shown in figure 4: 
$R_c=5$ and 
$R=10$ are the values used for the results in section § 4; and 
$R_c=15$ and 
$R=30$ give an extended domain so that the effects of damping and resonance can be seen more clearly. The damping clearly influences the eigenfunction shape, and for 
$\bar {\xi }\in \{0.1,1\}$ a clear standing wave shape is seen in figure 4. For 
$\bar {\xi }\in \{5,10\}$, the eigenfunctions are practically identical in the physical domain 
$0 < r < R_c$, and only differ in the damping layer 
$R_c\leq r \leq R$. For 
$\bar {\xi }=50$, however, the damping is too strong, and little oscillations can be seen for 
$r< R_c$, suggesting wave reflection by the edge of the damping layer. This is also supported by the eigenvalues in table 1, which show the sensitivity of 
$\textrm {Im}(\omega )$ to variations in damping strength, as expected.
Eigenfunctions for 
$m = 7$, 
$F = 0.3$ and different values of 
$\bar {\xi }$. A vertical red dashed line indicates the radial point 
$R_c$ where the damping effect begins. (a) A small computational domain 
$R_{c} = 5$, 
$R = 10$. (b) A large computational domain 
$R_c=15$, 
$R=30$.
Eigenvalues as function of the amount of damping for 
$m=7$ and 
$F = 0.3$, 
$R_{c} = 5$ and 
$R = 10$.
The second convergence study involves varying the width of the damping layer whilst maintaining the same size of the computational domain and a fixed damping strength 
$\bar {\xi } = 5$. The eigenfunctions are displayed in figure 5, and the corresponding eigenvalues in table 2. For most results in figure 5 the eigenfunctions can be seen to be very similar in the physical domain 
$0< r< R_c$ and to decay smoothly in 
$r$ in the damping layer, while for a damping layer of width 
$R-R_c=1$ a standing wave pattern can be seen for both domain sizes, implying significant wave reflection from the domain boundary at 
$r=R$. Again, this is also seen for the variations in the eigenvalue in table 2, with again 
$\textrm {Im}(\omega )$ being particularly sensitive.
Eigenfunctions for 
$m = 7$, 
$F = 0.3$, 
$\bar {\xi } = 5$, and different values of 
$R_{c}$, computed using a small domain 
$R=10$ (a) and a larger domain 
$R=30$ (b).
Eigenvalues as function of the initial position of the absorbing layer for 
$m=7$ and 
$F = 0.3$, 
$\bar {\xi } = 5$ and 
$R = 10$.
Finally, in figure 6 we compare the radiating mode for 
$m = 7$, 
$F = 0.3$ computed using two different sizes of the domain: 
$R_c = 5$, 
$R = 10$ on one hand and 
$R_c = 15$, 
$R = 30$ on the other. The two modes match very well over the physical interval 
$r \in [0, 5]$. This good match demonstrates that the damping layer successfully allows our numerical simulation on a small interval to reproduce results that would have been obtained using a larger interval, thus allowing our numerics to emulate an infinite unbounded domain using a finite computational domain.
Eigenfunction for 
$m = 7$, 
$F = 0.3$, 
$\bar {\xi } = 10$ computed using two different domain sizes: 
$R_c = 5$ and 
$R = 10$ for the blue solid curve; and 
$R_{c} = 15$, 
$R = 30$ for the red dashed curve.
3.5. Choice of numerical parameters
Based on numerical convergence studies, for the results that follow we take 
$R = 10$ and 
$R_{c} = 5$, with 
$\bar {\xi } = 5$ (to which results are compared with 
$\bar {\xi }=10$). This choice is motivated by the need for sufficiently high resolution to resolve all modes of interest in the range of Froude numbers and azimuthal wavenumbers considered. The eigenvalue tolerance and eigenfunction resolvedness tolerance are taken to be 
$\mathrm {tol} = 10^{-2}$ and 
$10^{-1}$, respectively, whereas other numerical parameters are taken to be 
$N_{x} = 50$, 
$N_{y} = 20$, 
$b_{x} = 12$ and 
$b_{y} = 4$, as explained in more detail in Appendix A. For the majority of results presented here, we use a fluid depth of 
$h_\infty = 5$. This is the depth shown in figure 2. This value is large enough to allow for a wide range of Froude numbers (
$0 \leq F < \sqrt {h_{\infty }/\log (2)} \simeq 2.68$) without forming a dry region near 
$r=0$, while small enough to exhibit finite-depth effects.
4. Results
We compute leading (i.e. least damped) eigenmodes for a range of Froude numbers and for a range of azimuthal wavenumbers as large as 
$m=20$, and find that above 
$m = 6$ surface gravity waves dominate. We shall first present the surface waves for a representative case, wavenumber 
$m=7$, and discuss in detail how the modes and eigenvalues depend on Froude number for this case. Following this, we consider the dependence of modes and eigenvalues on the azimuthal wavenumber.
4.1. Representative case 
$m=7$
We begin with a detailed description of the case azimuthal wavenumber 
$m=7$. The reason for taking such value of 
$m$ is dictated both by the explanatory picture in Patrick et al. (Reference Patrick, Coutant and Richartz2018, p. 6), who show a picture of the spiral structure of a normal mode solution for sufficiently high azimuthal wavenumber, as well as by the reasonable amount of computational time needed to get the leading surface waves eigenmodes. Indeed, we find numerically that the spatial structure of surface gravity waves becomes thinner and more localized close to the free surface as 
$m$ increases. Following particularly the latest argument, the case of 
$m=7$ has been taken as a reference study case and in the following we are going to show most of the interesting features applied to this case.
One of the key results of our study is the change in character of modes as functions of the rotation rate of the vortex. Specifically, we characterize two extremes of modes: radiating modes and trapped modes. The spatial structure of the radiating modes has a characteristic spiral shape, as shown in figure 7(a), and a long tail in the radial direction, as shown in figure 8(a). Such modes behave as travelling waves along the radial direction, and radiation to infinity is responsible for the dissipation of the initial energy put into the system to excite the mode. (An animation of the radiating mode in figures 7(a) and 8(a) is given as Movie 1 in the online Supplementary material available at https://doi.org/10.1017/jfm.2024.645.) In contrast, trapped modes are confined to a bounded region in the radial direction, and without the spiralling structure seen for the radiating modes. Those modes behave as standing waves along the radial direction and dissipate extremely slowly in comparison. Examples can be seen in figures 7(e) and 8(e). (An animation of the the trapped mode in figures 7(e) and 8(e) is given as Movie 2 in the online Supplementary material.) The distinction between radiating and trapped modes can further be seen in the eigenvalues. Indeed, radiating modes have eigenvalues with a significant negative imaginary part, indicating exponential decay in time, whereas trapped modes have a negligible imaginary part and so do not appreciably decay in time.
Plots of the free surface height 
$h(r, \theta, t=0) = \textrm {Re}[\phi (r, h_{0}(r))\exp \{\mathrm {i} m\theta \}]$ for 
$m = 7$: (a) 
$F = 0.3$; (b) 
$F = 0.32$; (c) 
$F = 0.33$; (d) 
$F = 0.34$; (e) 
$F = 0.36$. All modes displayed rotate clockwise, i.e. against the vortex flow. Animations of (a) and (e) are given in Movies 1 and 2 in the online Supplementary material.
Plots of the pressure distribution 
$|\phi _{m}(r, z)|$ for 
$m = 7$: (a) 
$F = 0.3$; (b) 
$F = 0.32$; (c) 
$F = 0.33$; (d) 
$F = 0.34$; (e) 
$F = 0.36$. Animations of (a) and (e) are given in Movies 1 and 2 in the online Supplementary material.
Figures 7 and 8 show representative examples of the continuous, but rapid, transition from radiating modes to trapped modes as the Froude number (the dimensionless rotation rate of the vortex), increases. This transition between radiating and trapped modes is also seen in figure 9 where the modulus of the free surface height is plotted as a function of 
$r$ for different Froude numbers.
Modulus of the free surface height along the radius for 
$m = 7$ in the transition regime 
$F \in [0.3, 0.36]$ from radiating to trapped modes. The trapped mode at 
$F = 0.5$ is also shown for comparison.
To look more closely at the trend in the eigensolutions, let us first denote by 
$n$ an integer representing the number of peaks in radial direction of the modulus of the pressure eigenfunctions, as displayed in figure 10. In this way each eigensolution will be indexed by both 
$m$ and 
$n$, with corresponding eigenvalues 
$\omega = \omega _{mn}$. In figure 11(a), the two types of modes can be seen: radiating modes with a significant negative 
$\textrm {Im}(\omega )$, and trapped modes having an almost null 
$\textrm {Im}(\omega )$. There is no sharp transition between radiating and trapped modes. Accordingly, here we set an arbitrary threshold to separate the two sets of modes by considering a mode to be trapped when its eigenvalue has an imaginary part smaller than 
$10^{-5}$ in modulus (although our results are relatively insensitive to this threshold; see § 4.3). Hence, while the trapped modes considered here are almost neutrally stable, they have eigenvalues with small negative imaginary part and so radiate very slightly. Interestingly, none of the modes computed here are observed to become linearly unstable, and we hypothesize that the base solution is at most marginally stable, but not unstable, to linear perturbations of this type at large Froude number. This stability can be viewed as a manifestation of the robustness of the Lamb–Oseen (Fabre et al. Reference Fabre, Sipp and Jacquin2006), even in the presence of a free surface.
Modulus of the pressure 
$|\phi _{m}(r, z)|$ for 
$m = 7$ and 
$F = 0.5$. We show the four least-stable eigenmodes for such parameters, having one, two, three and four peaks, respectively.
Imaginary (a) and real (b) parts of the eigenvalue 
$\omega _{mn}$ as function of the Froude number 
$F$ for 
$m = 7$ and varying 
$n$.
Another interesting phenomenon concerns the propagation direction of the surface waves with respect to the rotation of the base vortex flow as its rotation rate is varied. This is shown in figure 11(b), which shows how the real part of the eigenvalues varies with the Froude number, again for 
$m = 7$. For 
$\textrm {Re}(\omega ) < 0$, waves rotate opposite to the base flow (counter-rotating waves), while for 
$\textrm {Re}(\omega ) > 0$, waves rotate in the same direction as the base flow (corotating waves). It is clear that for each 
$n$, there is a value of Froude number separating counter-rotating from corotating surface waves.
4.2. Far-field behaviour and causality
In this section, we investigate the behaviour of our numerical solutions in the far field, and show that they have the correct behaviour expected of causal outgoing waves. Away from the vortex core, in the far field, the behaviour of surface waves is well understood. We first consider the behaviour of deep-water waves in the far-field (see, e.g. Batchelor Reference Batchelor2000, p. 371)
\begin{equation} h \propto r^{{-}1/2}\exp\{{\pm}\mathrm{i} kr\} \quad \text{with the dispersion relation}\ \omega^2=k. \end{equation}
For a complex frequency 
$\omega = \omega _r + \mathrm {i}\omega _i$, the dispersion relation gives 
$\textrm {Im}(k)=2\omega _r\omega _i$. Causality (Briggs Reference Briggs1964; Bers Reference Bers1983; Berti et al. Reference Berti, Cardoso and Starinets2009, pp. 17–18; Le Dizès & Billant Reference Le Dizès and Billant2009; Riedinger et al. Reference Riedinger, Le Dizès and Meunier2010) implies that the outgoing wave must decay when 
$\omega _i>0$, so we must take 
$\pm = \mathrm {sgn}(\omega _r)$. Thus, for temporally damped modes when 
$\omega _i<0$ we expect spatially growing behaviour in the far field. This growth in the far field is found in our numerical solutions, as shown in figure 12(a). This behaviour is also seen as 
$\omega _i \to 0$, but only becomes apparent at extremely large radii (e.g. 
$r= 10^{10}$ gives 
$|h(r)| = 0.01$ in figure 12b), hence our assertion that eigenmodes with sufficiently small 
$\omega _i$ can be deemed trapped.
Comparison of asymptotic and numerical solutions with parameters 
$m = 7$ and 
$h_{\infty } = 5$, with 
$R_c = 40$ and 
$\bar {\xi } = 10$. (a) A damped radiating mode with 
$F = 0.3$, giving 
$\omega =-1.38 - 0.013\mathrm {i}$. (b) A trapped mode with 
$F = 0.5$, giving 
$\omega = -0.54 - 6\times 10^{-10}\mathrm {i}$. Solid black lines denote the numerical eigenmode, and blue dashed lines denote the far-field asymptotic solution obtained from the deep-water dispersion relation (4.1). The vertical red dashed lines indicate the edge of the numerical damping region.
The growth of eigenmodes in the far field explains why our numerical method breaks down for radiating modes with large temporal decay rate observed at small 
$F$. When the temporal decay rate 
$-\omega _i$ is sufficiently large, the spatial growth of the eigenmode overcomes the artificial numerical damping in the damping region, and our numerical NRBC breaks down. This results in 
$\omega$ becoming dependent on the amount of damping in the damping region, and the numerical eigenvalue is automatically removed as being unresolved as previously described in § 3.3. A similar result has been obtained in Oliveira, Cardoso & Crispino (Reference Oliveira, Cardoso and Crispino2014) studying a model wave equation describing the vortex–waves interaction in the shallow-water limit. We leave to future studies the possibility of tracking the trend of the spectrum curve as the Froude number goes to zero. It should be noted, that modes having a large negative decay rate in time do not play a relevant role in the dynamics of the system.
4.3. Extension to other azimuthal wavenumbers.
The results described above for 
$m = 7$ are found to be typical at larger azimuthal wavenumbers. Figure 13 shows the trend of the imaginary part of the eigenvalues as a function of the Froude number. The results are qualitatively similar to the 
$m=7$ case. Eigenvalue branches shift to lower 
$F$ with increasing 
$m$, thus shifting to lower 
$F$ the value of 
$F$ separating radiating and trapped modes. This suggests that for very high azimuthal wavenumber perturbations we expect to get radiating eigenmodes at lower and lower Froude numbers, as shown in figure 14 for the example cases of 
$m = 15$ and 
$m = 20$.
Imaginary part of the eigenvalues as functions of 
$F$ for 
$m = 2$, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, with the arrow indicating the direction of increasing 
$m$.
Shape of the free surface height 
$h(r, \theta, t=0) = \textrm {Re}[\phi (r, h_{0}(r))\exp \{\mathrm {i} m\theta \}]$ in two high azimuthal wavenumber cases: (a) 
$m = 15$, 
$F = 0.19$; (b) 
$m = 20$, 
$F = 0.16$.
We have computed the range of Froude numbers in which modes radiate and transition towards a neutrally stable state for different values of 
$m$. This is summarized in figure 15 showing the classes of solutions previously described over a wide range of azimuthal orders 
$m$ and Froude numbers 
$F$. The contour lines plotted separate the regions of parameter space where waves are radiating and where waves are trapped, with the arrow indicating the direction of better trapping. While the individual contours range from 
$\textrm {Im}(\omega )=10^{-2}$ to 
$\textrm {Im}(\omega ) = 10^{-5}$, their close spacing shows that the exact threshold value of 
$\textrm {Im}(\omega )$ separating trapped and radiating behaviour is not that important, with all contours giving a similar boundary. Also plotted in figure 15 is a dashed line showing 
$\textrm {Re}(\omega )=0$, which is the boundary between counter-rotating and corotating trapped waves. Two notable features of the plot as 
$m$ is decreased are the shift of contours to larger 
$F$ and the widening of the separation between contours. Both of these features are consistent with figure 13, where one sees not only a shift in the eigenvalue curves with 
$m$, but also a steepening of the transition between radiating and trapped modes with increasing 
$m$.
Regimes of radiating (large 
$|\textrm {Im}(\omega )|$) and trapped (small 
$|\textrm {Im}(\omega )|$) modes for different azimuthal wavenumbers 
$m$ and Froude numbers 
$F$. Solid lines are contours of 
$|\textrm {Im}(\omega )|$, with the arrow indicating the direction of more perfectly trapped behaviour. The dash–dot line separates counter-rotating modes (
$T^{-}$) with 
$\textrm {Re}(\omega )<0$ from corotating modes (
$T^{+}$) with 
$\textrm {Re}(\omega )>0$.
While this section considers a wide range of values of 
$m \geq 3$, we find 
$m=0, 1$ and 
$2$ to be dominated by inertial waves rather than surface waves, which are discussed further in § 4.5.
4.4. Effect of the free-surface height on the eigenmodes
In this subsection we investigated the effects associated with a variation in the fluid depth, exploring the regime between shallow and deep water. We again focus on modes with azimuthal mode number 
$m = 7$ and show how the structure of eigenmodes and corresponding eigenvalues vary with fluid depth, and hence how the transition between radiating and trapped modes changes with depth.
The variation in the structure of the surface gravity eigenmodes with fluid depth 
$h_{\infty }$ is displayed in figure 16 for a typical trapped mode computed at Froude number 
$F = 0.5$. The same type of analysis is carried out for a radiating mode at 
$F = 0.3$, whose spatial structures are shown in figure 17. Finally, in figure 18, the transition between radiating and trapped is shown as the fluid depth is varied between 
$h_{\infty } = 0.25$, 0.5, 1. Figure 19 shows the variation in eigenvalue as the Froude number is varied, for different fluid depths. The trend in 
$\textrm {Im}(\omega )$ shown suggests that more shallow systems require a faster rotating vortex in order to get the same wave trapping as a deeper water configuration. The variation with depth of 
$\textrm {Re}(\omega )$ is less significant. For 
$h_{\infty } = 1$, the eigenvalues are essentially identical to those computed for even higher depths, for example, in this case 
$h_{\infty } = 5$. Hence, in this regard, for a trapped mode a fluid with 
$h_{\infty } \in [0.5, 1]$ can already be considered deep water.
Effect of the fluid height variation on a single peak (
$n = 1$) trapped pressure mode computed for 
$m = 7$, and 
$F = 0.5$: (a) 
$\omega = -0.3058$, 
$h_{\infty } = 0.25$; (b) 
$\omega = -0.5316$, 
$h_{\infty } = 0.5$; (c) 
$\omega = -0.5366$, 
$h_{\infty } = 1$; (d) 
$\omega = -0.5370$, 
$h_{\infty } = 5$. The imaginary part of the eigenvalues reported here is at least of order 
$10^{-9}$, hence can be considered as purely real.
Effect of the fluid height variation on a radiating pressure mode computed for 
$m = 7$, and 
$F = 0.3$: (a) 
$\omega = -1.3371 - 0.0530\mathrm {i}$, 
$h_{\infty } = 0.4$; (b) 
$\omega = -1.3633 - 0.0320\mathrm {i}$, 
$h_{\infty } = 0.5$; (c) 
$\omega = -1.3772 - 0.0140\mathrm {i}$, 
$h_{\infty } = 1$; (d) 
$\omega = -1.3770 - 0.0138\mathrm {i}$, 
$h_{\infty } = 5$.
Free surface height for 
$m = 7$, 
$F = 0.4$ and different fluid depths: (a) 
$h_{\infty } = 0.25$; (b) 
$h_{\infty } = 0.5$; (c) 
$h_{\infty } = 1$.
Effect of the fluid height variation on the spectrum curves for the least damped set of eigenmodes obtained for 
$m = 7$.
4.5. Inertial modes
While the focus of this study is surface gravity waves, our numerics finds modal solutions to the governing equations (2.11) without any assumption about the modes being surface gravity waves. By way of contrast, therefore, in this section we briefly discuss another type of mode, namely inertial modes. These modes persist even at low rotation rates, and are characterized by being neutrally stable as well as being concentrated within the vortex core and not being localized close to the free surface. Our results are therefore similar to the numerical results presented by Mougel et al. (Reference Mougel, Fabre and Lacaze2015) for a solid-body rotation, since the Lamb–Oseen vortex is very close to solid-body rotation close to the centre. Figure 20 shows two inertial modes computed for 
$m = 2$, 
$h_{\infty } = 1$, and 
$F = 0.2, 0.8$. The corresponding eigenvalues are purely real, and are also plotted as function of the Froude number in figure 20.
Inertial modes for 
$m = 2$, 
$h_{\infty } = 1$ and 
$F = 0.2$ (a); 
$F = 0.8$ (b). Eigenvalues curve as a function of the Froude number for the same parameters (c).
5. Conclusions and future developments
We have considered the linear response to small perturbations of a free surface Lamb–Oseen vortex flow in a laterally unbounded domain. Throughout this work, we have assumed a harmonic dependence 
$\exp \{-\mathrm {i}\omega t + \mathrm {i} m\theta \}$ (2.8) and looked for modal solutions of this form. We restricted our focus to surface waves, as it was these waves that were seen in the motivating experiments; we note that the numerical procedure described here also calculates other types of modes which we have not investigated in depth here, and which may in some circumstances be the least damped modes, especially at low Froude numbers.
Our study classified surface waves into two distinct classes: radiating modes and trapped modes. While these have previously been seen for shallow water vortices with inflow (such as a bathtub vortex), here we have shown that they arise in a fully three-dimensional non-shallow-water problem without the need for an inflow to help trap the modes. Radiating modes are dissipative, and have a spatial structure which extends in the horizontal plane, therefore behaving as radially travelling waves. By contrast, trapped modes, by virtue of their near zero growth rate, persist for long times with little dissipation and remain localized within the core region of the vortex and hence could be expected to be observed on the surface of a rotating vortex after all damped modes have dissipated. The spatiotemporal behaviour of these modes resembles that of a radially standing wave instead of a radially travelling wave. Our numerical predictions show the appearance of trapped modes provided the Froude number 
$F$ (the dimensionless rotation rate) is above a threshold value, with the threshold depending on the azimuthal mode number 
$m$ (as summarized in figure 15). Our computations also suggest that trapped modes will asymptotically approach a neutrally stable state in the limit of large Froude number without becoming unstable, although the exact asymptotic scaling of 
$\textrm {Im}(\omega )$ for large 
$F$ was not obvious from the numerical analysis carried out here. Interestingly, our numerics finds no unstable wave modes with 
$\textrm {Im}(\omega )>0$, unlike in either bounded vortex flows (Tophøj et al. Reference Tophøj, Mougel, Bohr and Fabre2013; Mougel et al. Reference Mougel, Fabre and Lacaze2014, Reference Mougel, Fabre, Lacaze and Bohr2017) or unbounded shallow water vortex flows (Ford Reference Ford1994), both of which possess unstable modes with 
$\textrm {Im}(\omega )>0$. (As shown in Appendix C, our numerics reproduces such instabilities, and thus would be expected to find instabilities in the present case if they existed.) Finally, we have confirmed what has been seen in the initial motivating experiment of the pool, namely that stable trapped surface wave modes exist and can propagate with or against the base swirling flow (as also summarized in figure 15).
Our assumed harmonic dependence 
$\exp \{-\mathrm {i}\omega t + \mathrm {i} m\theta \}$ (2.8) means modes with 
$\textrm {Im}(\omega )>0$ are unambiguously unstable and grow exponentially in time, while modes with 
$\textrm {Im}(\omega )<0$ are unambiguously stable and decay exponentially in time. In the edge case with 
$\textrm {Im}(\omega )=0$, a critical layer may be present, defined as a radial location 
$r=r_c$ for which 
$D_t = -\mathrm {i}\omega + \mathrm {i} mF\varOmega _0(r_c) = 0$. This possibility has been neglected here. The critical layer is related to the continuous spectrum of the differential operator, and analysing the critical layer requires a different analysis and numerical technique: the temporal Fourier transform implicit in the harmonic dependence must be inverted, and the singular point 
$r=r_c$ must be avoided by complex deformation of the solution into the complex 
$r$-plane, with the correct choice of deformation to maintain causality (such as has been done for aeroacoustic waves; e.g. King et al. (Reference King, Brambley, Liupekevicius, Radia, Lafourcade and Shah2022)). Despite the critical layer being present only for 
$\textrm {Im}(\omega )=0$, it may result in a solution with algebraic growth in time. Consequently, the lack of solutions with 
$\textrm {Im}(\omega )>0$ does not prove the stability of perturbations to a Lamb–Oseen vortex with a free surface, but only the lack of unstable solutions with exponential growth in time. However, we note that Fabre et al. (Reference Fabre, Sipp and Jacquin2006) found that the Lamb–Oseen vortex without a free surface is stable to axial-invariant perturbations as disturbances related to the critical layer are necessarily damped in that case.
Based on the results presented here, we may speculate on the underlying physical mechanism causing this trapping behaviour. One likely candidate might have been the non-uniform shape of the free surface, with trapped modes appearing to be localized at a radius close to the largest gradient of the undeformed free surface. However, in Appendix D results are shown for numerical simulations with an (artificial) flat free surface, and almost identical results are obtained, suggesting that, at least for the parameters investigated here, the free surface shape is less important than the base flow. We hypothesize, therefore, that it is the refraction of the surface gravity waves by the base flow, causing the waves to curve around the vortex with exactly the curvature needed to remain trapped, that is the underlying physical mechanism. This would explain why the non-shallow-water assumption is important here, as the dispersive nature of the waves in this case means a frequency can be chosen to obtain the particular wave group velocities needed for wave trapping, whereas in the shallow water case the waves are non-dispersive and have a prescribed group velocity independent of their frequency.
To numerically simulate a (horizontally) unbounded fluid on a bounded numerical domain, a far-field NRBC or buffer region is needed. Here, a novel additional term is introduced into the governing equations, to provide damping of the surface waves in the buffer region only. This method has proved more accurate than any NRBC we implemented while remaining computationally viable. Indeed, introduction of additional unknowns in the mathematical formulation was not needed, thus overcoming the main drawback of ‘perfectly matched layer’ methods. Furthermore, given that the background vortex flow vanishes in the far-field, the same absorbing layer formulation can be employed with other similar vortex distributions.
The numerical expense of our eigenvalue problem (which is a two-dimensional spatial problem involving both 
$r$ and 
$z$ coordinates) might be reduced by investigating a one-dimensional approximation along the radial coordinate only; for example in either the deep- or shallow-water limits. Moreover, new methods and formulations for the imposition of a NRBC in the far-field would certainly lead to a saving in the computational time as it will avoid, for example, the need to numerically resolve the unphysical buffer region.
Due to the computational expense, we have left to future studies the extension of our parametric study of waves to extreme Froude numbers. In particular it would be interesting to investigate the trend of modes for very small and very large Froude numbers. Also, consideration of lower or higher aspect ratios 
$h_{\infty }$ could be investigated. The assumption of modal solutions with a harmonic dependence 
$\exp \{-\mathrm {i}\omega t + \mathrm {i} m\theta \}$ (2.8) could also be relaxed, for example through time-domain simulation, which would allow investigation into the critical layer (e.g. as in Riedinger et al. (Reference Riedinger, Le Dizès and Meunier2010) and King et al. (Reference King, Brambley, Liupekevicius, Radia, Lafourcade and Shah2022)). Even with the assumption of harmonic dependence 
$\exp \{-\mathrm {i}\omega t + \mathrm {i} m\theta \}$, a different but related problem would be to investigate the scattering of an incoming wave encountering the vortex, which would require a different far-field boundary condition to introduce the incoming wave as well as to allow outgoing waves to propagate through the far-field boundary without reflection. Finally, our model also neglects both nonlinear and surface tension effects; while this is justified for the swimming pool application we model here, these assumptions break down for either large amplitudes or short wavelengths. Therefore, it would be interesting to investigate whether their inclusion could lead to an instability of the base vortex flow.
Supplementary material
Supplementary movies are available at https://doi.org/10.1017/jfm.2024.645.
Acknowledgements
The authors are grateful to J.M. Skipp for the initial motivation and videos motivating this project. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.
Funding
E.Z. was funded through the Warwick Mathematics Institute Centre for Doctoral Training, and gratefully acknowledges the support of the University of Warwick and the UK Engineering and Physical Sciences Research Council (EPSRC grant EP/W523793/1). E.J.B. gratefully acknowledges the support of the UK Engineering and Physical Sciences Research Council (EPSRC grant EP/V002929/1). D.B. gratefully acknowledges supported from the Simons Foundation (grant no. 662985).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Detailed derivation and analysis of the numerical solution
A.1. Numerical solution of the eigenvalue problem
In order to solve (2.11, 2.12, 3.1), we used a Galerkin spectral method and expand with a combination of Legendre polynomials as basis functions in both the radial and axial coordinate in order to satisfy the Dirichlet boundary conditions. First of all, to account for the shape of the computational domain with a variable surface height 
$h_0(r)$, we remap the domain from the physical domain 
$D = [0, R]\times [0, h_{0}(r)]$ to the computational square domain 
$S = [-1, 1]\times [-1, 1]$ using the non-orthogonal transformation
\begin{equation} x = \frac{2}{R}r - 1, \quad y = 2\frac{z}{h_{0}(r)} - 1. \end{equation}
In this new domain 
$S$, the problem becomes
$$\begin{gather} \mathrm{i} m F \varOmega_{0}(x)u - 2 F \varOmega_{0}(x) v + \frac{2}{R}\left[\frac{\partial\phi}{\partial x} - \frac{h'_{0}(x)}{h_0(x)}(y + 1)\frac{\partial\phi}{\partial y}\right] = \mathrm{i}\omega u, \end{gather}$$
$$\begin{gather}\mathrm{i} m F \varOmega_{0}(x)v + \frac{2F}{R(x+1)}\frac{{\rm d}}{{\rm d}\kern0.07em x}[ (x+1)U_{0}(x) ] u + \frac{2\mathrm{i} m}{R(x + 1)}\phi = \mathrm{i}\omega v, \end{gather}$$
$$\begin{gather}\mathrm{i} m F \varOmega_{0}(x)w + \frac{2}{h_0(x)}\frac{\partial \phi}{\partial y} = \mathrm{i}\omega w, \end{gather}$$
\begin{gather} \frac{2}{R(x+1)}\frac{\partial}{\partial x}[(x+1)u] + \frac{2}{R}\frac{h'_{0}(x)}{h_{0}(x)}u + \frac{2\mathrm{i} m}{R(x + 1)}v \nonumber\\ +\,\frac{2}{h_0(x)}\frac{\partial}{\partial y}\!\left[ w - \frac{h'_{0}(x)}{R}(y+1)u \right] + \xi(x)\phi = 0, \end{gather}
$$\begin{gather}\mathrm{i} m F \varOmega_{0}(x)\phi(x, 1) + \frac{2h_0'(x)}{R}u(x, 1) - w(x,1) = \mathrm{i}\omega\phi(x,1), \end{gather}$$
$$\begin{gather}w(x, -1) = 0, \quad \phi(1,y) = 0, \quad |\boldsymbol{u}({-}1, y)| < \infty, \quad |\phi({-}1, y)| < \infty. \end{gather}$$
For 
$R_c \leq r \leq R$, we include a damping layer with a non-zero 
$\xi (r)$. Here, we choose a simple form for 
$\xi (r)$, which we demonstrate in § 3.4 works well for the parameters we consider here:
\begin{equation} \xi(r) =\begin{cases} 0, & r < R_{c}, \\ \bar{\xi}\left(\dfrac{r - R_{c}}{R - R_{c}}\right)^{2}, & R_{c} \le r \le R. \end{cases} \end{equation}
Since 
$\xi (r)$ is continuously differentiable at 
$r=R_{c}$, no special transmission conditions are needed there (Sim Reference Sim2010, p. 46).
We next obtain the weak discretized formulation of this problem. The basis functions are chosen such that they identically satisfy the homogeneous Dirichlet boundary data. Also, the specific expansion of each variable in terms of basis functions will depend on the value of the azimuthal wavenumber, as this is associated with a singularity at the origin for the pressure and the azimuthal velocity component. In the following we make the distinction between the case of axisymmetric perturbations (
$m = 0$), and generic non-axisymmetric perturbations (
$m \neq 0$).
A.1.1. Discrete problem for non-axisymmetric perturbations
In the general case of a non-zero value of 
$m$, the weak formulation of the eigenvalue problem can be obtained in the following way. First of all, we notice the presence of singular terms at 
$r=0$ for 
$u, v, \phi$ so we will require these functions to be null at the origin. Consequently, even 
$w$ will be so. Thus, we look at the unknowns in the following spaces:
$$\begin{gather} [u, v, \phi] \in V_{H}(S) = \{ (u, v, \phi) \in H^{1}(S): (u, v, \phi) = 0, \ \text{at} \ r = 0 \}, \end{gather}$$
$$\begin{gather}w \in V_{v}(S) = \{ w \in H^{1}(S): w = 0, \ \text{at} \ r = 0, \ z = 0 \}, \end{gather}$$
where 
$H^{1}(S)$ is the usual Sobolev space (Quarteroni Reference Quarteroni2009, ch. 2). We multiply the azimuthal component of the momentum equation and the continuity equation by 
$(x+1)$. Then, by multiplying each of the equations by suitable test functions 
$(v_{x}, v_{t}, v_{y}, q)$ in the same space as the corresponding unknowns, after integrating over the square 
$S$ and exploiting the boundary conditions, the weak formulation of the differential problem reads – find 
$([u, v, \phi ], w) \in V_{H}(S) \times V_{v}(S)$ such that 
$\forall ([v_{x}, v_{t}, q], v_{y}) \in V_{H}(S) \times V_{v}(S)$ the following holds:
$$\begin{gather}\int_{S} v_{x}\textrm{i} m F \varOmega_{0}(x)u - \int_{S} v_{x} 2 F \varOmega_{0}(x) v - \int_{S} \frac{2}{R}\frac{\partial v_{x}}{\partial x}\phi - \int_{S} v_{x} \frac{2}{R}\frac{h'_{0}(x)}{h_0(x)}(y + 1)\frac{\partial\phi}{\partial y}\nonumber\\ = \textrm{i}\omega \int_{S} v_{x} u, \end{gather}$$
$$\begin{gather}\int_{S} v_{t}\left[\textrm{i} m F \varOmega_{0}(x)(x+1)v + \frac{2F}{R}\!\left[(x+1)U_{0}(x)\right]\!u + \frac{2\textrm{i} m}{R}\phi\right] = \textrm{i}\omega\int_{S} v_{t} (x+1) v, \end{gather}$$
$$\begin{gather}\int_{S} v_{y}\left[\textrm{i} m F \varOmega_{0}(x)w + \frac{2}{h_0(x)}\frac{\partial \phi}{\partial y}\right] = \textrm{i}\omega \int_{S} v_{y} w, \end{gather}$$
$$\begin{gather} \int_{S} \frac{2}{R} q \frac{\partial}{\partial x}[(x+1)u] + \int_{S} \frac{2}{R}\frac{h'_{0}(x)}{h_{0}(x)}(x+1) \left[q u + \frac{\partial q}{\partial y} u\right] + \int_{S} \frac{2\textrm{i} m}{R}q v\nonumber\\ -\int_{S} \frac{2}{h_{0}(x)}(x+1)\frac{\partial q}{\partial y} w + \int_{x={-}1}^{x=1} \frac{2\textrm{i} m F (x+1)\varOmega_{0}(x)}{h_{0}(x)}q(x, 1)\phi(x, 1)\nonumber\\ =\textrm{i} \omega \int_{x={-}1}^{x=1} \frac{2(x+1)}{h_{0}(x)} q(x, 1)\phi(x, 1). \end{gather}$$
Let us define the bilinear forms 
$\mathcal {A}: V_{H} \times V_{v} \rightarrow \mathbb {R}$ and 
$\mathcal {B}: V_{H} \times V_{v} \rightarrow \mathbb {R}$ such that the generalized eigenvalue problem above can be compactly written as – find 
$([u, v, \phi ], w) \in V_{H} \times V_{v}$ such that
\begin{align} &\mathcal{A}([u, v, w, \phi], [v_{x}, v_{t}, v_{y}, q]) = \omega \mathcal{B}([u, v, w, \phi], [v_{x}, v_{t}, v_{y}, q]), \nonumber\\ &\qquad \forall ([v_{x}, v_{t}, q], v_{y}) \in V_{H} \times V_{v}. \end{align}
At this point, in order to discretize the problem, we need to expand the unknowns in terms of proper basis functions. Such basis functions are taken in such a way that the homogeneous Dirichlet boundary conditions are automatically satisfied in both the axial and radial coordinates. In particular, we define the set of polynomials, 
$P_{n}^{\star }(x)$, as follows:
\begin{equation} P_{n}^{{\star}}(x) = P_{n}(x) + P_{n-1}(x), \quad n\ge 1, \end{equation}
where 
$P_{n}(x)$ are standard Legendre polynomials.
Thanks to the above definition, a Dirichlet boundary condition at 
$x = -1$ is automatically satisfied. Therefore, we expand the velocity components and the pressure as
\begin{equation} \left.\begin{gathered}
{[u, v, \phi]}(x, y) =
\sum_{i=1}^{N_x}\sum_{j=1}^{N_y}[u_{ij},
v_{ij},\phi_{ij}]P_{i}^{{\star}}(x)P_j(y), \\ w(x, y) =
\sum_{i=1}^{N_x}\sum_{j=1}^{N_y}w_{ij}P_{i}^{{\star}}(x)P_{j}^{{\star}}(y).
\end{gathered}\right\}
\end{equation}After substituting the expansion into the weak formulation (A6), we obtain the final discrete eigenvalue problem, which reads
\begin{equation} \boldsymbol{\mathsf{A}}\boldsymbol{w} = \omega \boldsymbol{\mathsf{B}}\boldsymbol{w}, \end{equation}
where 
$\boldsymbol {w} = (u_{ij}, v_{ij}, w_{ij}, \phi _{ij})$ is our array containing the spectral coefficients previously defined, 
$\boldsymbol{\mathsf{A}}$ and 
$\boldsymbol{\mathsf{B}}$ are the matrices of order 
$4N_{x}N_{y} \times 4N_{x}N_{y}$ coming from the evaluation of the integrals appearing in the weak formulation for each equation.
A.1.2. Discrete problem for axisymmetric perturbations
In the axisymmetric case (
$m = 0$), the differential problem (A2) over the square 
$S$ simplifies, and it can be noted that only the radial component of the velocity must go to zero as 
$r\rightarrow 0$, with the other unknowns allowed to take any finite value. Thus, studying the axisymmetric perturbation problem means to look for the eigensolutions in the following spaces: 
$u \in V_{H}(S)$, 
$(v, \phi ) \in H^{1}(S)$ and 
$w \in V_{v0} = \{ w \in H^{1}(S): w = 0, \textrm {at} \ z = 0 \}$. The weak formulation is then obtained as shown previously. Regarding the discretization process, we express the four unknowns (and corresponding test functions) as follows:
\begin{equation} \left.\begin{gathered} u(x, y) = \sum_{i=1}^{N_x}\sum_{j=1}^{N_y}u_{ij}P_{i}^{{\star}}(x)P_j(y), \\ [v,\phi](x, y) = \sum_{i=1}^{N_x}\sum_{j=1}^{N_y}[v_{ij}, \phi_{ij}]P_{i}(x)P_j(y), \\ w(x, y) = \sum_{i=1}^{N_x}\sum_{j=1}^{N_y}w_{ij}P_{i}(x)P_{j}^{{\star}}(y). \end{gathered}\right\} \end{equation}
After substitution of the expansion above into the weak formulation, the resulting algebraic generalized eigenproblem (A9) is recalled in the case 
$m = 0$.
A.2. Spurious numerical modes and resolvedness conditions
Solving the discretized eigenvalue problem (A9) is known to give rise to spurious eigenvalues and eigenfunctions. These are numerical artefacts caused by numerical under-resolution of highly oscillatory modes; that is, eigenfunctions do satisfy the discretized problem but do not approximate solutions of the continuous eigenvalue problem. Following Brambley & Peake (Reference Brambley and Peake2008) (see also Brambley Reference Brambley2007, pp. 78–81), we implement two tests to remove such spurious numerical modes.
The first condition requires that the eigenvalues do not move significantly when the numerical discretization is changed. Suppose eigenvalues 
$\{\omega _{p}\}$ have been computed using a discretization 
$(N_x, N_y)$, and suppose 
$\{\hat {\omega }_{q}\}$ have been computed with a slightly increased resolution, such as 
$(N_{x}+1, N_{y}+1)$. The first condition is then given by requiring that, for each eigenvalue 
$\omega _p$,
\begin{equation} \min_q|\omega_{p} - \hat{\omega}_{q}| < \text{tol}, \end{equation}
where 
$\textrm {tol}$ is a prescribed tolerance. For the results computed in this paper 
$\textrm {tol} = 10^{-2}$ was used.
The second resolvedness condition involves the spectral coefficients of the eigenfunctions. In particular, we require the modulus of the spectral coefficients to decay smoothly as the number of polynomials used gets larger. Both the aforementioned conditions have been applied to the pressure eigenfunctions as this is the unknown of most interest because it contains information on the shape of the free surface height. So, let 
$\phi _{ij}$ be a generic spectral coefficient, the following condition has been implemented:
\begin{equation} \frac{\sup_{(i, j) \in \mathcal{B}} |\phi_{ij}|}{\sup_{{all}\,i, j} |\phi_{ij}|} < \text{tol}, \end{equation}
where 
$\mathcal {B}$ is a part of the domain in the 
$N_{x}-N_{y}$ plane, defined as
\begin{equation} \mathcal{B} = \{ (i, j): N_{x}+1 - b_{x} \le i \le N_{x}+1 \text{ or } N_{y}+1 - b_{y} \le j \le N_{y}+1 \}, \end{equation}
$\textrm {tol}$ is a prescribed tolerance on the magnitude of the coefficients and 
$b_x, b_y$ are two prescribed borders widths in the 
$N_{x}-N_{y}$ plane. Figure 21 compares the magnitude of the coefficients for a well resolved mode and a poorly resolved mode for 
$N_{x} = 50$, 
$N_{y} = 20$, 
$b_x = 12$ and 
$b_y = 4$. The well-resolved mode can be seen to have spectral coefficient 
$|\phi _{ij}|$ decreasing exponentially quickly. By contrast, for the poorly resolved mode, the spectral coefficients 
$|\phi _{ij}|$ are oscillatory increasing as function of the number of polynomials used.
Example of a well-resolved mode (a) and a poorly resolved mode (b).
Appendix B. Damping layer formulation
In § 3.1, the continuity equation was modified by the introduction of the ‘damped compressibility’ term 
$\xi (r)\phi$ in (3.1) in order to emulate an infinite domain using our finite computational domain, by damping out disturbances as they reach the unphysical computational boundary at 
$r=R$. Here, we justify the inclusion of that term by analogy with the equations of acoustics.
Small (linear) homentropic perturbations 
$(\boldsymbol {u}, p, \rho )$ to a static fluid 
$(\boldsymbol {0}, p_0, \rho _0)$ are governed by the linearized Euler equations,
\begin{equation} \rho_{0}\frac{\partial\boldsymbol{u}}{\partial t} + \boldsymbol{\nabla} p = 0, \quad \frac{1}{c^2}\frac{\partial p}{\partial t} + \rho_0 \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \quad p = c^{2}\rho, \end{equation}
where 
$c$ is the speed of sound. These can be combined into a single wave equation for the pressure:
\begin{equation} \frac{1}{c^2}\frac{\partial^2 p}{\partial t^2} - \nabla^{2}p = 0. \end{equation}
Following Gao et al. (Reference Gao, Song, Zhang and Yao2017, pp. 81–82), in order to introduce a sponge layer to damp outgoing waves, we add a damping term of the form 
$\xi (r)\partial _{t}p$ into the wave equation (B2),
\begin{equation} \frac{1}{c^2}\frac{\partial^2 p}{\partial t^2} + \xi(r)\frac{\partial p}{\partial t} - \nabla^{2}p = 0. \end{equation}
This wave equation can be split back into the original physical mass- and momentum-equations in the original physical variables 
$(\boldsymbol {u}, p, \rho )$ as
\begin{equation} \rho_{0}\frac{\partial\boldsymbol{u}}{\partial t} + \boldsymbol{\nabla} p = 0, \quad \frac{1}{c^2}\frac{\partial p}{\partial t} + \xi(r)p + \rho_0 \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \quad p = c^{2}\rho. \end{equation}
If we now take the incompressible limit 
$1/c^2 \to 0$ in system (B4), we obtain the modified continuity equation (see (3.1)) introduced in § 3.1.
Figure 22 and table 3 show the effect of the compressibility 
$\alpha = 1/c^2$ on the structure of the radiating eigenmode and the eigenvalue, respectively, for parameters 
$m = 7$ and 
$F = 0.3$. From both, it is clear that the eigenmode computed does not vary significantly as 
$\alpha \to 0$, and so we conclude that the damping properties of the 
$\xi (r)p$ term are carried over in the limit 
$\alpha \to 0$.
Comparison of the structure of the eigenmode 
$|h(r, 0, 0)|$ for 
$m = 7$, 
$F = 0.3$ and different value of the compressibility 
$\alpha = 1/c^2$.
Eigenvalues as function of the compressibility 
$\alpha = {1}/{c^2}$ for a radiating mode, for 
$m = 7$, and 
$F = 0.3$.
Appendix C. Code validation and comparison with results for other background flow
Here we provide evidence of the correctness of our numerical methods. As a first test, we have validated our code against the results of Mougel et al. for surface waves on a Rankine vortex in a bounded domain (Mougel et al. Reference Mougel, Fabre and Lacaze2014, figure 5). We have adapted our code to this base flow and geometry by including a rigid wall at finite radius, and therefore not using a damping region to model an unbounded domain. With these changes, we find excellent agreement between our computations and theirs. For example, at 
$Fr = 1.502$ they find an unstable surface wave with eigenvalue 
$2.6101 + 0.0025313\mathrm {i}$ (J. Mougel, private communication). We confirm this instability and find 
$\omega = 2.6111+0.004556\mathrm {i}$. The moduli of these eigenvalues agree to better than 0.04 %. Mougel et al. include a small amount of viscous damping in their computations, while we do not, and this can account for the small absolute difference in the imaginary parts of the respective eigenvalues. By including an artificial viscous term (with artificial viscosity 
$\nu$) of the form 
$\nu \partial ^{2}_{x} u$ in the right-hand side of (A2a), we have confirmed the sensitivity of the small imaginary part of the eigenvalue to damping, but we have not attempted to reproduce the exact damping used by Mougel et al..
As a second test case, we have validated our code against results of Mougel et al. (Reference Mougel, Fabre and Lacaze2015) for Newton's bucket problem: solid body rotation with a rigid walls at 
$r=R$. We have adapted our code to this case. Results from our numerics are plotted in figure 23 and can be compared with figure 3 of Mougel et al. (Reference Mougel, Fabre and Lacaze2015, p. 223). For this case, we adopt the same non-dimensionalization as Mougel et al. (Reference Mougel, Fabre and Lacaze2015), so that we can directly compare our numerical solutions with theirs.
Eigenmodes comparison with Newton's bucket problem (Mougel et al. Reference Mougel, Fabre and Lacaze2015) using their same parameters and our artificial viscosity 
$\nu = 0.01$: (a) 
$\lambda = 4.2897 - 0.008\mathrm {i}$; (b) 
$\lambda = 7.4359 - 0.0123\mathrm {i}$; (c) 
$\lambda = 9.0436 - 0.0368\mathrm {i}$; (d) 
$\lambda = -5.2014 - 0.0029\mathrm {i}$; (e) 
$\lambda = -7.4580 - 0.0156\mathrm {i}$; ( f) 
$\lambda = -9.0349 - 0.0408\mathrm {i}$; (g) 
$\lambda = 1.4077 - 0.0112\mathrm {i}$; (h) 
$\lambda = 1.0974 - 0.0254\mathrm {i}$; (i) 
$\lambda = 1.7833 - 0.0180\mathrm {i}$; ( j) 
$\lambda = -1.6985 - 0.0111\mathrm {i}$; (k) 
$\lambda = -1.3292 - 0.0304\mathrm {i}$; (l) 
$\lambda = -1.9142 - 0.0122\mathrm {i}$; (m) 
$\lambda = -0.0346 - 0.0086\mathrm {i}$; (n) 
$\lambda = -0.0124 - 0.0139\mathrm {i}$; (o) 
$\lambda = -0.0077 - 0.0166\mathrm {i}$. Panels (a–f) represent surface gravity modes. Panels (g–l) represent row inertial modes. Panels (m–o) represent Rossby modes.
To obtain smooth and well-resolved eigenfunctions for Newton's bucket, we have again added the artificial viscous term explained just above for the Rankine vortex. Excellent agreement is found in this case as well, giving further confidence in the correctness of our numerics. For example, by setting our artificial viscosity 
$\nu = 0.01$, the eigenvalue of the first gravity mode from our computation is 
$\lambda = 4.2897 - 0.008\mathrm {i}$, whereas they found 
$\lambda = 4.290 - 0.0014 \mathrm {i}$. For the first inertial mode we find 
$\lambda = 1.4077 - 0.0112\mathrm {i}$, against their result 
$\lambda = 1.408 - 0.0083\mathrm {i}$. Finally for the first Rossby mode we have 
$\lambda = -0.0346 - 0.0086\mathrm {i}$, whereas they found 
$\lambda = -0.0372 - 0.0021\mathrm {i}$.
As a third and final test case, we have computed the stability the steady Lamb–Oseen vortex for the shallow-water equations presented in Ford (Reference Ford1994). Using the same convention as Ford, we obtain an unstable linear mode at the following set of parameters: 
$F = 1$, 
$m = 2$ and 
$f = 0$, where 
$f$ is the Coriolis parameter. Specifically, for these parameter values our code returns the most unstable eigenvalue 
$\omega = 0.2606 + 0.0017\mathrm {i}$, and this very agrees well with the eigenvalues plotted in Ford (Reference Ford1994, figure 4), even though the vorticity profile used here and by Ford are not identical. We again confirm that our code finds instability where it is known to exist.
Appendix D. Results without the base free surface deformation
In this appendix, we demonstrate that the dominant contribution to the trapping of modes by the vortex shown in § 4 is the swirl of the base flow, and that the deformation of the base free surface has little effect itself. In the following we artificially assume the base free surface to be flat, meaning that 
$h_{0}(r) = h_{\infty }$ and 
$h'_{0}(r) = 0$. All the remaining terms in the governing equations and boundary conditions do not change. Solving the eigenvalue problem in this case, we compare the eigenvalues with the complete case where the base free surface deformation has been taken into account in figure 24. Little difference is seen between the two. Indeed, figure 25 plots the difference in eigenvalues as function of the Froude number; i.e. 
$|\omega - \omega _{{flat}}|$, where 
$\omega _{{flat}}$ are the eigenvalues without the free surface variation. It is clear that for the range of Froude numbers considered, the two agree remarkably closely.
Imaginary and real parts of the eigenvalues as functions of the Froude number for modes with 
$m = 7$ and 
$n = 1$. Blue-circle curve: full numerics; red-asterisk curve: without the free surface contribution.
Absolute difference between eigenvalues computed with and without the free surface contribution, 
$|\omega - \omega _{{flat}}|$, for the first radial mode 
$n=1$ and 
$m = 7$.
In terms of the eigenfunctions, the free surface perturbation is displayed in figure 26, which should be compared with the free surface perturbation with a base free surface height variation plotted in figure 7; the two figures can be seen to be almost identical. The corresponding structure in the 
$r$–
$z$ plane is shown in figure 27, which should be compared with figure 8 in the results section; again, the two figures are almost identical apart from the base flow height.
Plots of the perturbation free surface height 
$h(r, \theta, t=0) = \textrm {Re}[\phi (r, h_{0}(r))\exp \{\mathrm {i} m\theta \}]$ for 
$m = 7$ without the base free surface deformation, to be compared with figure 7: (a) 
$F = 0.3$; (b) 
$F = 0.32$; (c) 
$F = 0.33$; (d) 
$F = 0.34$; (e) 
$F = 0.36$. All modes displayed rotate clockwise, i.e. against the vortex flow.
Plots of the pressure distribution 
$|\phi _{m}(r, z)|$ for 
$m = 7$ without the base free surface deformation, to be compared with figure 8: (a) 
$F = 0.3$; (b) 
$F = 0.32$; (c) 
$F = 0.33$; (d) 
$F = 0.34$; (e) 
$F = 0.36$.
We therefore conclude that the driving mechanism behind the trapping of modes by the vortex is the rotation of the base flow, and not the free surface deformation.
Open access
