# A primer on the analogue black hole bomb with capillary-gravity waves

Sam Patrick<EMAIL_ADDRESS>Department of Physics, King’s College
London, University of London, Strand, London, WC2R 2LS, UK  Théo Torres  Aix
Marseille Université, CNRS, Centrale Marseille, LMA UMR 7031, Marseille,
France

###### Abstract

Draining vortices with a free surface are frequently employed as rotating
black hole simulators, both in theory and experiments. However, most
theoretical work is restricted to the idealised regime, where wave dispersion
and dissipation are neglected. We investigate the role of these effects on the
analogue black hole bomb, an instability resulting from rotational
superradiant amplification in confined systems. We reveal that the dispersion
of deep water capillary-gravity waves significantly modifies the unstable mode
eigenfrequencies, whereas viscosity only affects those with high frequencies.
Furthermore, if the circulation is less than an order 1 multiple of the drain
rate, superradiance does not occur and the vortex is stable. The instability
is maximised in small systems with high flow velocities, provided there is
sufficient space between the vortex and the outer boundary for the first
excited state to lie inside the superradiant bandwidth. Implications for
experiments on analogue black holes and free surface vortices are discussed.

##  I Introduction

Draining vortices are the archetypal example of an analogue rotating black
hole [1]. The basic idea is that the draining fluid speeds up as it impinges
on the central outlet, and if the fluid velocity manages to surpass the speed
of long wavelength excitations of the fluid, the system possess an analogue
horizon. The presence of a horizon leads to the possibility of simulating
Hawking radiation from a black hole [2, 3]. Extensive work on the analogue
Hawking effect has been performed in one dimensional (1D) flows where the
fluid accelerates linearly and the analysis is simpler. Detailed numerical [4,
5] and analytic [6, 7, 8] analyses have revealed the persistence of the effect
in fluid-like systems and are supported by experimental measurements in both
classical [9, 10] and quantum systems [11, 12]. These results are generally
taken as an indication of the universality of the Hawking effect, that is, a
robustness against modifications to the short wavelength part of the governing
theory [13].

In two dimensions (2D), the analogy with black holes becomes much richer. In
particular, rotating systems like the draining vortex can also possess an
analogue ergosphere, the region of spacetime around a black hole where an
observer would need to travel faster than the speed of light to remain at rest
[14]. The existence of such a region implies the amplification of waves in a
certain frequency band, i.e. rotational superradiance [15, 16], which causes
an evaporating black hole to shed its angular momentum more rapidly than its
mass [17, 18]. Many theoretical studies demonstrate the occurrence of analogue
superradiance using a simple, shallow water model of a draining vortex [19,
20, 21, 22]. The effect was first measured in a water tank [23] and
subsequently in a (non-draining) photon superfluid [24].

Most experimental studies of the Hawking effect and superradiance to date have
employed one of two types of simulator. The first involves excitations
analogous to sound waves (density perturbations) in a Bose-Einstein condensate
(BEC), e.g. [11, 24]. The second employs a classical liquid, e.g. room
temperature water, where the excitations are surface waves [9, 10, 23]. Since
gravity and surface tension act as a restoring force, these are also known as
capillary-gravity waves (where capillary is synonymous with surface tension).
In both types of simulator, the dispersive nature of these waves complicates
the analogy with black holes. In relativity, all massless particles propagate
at the speed of light $c$. Similarly, in fluids, it is often the case that
long wavelength modes propagate at a constant (i.e. frequency independent)
speed $c$ – the “analogue speed of light”. However, at short wavelengths, the
waves probe the microscopic details of the fluid and the propagation speed
becomes frequency dependent. If short wavelengths propagate slower than $c$
(e.g. deep water gravity waves) the theory is called sub-luminally dispersive,
in analogy with relativity, whereas if they propagate faster than $c$ (e.g.
capillary waves and Bogoliubov modes in a BEC) the theory is called super-
luminally dispersive. In the case of capillary-gravity waves, the theory is a
mixed type, containing a sub-luminal part at intermediate wavelengths (the
deep water regime) and a super-luminal part as the wavelength goes to zero
(the capillary regime). The effect of each type of dispersion relation on
superradiance around a draining vortex was investigated in [25, 26]. The
effect of mixed dispersion in rotating systems that superradiate has yet to be
understood.

Through the analogy, many dynamical processes in supersonic fluid flows can be
interpreted in terms of their black hole scattering counterparts. For
instance, the Hawking effect in a trans-sonic BEC gives rise to the black hole
laser instability [27, 28, 29]. The Hawking effect is also responsible for the
phenomenon of the zero frequency undulation that appears when a fluid
accelerates over a submerged obstacle [30]. In rotating systems, the
relativistic phenomenon of rotational superradiance is analogous to the fluid
mechanical notion of over-reflection in shear flows [31]. In quantum fluids,
superradiance causes large vortices to decay into clusters of single vortex
quanta [32, 33, 34] whilst in classical fluids, it is related to the rotating
polygon instability of hollow core, free surface vortices [35].

In this work, we investigate the occurrence of another well-known black hole
instability in a draining water tank: the black hole bomb (BHB) [36]. This was
first studied by Press and Teukolsky [37] and involved a black hole surrounded
by an idealised spherical mirror. The essence of the BHB is that waves in a
certain frequency range get (superradiantly) scattered outward near the
boundary of the ergosphere then get partially or completely reflected back in
at some location far away from the black hole. These two reflecting locations
form a cavity which can trap the waves, and specific frequencies which fit
inside the cavity (by satisfying the appropriate boundary conditions) will be
resonances/bound states of the system. If the resonance lies within the
superradiant bandwidth, it will be unstable and grow exponentially in time.
Hence, the BHB is simply a superradiant bound state which occupies a region
outside the black hole’s ergosphere.

In gravity, a natural reflection mechanism arises for light massive fields,
which are prevented from propagating at infinity by the effect of their own
mass. This has led to intense interest on ultralight axions, which are
speculated to arise in various standard model extensions [38]. These fields
can form clouds around black holes through the superradiant instability,
sapping rotational energy from the system and dissipating into gravitational
waves. This process leads to distinct observational signatures, raising the
intriguing possibility of using black holes as observatories for axions [39,
40, 41, 42]. However, it has been argued that nonlinear interactions could
saturate growth of the cloud before it manages to spin down the black hole
[43]. It is therefore pertinent to develop a more thorough understanding of
nonlinear effects on BHB–type instabilities. In this regard, analogue models
of the BHB in fluids (which possess their own unique nonlinearities) could
help shed light on the influence of nonlinearities in a general sense, and can
be implemented experimentally to test predictions.

BHBs in analogue systems have already been discussed, e.g. in [21, 44].
However, these works do not take the inherent dispersion or dissipation of the
fluid into account, which will be present in any real experiment. For
instance, the detection of superradiance in [23] was performed using deep
water gravity waves, which are strongly dispersive. Subsequent theoretical
analysis [25] demonstrated that the character of superradiance (i.e. the
amount of amplification and the bandwidth of amplified frequencies) can differ
significantly from the simple long wavelength systems studied by many authors.
Furthermore, damping (e.g. due to viscosity) could significantly alter the
instability timescales or, in the most extreme cases, fully quench the
instability and stabilise the system. In view of this, our goal in this work
will be to assess the impact of dispersion and dissipation on the BHB
instability in analogue fluid systems. Using a simplified model of a draining
vortex, we will develop a calculational tool for computing unstable mode
frequencies in strongly dispersive settings, laying the theoretical ground
work for future experiments.

###  I.1 Method

To achieve this goal, the most suitable existing theory is the deep water
theory of superradiance outlined in [25]. The theory is based on a Wentzel-
Kramers-Brillouin (WKB) approximation, which considers wavelengths shorter
than the scale over which the background varies. The method reduces a
complicated wave scattering problem governed by partial differential equations
(PDEs) to a far simpler algebraic equation which governs the ray trajectories
through phase space. The rays are the paths traced out by a collection of
points along a wavefront of short wavelength (high frequency) waves and are
analogous to spacetime geodesics in relativity. Since conventional
relativistic theories are second order in space and time derivatives, the
solution of the linearised theory at a given frequency will be a superposition
of two waves, representing the two allowed directions of propagation (in-going
and out-going). Hence, in the WKB theory, there will be two rays and
intersections of the rays (caustics) are the locations where the wave scatters
off the inhomogeneous background medium. Similarly, the long wavelength (or
shallow water) limit of linear fluid dynamics is second order in space and
time derivatives and the two solutions correspond to the in/out-going waves in
the regions where the flow velocity is weak. For shorter wavelengths, however
(which is the natural regime for WKB techniques), higher order spatial
derivatives in the linear fluid equations become important and the general
solution becomes a superposition of more than two waves. Correspondingly, the
WKB theory contains more than two rays. One then has to consider all possible
intersections between neighbouring rays to fully understand the scattering of
waves off the inhomogeneous background. This was the approach taken in [25]
(and subsequently [45, 33]) which introduced a method of schematically
sketching the interactions between different rays to facilitate the
computation of scattering coefficients. These scattering diagrams will be the
principle tool used in this work as we attempt to understand the dispersive
BHB around a draining vortex.

However, the deep water theory in [25] contains unphysical negative energy
modes (with divergent wavenumber $k$) asymptotically far from the vortex core
and is therefore not suitable for an analysis of the BHB, which is dependent
on the nature of reflections at the outer boundary. These unphysical modes
arise when the (sub-luminal) deep water theory is applied beyond its regime of
validity. Surface tension regulates the theory by making it super-luminal at
large $k$, so that in the weak field (i.e. slow flow) limit there are only two
propagating modes outside the vortex corresponding to the usual in/out-going
ones. It is these modes which get reflected at the (mirror) boundary and give
rise to the BHB. Therefore, the primary aim of this work will be to establish
the WKB theory of capillary-gravity waves around a draining vortex, and derive
all the allowed scattering diagrams. Using the diagrams, we will identify the
necessary conditions for superradiant amplification then apply these to the
calculation of oscillation frequencies and growth rates of the BHB modes.

###  I.2 Outline and summary of results

In Section II, we start by reviewing the equations of fluid mechanics in the
context of the draining vortex, before deriving the equations of motion for
capillary-gravity waves. Then in Section III, we apply the WKB method to these
equations and introduce the necessary concepts to study the scattering of rays
in the inhomogeneous draining vortex background. Next, we illustrate how the
framework can be applied to the simpler case of shallow water waves in Section
IV. In this regime, the equations of motion can be cast in terms of an
effective rotating black hole metric, allowing us to identify the superradiant
instability as a BHB instability. The main result is that the oscillation
frequency $\omega_{n}$ is determined by a resonance condition of the form,

$\mathrm{Re}[S_{\mathrm{path}}(\omega_{n})]=\pi\left(2n+\frac{1}{2}\right),$
(1)

where $S_{\mathrm{path}}$ is the phase integral for the radial part of the
mode (which has to be evaluated numerically) along a trapped path consisting
of a pair of rays, and $n=0,1,...$ indexes the different possible bound
states. We find an expression for the growth rate $\Gamma_{n}$ which has the
intuitive interpretation that the exponential growth of the mode is equal to
the amount of superradiant amplification i.e.

$e^{\Gamma_{n}T_{\mathrm{path}}(\omega_{n})}=|\mathcal{R}(\omega_{n})|,$ (2)

where $T_{\mathrm{path}}$ is the time elapsed along the rays in the trapping
region, and $|\mathcal{R}|$ is the magnitude of the reflection coefficient
characterising the scattering of an in-going ray with the inhomogeneous
draining vortex flow. When rotational superradiance occurs, we have
$|\mathcal{R}|>1$ (i.e. amplification) and can immediately see that
$\Gamma_{n}>0$ so that the eigenmode is unstable. Since the destabilising
mechanism is superradiant scattering, we identify these modes as the BHB
modes. We calculate the complex eigenfrequencies
$\omega^{\mathbb{C}}_{n}=\omega_{n}+i\Gamma_{n}$ of the BHB modes for a range
of parameters, comparing our WKB predictions with exact numerical calculations
and finding excellent agreement. We then identify the key features of the
instability, allowing us to generalise the BHB to the dispersive regime.

In Section V, we analyse the scattering of capillary-gravity waves with the
draining vortex, identifying all possible scattering diagrams in the same vein
as [25]. From the diagrams, we infer a condition for superradiant
amplification of the form,

$0<\omega<\omega_{\star},$ (3)

where $\omega_{\star}$ is a characteristic location on the dispersion relation
(to be discussed in detail). We find in Fig. 12 that this band opens at a
critical rotation, hence, there is a minimum required speed of rotation for
superradiance to occur. With this information, we obtain a condition for the
existence of BHB modes in the system, accounting for both the dispersion and
dissipation of the fluid. The oscillation frequency is determined by a
resonance condition which has the same form as the shallow water formula (1).
Intuitively, one would expect the formula (2) for the growth rate to be
modified by the amount of amplification on the right-hand side becoming
suppressed. This can happen in two ways: firstly, amplification can be
weakened by dissipation and secondly, scattering with the extra dispersive
modes can lead to energy escaping the trapping region. We find that this
expectation is indeed correct and that $\Gamma_{n}$ is governed by the
relation,

$e^{\Gamma_{n}T_{\mathrm{path}}(\omega_{n})}=e^{-\nu\mathrm{Re}[\Sigma_{\mathrm{path}}(\omega_{n})]}\prod_{a}|X_{a}(\omega_{n})|,$
(4)

where the $X_{a}$ are the different scattering coefficients along the trapped
path consisting of a pair of two or more rays, and $\nu\Sigma_{\mathrm{path}}$
is the integral of the damping factor over the trapping region. The different
terms in this expression will be explained in detail in the main text and the
origin of the formula in Appendix F. We then calculate
$\omega^{\mathbb{C}}_{n}=\omega_{n}+i\Gamma_{n}$ for some exemplary system
parameters, finding in Fig. 15 that generally the instability is strongest in
small systems (i.e. when the encircling mirror has a small radius) when the
rotation speed is large and the drain is near the maximum allowed value that
permits a trapping region outside the vortex. We conclude in Section VII by
discussing future directions for our work.

##  II The system

In this section, we present the stationary background flow for our draining
vortex before deriving a set of linearised equations governing small amplitude
capillary-gravity waves. The key equations are the draining vortex flow
profile (11), the linear equations of motion (19) and statement of norm
conservation in (30).

###  II.1 Basic equations

We consider a cylindrical tank of water with a hole of radius $r=d$ in the
centre of the plane at $z=0$. Water drains toward the centre and is resupplied
at the edge $r=r_{B}$ of the cylindrical tank. The fluid has a free surface
located at $z=h(x,y)$. A schematic of the set-up is depicted on Fig. 1.

The incompressible Navier-Stokes equations in three spatial dimensions are,

$(\partial_{t}+\mathbf{v}\cdot\bm{\nabla})\mathbf{v}-\nu\nabla^{2}\mathbf{v}=-\frac{\bm{\nabla}P}{\rho}+\mathbf{g},$
(5)

where $\mathbf{v}$ is the fluid velocity, $P$ is the pressure and (for water)
the density $\rho\simeq 997\mathrm{kgm^{-3}}$ (which is taken to be constant
at room temperature) and the kinematic viscosity $\nu\simeq 8.93\times
10^{-7}\mathrm{m^{2}s^{-1}}$. The acceleration due to gravity is
$\mathbf{g}=-g\hat{\mathbf{e}}_{z}$ (with $g\simeq 9.81\mathrm{ms^{-2}}$ at
the surface of the earth) and we assume no other external forces. Since the
fluid is incompressible, the velocity field satisfies,

$\bm{\nabla}\cdot\mathbf{v}=0.$ (6)

We further assume an irrotational fluid satisfying,

$\bm{\nabla}\times\mathbf{v}=0\quad\Rightarrow\quad\mathbf{v}=\bm{\nabla}\Phi,$
(7)

where $\Phi$ is the velocity potential. Under these conditions, (5) and (6)
can be reformulated as the Bernoulli and Laplace equations,

$\partial_{t}\Phi+\frac{1}{2}\mathbf{v}^{2}+\frac{P}{\rho}+gz=\mathrm{const},\qquad\nabla^{2}\Phi=0,$
(8)

where the relation
$\nabla^{2}\mathbf{v}=\bm{\nabla}(\bm{\nabla}\cdot\mathbf{v})-\bm{\nabla}\times\bm{\nabla}\times\mathbf{v}$
has been used to drop the term proportional to $\nu$ and
$\mathbf{v}\cdot\bm{\nabla}\mathbf{v}=\frac{1}{2}\bm{\nabla}(\mathbf{v}^{2})-\mathbf{v}\times\bm{\nabla}\times\mathbf{v}$
to re-express the advective term.

Figure 1: A schematic illustration of our set-up. We display only a cross
section in the $(r,z)$–plane. There is a hard wall boundary situated at
$r=r_{B}$. Water enters through a cylindrical inlet at the bottom of the
container near the outer wall, then moves inward with velocity $v_{r}(r)$ and
exits through a hole of radius $d$ in the centre. The velocity field also has
a component in $\theta$–direction, not shown. The surface is given by the
function $h(r)$. In this work, we consider the dynamics of surface waves, also
schematically depicted.

The boundary condition at the free surface is [46, 47],

$P(z=h)=P_{a}-\sigma\bm{\nabla}_{\parallel}\cdot\left(\frac{\bm{\nabla}_{\parallel}h}{\sqrt{1+|\bm{\nabla}_{\parallel}h|^{2}}}\right)+2\rho\nu\partial_{z}^{2}\Phi$
(9)

where $P_{a}$ is atmospheric pressure which we take to be zero,
$\bm{\nabla}_{\parallel}$ is the gradient operator parallel to the $(x,y)$
plane and $\sigma\simeq 7.2\times 10^{-2}\mathrm{kgs^{-2}}$ is the surface
tension. The second term on the right results from the curvature of the free
surface and the final term encodes viscous friction within the fluid.
Bernoulli’s equation at the free surface is then,

$\begin{split}&\Bigg{[}\partial_{t}\Phi+\frac{1}{2}\mathbf{v}_{\parallel}^{2}+\frac{1}{2}v_{z}^{2}+gh-2\nu\nabla^{2}_{\parallel}\Phi\\\
&\qquad\qquad-\frac{\sigma}{\rho}\bm{\nabla}_{\parallel}\cdot\left(\frac{\bm{\nabla}_{\parallel}h}{\sqrt{1+|\bm{\nabla}_{\parallel}h|^{2}}}\right)\Bigg{]}_{z=h}=0,\end{split}$
(10)

whilst Bernoulli’s equation in the bulk simply gives an expression for the
pressure. Note, the assumption that the fluid motion is irrotational means
that the boundary condition on the tangential stress has been dropped. This
leads to an underestimation of dissipation since we are not accounting for
losses into the rotational degrees of freedom in the fluid. A more thorough
treatment for 1D deep water waves in a weakly viscous fluid leads to the
replacement $\nu\to 2\nu$ in the results, see e.g. [48]. Since we are only
concerned with the generic effect of dissipation, we do not consider such
losses here.

If the basic flow is assumed stationary and axisymmetric, it is natural to
work in polar coordinates $(r,\theta)$. In the region where the fluid height
is approximately uniform, the basic flow is given by,

$\mathbf{v}=-\frac{D}{r}\hat{\mathbf{e}}_{r}+\frac{C}{r}\hat{\mathbf{e}}_{\theta},\qquad
h\simeq h_{0}=\mathrm{const},$ (11)

which is valid in the region $r\gg r_{a}\equiv\sqrt{(C^{2}+D^{2})/2gh_{0}}$ We
refer the reader to Appendix A for a thorough discussion of the approximations
that lead to (11). Even working with these simplified profiles, solving the
scattering problem over the full parameter space is a challenging task,
especially with the inclusion of dispersive effects. Hence, order of magnitude
estimates for the various parameters will be useful when it comes computing
bound state frequencies. We assume that $d$ and $r_{B}$ are respectively on
the order of centimeters and tens of centimeters. To get an estimate of the
flow parameters achievable in such a set-up, we assume that (in practice)
$r_{a}$ will be some significant fraction of the radius of the drain hole
radius, say $r_{a}\sim\mathcal{O}(d)$. Assuming $h_{0}$ in the range of
millimeters to tens of centimeters, we find that $C$ and $D$ should roughly
lie in the range $10$–$100\mathrm{cm^{2}s^{-1}}$, in accord with experimental
realisations [49, 23].

###  II.2 Linear equations

To study perturbations of the basic flow, we work in the approximation that
gradients of the free surface can be neglected, i.e. $h^{\prime}\simeq 0$.
This approximation is valid suitably far from the vortex axis, i.e. $r\gg
r_{a}$. We now perturb the fluid variables as,

$\Phi\to\Phi+\phi,\qquad h\to h+\eta,$ (12)

where $\phi$ and $\eta$ are the velocity potential and height perturbations
respectively. (10) gives the first equation,

$\left(\partial_{t}+\mathbf{v}_{\parallel}\cdot\bm{\nabla}_{\parallel}\right)\phi+g\eta-\frac{\sigma}{\rho}\nabla^{2}_{\parallel}\eta-2\nu\nabla^{2}_{\parallel}\phi\Big{|}_{z=h}=0.$
(13)

The second is obtained by solving Laplace’s equation $\nabla^{2}\phi=0$ in the
bulk. In the region $r>d$, we can write a solution in terms of Fourier modes,

$\phi=\int\frac{d^{2}\mathbf{k}}{(2\pi)^{2}}\phi_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{x}},\qquad\phi_{\mathbf{k}}=A_{\mathbf{k}}\cosh(kz),$
(14)

where $k=|\mathbf{k}|$. Using this, we can write,

$\partial_{z}\phi|_{z=h}=F(-i\bm{\nabla}_{\parallel})\phi|_{z=h},$ (15)

where we have defined the operator,

$F(\mathbf{k})=k\tanh(hk).$ (16)

Now, perturbing the kinematic boundary condition (see (A.5) in Appendix A) we
obtain,

$\partial_{z}\phi|_{z=h}+\eta\partial_{z}v_{z}=(\partial_{t}+\mathbf{v}_{\parallel}\cdot\bm{\nabla}_{\parallel})\eta,$
(17)

which can be reformulated when combining with (15) as,

$(\partial_{t}+\bm{\nabla}_{\parallel}\cdot\mathbf{v}_{\parallel})\eta-F(-i\bm{\nabla})\phi|_{z=h}=0,$
(18)

where the divergence operator in parentheses acts on both
$\mathbf{v}_{\parallel}$ and $\eta$ to the right. Finally, we drop the
subscript $\parallel$ and $|_{z=h}$ and understand that all quantities in the
linear theory live on a two dimensional plane at $z=h$. The linearised
equations of motion are,

$\begin{split}\left(\partial_{t}+\mathbf{v}\cdot\bm{\nabla}\right)\phi+g\eta-\gamma\nabla^{2}\eta-2\nu\nabla^{2}\phi=0,\\\
(\partial_{t}+\bm{\nabla}\cdot\mathbf{v})\eta-F(-i\bm{\nabla})\phi=0,\end{split}$
(19)

where we have introduced $\gamma=\sigma/\rho$. The term proportional to $\nu$
encodes the effect of dissipation. The conservative part of the dynamics (i.e.
excluding the term proportional to $\nu$) can be derived from the action,

$\begin{split}\mathcal{S}=&\ \int dtd^{2}\mathbf{x}\bigg{[}-\eta
D_{t}\phi-\frac{1}{2}g\eta^{2}\\\
&\qquad\qquad\qquad-\frac{1}{2}\gamma|\bm{\nabla}\eta|^{2}-\frac{1}{2}\phi
F(-i\bm{\nabla})\phi\bigg{]},\end{split}$ (20)

where $D_{t}=\partial_{t}+\mathbf{v}\cdot\bm{\nabla}$ is the material
(convective) derivative on the free surface.

Although it is not the case for our draining vortex, if the basic flow were
independent of the coordinates, (19) admits an exact solution in terms of
plane waves $e^{i\mathbf{k}\cdot\mathbf{x}-i\omega t}$ with dispersion
relation,

$\omega=-i\nu k^{2}+\mathbf{v}\cdot\mathbf{k}\pm\sqrt{(gk+\gamma
k^{3})\tanh(hk)-\nu^{2}k^{4}},$ (21)

and $k=|\mathbf{k}|$. In Section III, we will see how a local version of this
dispersion relation can be obtained when $\mathbf{v}$ is not constant. We will
then be concerned with two limits. The low $k$ limit of the theory is,

$\omega=\mathbf{v}\cdot\mathbf{k}\pm\sqrt{gh}k,$ (22)

that is, one can neglect the effects of surface tension, viscosity and deep
water effects arising from the hyperbolic tangent. We also consider the deep
water limit $k\gg h^{-1}$, which leads to,

$\omega=-i\nu k^{2}+\mathbf{v}\cdot\mathbf{k}\pm\sqrt{gk+\gamma
k^{3}-\nu^{2}k^{4}}.$ (23)

In this case, the depth of the fluid is effectively infinite (relative to the
wavelength). The price to pay is that, in the regime of short wavelengths, the
effects of surface tension and viscosity are important and should be properly
taken into consideration.

###  II.3 Norm conservation

Splitting the fluctuations into positive and negative frequency components,

$\phi\to\phi+\phi^{*},\qquad\eta\to\eta+\eta^{*},$ (24)

the action (20) becomes,

$\begin{split}\mathcal{S}=&\ \int
dtd^{2}\mathbf{x}\bigg{[}-\eta^{*}D_{t}\phi-\eta
D_{t}\phi^{*}-g|\eta|^{2}-\gamma|\bm{\nabla}\eta|^{2}\\\
&\quad-\frac{1}{2}\phi
F(-i\bm{\nabla})\phi^{*}-\frac{1}{2}\phi^{*}F(-i\bm{\nabla})\phi\bigg{]}.\end{split}$
(25)

To recover the equations of motion, we vary the action with respect to $\phi$
and $\eta$ independently. This operation involves the commutation of
derivative terms, hence, we define $F$ by its Taylor expansion about $k=0$,

$F(\mathbf{k})=\sum_{m=1}^{\infty}F_{m}k^{2m},$ (26)

where we have used the fact that $F$ is an even function of $k$. This allows
us to write,

$\begin{split}\phi
F(-i\bm{\nabla})&\phi^{*}=\left[F(-i\bm{\nabla})\phi\right]\phi^{*}-i\bm{\nabla}\cdot\Bigg{(}\sum_{m=1}^{\infty}F_{m}\times\\\
&\sum_{n=0}^{2m-1}(-1)^{n}(-i\bm{\nabla})^{n}\phi(-i\bm{\nabla})^{2m-1-n}\phi^{*}\Bigg{)}.\end{split}$
(27)

Varying the action then leads to,

$\begin{split}&\delta\mathcal{S}=\ \int
dtd^{2}\mathbf{x}\Bigg{\\{}\delta\eta^{*}\left[-D_{t}\phi-g\eta+\gamma\nabla^{2}\eta\right]+\\\
&\qquad\delta\phi^{*}\left[\partial_{t}\eta+\bm{\nabla}\cdot(\mathbf{v}\eta)-F(-i\bm{\nabla})\phi\right]-\partial_{t}(\eta\delta\phi^{*})\\\
&\qquad-\bm{\nabla}\cdot\bigg{[}\gamma\delta\eta^{*}\bm{\nabla}\eta-\mathbf{v}\delta\phi^{*}\eta+\frac{i}{2}\sum_{m=1}^{\infty}F_{m}\times\\\
&\sum_{n=0}^{2m-1}(-1)^{n}(-i\bm{\nabla})^{n}\phi(-i\bm{\nabla})^{2m-1-n}\delta\phi^{*}\bigg{]}+\mathrm{c.c.}\Bigg{\\}}.\end{split}$
(28)

The total derivatives can be rewritten as terms on the boundary where the
variations are taken to vanish. Setting $\delta\mathcal{S}=0$ then yields the
equations of motion. Alternatively, the total derivative terms can be used to
deduce conserved currents of the system. Indeed, the action is invariant under
continuous phase rotations of the form,

$\phi\to\phi e^{-i\alpha}\simeq\phi+\delta\phi,\qquad\delta\phi=-i\alpha\phi,$
(29)

where $\alpha$ is infinitesimal, and similarly for $\eta$. Stripping out
$\alpha$, we obtain a conservation equation for the norm current,

$\partial_{t}\rho_{n}+\bm{\nabla}\cdot\mathbf{J}_{n}=0,$ (30)

whose components are given by,

$\begin{split}\rho_{n}=&\ 2\mathrm{Im}[\eta\phi^{*}],\\\ \mathbf{J}_{n}=&\
2\mathrm{Im}\bigg{[}\gamma\eta^{*}\bm{\nabla}\eta+\mathbf{v}\phi^{*}\eta+\frac{i}{2}\sum_{m=1}^{\infty}F_{m}\times\\\
&\sum_{n=0}^{2m-1}(-1)^{n}(-i\bm{\nabla})^{n}\phi(-i\bm{\nabla})^{2m-1-n}\phi^{*}\bigg{]}.\end{split}$
(31)

This conserved current is related to the expression for the varying wave
amplitude that we will shortly derive.

##  III WKB method

In the previous section, we derived the set of equations obeyed by capillary
gravity waves (19) in the vortex flow field of (11). In this section, we
present an approximate method of solving these equations. Key concepts will be
the dispersion relation (34) and transport equation (41) which determine the
local wavelength and amplitude of a wave. A formula is given in (55) that
describes the local scattering of waves with the vortex and in (62), we argue
that norm conservation implies superradiant amplification for certain
frequencies.

###  III.1 Ray tracing

The equations of motion in (19) are complicated to solve in practice. However,
simplifications occur if we consider the limit of short wavelengths. In
particular, if the local wavelength of the fluctuations is much smaller than
the scale over which the basic flow varies, we can assume that the solution
looks locally like a plane wave, with an amplitude and wavevector which vary
gradually over the flow. Each point along a wavefront will then trace out a
path as it moves through the system. These paths are known as rays, and are
the fluid mechanical/optical analogue of geodesics in curved spacetime. This
level of approximation, where only the ray trajectories are considered, is
known as the geometric optics limit. The geometric optics result then forms
the leading order term in the WKB expansion, with the next-to-leading order
term encapsulating variations of the wave amplitude. More details about the
method can be found in e.g. [50, 51, 52].

To find the rays, we write,

$\begin{bmatrix}\phi\\\
\eta\end{bmatrix}=\begin{bmatrix}\mathcal{A}(\mathbf{x},t)\\\
\mathcal{B}(\mathbf{x},t)\end{bmatrix}e^{i\mathcal{S}(\mathbf{x},t)},$ (32)

where the amplitudes $\mathcal{A},\mathcal{B}$ vary smoothly with respect to
the phase $\mathcal{S}$. The wavevector $\mathbf{k}$ and frequency $\omega$ of
the wave are defined by,

$\omega=-\partial_{t}\mathcal{S},\qquad\mathbf{k}=\bm{\nabla}\mathcal{S}.$
(33)

Importantly, these can be functions of $\mathbf{x}$ and $t$ in inhomogeneous
media. Substituting (32) into (19) and neglecting derivatives of
$\mathcal{A}$, $\mathcal{B}$, $\mathbf{k}$ and $\omega$ (to capture the
smoothly varying assumption) we obtain a local dispersion relation (see
Appendix B for further details),

$\begin{split}\Omega^{2}=&\
\mathcal{D}(k),\qquad\Omega=\omega-\mathbf{v}\cdot\mathbf{k},\\\
\mathcal{D}(k)=&\ f(k)F(k),\qquad f(k)=g+\gamma k^{2},\end{split}$ (34)

where $\mathcal{D}(k)$ is the dispersion function and $k=|\mathbf{k}|$ is the
modulus of the wavevector. The inverse length scale $k_{c}=\sqrt{g/\gamma}$
determines the importance of surface tension relative to gravity. For room
temperature water, we find $k_{c}\simeq 370\mathrm{m^{-1}}$. Surface tension
therefore becomes important when the local wavelength is $\lesssim
1.7\mathrm{cm}$.

The rays are defined to be perpendicular to the wavefronts (i.e. surfaces of
constant phase) hence, the ray equations are obtained by extremising
$\mathcal{S}$. To ensure that (34) is satisfied, a constraint $\mathcal{H}=0$
is added to $\mathcal{S}$ to define an action for the rays,

$S=\int d\lambda\left(k_{\mu}\dot{x}^{\mu}-\xi\mathcal{H}\right)$ (35)

where $x^{\mu}=(t,\mathbf{x})$, $k_{\mu}=(-\omega,\mathbf{k})$, $\xi$ is a
Lagrange multiplier, $\lambda$ is an arbitrary parameter along the rays and
$\dot{}=d/d\lambda$. Rescaling $\lambda\to\lambda/\xi$ absorbs the Lagrange
multiplier. The effective Hamiltonian is,

$\begin{split}\mathcal{H}=&\
-\frac{1}{2}\left(\omega-\omega_{D}^{+}\right)\left(\omega-\omega_{D}^{-}\right),\\\
\omega_{D}^{\pm}=&\
\mathbf{v}\cdot\mathbf{k}\pm\sqrt{\mathcal{D}(k)},\end{split}$ (36)

where $\omega_{D}^{\pm}$ are the upper ($+$) and lower ($-$) branches of the
dispersion relation and the rays are the solutions of Hamilton’s equations,

$\dot{x}^{\mu}=\frac{\partial\mathcal{H}}{\partial
k_{\mu}},\qquad\dot{k}_{\mu}=-\frac{\partial\mathcal{H}}{\partial x^{\mu}}.$
(37)

When $\mathcal{D}(k)\sim k^{2}$, the rays coincide with geodesics of the
effective spacetime geometry seen by shallow water waves (discussed in Section
IV.1).

Figure 2: A schematic of the dispersion relation for capillary-gravity waves.
The curves $\omega_{D}^{\pm}$ in (36) at fixed $m$ and $r$ are shown as black
curves. The dispersion relation in (34) is equivalent to
$\omega=\omega_{D}^{\pm}$. Hence, for a given frequency (see the dashed line
for an example) the solutions of this equation define the different allowed
roots $p^{j}$ (red dots). For our system, there are a maximum of 6 real
(propagating) roots at a fixed $r$. As $r$ is varied, the form of the
$\omega^{\pm}_{D}(p)$ changes and the six roots will shift location. Note the
$\omega_{D}^{\pm}$ are exaggerated here to clearly illustrate the six
propagating roots. This schematic is representative of when surface tension is
smaller than a critical value $\gamma_{c}$, discussed in detail in Section V.

An important quantity in the following will be the group velocity, which we
will soon see determines the direction along which the norm is transported.
The group velocity is defined as,

$\mathbf{v}_{g}\equiv\partial_{\mathbf{k}}\omega=\partial_{\mathbf{k}}\omega_{D}^{\pm},$
(38)

where we have used the fact that the frequency satisfies
$\omega=\omega_{D}^{\pm}$. Hence, the group velocity is given by the gradient
of the branches in the $\mathbf{k}$ direction at the points where
$\omega=\omega_{D}^{\pm}$ has solutions. Two useful properties of the group
velocity are,

$\partial_{\mathbf{k}}\mathcal{H}=\mathbf{v}_{g}\Omega,\qquad\mathbf{v}_{g}=-\frac{\partial_{\mathbf{k}}\mathcal{H}}{\partial_{\omega}\mathcal{H}}=\frac{\dot{\mathbf{x}}}{\dot{t}}=\frac{\partial\mathbf{x}}{\partial
t}.$ (39)

where the latter follows from Hamilton’s equations (37) and demonstrates that
a point travelling through space along a ray parametrised by coordinate time
moves in the direction of $\mathbf{v}_{g}$ as expected.

Solving Hamilton’s equations yields the rays as function of the parameter
$\lambda$. There exists, however, a more direct approach when the system has
symmetries. In particular, (11) does not depend on $t$ or $\theta$. Working in
polar coordinates $k_{\mu}=(-\omega,p,m/r)$ and Hamilton’s equations tell us
that $\omega$ and $m$ are conserved along the rays. Then, instead of solving
$\dot{r}=\partial_{p}\mathcal{H}$ and $\dot{p}=-\partial_{r}\mathcal{H}$ for
$r(\lambda)$ and $p(\lambda)$, one can search for the solutions $p(r)$ to
$\mathcal{H}=0$, thereby solving the original dispersion relation in (34).
This amounts to searching for the intersections of $\omega_{D}^{\pm}(p)$ with
$\omega=\mathrm{const}$ at fixed $m$ and $r$ (see Fig. 2).

At each radius, $\mathcal{H}=0$ has solutions $p^{j}$ where $j=1...N$ and $N$
is an integer which depends on the functional form of $\mathcal{D}$. The
direction of travel along a ray is determined by the group velocity in the
radial direction,

$v^{j}_{g,r}\equiv\hat{\mathbf{e}}_{r}\cdot\mathbf{v}^{j}_{g}=\partial_{p}\omega^{\pm}_{D}|_{p=p^{j}},$
(40)

and one selects $\pm$ depending on which branch $p^{j}$ is on. For the
dispersion relation (34), there are a maximum of 6 propagating modes at a
given radius, $j\in\\{b,d,-,+,u,t\\}$. The $p^{\pm}$ modes are the ones
present in the low $k$ theory, whilst the others appear in pairs as a result
of dispersion: $p^{d,u}$ arise from deep water effects [25] whereas $p^{b,t}$
are due to surface tension 111The labels are a nod to quark pairs in particle
physics, i.e. up/down and top/bottom.. The modes are labelled in order of
increasing $\mathrm{Re}[p^{j}]$ so that they appear from left to right on the
dispersion relation when all 6 solutions are real, see Fig. 2. A subtlety of
this construction is that, in many cases, not all 6 modes are simultaneously
propagating at a given location. However, we will see in Section V.2 that
there is a natural, unambiguous way to label the modes.

For the flow profile (11) considered in this work, the solutions of the
dispersion relation have the following asymptotic properties (these will be
derived properly in Section V.2). Approaching the vortex axis, there are four
propagating modes: an in-going/out-going pair on $\omega_{D}^{+}$ and another
in-going/out-going pair on $\omega_{D}^{-}$. Far away from the vortex core,
there will be two propagating modes (one in-going and one out-going) which are
on $\omega_{D}^{+}$ ($\omega_{D}^{-}$) for $\omega>0$ ($\omega<0$). These are
the parts of the solution capable of forming a BHB instability, provided a
suitable condition for amplification is met.

###  III.2 Transport equation

In Appendix B, we derive the following conservation equation describing the
transport of the wave amplitude along the rays,

$\partial_{t}\left(\frac{\Omega\mathcal{A}^{2}}{f}\right)+\bm{\nabla}\cdot\left(\frac{\mathbf{v}_{g}\Omega\mathcal{A}^{2}}{f}\right)+2\nu
k^{2}\frac{\Omega\mathcal{A}^{2}}{f}=0.$ (41)

(A similar equation for $\mathcal{B}$ can be obtained using
$\mathcal{B}=i\Omega\mathcal{A}/f$). Note that the first two terms are
equivalent to the conservation equation for the norm in (30) whilst the final
term encodes dissipative losses. Since the basic flow is independent of
$(t,\theta)$, we can solve (41) exactly by writing it as,

$\partial_{r}\log\left(\frac{rv_{g,r}\Omega\mathcal{A}^{2}}{f}\right)=-\frac{2\nu
k^{2}}{v_{g,r}},$ (42)

and integrating,

$A=\frac{\alpha e^{i\varphi}}{|v_{g,r}\Omega
f^{-1}|^{1/2}}\exp\left(-\nu\int\frac{dr}{v_{g,r}}k^{2}\right),\quad\mathcal{A}=\frac{A}{\sqrt{r}},$
(43)

where $\alpha$ is a complex valued integration constant. The slowly varying
phase $\varphi(r)$ comes from the factor on the denominator whose modulus we
have taken. In shallow water $\varphi$ is a constant and can be absorbed into
$\alpha$. Note that from Hamilton’s equations (37), we have the relation
$dt=dr/v_{g,r}$. Hence, the exponent is simply the damping coefficient $-\nu
k^{2}$ integrated over the time elapsed along a given ray. The mode amplitude
therefore decays in the direction of $v_{g,r}$ as expected. In the following,
it will be easier to work with the amplitudes $A$, which are simply the full
mode amplitudes $\mathcal{A}$ with $\sqrt{r}$ term factored out. This term
encodes the fact that as a radial wave propagates towards small $r$, its
amplitude will increase by virtue of the fact that the norm density $\rho_{n}$
in (31) has been focused onto a disk of smaller radius. It is also useful to
note that in regions where the $j^{\mathrm{th}}$–mode propagates, i.e. where
$p^{j}\in\mathbb{R}$, the phase of $A$ can be entirely absorbed into the
integration constant and we can write,

$A^{j}=\alpha^{j}\left|\frac{f^{j}}{\partial_{p}\mathcal{H}^{j}}\right|^{\frac{1}{2}}\exp\left(-\nu\int\frac{dr}{v^{j}_{g,r}}{k^{j}}^{2}\right),$
(44)

where we have used (39) to write $v_{g,r}\Omega=\partial_{p}\mathcal{H}$ and
the superscript $j$ indicates that the amplitude is evaluated for a particular
$p^{j}$. In the evanescent region, the phase $\varphi(r)$ can acquire a non-
zero value, although we will argue shortly that it does not contribute to any
of the physics.

###  III.3 Turning points

The amplitude (43) clearly diverges at locations where $v_{g,r}$ vanishes. By
considering the branches of the dispersion relation in Fig. 2, it is apparent
at locations where $v_{g,r}\equiv\partial_{p}\omega_{D}^{\pm}=0$, two
neighbouring $p^{j}$ can become degenerate and the corresponding rays will
intersect. The corresponding conditions on $\mathcal{H}$ are,

$\mathcal{H}_{tp}=0,\quad\partial_{p}\mathcal{H}_{tp}=0,$ (45)

which, for fixed $(\omega,m)$ and a given set of background parameters, are
solved by the pair $(r_{tp},p_{tp})$. Hamilton’s equations (37) imply that the
points of intersection satisfy $\dot{r}=0$, i.e. they are turning points for
the rays. Physically, the ray is scattered (reflected) by the inhomogeneous
nature of the medium; a ray which is initially travelling in one direction
instantaneously comes to rest before being reflected as a different root of
the dispersion relation. During this reflection, the adiabatic assumption of
the WKB approximation (32) is violated since the local amplitude (43) varies
over a much shorter scale than the local wavenumber $p^{j}(r)$.

There exists a standard procedure to fix this pathological behaviour. At each
location where $v_{g,r}=0$, one can expand the equations of motion (19) in the
vicinity of the turning point and find an exact solution for the locally
degenerate roots. This solution can then be mapped onto the WKB solution
either side of the turning point (see Appendix C for a detailed description of
this procedure). The result is a matrix equation relating the WKB amplitudes
either side of the turning point, which we now state.

Let the intersecting rays be denoted $p^{\mathrm{R}}$ and $p^{\mathrm{L}}$,
such that $p^{\mathrm{R}}>p^{\mathrm{L}}$ away from the turning point
$r_{tp}$. The labels $\mathrm{R},\mathrm{L}$ are placeholders for a
neighbouring pair of $j\in\\{b,d,-,+,u,t\\}$ modes. On the opposite side of
$r_{tp}$, the two modes will be a complex conjugate pair which we call
$p^{\downarrow}$ (the mode which decays in the direction of increasing $r$)
and $p^{\uparrow}$ (the mode which grows in the direction of increasing $r$).
These satisfy $\mathrm{Re}[p^{\uparrow}]=\mathrm{Re}[p^{\downarrow}]$ and
$\mathrm{Im}[p^{\uparrow}]=-\mathrm{Im}[p^{\downarrow}]$. When the pair are
propagating for $r<r_{tp}$ and evanescent for $r>r_{tp}$, the local amplitudes
in (43) are related by,

$\begin{pmatrix}A^{\mathrm{R}}_{tp}\\\
A^{\mathrm{L}}_{tp}\end{pmatrix}=e^{\frac{i\pi}{4}}\begin{pmatrix}1&-\frac{i}{2}\\\
-i&\frac{1}{2}\end{pmatrix}\begin{pmatrix}A^{\downarrow}_{tp}\\\
A^{\uparrow}_{tp}\end{pmatrix},$ (46)

whilst in the reverse scenario (evanescent for $r<r_{tp}$ and propagating for
$r>r_{tp}$) we have,

$\begin{pmatrix}A^{\uparrow}_{tp}\\\
A^{\downarrow}_{tp}\end{pmatrix}=e^{\frac{i\pi}{4}}\begin{pmatrix}\frac{1}{2}&-\frac{i}{2}\\\
-i&1\end{pmatrix}\begin{pmatrix}A^{\mathrm{R}}_{tp}\\\
A^{\mathrm{L}}_{tp}\end{pmatrix}.$ (47)

Each matrix equation is defined so that it takes mode amplitudes on the side
of large $r$ and returns the same on the side of small $r$.

An important point is that whilst the prefactor in (43) diverges at the
turning points, the integral in the exponent containing a factor of
$1/v_{g,r}$ does not. This point is worth mentioning since the dissipation
integral contributes to our expressions for the growth/decay rates, hence it
is important to check that it is well-behaved. This can be seen directly by
expanding the effective Hamiltonian around the turning point,
$\mathcal{H}=\partial_{r}\mathcal{H}(r-r_{tp})+\frac{1}{2}\partial_{p}^{2}\mathcal{H}_{tp}(p-p_{tp})^{2}$,
and using (39) to write
$v_{g,r}=\Omega_{tp}^{-1}\partial_{p}^{2}\mathcal{H}_{tp}(p-p_{tp})=\Omega_{tp}^{-1}[-2\partial_{r}\mathcal{H}_{tp}\partial_{p}^{2}\mathcal{H}_{tp}(r-r_{tp})]^{1/2}$,
where the last equality follows from $\mathcal{H}=0$. We can then evaluate the
integral proportional to $\nu$ in (43) in a small region around the turning
point,

$\begin{split}&\frac{\nu
k_{tp}^{2}\Omega_{tp}}{\sqrt{-2\partial_{r}\mathcal{H}_{tp}\partial_{p}^{2}\mathcal{H}_{tp}}}\int^{r_{tp}+\epsilon}_{r_{tp}}\frac{dr}{\sqrt{r-r_{tp}}}\\\
&\qquad\qquad\qquad\qquad=\frac{\nu
k_{tp}^{2}\Omega_{tp}}{\sqrt{-\partial_{r}\mathcal{H}_{tp}\partial_{p}^{2}\mathcal{H}_{tp}}}\sqrt{2\epsilon},\end{split}$
(48)

which is manifestly regular as $\epsilon\to 0$.

In general, there can be multiple turning points over the range in $r$, and to
calculate the asymptotic scattering coefficients, one should take into account
the reflections from each. Our procedure to compute these amplitudes will
follow that outlined in [25, 26, 33] in that we only account for reflections
at turning points, which give the dominant contributions to scattering. Sub-
dominant contributions can be obtained e.g. around saddle points [54, 45] and
complex turning points [8].

There is a convenient method to identify the locations of the turning points
which we will make use of in the coming sections. This method is equivalent to
solving the conditions (45), although it has a more natural interpretation in
terms of the dispersion relation. For fixed $m$ and $r$ (and a given
background flow) one can identify the extrema $p_{\mathrm{ex}}$ of
$\omega_{D}^{\pm}$ by solving
$\partial_{p}\omega_{D}^{\pm}|_{p=p_{\mathrm{ex}}}$=0. Since there are a
maximum of three per branch (see Fig. 2), we will label these
$p^{a}_{\mathrm{ex}}$ with $a=1,2,3$, noting that the extrema on the lower
branch are located at equal and opposite $p$ values to those on the upper
branch. We label the extrema by their proximity to the $p=0$ axis, with
$p_{\mathrm{ex}}^{1}$ ($-p_{\mathrm{ex}}^{1}$) being the closest on the upper
(lower) branch. By solving for all $r$, we can promote these solutions to
functions of the radius, i.e. $p_{\mathrm{ex}}^{a}(r)$. Substituting back into
the appropriate branch of the dispersion relation, we obtain two sets of three
curves,

$\omega^{+}_{a}=\omega^{+}_{a}(r,p_{\mathrm{ex}}^{a}(r)),\qquad\omega^{-}_{a}=\omega^{-}_{a}(r,-p_{\mathrm{ex}}^{a}(r)),$
(49)

which give the $\omega$ values at the extrema of the dispersion relation. The
utility of these curves is that once we select a particular frequency, we can
solve,

$\omega=\omega^{\pm}_{a}(r_{tp}),$ (50)

to find all the turning points. We will see a simple example of this in the
next section on shallow water waves. However, it will be particularly useful
in the analysis of deep water waves, where there is no simple analytic
expression for the $r_{tp}$.

At each turning point, the amplitudes are related by (46) or (47) depending on
the orientation. Between turning points, we assume that the different $p^{j}$
do not interact with each other and evolve adiabatically according to the WKB
approximation. Let us write the full solution for velocity potential
perturbations $\phi$ at an initial location $r=r_{b}$ in the form,

$\begin{split}\phi_{b}=&\ \frac{e^{im\theta-i\omega t}}{\sqrt{r}}R_{b},\\\
R_{b}=&\
\sum_{j}A_{b}^{j}(r)e^{i\int_{\mathrm{ref}}p^{j}_{b}dr}\equiv\mathbf{P}^{\mathrm{ref}}_{b}\cdot\mathbf{A}_{b}\end{split}$
(51)

where $\phi_{b}=\phi(r_{b})$ and similarly for the other quantities. In the
second equality for $R_{b}$, $\mathbf{P}^{\mathrm{ref}}_{b}$ is a row vector
contain the reference phases of the $j$ modes at $r_{b}$ and $\mathbf{A}_{b}$
is a column vector containing the mode amplitudes. The column vector is
ordered so that modes with largest $\mathrm{Re}[p^{j}]$ appear toward the top.
We can construct the full solution at a new point $r=r_{a}$, provided the WKB
approximation is valid everywhere for $r_{a}\leq r\leq r_{b}$, through the
relation $R_{a}=\mathbf{P}^{\mathrm{ref}}_{b}\cdot\mathbf{A}_{a}$, where the
amplitudes are related by,

$\mathbf{A}_{a}=\mathrm{diag(\mathcal{F}^{j}_{ab})}\mathbf{A}_{b}.$ (52)

The factors $\mathcal{F}^{j}_{ab}$ shift the phase and amplitude of the mode
adiabatically according to (43) and are given by,

$\begin{split}\mathcal{F}^{j}_{ab}=&\
\left|\frac{f^{j}_{a}\partial_{p}\mathcal{H}^{j}_{b}}{f^{j}_{b}\partial_{p}\mathcal{H}^{j}_{a}}\right|^{\frac{1}{2}}e^{i\left(\varphi^{j}_{a}-\varphi^{j}_{b}\right)}\exp\left(-i\int^{r_{b}}_{r_{a}}p^{j}dr\right)\times\\\
&\qquad\qquad\qquad\qquad\exp\left(\nu\int^{r_{b}}_{r_{a}}\frac{dr}{v^{j}_{g,r}}{k^{j}}^{2}\right).\end{split}$
(53)

The phases $\varphi^{j}$ are defined so that they go to zero on the turning
points and thus, when $r_{a}$ and $r_{b}$ are turning points, the factor
$\exp[i(\varphi^{j}_{a}-\varphi^{j}_{b})]$ drops out of the expression above.
By (44), the $\varphi^{j}$ are also zero in regions where
$p^{j}\in\mathbb{R}$. In this way, all information about the phase shifts at
the turning points are contained in the integration constants $\alpha$ in
(43). It is worth noting that one case where the $\varphi^{j}$ vanish
identically is in the relativistic (non-dispersive) regime discussed in
Section IV.

In general, we can describe the asymptotic scattering properties of the system
through the transfer matrix $\mathcal{M}$,

$\mathbf{A}_{0}=\mathcal{M}\mathbf{A}_{\mathrm{as}},$ (54)

where $r_{0}$ ($r_{\mathrm{as}}$) is a point such that there are no turning
points for $r<r_{0}$ ($r>r_{\mathrm{as}}$), ensuring that the scattering
amplitudes computed at $r_{0}$ ($r_{\mathrm{as}}$) will be invariant moving
toward smaller (larger) $r$. The matrix $\mathcal{M}$ is built from
combinations of the matrices in (46), (47) and (52). Its specific structure is
determined by the precise number and location of turning points for a given
$\omega,m$–mode. We will see a specific example of how to construct
$\mathcal{M}$ later on in Section V.2. The key point is that knowledge of the
turning points for a given set of parameters uniquely determines
$\mathcal{M}$, from which we can infer the scattering amplitudes. Hence, to
completely characterise the scattering of waves in our draining water tank, we
should find all radial turning points $(r_{tp},p_{tp})$ by solving (45) over
the whole parameter space of the equations in (19). This procedure was
previously carried out for shallow water ($\mathcal{D}=ghk^{2}$) [26], weakly
dispersive ($\mathcal{D}=ghk^{2}-gh^{3}k^{2}/3$) [45] and deep water
($\mathcal{D}=gk$) [25] gravity waves, where the focus was on unbounded
systems with $\nu=\gamma=0$. The extension in this work to the case of bounded
systems, with surface tension and dissipation properly accounted for, will
enable a computation of unstable BHB growth rates.

###  III.4 Local reflection coefficient

An important case is when there are two turning points (say $r_{a}<r_{b}$)
such that interacting modes are evanescent for $r_{a}<r<r_{b}$. In this case,
we can combine the matrices (46) and (47) with the shift factors (52) to
write,

$\begin{split}\begin{pmatrix}A^{\mathrm{R}}_{a}\\\
A^{\mathrm{L}}_{a}\end{pmatrix}=&\
\mathcal{N}_{ab}\begin{pmatrix}A^{\mathrm{R}}_{b}\\\
A^{\mathrm{L}}_{b}\end{pmatrix},\\\ \mathcal{N}_{ab}=&\
\mathcal{F}^{\downarrow}_{ab}\begin{bmatrix}1+f_{ab}^{2}/4&i(1-f_{ab}^{2}/4)\\\
-i(1-f_{ab}^{2}/4)&1+f_{ab}^{2}/4\end{bmatrix},\end{split}$ (55)

where $\mathrm{Re}[p^{\mathrm{R}}]>\mathrm{Re}[p^{\mathrm{L}}]$ and we have
defined,

$\begin{split}f_{ab}=&\
\exp\left(-\int^{r_{b}}_{r_{a}}\mathrm{Im}[p^{\downarrow}]dr\right)\times\\\
&\qquad\qquad\quad\exp\left(-i\nu\int^{r_{b}}_{r_{a}}\mathrm{Im}\left[\frac{{k^{\downarrow}}^{2}}{v_{g,r}^{\downarrow}}\right]dr\right).\end{split}$
(56)

We will be particularly interested in the scattering of the long wavelength
modes $p^{\pm}$. For large $r$ and $\omega>0$, $p^{+}$ ($p^{-}$) is out-going
(in-going) since $v_{g,r}^{+}>0$ ($v_{g,r}^{-}<0$) as $r\to r_{\mathrm{as}}$.
At small $r$, the same is true if $p^{\pm}$ are on the upper branch of the
dispersion relation $\omega_{D}^{+}$ as $r\to r_{0}$. However, for certain
parameters ranges (to be deduced later) $p^{\pm}$ are solutions of
$\omega=\omega_{D}^{-}$, in which case $p^{+}$ ($p^{-}$) is the in-going (out-
going) mode as seen from $v_{g,r}^{+}<0$ ($v_{g,r}^{-}>0$) as $r\to r_{0}$.

When $v^{\mathrm{R}}_{g,r}(r>r_{b})>0$ and $v^{\mathrm{L}}_{g,r}(r>r_{b})<0$,
the local reflection coefficient
$\mathcal{R}_{b}=A^{\mathrm{R}}_{b}/A^{\mathrm{L}}_{b}$ is given by,

$\mathcal{R}_{b}=-i\left(\frac{1-f_{ab}^{2}/4}{1+f_{ab}^{2}/4}\right)^{\pm
1},$ (57)

where the $\pm$ corresponds to the sign of $v_{g,r}^{\mathrm{R}}(r<r_{a})$.
Note that for $\nu=0$, the phase of $\mathcal{R}_{b}$ is entirely determined
by the prefactor $-i$. Treating $\nu$ perturbatively, we can expand $f_{ab}$
to show that,

$\mathcal{R}_{b}=\mathcal{R}_{b}^{\nu=0}\left[1\pm\mathcal{O}\left(\nu|f_{ab}|^{2}\right)\right],$
(58)

where we have assumed that $f_{ab}$ is itself small, which is true in our case
when $r_{a,b}$ are well separated. Later, we present results to leading order
in $\nu$ and $|f_{ab}|^{2}$, hence, reflection coefficients can be evaluated
for $\nu=0$.

###  III.5 Energy conservation

The expression for the reflection coefficient in (57) can exceed unity in
certain cases, implying amplification of scattered waves. This is precisely
the phenomenon of rotational superradiance. We can understand the origin of
this amplification using energy considerations in the limit $\nu=0$. First, we
note that in stationary backgrounds (which is our case) the conserved energy
current is related to the norm current in (30) by a factor of $\omega$. For
WKB modes, we show in Appendix B that the components of the norm current for
each WKB mode (labelled by $j$) is,

$\displaystyle\rho^{j}_{n}=2\frac{\Omega^{j}}{f^{j}}|A^{j}|^{2},\qquad\mathbf{J}_{n}^{j}=\mathbf{v}^{j}_{g}\rho^{j}_{n}.$
(59)

The second equality implies that the norm is carried in the direction of the
group velocity. The sign of the norm is determined by the comoving frequency
$\Omega^{j}$, which is positive (negative) on $\omega_{D}^{+}$ and
$\omega_{D}^{-}$. Therefore, modes on $\omega_{D}^{-}$ have negative norm, and
if they also have $\omega>0$ they have negative energy densities
$\omega\rho_{n}$. The region where this occurs is the dispersive analogue of
the ergosphere around rotating black holes (discussed in Section IV). In
particular, its presence implies that there is an ingoing mode
($v_{g,r}^{j}<0$) with negative energy density, which contributes positively
to the outgoing energy current. This is ultimately the phenomenon responsible
for superradiance.

To see this explicitly, we can rewrite (30) as a statement of energy
conservation relating the various WKB modes evaluated near $r=0$ and $r\to
r_{\mathrm{as}}$,

$\sum_{j=1}^{4}\frac{\mathcal{H}^{\prime}(p^{j}_{0})}{f(p^{j}_{0})}|A^{j}_{0}|^{2}=\sum_{j=1}^{2}\frac{\mathcal{H}^{\prime}(p^{j}_{\mathrm{as}})}{f(p^{j}_{\mathrm{as}})}|A^{j}_{\mathrm{as}}|^{2},$
(60)

where we have used the fact that there are 4 propagating modes approaching the
vortex axis and 2 propagating modes far from the vortex, as discussed at the
end of Section III.1. Now, let us consider $\omega>0$ so that the modes far
from the vortex are on $\omega_{D}^{+}$. Using the expression for the
amplitudes in (44) we can write the equation above as,

$|\mathcal{R}|^{2}-1=\left|\frac{\alpha_{0}^{\mathrm{out}+}}{\alpha_{\mathrm{as}}^{\mathrm{in}+}}\right|^{2}-\left|\frac{\alpha_{0}^{\mathrm{in}+}}{\alpha_{\mathrm{as}}^{\mathrm{in+}}}\right|^{2}-\left|\frac{\alpha_{0}^{\mathrm{out}-}}{\alpha_{\mathrm{as}}^{\mathrm{in}+}}\right|^{2}+\left|\frac{\alpha_{0}^{\mathrm{in}-}}{\alpha_{\mathrm{as}}^{\mathrm{in+}}}\right|^{2},$
(61)

where we use the suffixes $\mathrm{in}\pm$ and $\mathrm{out}\pm$ to label the
in-going and out-going modes on $\omega_{D}^{\pm}$ respectively (note that
these labels are placeholders for the mode labels $j\in\\{b,d,-,+,u,t\\}$
which we will eventually solve for by drawing the scattering diagrams). We
have also introduced the asymptotic reflection coefficient
$\mathcal{R}=\alpha^{\mathrm{out}+}_{\mathrm{as}}/\alpha^{\mathrm{in}+}_{\mathrm{as}}$,
whose modulus coincides with the one introduced in (57). Our boundary
conditions in the vortex core (discussed in detail in Appendix D) amount to
assuming the vortex does not emit. The condition above then becomes,

$|\mathcal{R}|^{2}-1=-\left|\frac{\alpha_{0}^{\mathrm{in}+}}{\alpha_{\mathrm{as}}^{\mathrm{in+}}}\right|^{2}+\left|\frac{\alpha_{0}^{\mathrm{in}-}}{\alpha_{\mathrm{as}}^{\mathrm{in+}}}\right|^{2}.$
(62)

Depending on the background and wave parameters, the constraint that the
vortex does not emit fixes one (but not both) of the terms on the right-hand-
side to be zero. The condition on the parameters can then be expressed as a
condition on $\omega$ such that when (3) is satisfied, $|\mathcal{R}|>1$ and
we have superradiance. When this occurs, (62) tells us that the vortex absorbs
a negative energy, and energy conservation results in the amplification of the
wave reflected to infinity. It will be the challenge of Section V to identify
the characteristic frequency $\omega_{\star}$ in (3) that governs when this
absorption of negative energy occurs.

##  IV Shallow water waves

So far, we have established a set of equations governing capillary-gravity
waves around our draining vortex and presented a general method capable of
handling effects of dispersion and dissipation. In this section, we will apply
our method in the shallow water regime where the analogy with fields around a
black hole is exact. In this way, we demonstrate the manner in which the
instability of a draining vortex in a cylindrical container is equivalent to
the BHB instability, before generalising this effect to situations where
dispersion and dissipation are important. The key results will be presented in
Figs. 5 and 6, where we validate our procedure against standard numerical
methods available to us in this regime.

In the low $k$ limit of the theory, $k^{-1}$ is assumed much larger than all
other length scales in the problem (except $r_{B}$). In this regime (16)
becomes $F(k)\simeq hk^{2}$, and the effects of surface tension and viscosity
can be neglected. The effective Hamiltonian (36) becomes,

$\mathcal{H}=-\frac{1}{2}(\omega-\mathbf{v}\cdot\mathbf{k})^{2}+\frac{1}{2}c^{2}k^{2},\qquad
c=\sqrt{gh},$ (63)

where $c$ is the shallow water wave speed. The corresponding wave equation for
velocity potential fluctuations is obtained from (19),

$(\partial_{t}+\bm{\nabla}\cdot\mathbf{v})(\partial_{t}+\mathbf{v}\cdot\bm{\nabla})\phi-\bm{\nabla}\cdot(c^{2}\bm{\nabla}\phi)=0.$
(64)

The magnitude of the group velocity in a frame comoving with the background
fluid flow is $|\partial_{\mathbf{k}}\Omega|=c$, that is, all frequencies
propagate at the same speed $c$; the system is non-dispersive. This greatly
simplifies the analysis of (64) and provides a good starting point to
illustrate the physics involved in the BHB instability.

###  IV.1 Analogy to black holes

In this subsection, we review the standard analysis of wave scattering in the
analogue metric (e.g. [1, 20]). It is included for completeness so that we may
(a) map the system at low $k$ to an analogue black hole and (b) identify the
key features of the instability that should be present for the instability to
be called a BHB. This is an important point since there is no strict notion of
a horizon in the dispersive regime, however, there is still an instability
which coincides with the BHB in the shallow water limit.

The analogy to the dynamics of fields in curved spacetime arises due to the
equivalence of (64) with the Klein-Gordon equation,

$\frac{1}{\sqrt{-g}}\partial_{\mu}\left(\sqrt{-g}g^{\mu\nu}\partial_{\nu}\phi\right)=0,$
(65) 222Note, the $g$ here is the determinant of the metric, not the
gravitational acceleration used in (63) and elsewhere, which is conventionally
denoted with the same letter.

where $\phi$ is a massless scalar field and the effective metric tensor
$g_{\mu\nu}$ is given by,

$\begin{split}&g_{\mu\nu}dx^{\mu}dx^{\nu}\\\
&\quad=c^{2}\left[-c^{2}dt^{2}+\left(dr-v_{r}dt\right)^{2}+\left(rd\theta-
v_{\theta}dt\right)^{2}\right].\end{split}$ (66)

Note that when taking $h\simeq h_{0}$ (as we do in this work) the pre-factor
$c^{2}$ is a constant and can be dropped. Since (65) is invariant under
Lorentz transformations, we have an effective relativistic theory, with the
shallow water wave speed $c$ playing the role of the speed of light. Working
with the basic flow field in (11), this metric has a horizon $r_{h}$ and an
ergosphere $r_{e}$ at the locations,

$r_{h}=\frac{D}{c},\qquad r_{e}=\frac{\sqrt{C^{2}+D^{2}}}{c},$ (67)

where horizon is the boundary of the region where the radial velocity exceeds
the wave speed, $|v_{r}(r_{h})|=c$, and the ergosphere is the boundary of the
region where the magnitude of the total velocity field exceeds the wave speed,
$|\mathbf{v}(r_{e})|=c$. Therefore, inside $r<r_{h}$ it is not possible for a
wave satisfying (64) to move in the direction of increasing $r$. Similarly,
for $r<r_{e}$ waves satisfying (64) cannot move against the rotation of the
vortex. Inside the ergosphere, shallow water waves are permitted to have
negative energy densities $\omega\Omega^{j}|\mathcal{A}^{j}|^{2}$. In
particular if $\Omega<0$, the positive frequency mode crossing $r_{h}$ carries
a negative energy in that direction, which results in amplification of the
wave escaping to infinity.

It is possible to define the following dimensionless quantities,

$\frac{r}{r_{h}}\to r,\qquad\frac{ct}{r_{h}}\to t,\qquad\frac{C}{D}\to C,$
(68)

which amounts to setting $c=D=1$. In the text, we will work in these units to
lighten the notation whilst in all figures, quantities will be dimensionful.
This scaling will facilitate an efficient analysis of the parameter space
accessible to shallow water waves around a draining vortex.

To study the scattering of waves described by (64), the conventional method
involves first transforming the field,

$\phi=\sum_{\omega m}e^{im\theta-i\omega
t}\exp\left(i\int\frac{\tilde{\omega}(r)dr}{r\Delta(r)}\right)\frac{R(r)}{\sqrt{r}},$
(69)

where the radial functions $R(r)$ are specific to each $\omega,m$ mode and we
have defined,

$\tilde{\omega}=\omega-\frac{mC}{r^{2}},\qquad\Delta(r)=1-\frac{1}{r^{2}}.$
(70)

Next, the tortoise coordinate $r_{*}$ is defined as,

$r_{*}=r+\frac{1}{2}\log\left|\frac{r-1}{r+1}\right|,$ (71)

which maps $r\in(1,\infty)$ on to $r_{*}=(-\infty,\infty)$, i.e. the horizon
is pushed to negative infinity. The wave equation (64) then becomes a one-
dimensional scattering problem in the radial direction

$-\partial^{2}_{r_{*}}R+V(r)R=0,$ (72)

where the effective potential is,

$V(r)=-\tilde{\omega}^{2}+\Delta\left(\frac{m^{2}-1/4}{r^{2}}+\frac{5}{4r^{4}}\right).$
(73)

(72) has asymptotic solutions,

$R(r)\sim\begin{cases}A^{-}_{h}e^{-i\tilde{\omega}_{h}r_{*}},\qquad\quad\quad\
\qquad r_{*}\to-\infty\\\ A^{+}_{\infty}e^{i\omega
r_{*}}+A^{-}_{\infty}e^{-i\omega r_{*}},\quad\ \ r_{*}\to\infty\end{cases}$
(74)

where we have imposed the boundary condition that waves are purely ingoing at
the horizon and the subscript $h$ implies that $\tilde{\omega}$ is evaluated
at $r=1$ in our dimensionless units. Using the conservation of the Wronskian
$W=R^{*}\partial_{r_{*}}R-R\partial_{r_{*}}R^{*}$, i.e. $\partial_{r_{*}}W=0$
which is equivalent to conservation of the norm current $\mathbf{J}_{n}$ in
(31), the amplitudes can be shown to obey,

$|A_{\infty}^{+}|^{2}+\frac{\tilde{\omega}_{h}}{\omega}|A_{h}^{-}|^{2}=|A_{\infty}^{-}|^{2}.$
(75)

Hence positive frequency modes with $\omega<mC$ will be superradiantly
amplified, i.e. $|A^{+}_{\infty}|>|A^{-}_{\infty}|$. This relation is
equivalent to the one in (62) as we will discuss below.

When the system is enclosed by an impenetrable boundary, superradiance can
cause the vortex to become unstable: amplified waves are reflected from the
boundary back into the vortex where they are further amplified, leading to a
runaway process. We will now show how the WKB method can be used to predict
the frequencies of unstable modes formed in this way, before verifying the
method’s accuracy with by comparing with exact numerical solutions of (64).

###  IV.2 WKB method

The calculation of the unstable mode frequencies proceeds in three steps.
Firstly, we find the radial momentum $p$ for the two shallow water modes and
determine which of these represents out/in-going waves. Next, we determine the
reflection coefficient $\mathcal{R}$ for a wave incident on the vortex.
Finally, we use $p$ and $\mathcal{R}$ to deduce a resonance condition for
waves which are scattered between the vortex and the boundary of the system.
We assume here that $m\geq 0$ and allow $\omega$ to take on positive and
negative values. The results for $m<0$ can then be inferred via the symmetry
of the problem under $(\omega,m,p)\to(-\omega,-m,-p)$.

####  IV.2.1 Scattering

In polar coordinates, the shallow water effective Hamiltonian (63) becomes,

$\mathcal{H}=-\frac{1}{2}\left(\omega-\frac{mC}{r^{2}}+\frac{p}{r}\right)^{2}+\frac{1}{2}\left(p^{2}+\frac{m^{2}}{r^{2}}\right),$
(76)

where we have used the rescalings in (68) to set
$\mathcal{H}r_{h}^{2}/c^{2}\to\mathcal{H}$. Since this is quadratic in $p$,
$\mathcal{H}=0$ only has two solutions corresponding the modes labelled
$p^{\pm}$ on Fig. 2. These can be expressed concisely as,

$p^{\pm}=\frac{\tilde{\omega}/r\pm\sqrt{-V}}{\Delta},\qquad
V=-\tilde{\omega}^{2}+\Delta\frac{m^{2}}{r^{2}}.$ (77)

Note that the $V$ defined here is equivalent to that in (73) in the limit
$|m|\gg 1$, which is a formal requirement for the WKB approximation to hold
(although good agreement with more precise methods is still obtained for the
lowest $m$ modes as we shall soon see). With this identification, (51) and
(69) become equivalent for a single $\omega,m$ mode upon noting that the WKB
approximation applied to (72) gives the phase factor as $\exp(\pm
i\int\sqrt{-V}dr_{*})$, where $dr_{*}=dr/\Delta$.

The expression for the amplitude in (44) is easily calculated for
$\gamma=\nu=0$ using $\partial_{p}\mathcal{H}^{\pm}=\pm\sqrt{-V}$ which gives,

$A^{\pm}=\alpha^{\pm}|V|^{-1/4}.$ (78)

This diverges for $V=0$, hence the turning points $r_{tp}$ defined by (45) are
given by the zeros of the potential. Note that, in shallow water, the
adiabatic amplitude change in (78) for the two modes $p^{\pm}$ is the same,
something which won’t be true for dispersive systems in general. We also
mention that the expression for $p^{\pm}$ in (77) manifestly changes between
real and complex when $V$ changes sign.

Figure 3: The curves $\omega^{\pm}$ in (79) are illustrated for two different
cases. In shallow water, the curves terminate at $r_{h}$ since
$\omega_{D}^{\pm}$ have no extrema in the $p$ direction for $r<r_{h}$. Turning
points are solutions of $\omega=\omega^{\pm}(r_{tp})$ and modes are evanescent
for $\omega^{-}<\omega<\omega^{+}$. We display the $r_{tp}$ for an example
frequency $\omega_{0}$ in panel (a). The effect of increasing $C$ is to
increase the asymmetric skew between $\omega^{+}$ and $\omega^{-}$ through the
$mC/r^{2}$ term, whilst increasing $|m|$ results in an overall scaling by a
factor of $|m|$. The light ring frequencies $\omega_{lr}^{\pm}$ are the
extrema of $\omega^{\pm}(r)$.

To identify the turning points, we can write $V$ in the form,

$V=-(\omega-\omega^{+})(\omega-\omega^{-}),\qquad\omega^{\pm}=\frac{mC}{r^{2}}\pm\sqrt{\Delta\frac{m^{2}}{r^{2}}},$
(79)

where $\omega^{\pm}$ are defined in (49) (the label $a$ is redundant here
since in shallow water, $\omega_{D}^{\pm}$ each possess a maximum of one real
extremum). These curves are illustrated in Fig. 3. In shallow water, (50)
admits closed form solutions for the two turning points,

$r^{\pm}_{tp}=\left[\frac{2m(C^{2}+1)}{2\omega C+m\mp
2\sqrt{(\omega-\omega_{lr}^{+})(\omega_{lr}^{-}-\omega)}}\right]^{1/2},$ (80)

where the $\pm$ appearing on $r^{\pm}_{tp}$ labels the location of the turning
point such that $r^{-}_{tp}<r^{+}_{tp}$ and not the branch of the dispersion
relation. The extrema of $\omega^{\pm}$ in the $r$ direction are given by,

$\begin{split}r_{lr}^{\pm}=&\ \sqrt{\frac{2(C^{2}+1)}{C^{2}+1\pm
C\sqrt{C^{2}+1}}},\\\ \omega^{\pm}_{lr}=&\
\frac{m(C\pm\sqrt{C^{2}+1})}{2},\end{split}$ (81)

i.e. $\partial_{r}\omega^{\pm}|_{r=r^{\pm}_{lr}}=0$ and
$\omega_{lr}^{\pm}=\omega^{\pm}(r_{lr}^{\pm})$, where the subscript stands for
light-ring to make the connection with the analogous concept around black
holes [56]. It is manifest from the expression in (80) that the $r_{tp}^{\pm}$
are only real over the range $\omega_{lr}^{-}<\omega<\omega_{lr}^{+}$. Using
the relation $2m(C^{2}+1)={r_{lr}^{\pm}}^{2}(2\omega_{lr}^{\pm}C+m)$, the
expression also makes it apparent that
$r_{tp}^{\pm}(\omega=\omega_{lr}^{+})=r_{lr}^{+}$ and
$r_{tp}^{\pm}(\omega=\omega_{lr}^{-})=r_{lr}^{-}$.

We clarify here the use of $\pm$ in the expressions above since there are
three different meanings. Firstly, in $\omega^{\pm}_{lr}$ and $r_{lr}^{\pm}$,
the $\pm$ refers to the upper and lower branches of the dispersion relation
$\omega_{D}^{\pm}$. Secondly, in $p^{\pm}$, the $\pm$ is just a label for the
modes. For $\omega>\omega^{+}$, both $p^{\pm}$ are on the upper branch
$\omega_{D}^{+}$ whilst for $\omega<\omega^{-}$ they are both on the lower
branch $\omega_{D}^{-}$. Thirdly, in $r_{tp}^{\pm}$, the $\pm$ denotes
ordering of turning points. Depending on the value of $\omega$, $r_{tp}^{\pm}$
are either both on $\omega^{+}$, both on $\omega^{-}$, or $r_{tp}^{-}$ is on
$\omega^{-}$ and $r_{tp}^{+}$ is on $\omega^{+}$. The latter case is when
superradiance occurs. These definitions are readily deducible by inspecting
the left panel of Fig. 3.

In terms of the potential, we can use (79) to show that
$\partial_{r}V_{r=r_{lr}}=0$ and $V(r_{lr},\omega_{lr})=0$. In other words,
$r_{lr}$ is the location of the peak in the potential and $\omega_{lr}$ is the
frequency required for an $m$–mode to orbit the vortex on a circular
trajectory at $r=r_{lr}$. This last fact can be deduced from the condition on
the effective Hamiltonian (36) at the light-ring,

$\mathcal{H}_{lr}=0,\quad\partial_{p}\mathcal{H}_{lr}=0,\quad\partial_{r}\mathcal{H}_{lr}=0,$
(82)

which, using Hamilton’s equations (37), implies $\dot{r}=0$ and $\dot{p}=0$.
That is, the ray which solves these equations has vanishing velocity and
acceleration in the radial direction and therefore remains at $r=r_{lr}$
indefinitely. For our purpose, the light-ring is an important notion in
determining the frequency ranges in which there are turning points (and
therefore scattering within our approximation) and those without.
Specifically, for $\omega_{lr}^{-}<\omega<\omega_{lr}^{+}$ there is a coupling
between the two $p^{\pm}$ modes and in-going modes reflect will off the
vortex.

The added utility of the curves $\omega^{\pm}(r)$ are that they allow us to
infer the local energy density of the solutions. A ray with frequency $\omega$
at a particular location $r$ lies on the upper (lower) branch of the
dispersion relation $\omega_{D}^{+}$ ($\omega_{D}^{-}$) if the point
$r,\omega$ is above $\omega^{+}$ (below $\omega^{-}$). Using (59), one can
infer that the sign of $\Omega$, and hence the sign of the WKB norm $\rho_{n}$
is positive (negative) on $\omega_{D}^{+}$ ($\omega_{D}^{-}$). In stationary
systems, the local energy density is simply $\omega\rho_{n}$ and thus,
positive frequency modes on $\omega_{D}^{+}$ ($\omega_{D}^{-}$) have a
positive (negative) energy flux in the direction of $v_{g,r}$ (see (59)). By
inspecting Fig. 3 and noting that $\omega^{\pm}(r\to\infty)\to 0$, we see that
the condition for locally negative energy densities is fulfilled for modes at
low $r$ in the frequency range $0<\omega<\omega^{\pm}(r_{h})=mC/r_{h}^{2}$).
When the in-going wave at the horizon has a negative energy flux, conservation
of energy means the wave reflected to infinity must have greater energy than
the incident wave. Hence, we recover the criterion for superradiant
amplification $\omega<mC/r_{h}^{2}$ (for positive frequency waves).

####  IV.2.2 Reflection coefficient

The WKB approximation allows us to find a semi-analytic expression for the
reflection coefficient and therefore quantify the amount of amplification that
occurs. For this, we first need to know which of the modes $p^{\pm}$ is in-
going on the horizon. Evaluating (40) for the two modes, we find,

$v_{g,r}^{\pm}=\frac{\pm\Delta\sqrt{-V}}{\tilde{\omega}\pm\sqrt{-V}/r},$ (83)

which approaching the horizon from outside (i.e. $\Delta\to 0^{+}$) becomes,

$\begin{split}v^{+}_{g,r}&\ \overset{\Delta\to
0^{+}}{\sim}\begin{cases}\Delta/2,\qquad\ \ \,\tilde{\omega}_{h}>0,\\\
-\frac{2\tilde{\omega}_{h}^{2}}{\tilde{\omega}_{h}^{2}+m^{2}},\quad\tilde{\omega}_{h}<0,\end{cases}\\\
v^{-}_{g,r}&\ \overset{\Delta\to
0^{+}}{\sim}\begin{cases}-\frac{2\tilde{\omega}_{h}^{2}}{\tilde{\omega}_{h}^{2}+m^{2}},\quad\tilde{\omega}_{h}>0,\\\
\Delta/2,\qquad\ \ \,\tilde{\omega}_{h}<0.\\\ \end{cases}\end{split}$ (84)

For $\tilde{\omega}_{h}>0$, we see that $p^{+}$ is out-going and $p^{-}$ is
in-going just outside the horizon whilst for $\tilde{\omega}_{h}<0$, $p^{-}$
is out-going and $p^{+}$ is in-going. Consistency with the horizon boundary
condition in (74) then tells us that the amplitude of the out-going mode
should be set to zero. Using the expression in (57) for the reflection
coefficient, we obtain,

$\begin{split}\mathcal{R}=&\
-i\left(\frac{1-e^{-2S_{-+}}/4}{1+e^{-2S_{-+}}/4}\right)^{\mathrm{sgn}(\tilde{\omega}_{h})},\\\
S_{-+}=&\
\int^{r_{tp}^{+}}_{r_{tp}^{-}}\frac{|\sqrt{-V}|dr}{\Delta}.\end{split}$ (85)

Since $S_{-+}$ is real and positive, this approximate expression makes it
manifest that $|\mathcal{R}|<1$ for $\tilde{\omega}_{h}>0$ whilst
$|\mathcal{R}|>1$ for $\tilde{\omega}_{h}<0$. This is consistent with (62),
where for $\tilde{\omega}_{h}<0$ ($\tilde{\omega}_{h}>0$), the in-going mode
in the core is on $\omega_{D}^{-}$ ($\omega_{D}^{+}$) implying that the non-
zero term on the right-hand zero of (62) is positive (negative) and the
reflected wave is amplified (attenuated).

In Fig. 4, we explicitly compare the expression in (85) to the exact value
obtained by numerically integrating the wave equation over the bandwidth of
superradiant frequencies (see the Appendix A of [25] for details of the
integration). We observe excellent agreement at low circulation whilst for
higher values, our expression captures only the overall behaviour and order of
magnitude. The deviation occurs since for larger $C$, the tunnelling region
between $\omega^{+}$ and $\omega^{-}$ becomes smaller for higher frequency
(see Fig. 3 panel (b)) whilst WKB works best for well separated turning
points. Another feature that is not captured by the WKB approximation is that
the correct $|\mathcal{R}|$ smoothly passes through unity crossing the
superradiant threshold $\tilde{\omega}_{h}=0$. The reason is that at the
geometric optics level, the change from superradiant amplification to
absorption corresponds to a discontinuous transition where the two modes
$p^{\pm}$ transfer from the lower branch of the dispersion relation at small
$r$ to the upper branch. These differences between will be important in the
next section when comparing the WKB growth rates of unstable modes with those
obtained by a numerical computation. One feature that is well captured by WKB
is the decrease in the amount of amplification with $m$. This can be seen from
(85) since, for a given $\omega/m$, $S_{-+}$ increases linearly with $|m|$
which means deviations of $|\mathcal{R}|$ from unity decrease. Hence the most
amplified modes in the system will occur for $m=1$.

Figure 4: Comparison between the WKB expression (85) for the reflection
coefficient (solid lines) and the exact value obtained numerically (dashed
lines). Each panel represents a different dimensionless rotation parameter,
and the horizontal axis is scaled so that the superradiant threshold coincides
with 1. On each panel, we display the $m=1$ (black), $m=2$ (blue) and $m=3$
(red) results, which decrease in overall magnitude with $m$. The agreement is
best at low $C$. Furthermore, the WKB result should not be trusted at the
superradiant threshold since the geometric optics limit erroneously predicts a
discontinuous change from amplification to absorption.

####  IV.2.3 Bound states

Having solved for the radial momenta $p^{\pm}$ and the reflection coefficient
$\mathcal{R}$, we are poised to compute the resonance condition for trapped
waves between the vortex and the boundary. First we note that, in our
formalism, this trapping can only occur when we have real turning points. The
frequency of any trapped modes must lie between the light-ring frequencies
$\omega_{lr}^{-}<\omega<\omega_{lr}^{+}$. A further constraint can be imposed
by the requirement that the solution should be propagating approaching the
boundary at $r_{B}$, leading to the constraint
$\omega_{lr}^{-}<\omega<\omega^{-}(r_{B})$ for negative frequencies and
$\omega^{+}(r_{B})<\omega<\omega_{lr}^{+}$ for positive frequencies. The
resonance condition can then be found by transporting the amplitude of the in-
going mode at the boundary $A^{-}_{B}$ up to the turning point $r_{tp}^{+}$
(according to (78) with the appropriate phase factor) where it gets reflected
and transported back out to the boundary. A detailed description of the
computation can be found in Appendix D, where it is described as bound state
A. This leads to the condition for the real part of the frequency
$\omega_{n}$,

$S_{+B}(\omega_{n})\equiv\int^{r_{B}}_{r_{tp}^{+}(\omega_{n})}\frac{\sqrt{-V(\omega_{n})}dr}{\Delta}=\pi\left(n+\frac{1}{4}\right),$
(86)

where $n=0,1,...$ indexes the different modes which solve the condition. A
subtlety of the derivation is that we approximate the radial velocity at
$r_{B}$ as zero so that the no penetration boundary condition implies the
$\pm$ modes have equal amplitudes there. In reality, the resonance condition
will be modified by either the small value of $v_{r}(r_{B})$ or the inlet
region where $v_{r}$ smoothly goes to zero on the boundary. However, we will
soon see that the predictions of (86) agree exceptionally well with numerical
calculations, suggesting that any differences due to the flow field over the
inlet will be negligible.

Assuming the imaginary part of the complex eigenfrequency $\Gamma_{n}$ is a
small correction to the real part, we find in Appendix D the expression,

$\Gamma_{n}=\frac{\log|\mathcal{R}|}{2\partial_{\omega}S_{+B}}\Bigg{|}_{\omega=\omega_{n}},$
(87)

where $\mathcal{R}$ is given in (85). Hence, we see explicitly that if the
resonant frequency $\omega_{n}$ lies within the superradiant range where
$|\mathcal{R}|>1$, $\Gamma_{n}>0$ and the bound state will be unstable. Since
we are only interested in the unstable modes, we can place constraints on
values of the rotation parameter which result in an instability. A minimum
constraint for the existence of a bound state is that $r_{tp}^{+}<r_{B}$ in
(86) and, for superradiant frequencies, the minimum value of the outer turning
point occurs at the superradiant threshold $\tilde{\omega}_{h}=0$. We can then
invert the requirement $r_{B}>r_{tp}^{+}(\omega=mC)\equiv\sqrt{C^{2}+1}/C$ to
find the minimum $C$ value below which superradiant bound states cannot exist,

$C_{\mathrm{min}}=(r_{B}^{2}-1)^{-\frac{1}{2}}.$ (88)

For our analysis to apply, we must clearly have $r_{B}>1$ (i.e. the mirror is
outside the horizon). Thus, as we bring the reflective boundary to its lowest
allowed location, the rotation required for a superradiantly unstable mode
tends to diverges. In practice, this means that once we select a value of
$r_{B}$, only $C>C_{\mathrm{min}}$ values are relevant to the study of the
instability. Alternatively, if we fix $C$, the boundary should satisfy
$r_{B}>\sqrt{C^{2}+1}/C$.

Figure 5: We display the bound state frequencies
$\omega^{\mathbb{C}}_{n}=\omega_{n}+i\Gamma_{n}$ obtained by solving the WKB
conditions in (86) and (87) as solid points (or solid lines when there are
sufficiently many states that the points become indiscernible). These are
compared with numerical results of Section IV.3, shown as open circles/broken
lines. The real part $\omega_{n}$ is scaled by the dimensionless rotation so
that the horizontal axis is the same for all panels. The agreement of the real
part is excellent in all cases whilst the imaginary part begins to deviate for
large rotation. Results of the time domain simulations from Section IV.4 are
displayed in a few cases as crosses and match the frequency domain
simulations.

On Fig. 5, we display the complex frequencies
$\omega^{\mathbb{C}}_{n}=\omega_{n}+i\Gamma_{n}$ for the first three $m$–modes
for a variety of $C$ and $r_{B}$ values. We find that the $m=1$ mode is the
most unstable in all cases considered. This is because reflection coefficient
is largest for the lowest $m$–modes, as seen in Fig. 4. For a given $C$ and
$r_{B}$, the spacing between modes is roughly constant, and unstable modes
exist all the way up to the superradiant threshold. These two facts mean that,
on a given panel, there are roughly $m$ times as many bound states in the
$m>1$ channel compared with $m=1$ (at least in the limit of a large number of
states).

We observe that increasing both $r_{B}$ and $C$ leads to the existence of more
bound states. In the first case, increasing $r_{B}$ increases the density of
states and the modes pack closer together in the allowed $\omega$–range. In
the second case, increasing $C$ raises the superradiant threshold leaving more
room in the allowed $\omega$–range for the existence of trapped modes. This is
further exemplified by Fig. 6. At fixed $C$, as we increase $r_{B}$ the
oscillation frequency of each mode is $\sim 1/r_{B}$ asymptotically whilst at
fixed $r_{B}$, they increase linearly with $C$ asymptotically.

There are two competing factors that determine the imaginary part. The first
is that the factor $\log|\mathcal{R}|$ is larger for higher frequencies. This
leads to the increase of the growth rate with increasing $\omega_{n}$ at low
$C$ and high $r_{B}$ on Fig. 5. However, the time taken to orbit the cavity
$2\partial_{\omega}S_{+B}$ goes to zero as the oscillation frequency goes to
its lowest allowed value. If the frequency is small enough, this can outweigh
the effect of the $\log|\mathcal{R}|$ and make small frequencies grow faster,
e.g. in the top right panel of Fig. 5. This can also be inferred from the
lower left panel of Fig. 6 where, as $r_{B}$ decreases, the $n=0$ mode
eventually attains a larger growth rate than any of the higher $n$.

Figure 6: Variation of the bound state frequencies with the cylinder size
(left) and the dimensionless rotation (right). Only a few modes are displayed
for the imaginary part to improve their visibility (on the plots of the real
part, the modes are indexed $n=0,1,2,...$ as they appear from left to right).
Increasing both $C/D$ and $r_{B}c/D$ increases the number of bound states,
with new modes first appearing at the superradiant threshold. The WKB
calculations are shown as solid black lines whilst numerical results are
broken red. Agreement in the real part is excellent in all cases, whilst the
imaginary part captures the order of magnitude and deviates significantly for
modes near the superradiant threshold.

###  IV.3 Frequency-domain simulations

To validate the WKB used in the previous section, we compute numerically the
spectrum of unstable modes in the shallow water regime by integrating (72)
subject to in-going boundary conditions at the horizon up to $r=r_{B}$. It is
numerically more convenient to work with the radial coordinate $r$ rather than
$r_{*}$, therefore the equation to be integrated is,

$\displaystyle\Delta^{2}R^{\prime\prime}(r)+2r^{-3}\Delta
R^{\prime}(r)+V(r)R(r)=0.$ (89)

The initial conditions are applied just outside the horizon, i.e. at
$r=1+\epsilon$. To ensure that our solution satisfies the ingoing boundary
condition, we write $R(r)$ as,

$R(r)=e^{-i\tilde{\omega}_{h}r_{*}}\sum_{n=0}^{N}a_{n}(r-1)^{n}.$ (90)

The coefficients ${a_{n}}$ satisfy a 6–term recurrence relation given by,

$\displaystyle\alpha_{n}a_{n+1}$ $\displaystyle+$
$\displaystyle\beta_{n}a_{n}+\gamma_{n}a_{n-1}$ (91) $\displaystyle+$
$\displaystyle\delta_{n}a_{n-2}+\epsilon_{n}a_{n-3}+\lambda_{n}a_{n-4}=0,$

where the coefficients are explicitly,

$\displaystyle\alpha_{n}$ $\displaystyle=$ $\displaystyle
8(1+n)(1+n-i\tilde{\omega}_{h}),$ (92) $\displaystyle\beta_{n}$
$\displaystyle=$ $\displaystyle-4[1+m^{2}-5n^{2}+n(3+8i\tilde{\omega}_{h})$
$\displaystyle-2Cm\tilde{\omega}_{h}],$ $\displaystyle\gamma_{n}$
$\displaystyle=$ $\displaystyle 34+8(C^{2}-1)m^{2}-48n+16n^{2}$
$\displaystyle+\tilde{\omega}_{h}[48i(1-n)+32Cm],$ $\displaystyle\delta_{n}$
$\displaystyle=$ $\displaystyle 25+4(5C^{2}-1)m^{2}-20n+4n^{2}$
$\displaystyle+\tilde{\omega}_{h}[32i(2-n)+48Cm],$ $\displaystyle\epsilon_{n}$
$\displaystyle=$ $\displaystyle 16C^{2}m^{2}+8\tilde{\omega}_{h}[4Cm+i(3-n)],$
(96) $\displaystyle\lambda_{n}$ $\displaystyle=$ $\displaystyle
4Cm(Cm+2\tilde{\omega}_{h}).$ (97)

From the initial condition (90) with $N=100$, we numerically integrate (89)
from $r=1+\epsilon$, with $\epsilon=10^{-5}$, up to $r=r_{B}$. We then extract
the value of the derivative of the numerical solution at the wall
$R^{\prime}_{\mathrm{num}}(\omega,r=r_{B})$. The numerical integration is
performed using the Mathematica function NDSolve with the stiffness switching
option. Bound states frequencies are then defined by the condition,

$R^{\prime}_{\mathrm{num}}(\omega,r=r_{B})=0.$ (98)

This condition is solved numerically using the FindRoot function with the WKB
prediction as an initial guess.

The results are displayed on Figs. 5 and 6. In all cases, we observe excellent
agreement with the WKB prediction for the oscillation frequency. This is
because it is determined by the phase integral in a region of the flow where
the velocity field of the central vortex is weakest. Consequently, the
background varies smoothly and the WKB approximation works well. Conversely,
agreement with the growth rate works well for reasonably low circulation but
deviates notably for larger values. The reason for this discrepancy can be
traced back to the reflection coefficient in Fig. 4.

###  IV.4 Time-domain simulations

Figure 7: The time development of the signal is extracted at the location
$r/r_{h}=2$ and is plotted on a logarithmic scale in the upper panel. The
growth rate $\Gamma_{n}$ is extracted by fitting the signal at late times to a
straight line. The oscillation frequency $\omega_{n}$ is obtained by finding
the peak in Fourier space (lower panel). The simulation parameters were $m=1$,
$C=1$, $r_{B}=5$.

The 1D equation in (72) can be simulated in the time-domain by restoring the
time derivatives,

$\begin{split}\partial_{t}\begin{pmatrix}R\\\ \Pi\end{pmatrix}=&\
\begin{pmatrix}-imC/r^{2}&1\\\
\partial_{r_{*}}^{2}-U(r)&-imC/r^{2}\end{pmatrix}\begin{pmatrix}R\\\
\Pi\end{pmatrix},\\\ U(r)=&\
\Delta(r)\left(\frac{m^{2}-1/4}{r^{2}}+\frac{5}{4r^{4}}\right),\end{split}$
(99)

where $\Pi=(\partial_{t}+imC/r^{2})R$ and $r$ is understood as an implicit
function of $r_{*}$ via (71). This can be simulated using a
$4^{\mathrm{th}}$–order Runge-Kutta algorithm. At the inner and outer
boundaries, we impose absorbing and reflective (Neumann) boundary conditions
respectively,

$\left(\Pi-\partial_{r_{*}}R\right)|_{r_{*}=r_{*1}}=0,\qquad\partial_{r_{*}}R|_{r_{*}=r_{*B}}=0,$
(100)

where $r_{*1}$ is large and negative to approximate the location of the
horizon (which is negative infinity) and $r_{*B}=r_{*}(r_{B})$. The initial
condition is taken as a plane wave with wavenumber $\omega_{0}$ whose envelope
is a Gaussian of width $\sigma_{0}$ centred at a location $r_{*0}$ near the
outer boundary,

$\begin{split}R(t=0)=&\
\exp\left[-\frac{(r_{*}-r_{*0})^{2}}{2\sigma_{0}^{2}}\right]e^{i\omega_{0}r_{*}},\\\
\Pi(t=0)=&\ 0.\end{split}$ (101)

We use the parameters $\sigma_{0}=3$, $r_{*0}=3$ and $\omega_{0}$ is the
frequency of the bound state with the largest growth rate in the WKB
approximation.

In the time-domain, the fastest growing mode in the system will dominate the
signal and be measurable at late times. We analyse the signal by plotting the
time development at a single location, see Fig. 7. The oscillation frequency
of this mode is extracted by performing a Fourier transform, and the growth
rate is found by fitting a straight line to the logarithm of the signal
magnitude. On Fig. 5, we indicate the complex frequencies extracted in this
manner with crosses. When there are several nearby frequencies with comparable
growth rates, the signal can exhibit a strong beating effect and it is not
possible to reliably extract the growth rate with this method. Hence, results
are only presented for the cases where a single mode clearly dominates the
response. In each of these cases, we observe excellent agreement with results
obtained in our frequency domain simulations.

###  IV.5 Validity

In order for the shallow water analysis above to hold, the wavelength of the
bound states must be much larger than the height of the fluid. Since the
characteristic wavelength of the lowest $k$ modes in the system is set by
$r_{B}$, the relevant condition becomes $h_{0}\ll r_{B}$. A further
approximation was that the curvature of the free surface was assumed
negligible, i.e. $h^{\prime}\ll 1$. To see when this is satisfied, consider
the basic velocity profile in (11) and the height profile with the leading
order curvature correction (see Appendix A),

$h=h_{0}\left(1-\frac{r_{a}^{2}}{r^{2}}\right).$ (102)

The location of horizon will be shifted in this case but it will still be
$\mathcal{O}(D/\sqrt{gh_{0}})$ provided the circulation is not too large. The
constraint on $h^{\prime}$ then leads to the condition,

$r_{h}\gg r_{s}\equiv\left(\frac{D^{2}+C^{2}}{g}\right)^{\frac{1}{3}}.$ (103)

This condition improves for small $C$, and works best when $C=0$. In that
case, it can be reformulated as,

$D\gg\sqrt{gh_{0}^{3}}.$ (104)

Altogether, we end up with a hierarchy of scales,

$h_{0}\ll r_{s}\ll r_{h}<r_{B}.$ (105)

For $r_{B}$ on the order of tens of centimeters, the requirement above is
satisfied for $h_{0}$ on the order of millimeters. For these low water
heights, the viscous boundary layer is expected to spread throughout the whole
fluid volume, at least near the drain where the flow velocities are large.
Furthermore, the drain is powered by gravity itself, hence, there is a
complicated relation between $D$ and $h_{0}$. A basic argument is the
following. Starting from initial values $D_{i}$ and $h_{0i}$ and increasing
the flow rate, the system would reach a new equilibrium with final values
$D_{f}>D_{i}$ (to compensate the increased flow rate) and $h_{0f}>h_{0i}$ (due
to the influx of more water). In other words, one cannot increase $D$
arbitrarily without also increasing $h_{0}$, thereby making the criterion
above impractical to satisfy.

###  IV.6 Key features

We conclude this section by highlighting the key features of the BHB which
should be present when generalising to the dispersive regime. These are:

  * •

A region inside the vortex permitting negative energy waves.

  * •

A mechanism for absorbing negative energies.

The first criterion is realised by the ergosphere in (67) which leads to
positive frequency modes with negative norm, i.e. $\Omega<0$ in the WKB
framework. Hence, this is the feature to anticipate when we extend to more
general scenarios. The second criterion is provided in shallow water by the
analogue horizon, however, the practical requirement is the in-going boundary
condition into the centre of the vortex in (74). This could be achieved for a
more generic absorption mechanism without the strict notion of a horizon. In
our case, for example, the draining vortex has a hole in the middle which is a
sink for the wave norm. Furthermore, one could even relax the purely in-going
boundary condition and require only that the out-going waves are sufficiently
low amplitude that they can be treated as a higher order effect (see e.g.
related work in [57, 58]).

##  V Deep water waves

To recap, so far we have derived a set of equations which govern the
propagation of surface waves around a draining vortex (Section II), presented
a general general method for studying the scattering of such waves (Section
III) and demonstrated that these methods accurately predict the unstable mode
frequencies in the shallow water regime (Section IV).

However, the difficulty in realising shallow water conditions in practice
motivates the study of the opposite regime $hk\gg 1$, known as the deep water
regime. When the wavelength becomes suitably short, the effect of surface
tension and viscosity becomes non-negligible. The main difference compared
with the low $k$ theory is that in deep water, one has to keep track of all 6
modes shown in Fig. 2. This significantly enhances the number of possible
scattering processes and thereby greatly increases the complexity of the
problem.

Our approach in this section will be to systematically analyse all scattering
processes allowed by the deep water dispersion relation (which will be shown
in Fig. 9) before identifying those relevant to the BHB instability in a
physically motivated parameter range. The main result will be a compact
expression governing the unstable mode frequencies and growth rates in (117)
and (118), the properties of which are exemplified in the figures shown toward
the end of the section.

###  V.1 Simplifications

We again work with the basic flow field in (11) where the water height $h$ is
treated as a constant. To completely characterise the scattering properties of
the system, one must solve for the $p^{j}(r)$ at each point in the parameter
space of (36). Solving $\mathcal{H}=0$ with $\mathcal{D}(k)$ defined in (34)
is a daunting task, since one must find the roots of a transcendental equation
over the 7–dimensional parameter space $(\omega,m,C,D,\gamma,g,h)$. We can
rescale $\mathcal{H}$ to effectively reduce the dimensionality of parameter
space. Following [25], we choose to remove $D$ and $g$ axes. We define,

$l=(D^{2}/g)^{\frac{1}{3}},\qquad\tau=(D/g^{2})^{\frac{1}{3}},$ (106)

and perform the following rescalings,

$\begin{split}&\frac{r}{l}\to r,\quad\frac{h}{l}\to h\quad pl\to
p,\quad\omega\tau\to\omega,\\\
&\frac{\gamma}{gl^{2}}\to\gamma\quad\frac{C}{D}\to
C,\quad\mathcal{H}\tau^{2}\to\mathcal{H},\quad\frac{\nu}{D}\to\nu,\end{split}$
(107)

where, in anticipation of including the effects of viscosity at next to
leading order, we have also rescaled $\nu$. The result is equivalent to
setting $D=g=1$:

$\begin{split}\mathcal{H}=&\
-\frac{1}{2}\bigg{[}\left(\omega-\frac{mC}{r^{2}}+\frac{p}{r}\right)^{2}-\left(k+\gamma
k^{3}\right)\times\\\ &\qquad\quad\tanh(hk)\bigg{]},\qquad
k=\sqrt{p^{2}+\frac{m^{2}}{r^{2}}}.\end{split}$ (108)

Unlike the shallow water regime, the magnitude of the group velocity
$|\partial_{\mathbf{k}}\omega|$ is no longer independent of frequency. In
other words, different frequencies will travel at different speeds; the medium
is dispersive. This is the physical reason for the appearance of extra spatial
modes in the system compared with the shallow water case.

Considerable further simplification occurs if we work in the deep water limit
$hk\gg 1$. In this case, we can approximate $\tanh(hk)\to 1$ in the expression
above, which reduces the parameter space to $(\omega,m,C,\gamma)$.
Furthermore, we needn’t assume that the water height is constant in this
regime; all deep water waves have wavelength much longer than the local water
height, therefore, variations in $h$ become insignificant. Although this is a
simplification of the full problem, many of the scattering processes we will
encounter result from the large $k$ properties of $\mathcal{H}$. Intuitively,
this is because the large velocities achieved at small $r$ eventually blue-
shift in-coming waves into the regime where $k\gg h^{-1}$. Hence, qualitative
features of the system in the $hk\to\infty$ limit will also carry over to the
full dispersion relation.

###  V.2 Scattering

Since we are ultimately concerned with the BHB instability, which arises for
$m>0$, we focus our attention on $m\neq 0$. In this case, we can scale out $m$
from the effective Hamiltonian. First we notice that $\mathcal{H}$ is
symmetric under the transformation $(\omega,m,p)\to(-\omega,-m,-p)$.
Therefore, positive frequency waves with $m<0$ will scatter the same as the
negative frequency waves with the opposite radial and angular momentum.
Without loss of generality, we can therefore fix $m>0$ and consider the
frequency as ranging from negative to positive infinity. Now we perform the
following rescalings,

$\begin{split}&r/m^{\frac{1}{3}}\to
r,\quad\omega/m^{\frac{1}{3}}\to\omega\quad p/m^{\frac{2}{3}}\to p,\\\ &\gamma
m^{\frac{4}{3}}\to\gamma\quad\mathcal{H}/m^{\frac{2}{3}}\to\mathcal{H},\end{split}$
(109)

and write the effective Hamiltonian as,

$\begin{split}\mathcal{H}=&\
-\frac{1}{2}\left[\left(\omega-\frac{C}{r^{2}}+\frac{p}{r}\right)^{2}-k-\gamma
k^{3}\right],\\\ k=&\ \sqrt{p^{2}+\frac{1}{r^{2}}}.\end{split}$ (110)

The parameter space is now $(\omega,C,\gamma)$. As we proceed, our tactic will
be to understand how the ray trajectories vary throughout the parameter space
at fixed $\gamma$. We will then vary $\gamma$ from small to larger values to
see how the structure of the parameter space changes.

Figure 8: An exemplary interaction between the 6 capillary-gravity wave modes
for the parameters $C/D=15$, $\gamma m^{\frac{4}{3}}/gl^{2}=10^{-4}$ and
$\omega\tau=2.3$. Main panel: On the vertical axis (which has non-uniform
spacing), we show the real part of the radial wavevector of the 6 modes that
solve $\mathcal{H}=0$. As $r$ is varied, the value of $p^{j}$ changes. Solid
lines represent propagating modes with $p^{j}\in\mathbb{R}$. Neighbouring
$p^{j}$ coincide at turning points (black dots) and are converted into
evanescent modes with $p^{j}\in\mathbb{C}$ (dashed lines). Note that a single
dashed line represents a complex conjugate mode pair. In the classification
scheme of Fig. 9, this diagram corresponds to Type VI scattering. Small
panels: Propagating modes are shown as black points on the two branches of the
dispersion relation $\omega_{D}^{\pm}$ for 4 illustrative radii (indicated on
the right of the panels). The shape of the branches varies with $r$, changing
the locations of the $p^{j}$ roots and generating the trajectories shown on
the main panel.

In shallow water, there are two modes ($p^{\pm}$) which either superradiantly
scatter, non-superradiantly scatter or do not scatter [26]. These three
possibilities cover the entire parameter space of the low $k$ theory. In deep
water, there are 6 propagating modes that can scatter with each other. We show
an example of this scattering for a particular choice of parameters in Fig. 8.
The diagram shows the trajectories of the WKB modes through the $(r,p)$–phase
space. There are 6 functions $p^{j}(r)$ from which one can reconstruct the WKB
solution of the scattering problem using (51). Hence, drawing the scattering
diagram is equivalent to knowing the general solution to the linear equations
of motion for a given parameter set. The mode amplitudes on the far left and
right sides of the diagram are related by the full transfer matrix
$\mathcal{M}$ defined in (54). For example, the matrix corresponding to Fig. 8
is,

$\begin{split}\mathcal{M}=&\
\mathrm{diag}\big{(}\mathcal{F}^{t}_{01},\mathcal{F}^{u}_{01},\mathcal{F}^{+}_{01},\mathcal{F}^{\uparrow}_{01},\mathcal{F}^{\downarrow}_{01},\mathcal{F}^{b}_{01}\big{)}\tilde{t}_{1}\\\
\times&\
\mathrm{diag}\big{(}\mathcal{F}^{t}_{12},\mathcal{F}^{u}_{12},\mathcal{F}^{+}_{12},\mathcal{F}^{-}_{12},\mathcal{F}^{d}_{12},\mathcal{F}^{b}_{12}\big{)}\tilde{t}_{2}\\\
\times&\
\mathrm{diag}\big{(}\mathcal{F}^{t}_{23},\mathcal{F}^{u}_{23},\mathcal{F}^{\uparrow}_{23},\mathcal{F}^{\downarrow}_{23},\mathcal{F}^{d}_{23},\mathcal{F}^{b}_{23}\big{)}\tilde{t}_{3}\\\
\times&\
\mathrm{diag}\big{(}\mathcal{F}^{t}_{34},\mathcal{F}^{u}_{34},\mathcal{F}^{+}_{34},\mathcal{F}^{-}_{34},\mathcal{F}^{d}_{34},\mathcal{F}^{b}_{34}\big{)}\tilde{t}_{4}\\\
\times&\
\mathrm{diag}\big{(}\mathcal{F}^{t}_{45},\mathcal{F}^{\uparrow^{+}}_{45},\mathcal{F}^{\downarrow^{+}}_{45},\mathcal{F}^{-}_{45},\mathcal{F}^{d}_{45},\mathcal{F}^{b}_{45}\big{)}\tilde{t}_{5}\\\
\times&\
\mathrm{diag}\big{(}\mathcal{F}^{t}_{56},\mathcal{F}^{\uparrow^{+}}_{56},\mathcal{F}^{\downarrow^{+}}_{56},\mathcal{F}^{-}_{56},\mathcal{F}^{\uparrow^{-}}_{56},\mathcal{F}^{\downarrow^{-}}_{56}\big{)}\tilde{t}_{6}\\\
\times&\
\mathrm{diag}\big{(}\mathcal{F}^{t}_{67},\mathcal{F}^{u}_{67},\mathcal{F}^{+}_{67},\mathcal{F}^{-}_{67},\mathcal{F}^{\uparrow^{-}}_{67},\mathcal{F}^{\downarrow^{-}}_{67}\big{)}\tilde{t}_{7}\\\
\times&\
\mathrm{diag}\big{(}\mathcal{F}^{\uparrow^{+}}_{7B},\mathcal{F}^{\uparrow^{+}}_{7B},\mathcal{F}^{+}_{7B},\mathcal{F}^{-}_{7B},\mathcal{F}^{\uparrow^{-}}_{7B},\mathcal{F}^{\downarrow^{-}}_{7B}\big{)}\end{split}$
(111)

where the locations numbered 1 to 7 correspond to the turning points in Fig.
8. The matrices $\tilde{t}_{j}$ are constructed from the transfer matrices in
(46) and (47) and are given explicitly by,

$\begin{split}\tilde{t}_{1}=&\ \begin{pmatrix}1&0&0&0&0&0\\\ 0&1&0&0&0&0\\\
0&0&1&0&0&0\\\ 0&0&0&\beta&\beta^{*}&0\\\ 0&0&0&\alpha^{*}&\alpha&0\\\
0&0&0&0&0&1\\\ \end{pmatrix},\ \ \tilde{t}_{2}=\begin{pmatrix}1&0&0&0&0&0\\\
0&1&0&0&0&0\\\ 0&0&\beta^{*}&\alpha&0&0\\\ 0&0&\beta&\alpha^{*}&0&0\\\
0&0&0&0&1&0\\\ 0&0&0&0&0&1\\\ \end{pmatrix},\\\ \tilde{t}_{3}=&\
\begin{pmatrix}1&0&0&0&0&0\\\ 0&1&0&0&0&0\\\ 0&0&\beta&\beta^{*}&0&0\\\
0&0&\alpha^{*}&\alpha&0&0\\\ 0&0&0&0&1&0\\\ 0&0&0&0&0&1\\\ \end{pmatrix},\ \
\tilde{t}_{4}=\begin{pmatrix}1&0&0&0&0&0\\\ 0&\beta^{*}&\alpha&0&0&0\\\
0&\beta&\alpha^{*}&0&0&0\\\ 0&0&0&1&0&0\\\ 0&0&0&0&1&0\\\ 0&0&0&0&0&1\\\
\end{pmatrix}\\\ \tilde{t}_{5}=&\ \begin{pmatrix}1&0&0&0&0&0\\\ 0&1&0&0&0&0\\\
0&0&1&0&0&0\\\ 0&0&0&1&0&0\\\ 0&0&0&0&\beta^{*}&\alpha\\\
0&0&0&0&\beta&\alpha^{*}\\\ \end{pmatrix},\ \
\tilde{t}_{6}=\begin{pmatrix}1&0&0&0&0&0\\\ 0&\beta&\beta^{*}&0&0&0\\\
0&\alpha^{*}&\alpha&0&0&0\\\ 0&0&0&1&0&0\\\ 0&0&0&0&1&0\\\ 0&0&0&0&0&1\\\
\end{pmatrix},\\\ \tilde{t}_{7}=&\ \begin{pmatrix}\beta^{*}&\alpha&0&0&0&0\\\
\beta&\alpha^{*}&0&0&0&0\\\ 0&0&1&0&0&0\\\ 0&0&0&1&0&0\\\ 0&0&0&0&1&0\\\
0&0&0&0&0&1\\\ \end{pmatrix},\end{split}$

with $\alpha=e^{i\pi/4}$ and $\beta=\alpha/2$ defined for compactness. Mode
mixing results from the off-diagonal elements in the $\tilde{t}_{j}$. Imposing
the boundary conditions at $r_{0}$ and $r_{B}$, $\mathcal{M}$ determines a
resonance condition. However, it is evidently extremely impractical to write
down $\mathcal{M}$ like this and in practice, it is much simpler to identify
the relevant rays that form the bound state loops in phase space and transport
the modes around them (see Appendix F).

Figure 9: Schematics illustrating the 13 different scattering possibilities
allowed by the deep water system in (110). Solid black (red) lines represent
modes propagating on the upper (lower) branch of the dispersion relation in
the direction shown by the arrow. Dashed lines represent a pair of evanescent
modes. The horizontal axis represents $r$ with the left and right sides
representing the vortex core and spatial infinity respectively. The vertical
axis represents $\mathrm{Re}[p]$, or more specifically, the relative values of
$\mathrm{Re}[p]$ of propagating modes at a given $r$ (compare, e.g. the Type
VI schematic with Fig. 8). The aim of the diagrams is to clearly indicate the
locations of the turning points and which modes interact there.

Since there is a large range of mode interactions allowed by (110), it is
helpful to classify them according to the turning point structure, i.e. the
locations of the various $r_{tp}$ and the modes that interact at each. The
different mode interactions can then be schematically illustrated as one of
the 13 diagrams in Fig. 9, which comprise an exhaustive list of the mode
interactions allowed by (110). The names of the diagrams (Types I–XIII) are
chosen in the following way. The first 6 (I–VI) are modifications of the deep
water diagrams in [25] that regulate the theory by removing the propagating
modes at infinity with divergent $k$. The next three (VI–IV) are
characteristic of super-luminal dispersion relations, e.g. those studied in
[26]. The last 4 (X–XIII) result from the mixed character of the capillary-
gravity wave dispersion relation, which contains both sub- and super-luminal
sections.

The diagrams are schematic illustrations of the $(r,p)$–phase space
exemplified by Fig. 8. These diagrams represent the real part of $p$ (vertical
axis) as a function of $r$ (horizontal axis). Solid lines indicate
$p\in\mathbb{R}$ with the arrow indicating the direction of propagation as
determined by the radial component of the group velocity. Dashed lines
indicate $p\in\mathbb{C}$. Black (red) lines represent modes on
$\omega_{D}^{+}$ ($\omega_{D}^{-}$) which have positive (negative) norm
density $\rho_{n}$ (59). The ordering of lines is such that $\mathrm{Re}[p]$
increases from bottom to top for propagating modes (note, this is not
necessarily true for the dashed lines since evanescent modes can pass round
propagating modes in the complex plane, although this behaviour is omitted
since the aim of the diagrams is to clearly indicate the different mode
interactions). Bottom to top ordering here translates to the left to right
ordering of the real roots of $\omega=\omega_{D}^{\pm}(p)$ (see Fig. 2) and
since the labels $j=\\{b,d,-,+,u,t\\}$ go from left to right on the dispersion
relation, this translates to an unambiguous labelling convention for the modes
in Fig. 9.

Locations where solid lines connect to dashed lines, i.e. vertices, represent
turning points $r_{tp}$. The number of modes is preserved either side of
$r_{tp}$ since the dashed lines actually represent two complex conjugate modes
with the same $\mathrm{Re}[p]$. Each diagram is unique in that it has a fixed
structure of vertices 333Note that for vertices connected by more than one
line (or not connected) the relative locations of the $r_{tp}$ are not
necessarily captured by the diagrams. In other words, the diagrams can be
deformed within the same classification provided: 1) vertices are neither
created nor destroyed, 2) the direction (right or left) of the arrows is
unaltered and 3) the ordering of propagating modes from top to bottom is
respected. For example, the turning point between the $p^{\pm}$ modes on the
Type II scattering schematic in Fig. 9 is not necessarily at larger $r$
compared to the turning point involving the $p^{u,t}$ modes (as can be seen by
comparing with panel (c) of Fig. 10).. The scattering processes for $\omega<0$
are inverted versions (i.e. $p\to-p$) of those for $\omega>0$, which we denote
with a superscript $\downarrow$. Note, in the inversions Type VII↓ and VIII↓,
it is more natural to label the modes which propagate as $r\to\infty$ as
$p^{\pm}$. In this way, the mode labels are correct when transitioning between
different diagrams as the frequency is varied.

Each diagram corresponds to a particular structure in the transfer matrix (54)
and therefore a particular computation of the asymptotic scattering
coefficients or (in our case) bound state integrals. However, to evaluate
these integrals, we need quantitative values of the different $p^{j}$ between
the various turning points. The first step in this procedure is to find
accurate values for the $r_{tp}$. In deep water (unlike the low $k$ theory
where we had a simple analytic formula (80)) there is no closed form
expression for the $r_{tp}$. Instead, we consider the curves
$\omega^{\pm}_{a}(r)$ defined in (49) (see Appendix E.1 where we detail our
procedure to calculate these). Examples of the $\omega^{\pm}_{a}(r)$ are shown
in Fig. 10 for the case where the $\omega_{D}^{\pm}$ have three (low $\gamma$)
and one (high $\gamma$) extrema respectively. When there are three extrema
(and therefore three curves $a=1,2,3$) the two which are real as $r\to\infty$
are $\omega_{1}^{\pm}$ and describe the scattering of the $p^{\pm}$ modes on
the upper and lower branches of the dispersion relation respectively.
$\omega_{1}^{+}$ then joins onto $\omega_{2}^{+}$ (which determines how
$p^{+,u}$ modes scatter) and $\omega_{2}^{+}$ in turn joins onto
$\omega^{+}_{3}$ (which determines how $p^{u,t}$ modes scatter). Similarly,
$\omega_{1}^{-}$ joins onto $\omega_{2}^{-}$ (which determines how $p^{d,-}$
modes scatter) and $\omega_{2}^{-}$ in turn joins onto $\omega_{3}^{-}$(which
determines how $p^{b,d}$ modes scatter). When $\gamma$ is sufficiently large,
there is only a single real extremum on each of the $\omega^{\pm}_{D}$ and the
two corresponding curves are called $\omega^{\pm}_{1}(r)$. With these
functions in hand, we can find the $r_{tp}$ by selecting a frequency and
solving (50) graphically.

Figure 10: In panel (a), we display examples of the curves defined in (49)
when $\gamma m^{\frac{4}{3}}/gl^{2}$ is small. The specific parameters used
are $C/D=8$ and $\gamma m^{\frac{4}{3}}/gl^{2}=10^{-3}$. The characteristic
frequencies defined in (112) are also indicated. The black curves are the
$\omega^{+}_{a}$, with $a=1$ extending to radial infinity, $a=2$ connecting
the first curve at $\omega^{+}_{I1}$ and $a=3$ connecting the second curve at
$\omega^{+}_{I2}$. The red curves are the $\omega^{-}_{a}$ with the labelling
the same. In panels (b,c), we show the scattering of the 6 modes by plotting
the $(r,p)$ phase space for two example frequencies. The turning points are
shown as black points and are determined by the points of intersection between
the curves on panel (a) and a line of constant $\omega$. In panels (d-f), we
display the same but for high $\gamma$, with the specific parameters $C/D=4$
and $\gamma m^{\frac{4}{3}}/gl^{2}=0.1$. In this case, the characteristic
frequencies are those in (115). The black curve is $\omega^{+}_{1}$ and the
red curve is $\omega^{-}_{1}$. Using the classification scheme of Fig. 9,
panels (b,c,e,f) correspond respectively to Type V, II, IX, VII scattering.

To determine which of the 13 scattering possibilities occurs for a particular
frequency, we can calculate various characteristic frequencies
$\omega(C,\gamma)$ which divide up the parameter space into regions where
different scattering processes occur. We defer a thorough discussion of how
these frequencies arise to Appendix E.2 and simply quote the results here.
When $\gamma$ is smaller than the critical value $\gamma_{c}\simeq 0.0522$,
the relevant frequencies are,

$\omega_{lr}^{+},\quad\omega_{lr}^{-},\quad\omega_{I1}^{+},\quad\omega_{I1}^{-},\quad\omega_{I2}^{+},\quad\omega_{I2}^{-},$
(112)

where the light-ring frequencies $\omega_{lr}^{\pm}$ satisfy the same
conditions (82) as in shallow water. The other four frequencies are the
primary ($\omega_{I1}^{\pm}$) and secondary ($\omega_{I2}^{\pm}$) inflection
points, (where $\pm$ indicates whether these points are on $\omega_{D}^{\pm}$)
which all satisfy the relations,

$\mathcal{H}_{I}=0,\quad\partial_{p}\mathcal{H}_{I}=0,\quad\partial_{p}^{2}\mathcal{H}_{I}=0.$
(113)

These are the locations where the extrema of $\omega_{D}^{\pm}$ (in the $p$
direction) coalesce before moving off in the complex plane. They therefore
correspond to the locations where the $\omega^{\pm}_{a}$ meet and terminate.
They are important since the value of $\omega$ relative to the $\omega_{I}$
determines which mode is the propagating one moving in the direction of
decreasing $r$.

Above $\gamma_{c}$, each branch of the dispersion relation has a single real
extremum in the $p$ direction. In this regime, $\mathcal{H}$ has no real
inflection points and the transitions between different scattering processes
are governed by a different set of frequencies. In Appendix E.2, we argue that
the important frequencies are complex saddle points of $\mathcal{H}$ which
satisfy,

$\begin{split}\mathrm{Re}[\mathcal{H}_{S}]=&\
0,\quad\mathrm{Im}[\mathcal{H}_{S}]=0,\\\
\mathrm{Re}[\partial_{p}\mathcal{H}_{S}]=&\
0,\quad\mathrm{Im}[\partial_{p}\mathcal{H}_{S}]=0.\end{split}$ (114)

In this case there are four conditions as we need to solve for
$(r_{S},\mathrm{Re}[p_{S}],\mathrm{Im}[p_{S}],\omega_{S})$. There are two
different saddle points for the upper and lower branches of the dispersion
relation which we identify with $\pm$ below. Hence, for $\gamma>\gamma_{c}$,
the characteristic frequencies of the system are,

$\omega_{lr}^{+},\quad\omega_{lr}^{-},\quad\omega_{S}^{+},\quad\omega_{S}^{-}.$
(115)

On Fig. 10, we indicate the frequencies in (112) and (115) in panel (a) and
(d) respectively. Note that since $\omega_{S}^{\pm}$ (and the associated
locations in the $(r,p)$–phase space) correspond to saddle points in the
complex plane, they are not positioned directly on top of
$\omega^{\pm}_{1}(r)$.

Figure 11: $\gamma m^{\frac{4}{3}}/gl^{2}=\mathrm{const}$ slices of the full
parameter space for $\gamma m^{\frac{4}{3}}/gl^{2}=10^{-4}$ (a1), $5\times
10^{-3}$ (b1), $0.01$ (c1) and $0.08$ (d1). The plane is divided up by various
characteristic frequencies into regions where different scattering processes
occur. $\omega^{\pm}_{lr}$ are the curves which are not linear in $C/D$
whereas $\omega^{\pm}_{I1,2}$ are linear in $C/D$, with $\omega^{\pm}_{I1}$
being the steeper of the two pairs. Note that panel (d1) is the only case
where $\gamma>\gamma_{c}$. The colour axis indicates which of 13 processes
occur in each region. We only include region labels for $\omega>0$ in panels
(a1) and (d1) since the two central panels and the $\omega<0$ regions become
obscured by including them. A detailed discussion of the different parameter
space regions, with clearer diagrams, can be found in Appendix E.4. In panels
(a2,a3,b2,b3,…), we display the turning point curves $\omega^{\pm}_{a}$ at
$C/D=5,15$ indicated by the black arrows, where the colour convention and
solid points have the same meaning as in Fig. 10.

In Fig. 11, we give examples of a $(\omega,C,\gamma=\mathrm{const})$ slice of
the parameter space for various $\gamma$, and show how the characteristic
frequencies divide the space into regions where the different scattering
processes occur. At our level of approximation, there is a discontinuous
change in the mode scattering crossing the characteristic frequencies (this is
the same phenomenon that lead to the errorenous jump in the reflection
coefficient on Fig. 4 at the superradiant threshold). In reality, these
discontinuities will be smoothed over by mode scattering away from the turning
points, which we do not account for. Since Fig. 11 does not clearly indicate
all the different regions, we provide a detailed discussion in Appendix E.4 of
how varying $\gamma$ changes $(\omega,C)$ plane, giving the critical
parameters where the plane undergoes a structural change. Due to the symmetry
$(\omega,m,p)\to(-\omega,-m,-p)$, the parameter space for $m<0$ can be
obtained by inverting the $m>0$ parameter space through $\omega=0$, and also
inverting all the scattering diagrams through $\mathrm{Re}[p]=0$. Note that
the arrows point in the same direction upon inversion since
$v_{g,r}=\partial_{p}\omega$ does not change sign under this transformation.

To recap before moving on, in this section we characterized all possible
scattering processes in the deep water system governed by (110). These are
divided in 13 possible mode interactions which are summarized by the diagrams
in Fig. 9. In general, there are two propagating modes at spatial infinity
(right side of the diagrams) and four approaching the vortex axis (left side).
However, they differ in the modes which interact in the intermediate region
and which modes propagate in the core and at infinity. The particular
scattering configuration is then dictated by the dimensionless parameters
$(\omega,C,\gamma)$, as shown in Fig. 11.

###  V.3 Superradiance

Combining knowledge of the various scattering processes (Fig. 9) with where
they occur in the parameter space (Fig. 11), we can now deduce a superradiance
condition for deep water waves. Superradiant amplification occurs for
$\omega>0$ whenever the propagating modes on $\omega_{D}^{+}$ as $r\to\infty$
become propagating modes on $\omega_{D}^{-}$ at low $r$. This can be seen by
noting that when $p^{\pm}$ are on $\omega_{D}^{-}$, $p^{+}$ is the in-going
mode, meaning that one takes the $-$ sign in the exponent of (57) and
$|\mathcal{R}|>1$. In Fig. 9, the clearest signature is that the black part of
the diagram connects to the red part of the diagram, which occurs for the
following cases.

In Type V and IX scattering, all amplified radiation gets scattered to
infinity because the out-going $p^{+}$ mode does not interact with dispersive
$p^{u}$ mode. In Type VI scattering, most of the amplified radiation is
scattered back into the vortex and only a small amount tunnels out. This can
be seen because the out-going $p^{+}$ mode at infinity is connected to the
superradiant part of the diagram by an evanescent (dashed) line. For Type XI,
all amplified radiation is scattered back into the vortex, which is a
consequence of the out-going mode at infinity ($p^{t}$ in this case) not being
connected to the superradiant part of the diagram. For Type XII, most of the
amplified radiation escapes (since the out-going mode is connected to the
superradiant part by propagating lines) and only a small amount tunnels back
into the vortex through the evanescent line on the $p^{u,t}$ modes.

The criterion for amplification (of positive frequency modes) when
$\gamma<\gamma_{c}$ is $\omega<\omega_{I1}^{-}$, i.e. the superradiant
threshold is determined by an inflection point of $\omega_{D}^{-}(p)$.
However, the additional radiation is only guaranteed to reach infinity if
$\omega<\omega_{lr}^{+}$ is also satisfied (note that similar behaviour was
found in the absence of surface tension in [25]). For $\gamma>\gamma_{c}$,
superradiance occurs for $\omega<\omega_{S}^{-}$, i.e. a saddle point in the
complex $p$–plane plays the role of the threshold. Altogether, we can write,

$0<\omega<\omega_{\star}=\begin{cases}\mathrm{min}[\omega_{I1}^{-},\omega_{lr}^{+}],\qquad\gamma<\gamma_{c}\\\
\omega_{S}^{-},\qquad\qquad\qquad\,\gamma>\gamma_{c}\end{cases},$ (116)

as a sufficient requirement for the amplified radiation to escape the vortex.
An interesting consequence of the bound above is that, below a critical
rotation, the superradiant bandwidth closes up and amplification can no longer
occur at any frequency. This happens below (above) $\gamma_{c}$ when
$\omega_{I1}^{-}$ ($\omega_{S}^{-}$) passes through zero. In Fig. 12, we show
the critical rotation at which this happens as a function of the dimensionless
surface tension parameter. As $\gamma\to 0$, the critical value goes to
$2^{3/2}$, which agrees with the analytic prediction for deep water waves in
the absence of capillary effects [25]. In the limit of large $\gamma$, we find
numerically that the critical rotation tends to unity.

Figure 12: The critical rotation is shown as a function of the dimensionless
surface tension parameter. Superradiance cannot occur for rotations slower
than this since modes propagating on $\omega_{D}^{+}$ at infinity do not reach
$\omega_{D}^{-}$ in the vortex core.

###  V.4 Bound states

Superradiance is the key ingredient in the BHB instability. We therefore only
need to analyse the diagrams which satisfy (116) to see when we can have
superradiantly unstable bound states. We do this by looking in the parameter
space of Fig. 11 to find any superradiant diagrams for the chosen parameters,
then identifying the segment of the relevant diagram in Fig. 9 which forms a
bound state between the vortex and the boundary. This means that, starting
from the in-going $p^{-}$ mode on the right edge of the diagram, one should be
able to follow the arrows and end up back on the right edge as an out-going
mode. Once such a path has been identified, one then transports the amplitudes
along the path and imposes the boundary condition. In the present case, we
also account for the $\nu$ dependent attenuation in (43). This procedure is
detailed in full in Appendix F. The result is that the oscillation frequency
$\omega_{n}$ satisfies a Bohr-Sommerfeld resonance condition,

$\begin{split}\mathrm{Re}[S_{\mathrm{path}}(\omega_{n})]=&\
\pi\left(2n+\frac{1}{2}\right),\\\ S_{\mathrm{path}}=&\
\int_{\mathrm{path}}p(r)dr,\end{split}$ (117)

where $S_{\mathrm{path}}$ is the phase integral along the sequence of rays
comprising the path of the bound state through $(r,p)$–phase space. Assuming
$\Gamma_{n}\ll\omega_{n}$, we find that the growth rate is given by,

$\begin{split}&\Gamma_{n}=\frac{\sum_{a}\log|X_{a}|-\nu\mathrm{Re}[\Sigma_{\mathrm{path}}]}{\partial_{\omega}\mathrm{Re}[S_{\mathrm{path}}]}\Bigg{|}_{\omega=\omega_{n}},\\\
&\qquad\Sigma_{\mathrm{path}}=\int_{\mathrm{path}}\frac{k^{2}dr}{v_{g,r}},\end{split}$
(118)

where $\Sigma_{\mathrm{path}}$ is the attenuation due to viscosity (which
comes from the exponent containing $\nu$ in (44)) and the $X_{a}$ are
different scattering coefficients along the path. For Type V and IX scattering
(which exhibit bound state A in Appendix F) there is a single superradiant
reflection off the vortex so that $X_{a}=\mathcal{R}_{1}$. For Type VI
scattering (bound state B), there is a superradiant reflection followed by a
regular transmission so that $X_{a}\in\\{\mathcal{R}_{1},\mathcal{T}_{23}\\}$.
For Type XII scattering (bound state C), there is a superradiant reflection
followed by two regular reflections so that
$X_{a}\in\\{\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{3}\\}$. Upon noting
that $T_{\mathrm{path}}=\partial_{\omega}\mathrm{Re}[S_{\mathrm{path}}]>0$ is
the time taken for a point on the radial wavefront to complete a full orbit of
the trapping region, the expressions in (117) and (118) become equivalent to
those presented in the introduction in (1) and (4) respectively.

The result in (118) demonstrates that the effect of the additional
reflections/transmissions induced by dispersion is to reduce the magnitude of
the instability. Specifically, for bound states B and C we have
$\log|\mathcal{R}_{1}\mathcal{T}_{23}|$ and
$\log|\mathcal{R}_{1}\mathcal{R}_{2}\mathcal{R}_{3}|$ respectively, where
$|\mathcal{R}_{1}|>1$ and all the rest are less than unity. Hence, for an
instability to occur, the amount of energy gained through superradiance must
exceed the decrease in energy within the trapping at each additional
reflection/transmission.

Furthermore, it is clear from (118) that the effect of viscosity is to further
reduce the strength of the instability. We can define a critical value of
$\nu$ for which in the instability is quenched via,

$\nu_{c}(\omega_{n})=\frac{\sum_{a}\log|X_{a}|}{\mathrm{Re}[\Sigma_{\mathrm{path}}]}\bigg{|}_{\omega=\omega_{n}},$
(119)

where this expression only makes sense when $\Gamma_{n}^{\nu=0}>0$. The growth
rate can then be expressed in terms of its value in the absence of viscosity,

$\Gamma_{n}=\Gamma_{n}^{\nu=0}\left(1-\frac{\nu}{\nu_{c}}\right),$ (120)

where the factor in parentheses softens the instability for $\nu<\nu_{c}$ and
completely quenches it for $\nu\geq\nu_{c}$. This representation will be
useful for comparing the undamped BHB modes to the damped ones. Another useful
representation of the damping term is,

$\frac{\mathrm{Re}[\Sigma_{\mathrm{path}}]}{T_{\mathrm{path}}}=\frac{\int_{\mathrm{path}}k^{2}dt}{T_{\mathrm{path}}}=\langle\langle
k^{2}\rangle\rangle_{\mathrm{path}},$ (121)

where we have used $v_{g,r}=dr/dt$ which follows from (39) and the double
angular brackets denote the average value over the path. Thus, modes with
larger $k$ will be more heavily attentuated by viscous damping.

Figure 13: BHB eigenfrequencies for $D=10\,\mathrm{cm^{2}s^{-1}}$ and a range
of $r_{B}$ and $C/D$. We display results for $m=1,2,3$ in black, blue and red
respectively. The $\nu=0$ results are shown as solid points. For $\nu\neq 0$,
the eigenfrequencies have identical real parts to the undamped values but
different imaginary parts through (120). These are illustrated either as solid
lines, or open circles where there is only a single point. The lines terminate
with crosses when the instability is quenched by viscosity (i.e. $\Gamma_{n}$
goes to zero). No line is shown if every unstable mode in the undamped
spectrum is stabilised by damping. The trend is similar to the shallow water
results in Fig. 5 in the sense that larger $r_{B}$ and $C/D$ give rise to more
modes and the growth rates are increased with rotation. However, the
functional dependence on the system’s parameters is quite different from
shallow water.

##  VI Examples

Figure 14: Dependence of the BHB eigenfrequencies on $r_{B}$, $D$ and $C/D$
for $m=1$. Each column represents a different set of system parameters
indicated at the top of the column. The lines are colour coded so that the
larger frequencies appear lighter orange. The oscillation frequencies are
shown in panels (a1,b1,c1), undamped growth rates in panels (a2,b2,c2) and the
factor which multiplies the undamped growth rates when viscosity is non-zero
in panels (a3,b3,c3).

We now compute $\omega^{\mathbb{C}}_{n}=\omega_{n}+i\Gamma_{n}$ for exemplary
parameters to illustrate the effect of dispersion and dissipation on the BHB
instability. From hereon, we revert to using physical units so that our
discussion is relevant to experiments. For water at room temperature, the
viscosity and surface tension parameters are $\nu=8.93\times
10^{-7}\mathrm{m^{2}s^{-1}}$ and $\gamma=7.22\times
10^{-5}\mathrm{m^{3}s^{-2}}$ respectively. We restrict our attention to trap
sizes in the range $r_{B}\in[5,20]\,\mathrm{cm}$ and drain rates in the range
$D\in[10,100]\,\mathrm{cm^{2}s^{-1}}$. For our third variable, we take
$C/D\in[1,10]$ since Fig. 12 shows that this ratio must be at least greater
than unity. Using these variables, we analyse the parameter space in Fig. 11,
finding that the relevant superradiant diagrams are Type V and Type IX. This
restriction comes from the chosen range for $C/D$, since the more exotic
superradiant diagrams occur only at larger $C/D$ ratios. Both diagrams for
Type V and IX are concerned with the same class of bound state (bound state A
in Appendix D) which scatters a single time with the vortex before reflecting
from the boundary. This greatly simplifies the analysis since we only need to
locate the two turning points $r_{{tp}1,2}$ where the $p^{\pm}$ modes scatter
on $\omega_{D}^{-}$ and $\omega_{D}^{+}$ respectively. Solving $\mathcal{H}=0$
for $p^{\pm}(r)$, we then have sufficient information to evaluate all the
phase integrals in (117) and (118). At a given selection of parameters
$(r_{B},D,C/D)$, these integrals are evaluated over the range of superradiant
frequencies inside the trapping region. The upper limit of this range is
$\omega_{\star}$ defined in (116) whilst the lower bound is value of the
turning point curve at the boundary, i.e. $\omega_{1}^{+}(r_{B})$, in Fig. 10.

In the deep water regime, one has to solve (45) numerically for each frequency
in the range (this step was trivial in shallow water where turning points are
given analytically by (80)). To do this, we first find an approximate solution
to $\omega=\omega^{-}_{1}(r_{{tp}1})$ and $\omega=\omega^{+}_{1}(r_{{tp}2})$
graphically, where we can use the fact that
$\omega_{1}^{\pm}(r\to\infty)\sim\pm 1/\sqrt{r}$ to identify which of the
curves we need for $\gamma<\gamma_{c}$ (note that this implies the relation
$(r_{tp},p_{tp})\sim(\omega^{-2},2\omega^{5})$ in the $\omega\to 0$ limit). We
then take these solutions as initial guesses and solve (45) using Matlab’s
inbuilt fsolve function. Once we have the turning points, $\mathcal{H}=0$