# Quantum vortex instability and black hole superradiance

Sam Patrick  Department of Physics and Astronomy, University of British
Columbia, Vancouver, British Columbia, V6T 1Z1, Canada  Institute of Quantum
Science and Engineering, Texas A&M University, College Station, Texas, 77840,
US  August Geelmuyden  School of Mathematical Sciences, University of
Nottingham, University Park, Nottingham, NG7 2RD, UK  Sebastian Erne  Vienna
Center for Quantum Science and Technology, Atominstitut, TU Wien, Stadionallee
2, 1020 Vienna, Austria  Carlo F. Barenghi  Joint Quantum Centre Durham-
Newcastle, School of Mathematics, Statistics and Physics, Newcastle
University, Newcastle upon Tyne, NE1 7RU, UK  Silke Weinfurtner  School of
Mathematical Sciences, University of Nottingham, University Park, Nottingham,
NG7 2RD, UK  Centre for the Mathematics and Theoretical Physics of Quantum
Non-Equilibrium Systems, University of Nottingham, Nottingham, NG7 2RD, UK

###### Abstract

Vortices and black holes set the scene for many interesting dynamical
processes in physics. Here, we study the dynamical instability of quantised
vortices and rotational superradiance around rotating black holes,
illustrating in the process that the same physics is at play in these two
seemingly disparate phenomena. We also compare the instability of the vortex
to the black hole bomb instability, which occurs for massive scalar fields in
the Kerr spacetime. Taking inspiration from the analogy between black hole
bomb modes and the hydrogen spectrum, the vortex instability is compared with
nuclear resonances involved in $\alpha$-decay.

##  I Introduction

It is a remarkable fact that phenomena in completely separate areas of physics
can admit the same theoretical description. Historically, the pervasiveness of
certain descriptions has served not only as a tool to aid our understanding of
known phenomena, but also as a driving force in the development of new
theories. A prominent example is Maxwell’s use of the fluid equations as a
“vehicle of mathematical reasoning” in his development of classical
electrodynamics [1]. Another example is Bekenstein’s identification of the
similarity between the laws of thermodynamics and black hole (BH) physics [2,
3], which ultimately led to Hawking’s discovery of black hole evaporation [4].
Theoretical similarities are often indicative that the same conceptual
mechanism underpins different physical processes. In this paper, we
investigate such a similarity between the instability of quantum vortices in
superfluids and the superradiant amplification of waves around rotating BHs.

In recent years, there has been a growing body of work in the literature
drawing comparisons between hydrodynamical vortices and rotating BHs [5, 6, 7,
8, 9, 9, 10, 11, 12, 13]. This is largely due to the success of the analogue
gravity programme, which aims to shed light on processes that occur for
(quantum) fields in curved spacetime by studying their counterparts in
condensed matter physics [14]. The basic notion at the heart of analogue
gravity is that the excitations of many fluid-like systems obey (in the non-
dispersive regime) an equation which is formally equivalent to that of a
Klein-Gordon field in a curved spacetime [15]. This correspondence has led to
the detection of analogue Hawking radiation [16, 17, 18, 19, 20] and
superradiance [21] in the laboratory. Here, we use the term superradiance to
refer to rotational superradiance (commonly studied in BH physics [22]) as
distinct from Dicke superradiance in quantum optics [23]. In view of the
experiments of [21], which were performed in the dispersive regime, two of the
authors have theoretically demonstrated that superradiance persists well-
beyond the regime of the correspondence with the Klein-Gordon equation [24,
25], highlighting the possibility of studying BH processes in a much wider
range of scenarios than previously thought possible. The mechanism which
enables superradiance is the existence of a region where excitations of the
system can have negative energies. Such a region is commonly called an
ergoregion, as it is a region where excitations do work on the system to
extract energy from it. The negative energy of such excitations then simply
captures the fact that the energy of the system is being lowered.

A promising system within which to pursue this analogy involves flowing Bose-
Einstein condensates (BECs). One of the most striking features of BECs (and
superfluids in general) is that circulation is quantised [26]. This has been
known since the predictions of Onsager and Feynman in the mid-20th century
[27, 28] and has been researched heavily in the decades since, see e.g. [29]
and references therein. From the BH perspective, these systems are interesting
since they provide a way to simulate effective spacetimes with quantum
features. From the condensed matter point of view, the analogy with BHs is
useful since it can offer new perspectives on known BEC phenomena.

It has long been recognised that multiply-quantised (also described as highly-
charged or large) vortices are unstable and will decay into clusters of
singly-quantised vortices. This decay occurs in dissipative systems since
splitting into a cluster lowers the total energy [30], however, it can also
occur in non-dissipative systems if the vortex is able to transfer energy into
sound waves [31]. This instability has been confirmed both numerically and
experimentally [32, 33, 34]. Interestingly, it has been argued that multiply-
quantised vortices in condensates can be stabilised by a draining component in
the flow field [35, 36]. In general, a deeper understanding of vortex dynamics
is desirable due to the implications for on-going research on (e.g.) wave-
vortex interactions [37] and quantum turbulence [38].

The connection between vortex instabilities and superradiance (or more
precisely ergoregions) has been alluded to on several occasions in the
literature. Early studies focussed on the non-dispersive regime and ad hoc
boundary conditions (BCs) were imposed near the vortex axis, where the fluid
density drops to zero, to circumvent the breakdown of the hydrodynamic
approximation there [39, 40]. These works described this instability as an
ergoregion instability, a term coined in the 1970s for describing
instabilities of rapidly rotating neutron stars [41]. The ergoregion
instability has also been referred to as an “inverted BH bomb” [42], where the
BH bomb is a similar superradiant instability which results from the trapping
of amplified modes outside the ergoregion. This occurs, for example, for
massive scalar fields in the Kerr spacetime [43]. Recently, it has been shown
in [44] that the link between the ergoregion and vortex instabilities also
carries over to the dispersive regime. Specifically, [44] demonstrated that
the vortex instability originates in the vortex core where the mode has
negative energy, thereby identifying the phenomenon as an ergoregion
instability and not a BH bomb. In this work, we provide further insight into
the link between vortex instabilities, ergoregions and superradiance by
illustrating how each can be analysed within the same framework (introduced
below). We also perform an explicit comparison between the ergoregion and BH
bomb instabilities.

##  II Methods

Our main tool to study the physics underlying vortex instabilities and
superradiance will be based on the dispersion relation. To obtain the
dispersion relation, one assumes that the solutions to the linear equations of
motion look locally like plane waves, which amounts to a Wenzel-Kramers-
Brillouin (WKB) approximation. The essence of the WKB approximation is that
the waves vary rapidly when compared with the changes in the background they
move through, allowing one (to leading order) to obtain a relation between the
local wavelength and the frequency [45, 46]. At next to leading order, one
obtains a correction to this plane wave approximation in the form of a
spatially varying amplitude. At certain locations in the system, the amplitude
diverges and the approximation breaks down, which signals the interaction of
the different modes of the dispersion relation. By solving the wave equation
locally and matching asymptotically onto the WKB modes, one obtains a transfer
matrix which relates the mode amplitudes in different regions of the system.
Armed with the transfer matrices, it is possible to estimate scattering
coefficients and resonant frequencies.

The benefit of this method is that (1) one has precise control of the
locations in the system where scattering and trapping of the waves take place,
(2) it can be applied with relative ease even when conventional exact methods
fail and (3) it converges to the exact solution in the limit of vanishingly
small wavelengths (which corresponds in this work to high frequencies). The
first point is appealing to us here since the WKB method will allow us to
pinpoint precisely where the trapping of unstable modes occurs. The second is
also important since BH instabilities are usually computed using a continued
fraction algorithm [43] which, as of yet, has not been reformulated to apply
to dispersive systems. The WKB method, however, can easily be adapted to
handle dispersion, e.g. [47, 48, 24, 25]. The final point is relevant since
the convergence of WKB to the exact solution often happens so quickly that
approximation yields acceptable results even for relatively long wavelength
modes. Indeed, we will see examples in this work where WKB results agree with
more accurate methods to within a couple of percent, even for low frequencies.

##  III Outline and main results

We now present an outline of the paper along with a summary of our main
results. In Section IV, we study the instability of multiply-quantised
vortices. After establishing the preliminaries, we show how the WKB method can
be applied to study wave scattering in Section IV.5. Then, in Section IV.8, we
compute the unstable frequencies of a vortex in an infinite system using the
WKB method. This analysis reveals that the instability is a superradiant bound
state in the ergoregion, i.e. a trapped negative energy excitation.
Equivalently, one can think of this as a trapped positive frequency mode with
negative norm (the norm is the particle number) which turns out to be a more
useful notion for our discussion. The validity of these results is confirmed
by a comparison with full numerical simulations of the linearised equations.
Since real experiments are performed in finite sized systems, we show in
Section IV.9 how our WKB formula for the instabilities is modified when a
reflective boundary is placed far from the vortex core. An important result
here is that the WKB instability condition reduces exactly to the open system
condition when taking the outer boundary to infinity, thereby validating
existing approaches in the literature [31, 44] which inferred stability
properties of open systems by considering limits of closed systems. In our
companion paper [49], we study the instability in the non-linear regime and
observe the subsequent decay of a doubly-quantised vortex into a pair of phase
singularities. There, we also uncover novel behaviour of the system at late
times: the system does not completely decay into two well-separated vortices
but instead forms a bound state of phase defects orbiting in close proximity,
whilst the separation between them undergoes a series of modulations.

Having identified the mechanism underpinning the vortex instabilities, next we
study a non-dispersive draining vortex in Section V.1, which is an analogue
rotating BH. This serves three main purposes. Firstly, we demonstrate in
Section V.1.2 that BH superradiance is also the result of positive frequency
modes tunnelling to a region where they have negative norm, i.e. the same
process that causes the vortex instabilities. Secondly, the analogy to BHs
allows us to identify the negative norm region of the system as the
ergoregion, thereby justifying the name ergoregion instability for the
unstable vortex modes. Thirdly, we show in Section V.1.3 that the addition of
a drain to a vortex in an infinite system can stabilise the system (at least
in the non-dispersive regime).

Finally, we study the BH bomb instability of a massive scalar field in the
Kerr spacetime. After introducing the system in Section V.2, in Section V.3 we
apply the WKB method to demonstrate the existence of instabilities, showing in
the process how our method can be used to recover the well-known hydrogenic
approximation for the instability spectrum [43]. We then explicitly compare
the two confinement mechanisms for the BH bomb and vortex/ergoregion
instabilities. Taking inspiration from the hydrogen analogy, we conclude by
arguing that the vortex instability shares common features with the nuclear
resonances involved in $\alpha$-decay.

##  IV Quantised vortices

###  IV.1 Gross-Pitaevskii equation

Consider a two-dimensional condensate of bosons with nearly zero temperature.
In the limit of many bosons, the system is well-described by the mean-field
condensate wavefunction $\Psi(t,\mathbf{x})$, where $t$ is time and
$\mathbf{x}$ are spatial coordinates on a 2D plane [26]. The action for the
system can be written as,

$\begin{split}\mathcal{S}=&\ \int
dtd^{2}\mathbf{x}\bigg{[}\frac{i\hbar}{2}\left(\Psi^{*}\partial_{t}\Psi-\Psi\partial_{t}\Psi^{*}\right)\\\
&\ \
-\frac{\hbar^{2}}{2M}\bm{\nabla}\Psi\cdot\bm{\nabla}\Psi^{*}-U(\mathbf{x})|\Psi|^{2}-\frac{1}{2}g|\Psi|^{4}\bigg{]},\end{split}$
(1)

where $g$ is the 2D interaction parameter, $M$ is the mass of bosons in the
condensate and $U$ is the trapping potential which depends on the spatial
coordinate. The square of the wavefunction has the interpretation of the
particle density $\rho=|\Psi|^{2}$ and its integral gives the total number of
atoms in the condensate, i.e. $N=\int d^{2}\mathbf{x}|\Psi|^{2}$. In the
ground state, the condensate oscillates at the chemical potential $\mu$, i.e.
$\Psi\sim e^{-i\mu t}$, where $\mu$ sets the typical energy scale of the
system. To lighten notation, we can define the following parameters,

$\xi=\hbar/\sqrt{M\mu},\qquad\tau=\hbar/\mu,\qquad\rho_{c}=\mu/g,$ (2)

where $\xi$ is the healing length, $\tau$ is the healing time and $\rho_{c}$
is a reference density. We then rescale all quantities in (1) by these
parameters. Setting $\mathbf{x}/\xi\to\mathbf{x}$, $t/\tau\to t$,
$\Psi/\sqrt{\rho_{c}}\to\Psi$, $U/\mu\to U$ and
$\mathcal{S}/\hbar\to\mathcal{S}$, which amounts to setting $\hbar=M=\mu=g=1$,
we obtain the action in dimensionless form,

$\begin{split}\mathcal{S}=&\ \int
dtd^{2}\mathbf{x}\bigg{[}\frac{i}{2}\left(\Psi^{*}\partial_{t}\Psi-\Psi\partial_{t}\Psi^{*}\right)\\\
&\qquad-\frac{1}{2}\bm{\nabla}\Psi\cdot\bm{\nabla}\Psi^{*}-U(\mathbf{x})|\Psi|^{2}-\frac{1}{2}|\Psi|^{4}\bigg{]}.\end{split}$
(3)

Variation of $\mathcal{S}$ with respect to $\Psi$ yields the dimensionless
Gross-Pitaevskii equation (GPE),

$i\partial_{t}\Psi=-\frac{1}{2}\nabla^{2}\Psi+\left[U(\mathbf{x})+|\Psi|^{2}\right]\Psi,$
(4)

which is a non-linear Schrödinger equation. Instead varying (3) with respect
to $\Psi^{*}$, i.e. the complex conjugate of the wavefunction, yields the
complex conjugate of (4).

###  IV.2 Hydrodynamic variables

The GPE can be put into hydrodynamic form by writing $\Psi$ in the Madelung
representation,

$\Psi=\sqrt{\rho}e^{i\Phi}.$ (5)

The GPE then splits into two separate equations,

$\displaystyle\partial_{t}\Phi+\frac{1}{2}(\bm{\nabla}\Phi)^{2}+\rho-\frac{1}{2}\frac{\nabla^{2}\sqrt{\rho}}{\sqrt{\rho}}=$
$\displaystyle\ 0,$ (6)
$\displaystyle\partial_{t}\rho+\bm{\nabla}\cdot(\rho\bm{\nabla}\Phi)=$
$\displaystyle\ 0,$ (7)

and $\mathbf{v}=\bm{\nabla}\Phi$ is identified as the velocity field of the
flowing condensate. The first of these is similar to the Bernoulli equation
from classical fluid dynamics, except for the final term which referred to as
the quantum pressure. This term acts on scales on the order of the healing
length. The second equation is simply a continuity equation for boson density.
Together, (6) and (7) give a hydrodynamic interpretation to the flowing
condensate.

###  IV.3 Vortex solutions

Let us now specialise to a stationary condensate in polar coordinates
$(r,\theta)$ with the wavefunction,

$\Psi(r,\theta,t)=Y(r)e^{i\ell\theta-it},$ (8)

where the $e^{-it}$ dependence is the oscillation due to the chemical
potential (which is unity in our dimensionless variables). Since
$\Psi(\theta=0)=\Psi(\theta=2\pi)$ by continuity of the wavefunction, the
winding number $\ell$ is constrained to be an integer. Assuming no flow in the
radial direction (which amounts to taking $Y$ to be real) the velocity field
is simply,

$\mathbf{v}=\frac{\ell}{r}\vec{\mathbf{e}}_{\theta}.$ (9)

This velocity profile is known as a vortex solution and, since $\ell$ can only
take on integer values, the vortex circulation is said to be quantised. In
this work, we will take $\ell>0$ so that our vortex rotates in the direction
of positive $\theta$, i.e. anti-clockwise.

The corresponding ground state density profile in an unbounded condensate,
i.e. with $U=0$, is obtained by solving,

$\frac{1}{2r}\partial_{r}\left(r\partial_{r}Y\right)+\left(1-\frac{\ell^{2}}{2r^{2}}\right)Y=Y^{3},$
(10)

subject to the BCs $Y(r=0)=0$ and $Y(r\to\infty)=1$. Then the density is
obtained through $\rho=Y^{2}$ and the parameter $\rho_{c}$ introduced in (2)
is identified as the density of the vortex at spatial infinity. In Fig. 1, we
show the numerically obtained density profiles corresponding to the first six
values of $\ell$. The details of the numerical solver can be found in Appendix
A. By noting that (10) asymptotes to Bessel’s equation in the limit $r\to 0$
and a polynomial equation for $r\to\infty$, we can write down the asymptotics
for the density,

$\rho(r)\sim\begin{cases}\rho_{0}\,J^{2}_{\ell}(\sqrt{2}r),\qquad r\to 0,\\\
1-\ell^{2}/2r^{2},\ \qquad r\to\infty.\end{cases}$ (11)

where $\rho_{0}$ is a constant amplitude.

Figure 1: The density profiles $\rho(r)$ for the first 6 non-zero $\ell$. It
is apparent that, the larger the winding number $\ell$ is, the larger the
vortex core is at a prescribed level of density.

###  IV.4 Fluctuations

Now we consider fluctuations of the condensate by writing $\Psi\to\Psi+\psi$
with $\psi\ll\Psi$. Linearising the GPE (4) as well as its complex conjugate,
the linear equations can be collected into a single matrix equation called the
Bogoliubov-de Gennes (BdG) equation [44],

$i\partial_{t}|\psi\rangle=\widehat{L}|\psi\rangle,$ (12)

where $\psi$ has been split according to,

$\begin{bmatrix}\psi(t,r,\theta)\\\
\psi^{*}(t,r,\theta)\end{bmatrix}=\sum_{m=-\infty}^{\infty}e^{im\theta}\begin{bmatrix}u_{+}(t,r)e^{+i\ell\theta-
it}\\\ u_{-}(t,r)e^{-i\ell\theta+it}\end{bmatrix},$ (13)

and the operator $\widehat{L}$ and the state vector $|\psi\rangle$ are defined
as,

$\widehat{L}=\begin{bmatrix}\mathcal{D}_{+}&\rho\\\
-\rho&-\mathcal{D}_{-}\end{bmatrix},\qquad|\psi\rangle=\begin{bmatrix}u_{+}\\\
u_{-}\end{bmatrix}$ (14)

with,

$\mathcal{D}_{\pm}=-\frac{1}{2}\left[\partial_{r}^{2}+\frac{1}{r}\partial_{r}-\frac{(m\pm\ell)^{2}}{r^{2}}\right]+2\rho+U-1.$
(15)

The BdG conserves the norm, which measures the particle content of a given
mode and is defined as,

$\mathcal{N}=\langle\psi|\sigma_{3}|\psi\rangle=\int
d^{2}\mathbf{x}\left(|u_{+}|^{2}-|u_{-}|^{2}\right),$ (16)

where $\sigma_{3}$ is the third Pauli matrix. $\partial_{t}\mathcal{N}=0$ then
follows by inserting the BdG relation in (12) and invoking $u_{\pm}=0$ at
either the system’s outer boundary or as $r\to\infty$.

Alternatively, we can instead write down the linear equations in terms of
fluctuations of the hydrodynamic variables defined in Section IV.2 by
perturbing $\Phi\to\Phi+\phi$ and $\rho\to\rho\,(1+\eta)$. These variables are
related to the perturbed wavefunction via,

$\psi=\Psi\left(\tfrac{1}{2}\eta+i\phi\right),$ (17)

and obey the following set of coupled equations,

$\displaystyle D_{t}\eta+\rho^{-1}\bm{\nabla}\cdot(\rho\bm{\nabla}\phi)$
$\displaystyle\ =0,$ (18a) $\displaystyle
D_{t}\phi+\rho\eta-\tfrac{1}{4}\rho^{-1}\bm{\nabla}\cdot\left(\rho\bm{\nabla}\eta\right)$
$\displaystyle\ =0,$ (18b)

where $D_{t}=\partial_{t}+\mathbf{v}\cdot\bm{\nabla}$ is the material
derivative. Of course, (12) and (18) are completely equivalent to one another,
however, both representations will be used at different points in this work
depending on which best suits our needs.

It is useful to separate the two real fields $\phi$ and $\eta$ into positive
and negative frequency components,

$\displaystyle\phi=$ $\displaystyle\
\sum_{\lambda}\nolimits\left(\alpha_{\lambda}\varphi_{\lambda}+\alpha_{\lambda}^{*}\varphi_{\lambda}^{*}\right),$
(19a) $\displaystyle\eta=$ $\displaystyle\
\sum_{\lambda}\nolimits\left(\alpha_{\lambda}n_{\lambda}+\alpha_{\lambda}^{*}n^{*}_{\lambda}\right),$
(19b)

where the $\alpha_{\lambda}$ are constants and we have used the $t$ and
$\theta$ independence of the background to write,

$\begin{bmatrix}\varphi_{\lambda}\\\
n_{\lambda}\end{bmatrix}=\begin{bmatrix}\tilde{\varphi}(r)\\\
\tilde{n}(r)\end{bmatrix}e^{im\theta-i\omega t}.$ (20)

The sum over $\lambda$ is shorthand for,

$\sum_{\lambda}\nolimits=\int^{\infty}_{0}d\omega\sum_{m=-\infty}^{\infty}.$
(21)

Since the equations are linear by definition, each mode satisfies (18)
independently. The mode equations can be derived from the following
Lagrangian,

$\begin{split}\mathcal{L}=&\
\rho\bigg{[}\frac{1}{2}\left(\varphi_{\lambda}^{*}D_{t}n_{\lambda}+\varphi_{\lambda}D_{t}n^{*}_{\lambda}-n^{*}_{\lambda}D_{t}\varphi_{\lambda}-n_{\lambda}D_{t}\varphi^{*}_{\lambda}\right)\\\
&\qquad-\rho\,n_{\lambda}n_{\lambda}^{*}-\bm{\nabla}\varphi_{\lambda}\cdot\bm{\nabla}\varphi^{*}_{\lambda}-\frac{1}{4}\bm{\nabla}n_{\lambda}\cdot\bm{\nabla}n^{*}_{\lambda}\bigg{]}\\\
&\qquad-\frac{1}{2}\bm{\nabla}\cdot\left[\bm{\nabla}\rho\left(\varphi_{\lambda}\varphi_{\lambda}^{*}+\frac{1}{4}n_{\lambda}n_{\lambda}^{*}\right)\right].\end{split}$
(22)

This is invariant under phase rotations, which results in a conservation
equation for the norm $\mathcal{N}$, whose components are,

$\begin{split}\rho_{n}[\varphi]=&\
i\rho\left(\varphi_{\lambda}n_{\lambda}^{*}-\varphi_{\lambda}^{*}n_{\lambda}\right),\\\
\mathbf{j}_{n}[\varphi]=&\
i\rho\big{[}\mathbf{v}\left(\varphi_{\lambda}n_{\lambda}^{*}-\varphi_{\lambda}^{*}n_{\lambda}\right)+\varphi_{\lambda}\bm{\nabla}\varphi_{\lambda}^{*}\\\
&-\varphi_{\lambda}^{*}\bm{\nabla}\varphi_{\lambda}+\tfrac{1}{4}\left(n_{\lambda}\bm{\nabla}n_{\lambda}^{*}-n_{\lambda}^{*}\bm{\nabla}n_{\lambda}\right)\big{]},\end{split}$
(23)

and satisfies,

$\partial_{t}\rho_{n}[\varphi]+\bm{\nabla}\cdot\mathbf{j}_{n}[\varphi]=0.$
(24)

Using the relations $u_{\pm}=\rho^{-\frac{1}{2}}(\frac{1}{2}n_{\lambda}\pm
i\varphi_{\lambda})\alpha_{\lambda}e^{-i\omega t}$, one recovers the
expression for the total norm in (16) by integrating the norm density over
space, i.e. $\mathcal{N}=\int\rho_{n}\,d^{2}\mathbf{x}$. Time translational
invariance of (22) gives rise to the conservation of mode energy. The
conservation equation then contains the components $\rho_{e}$ and
$\mathbf{j}_{e}$ where (for stationary set-ups) the relation to the components
of the norm are $\rho_{e}=\omega\rho_{n}$ and
$\mathbf{j}_{e}=\omega\mathbf{j}_{n}$. This reinforces the idea of the norm as
the particle number, since the each particle in the mode contributes an energy
proportional to $\omega$. For our purposes, the norm density $\rho_{n}$ will
be particularly useful in identifying the mechanism behind vortex
instabilities and how they relate to superradiance.

###  IV.5 WKB method

Since the equations in (18) lack a closed form solution, we study them
approximately using the WKB method. The fundamental principle underlying the
WKB method is the assumption that $\phi$ and $\eta$ look locally like plane
waves. That is, the amplitude of the waves varies slowly when compared with
the phase or, equivalently, the background varies over much longer scales that
the local wavelength. In a homogeneous medium, plane waves are exact solutions
and the WKB method therefore becomes exact. The equations of motion then
determine a relation between the wavevector $\mathbf{k}$ and the frequency
$\omega$ of the waves called the dispersion relation. In an inhomogeneous
medium, the dispersion relation appears as the leading order solution in the
WKB expansion and gives the relation between the local values of $\mathbf{k}$
and $\omega$. Therefore, since the vortex profile in (9) depends only on the
radial coordinate, the dispersion relation will determine the local value of
the radial component of the wavevector for a given choice of $\omega$ and $m$.
As shown in e.g. [24, 25], one can then use the dispersion relation to study
the scattering of the different modes in the system. This will be the aim of
the present section.

Before applying the WKB method, let us first write down the equations
satisfied by the radial eigenfunctions $\tilde{\varphi}$ and $\tilde{n}$.
Plugging (20) into (18) results in the following coupled set of equations,

$\displaystyle-i\Omega\tilde{n}+\frac{1}{\rho r}\partial_{r}\left(\rho
r\partial_{r}\tilde{\varphi}\right)-\frac{m^{2}}{r^{2}}\tilde{\varphi}=$
$\displaystyle\ 0,$ (25a)
$\displaystyle-4i\Omega\tilde{\varphi}+4\rho\tilde{n}-\frac{1}{\rho
r}\partial_{r}\left(\rho
r\partial_{r}\tilde{n}\right)+\frac{m^{2}}{r^{2}}\tilde{n}=$ $\displaystyle\
0,$ (25b)

where $\Omega$ is the frequency in the fluid frame,

$\Omega=\omega-\frac{m\ell}{r^{2}}.$ (26)

However, the naive application of the WKB approximation to the equations in
(25) does not yield the correct mode behaviour for the modes at the origin.
The simplest way to see what goes wrong is to define a new coordinate $y=\log
r$ and rewrite (25) as,

$\displaystyle\tilde{\varphi}^{\prime\prime}+(\log\rho)^{\prime}\tilde{\varphi}^{\prime}-m^{2}\tilde{\varphi}=$
$\displaystyle\ i\Omega r^{2}\tilde{n},$ (27a)
$\displaystyle\tilde{n}^{\prime\prime}+(\log\rho)^{\prime}\tilde{n}^{\prime}-m^{2}\tilde{n}=$
$\displaystyle\ 4\rho r^{2}\tilde{n}-4i\Omega r^{2}\tilde{\varphi},$ (27b)
where ${}^{\prime}=\partial_{y}$.

Then using the asymptotics of $\rho$, we recognise that
$(\log\rho)^{\prime}\to 2\ell$ approaching the origin which is comparable to
$m$. Since $m$ sets the scale on which the phase of the waves changes, the
assumption of a separation of scales between the waves and the background
(which is a requirement of WKB) is badly violated. To proceed, we can factor
out the offending term by writing,

$\tilde{\varphi}=\rho^{-\frac{1}{2}}f,\qquad\tilde{n}=\rho^{-\frac{1}{2}}g,$
(28)

which reduces the equations in (27) to,

$\displaystyle f^{\prime\prime}-\widetilde{m}^{2}f=$ $\displaystyle\ i\Omega
r^{2}g,$ (29a) $\displaystyle g^{\prime\prime}-\widetilde{m}^{2}g=$
$\displaystyle\ 4\rho r^{2}g-4i\Omega r^{2}f,$ (29b)

where we have used the equation of motion for the density in (10) to write,

$\widetilde{m}^{2}=m^{2}+\ell^{2}+2r^{2}(\rho-1).$ (30)

Since the factor of $\rho^{-\frac{1}{2}}$ in (28) cancels that contained in
$\Psi$ in (17), the new equations in (29) can be viewed as an intermediate
step in converting the BdG equation (12) into the hydrodynamic equations (18).
The benefit of this rewriting is that there are no derivatives of any
background quantities which become large approaching the origin. Now we can
write a WKB ansatz for the waves using the $y$ coordinate,

$f(y)=\mathcal{A}(y)e^{i\int q(y)dy},\quad g(y)=\mathcal{B}(y)e^{i\int
q(y)dy},$ (31)

where $\mathcal{A}^{\prime}\ll q\mathcal{A}$, $\mathcal{B}^{\prime}\ll
q\mathcal{B}$ and $q^{\prime}\ll q^{2}$. Inserting into (29), the leading
order contributions can be combined into a single equation,

$4\Omega^{2}r^{4}=4\rho(q^{2}+\widetilde{m}^{2})+(q^{2}+\widetilde{m}^{2})^{2}.$
(32)

Then defining $q=pr$ and $k^{2}=p^{2}+\widetilde{m}^{2}/r^{2}$, this equation
can be rewritten as a dispersion relation,

$\Omega^{2}=F(k)k^{2},\qquad F(k)=\rho+\tfrac{1}{4}k^{2},$ (33)

which has two separate branches,

$\omega_{D}^{\pm}=\frac{m\ell}{r^{2}}\pm\sqrt{F(k)k^{2}},$ (34)

where modes with $\omega=\omega_{D}^{+}$ are said to be on the upper branch
and those with $\omega=\omega_{D}^{-}$ are on the lower branch. We will
shortly see that modes on the upper branch are the positive norm solutions
whereas those on the lower branch have negative norm.

Since (33) is a quartic polynomial in $p$, it has four different solutions for
$p$ which we call,

$\displaystyle p^{\pm}(r)=$ $\displaystyle\ \pm\sqrt{-W_{+}},\qquad
p^{t,b}(r)=\mp i\sqrt{W_{-}},$ (35a) $\displaystyle W_{\pm}=$ $\displaystyle\
\mp 2\sqrt{\rho^{2}+\Omega^{2}}+2\rho+\widetilde{m}^{2}/r^{2},$ (35b)

The solutions $p^{t,b}$ are evanescent for all $r$ with $ip^{t}>0$
representing the mode which grows with increasing $r$ and $ip^{b}<0$ the one
which decays. The $p^{\pm}$ solutions correspond to in and out-going plane
waves at large $r$. The direction of propagation of these modes is determined
by the radial component of the group velocity,

$v^{r}_{g}=\partial_{p}\omega_{D}^{\pm}=\frac{\rho+k^{2}/2}{\Omega}p.$ (36)

Now that the phase has been determined, it remains to find an expression for
the amplitudes. In Appendix B, we show that the next to leading order equation
gives an equation for the amplitude which can be put in the form,

$\partial_{r}\left(rQ\mathcal{A}^{2}\right)=0,$ (37)

with,

$Q=F^{-1}\Omega v_{g}^{r}=\frac{\rho+k^{2}/2}{\rho+k^{2}/4}\,p,$ (38)

which we can solve for,

$\mathcal{A}\propto\left|rQ\right|^{-\frac{1}{2}}.$ (39)

The amplitude $\mathcal{B}$ can then be found from the relation
$F\mathcal{B}=i\Omega\mathcal{A}$, which is the leading order equation from
(18b). It is also useful to define two new amplitudes $A$ and $B$ which
satisfy,

$A=\sqrt{r}\mathcal{A},\qquad B=\sqrt{r}\mathcal{B},$ (40)

where the $\sqrt{r}$ factors out the part of the amplitude which grows when an
in-going mode gets focussed onto a smaller disk.

Putting everything together, the original $\varphi_{\lambda}$ and
$n_{\lambda}$ modes are given by,

$\displaystyle\varphi_{\lambda}\sim$ $\displaystyle\ \left|r\rho F^{-1}\Omega
v_{g}^{r}\right|^{-\frac{1}{2}}e^{i\int p(r)dr+im\theta-i\omega t},$ (41a)
$\displaystyle n_{\lambda}\sim$ $\displaystyle\ i\left|r\rho
F\Omega^{-1}v_{g}^{r}\right|^{-\frac{1}{2}}e^{i\int p(r)dr+im\theta-i\omega
t}.$ (41b)

Plugging these modes into the expression for the norm in (23), we find,

$\rho_{n}[\varphi_{\lambda}]=2\rho F^{-1}\Omega|\varphi_{\lambda}|^{2}.$ (42)

Recognising that $\Omega>0$ ($\Omega<0$) on the upper (lower) branch of the
dispersion relation by definition, and also that $F$ is everywhere positive
for propagating modes, we confirm the earlier statement that modes of the
upper branch have positive norm whilst those on the lower branch have negative
norm. Thus, $\omega_{D}^{\pm}$ will henceforth be referred to as the
positive/negative norm branches. When a positive frequency mode is on the
negative norm branch, the mode will have negative energy. We will see later on
(Sec. V.1.2) how the region where this occurs can be naturally identified
using our framework with the ergoregion, a notion commonly invoked BH physics
[22].

###  IV.6 Scattering

Scattering occurs when one of the modes in (35a) gets converted into one of
the others, which results in waves being (partially) reflected due to the
inhomogeneity of the medium they propagate through. The WKB approximation
encodes no information about scattering, since it assumes that each mode
evolves adiabatically as it moves through the flow. Hence, the dominant
contributions to scattering come from the locations in the flow where the WKB
approximation fails. An important example of where this occurs is at turning
points in the radial direction, $r_{tp}$, which are the locations where a wave
radially stagnates, i.e. $v_{g}^{r}=0$, before reversing its direction. From
(36), we can see that this occurs when the local value of $p$ vanishes thereby
causing the WKB amplitude to suddenly diverge, violating the slowly varying
assumption.

From the expressions for the four $p$ in (35a), we can see immediately that
nowhere in the flow will $p^{t,b}=0$ be satisfied, since $W_{-}>0$ everywhere.
Hence, the $p^{t,b}$ modes do not scatter. By contrast, since $W_{+}$ can be
either positive or negative, it is possible to have $p^{\pm}=0$ and so these
mode can scatter. On one side of $r_{tp}$, $W_{+}<0$ and the modes correspond
to propagating solutions whereas on the other side, $W_{+}>0$ and the modes
will be evanescent. To avoid a naming ambiguity, let the $p^{\pm}$ be labelled
$p^{\uparrow\downarrow}$ when they are evanescent: $p^{\uparrow}$ is the one
which grows in the direction of increasing $r$ whereas $p^{\downarrow}$ is the
one which decays. The relation between the different mode amplitudes on either
side of $r_{tp}$ can be found from an asymptotic matching procedure, see e.g.
[24, 48] for relevant examples, and the results concisely expressed using the
following transfer matrices,

$T=e^{\frac{i\pi}{4}}\begin{bmatrix}1&-\frac{i}{2}\\\
-i&\frac{1}{2}\end{bmatrix},\qquad\widetilde{T}=e^{\frac{i\pi}{4}}\begin{bmatrix}\frac{1}{2}&-\frac{i}{2}\\\
-i&1\end{bmatrix}.$ (43)

These matrices are defined so that they act on the mode amplitudes at
$r>r_{tp}$ and return the amplitudes at $r<r_{tp}$ according to the following
relations,

$\begin{bmatrix}A^{+}_{tp}\\\
A^{-}_{tp}\end{bmatrix}=T\begin{bmatrix}A^{\downarrow}_{tp}\\\
A^{\uparrow}_{tp}\end{bmatrix},\qquad\begin{bmatrix}A^{\uparrow}_{tp}\\\
A^{\downarrow}_{tp}\end{bmatrix}=\widetilde{T}\begin{bmatrix}A^{+}_{tp}\\\
A^{-}_{tp}\end{bmatrix},$ (44)

where the labels $+,-,\uparrow,\downarrow$ correspond to the ones on the
different $p$. One can then construct a global solution to the equations of
motion by using the WKB solutions in regions of the flow separated by the
$r_{tp}$, then patching these solutions together using the transfer matrices.

There exists a useful representation that we can use to see where the turning
points are located for different frequencies. Plugging $p=0$ into (33) defines
two curves,

$\omega^{\pm}(r)=\frac{m\ell}{r^{2}}\pm\sqrt{\rho\frac{\widetilde{m}^{2}}{r^{2}}+\frac{\widetilde{m}^{4}}{4r^{4}}},$
(45)

which give the frequency required for a mode to be on a turning point at a
given radius. Examples of these curves are displayed in Fig. 2. Thus, we can
find the locations of the turning points by looking for the intersections of a
line with $\omega=\mathrm{const}$ with the functions $\omega^{\pm}$. In
particular, when we have $\omega=\omega^{+}$ there will be a turning point on
the positive norm branch of the dispersion relation, whereas for
$\omega=\omega^{-}$ we have a turning point on the negative norm branch. That
means that if a particular frequency intersects with both the $\omega^{\pm}$
curves, there will be tunnelling between positive and negative norm branches
of the dispersion relation. This is the mechanism responsible for all the
phenomena to be discussed in this work.

Figure 2: The functions $\omega^{+}$ (solid black line) and $\omega^{-}$
(dashed black line) with $\ell=2$ for the first four $m>0$ values. Modes above
$\omega^{+}$ have positive norm whereas those below $\omega^{-}$ have negative
norm. Between the two curves (in the pink region) modes are evanescent, hence
do not propagate. In the frequency band between the horizontal red lines,
waves can be trapped in a cavity in the vortex core (shaded in turquoise) and
become (quasi-)bound states. The possibility of there being modes with
positive norm at infinity but negative norm in the cavity is what leads to the
vortex instability. Inverting the vertical axis about $\omega=0$ and swapping
the labels $\omega^{+}\leftrightarrow\omega^{-}$ gives the corresponding plots
for the $m<0$ modes.

Lastly, let us note in passing that,

$W_{+}W_{-}=-4(\omega-\omega^{+})(\omega-\omega^{-}),$ (46)

which means that, since $W_{-}>0$ everywhere, the zeros of $W_{+}$ coincide
with the turning points, which can also be seen directly from (35a). Therefore
one can think of $W_{+}$ as an effective potential governing the scattering of
the $p^{\pm}$ modes.

###  IV.7 Mode asymptotics

Having found a well-behaved set of WKB modes and looked at their scattering,
it remains be to shown that these modes indeed incorporate the correct
asymptotic behaviour. Let us start with the limit that $r\to\infty$. Using the
limit for the density in (11), we find that $\widetilde{m}^{2}\to m^{2}$, see
(30). Therefore the dispersion relation in (33) reduces to the well-known one
for plane waves in uniform condensate (albeit in polar coordinates) as we
would expect, see e.g. [50, 14].

In the opposite limit $r\to 0$, we need to do a bit more work. The first thing
to notice is that when the density drops to zero, the BdG equation in (12)
reduces to the usual Schrödinger equation for the decoupled $u_{\pm}$,

$\frac{1}{2r}\partial_{r}\left(r\partial_{r}u_{\pm}\right)-\frac{(m\pm\ell)}{2r^{2}}u_{\pm}+(1\pm\omega)u_{\pm}=0,$
(47)

which has the following general solution,

$\begin{split}u_{\pm}(r\to 0)=&\
\alpha_{\pm}J_{m\pm\ell}\big{(}\sqrt{2(1\pm\omega)}r\big{)}\\\
&\qquad+\beta_{\pm}Y_{m\pm\ell}\big{(}\sqrt{2(1\pm\omega)}r\big{)},\end{split}$
(48)

where $\alpha_{\pm}$ and $\beta_{\pm}$ are constants. If our WKB modes have
the correct limiting behaviour, they should agree with that in (48) when
taking $r\to 0$. For this, we first need to know the asymptotic behaviour of
the different $p$ in (35a), which is found to be,

$p^{2}\overset{r\to
0}{\sim}-\frac{(\ell\mp|m|)^{2}}{r^{2}}+2\left[1\mp\omega\,\mathrm{sgn}(m)\right]+\mathcal{O}(r^{2}),$
(49)

with the upper sign corresponding to $p^{\pm}$ and the lower sign for
$p^{t,b}$. The amplitudes are obtained by plugging the value of $p$ into (39).
Then to compare the WKB modes to the those in (48), we can use
$u_{\pm}\sim\tfrac{1}{2}g\pm if$. From (49), it is readily apparent that there
will be two distinctly different behaviours for $|m|\neq\ell$ and $|m|=\ell$,
which we address now separately.

####  IV.7.1 Asymptotics for $|m|\neq\ell$

In this case, the leading term in (49) as $r\to 0$ is the one proportional to
$r^{-2}$. The phases of the different modes are,

$e^{i\int p^{\uparrow\downarrow}dr}\sim r^{\pm||m|-\ell|},\qquad e^{i\int
p^{t,b}dr}\sim r^{\pm||m|+\ell|},$ (50)

and the amplitudes are constants to leading order. We then recognise that the
Bessel functions in (48) behave as,

$\begin{split}J_{m+\ell}\sim r^{|m+\ell|},&\qquad Y_{m+\ell}\sim
r^{-|m+\ell|},\\\ J_{m-\ell}\sim r^{|m-\ell|},&\qquad Y_{m-\ell}\sim
r^{-|m-\ell|}.\end{split}$ (51)

Hence, we find that the WKB modes do indeed incorporate the correct behaviour
at the origin for $|m|\neq\ell$.

####  IV.7.2 Asymptotics for $|m|=\ell$

In this case, the asymptotic behaviour of the $p^{t,b}$ modes will be the same
as the previous case with the $r^{-2}$ term dominating in (49).
Correspondingly, either $u_{+}$ or $u_{-}$ (depending on the sign of $m$) in
(48) will match this behaviour. For the $p^{\pm}$ modes, the leading term in
(49) is the one which does not depend on $r$. The behaviour of the modes now
depends on both the sign of $m$ and the frequency. For all $m<0$ modes and
also for $m>0$ modes with $\omega<1$, the modes are oscillatory,

$e^{i\int p^{\pm}dr}\sim e^{\pm i|\varepsilon|r},$ (52)

whereas for $m>0$ with $\omega>1$ the modes are evanescent,

$e^{i\int p^{\uparrow\downarrow}dr}\sim e^{\pm|\varepsilon|r},$ (53)

where we have written $\varepsilon=\sqrt{2-2\,\mathrm{sgn}(m)\omega}$ for
brevity. The amplitudes in each of these cases scale with $r^{-\frac{1}{2}}$.
Now we compare to the solutions in (48). The remaining solution, i.e. the one
which does not match the $p^{t,b}$ modes, contain the functions,

$J_{0}\sim 1-\tfrac{1}{4}\varepsilon^{2}r^{2},\qquad Y_{0}\sim\log r,$ (54)

which evidently do not match the behaviour of the WKB modes. In fact, we could
have already anticipated this since the WKB amplitude scaled with
$r^{-\frac{1}{2}}$ which diverges at the origin, thereby violating the
approximation. This is not a serious problem since this is precisely the type
of divergent behaviour that led to the definition of the transfer matrices
around turning points in (43). Thus we can follow the same procedure here to
remedy the situation. To recap, the goal is to take the exact solution at the
offending location (in this case a superposition of zeroth order Bessel
functions), then match the asymptotics of these solutions onto the WKB modes.
For large arguments, the Bessel functions asymptote to,

$\begin{split}J_{0}(\varepsilon r)&\,\sim
r^{-\frac{1}{2}}\cos\left(\varepsilon r-\pi/4\right),\\\ Y_{0}(\varepsilon
r)&\,\sim r^{-\frac{1}{2}}\sin\left(\varepsilon r-\pi/4\right).\end{split}$
(55)

For the case of propagating modes at $r=0$ (which is the case we shall shortly
be concerned with), the $J_{0}$ solution asymptotes to a combination of the
modes in (52) with the amplitudes related by $A^{-}_{0}=iA^{+}_{0}$ whereas
the combination of modes in $Y_{0}$ satisfy instead $A^{-}_{0}=-iA_{0}^{+}$.
Hence, although the WKB modes do not automatically contain the correct
behaviour approaching the origin in this case, they can nonetheless be mapped
onto the correct solutions using the same method we have already applied to
deal with the turning points.

To summarize, using the WKB modes defined in (41), we can map our solution
onto the correct asymptotic mode behaviour contained in the full equations of
motion, which will ultimately allow us to implement the correct BCs. This
analysis also offers an intuitive perspective on the dispersion relation we
are using in (33), namely, it interpolates between the usual Bogoliubov-type
dispersion for the excitations of a uniform condensate in limit $r\to\infty$
and the dispersion of the Schrödinger equation when the density drops to zero
in the limit $r\to 0$.

###  IV.8 Open systems

Having established the relevant preliminaries, we can now compute the
instabilities of an $\ell$-charged vortex. We consider first the case of an
open system, i.e. one with no boundary at large $r$.

####  IV.8.1 Boundary conditions

Since we have two coupled equations in (18), each of which being second order
in spatial derivatives, we need to supply a total of four BCs. Practically,
this will allow us to eliminate all the mode amplitudes from the scattering
computation, resulting in in a condition to be solved for the eigenfrequencies
(as we will see in the next section). The first BC is given by the fact that
the natural modes of vibration in an open system are determined by the
condition that there be no in-coming modes from infinity, i.e.
$A^{-}_{\infty}=0$. Furthermore, the amplitude of the $p^{t}$ mode diverges at
infinity, so we need to take $A^{t}_{\infty}=0$. Similarly, the amplitude of
the $p^{b}$ mode di verges at the origin, hence we should also take
$A^{t}_{0}=0$. Since the $p^{t,b}$ modes do not scatter, these last two
conditions fix the amplitudes of these modes to be zero everywhere, hence, we
need not consider them any further in our calculation.

Before proceeding to state the final BC, consider the following. If the system
only contains either one or no turning points for a mode of a particular
frequency, an in-coming mode will be completely reflected either at the
turning point or at the origin. We see then that setting the amplitude of the
in-coming mode to zero trivially fixes the amplitude of the fourth (and final)
mode in the system to also be zero. Therefore the system has no natural
vibrational modes in the frequency band containing containing one or no
turning points. What is needed then is a cavity capable of trapping modes such
that the system can oscillate under no external influence. Such a cavity
exists when there are either two (for $|m|=\ell$) or three (for $|m|\neq\ell$)
turning points, which can be seen, for example, in Fig. 2.

For $|m|\neq\ell$, let the turning points be labelled $r_{0}<r_{1}<r_{2}$.
Then, by using the $\widetilde{T}$ matrix in (44), we see that setting the
amplitude of the mode which diverges at the origin to zero, i.e.
$A^{\downarrow}_{0}=0$, gives a relation $A^{-}_{0}=iA^{+}_{0}$ between the
propagating modes immediately outside of $r_{0}$.

Similarly, for $|m|=\ell$ we have the two turning points $r_{1}<r_{2}$ and
also $r_{0}=0$. From the analysis in Section IV.7.2, we see that the
combination of modes which yields a regular solution at $r=0$ is also
$A^{-}_{0}=iA^{+}_{0}$. Hence this labelling has allowed us to deal with the
cases possessing both two and three turning points simultaneously.

In summary, the BCs for the excitations of a vortex in an open condensate are
the following,

$A^{-}_{\infty}=A^{t}_{\infty}=A^{b}_{0}=0,\qquad A^{-}_{0}=iA^{+}_{0}.$ (56)

####  IV.8.2 The resonance condition

To find the condition for the natural frequencies of the vortex, we first need
to relate the amplitudes of the $p^{\pm}$ modes at $r_{0}$ to those as
$r\to\infty$, which we achieve through,

$\begin{pmatrix}A^{+}_{0}\\\
A^{-}_{0}\end{pmatrix}=\left|\frac{Q(\infty)}{Q(r_{0})}\right|^{\frac{1}{2}}\begin{pmatrix}e^{-iS_{01}}&0\\\
0&e^{iS_{01}}\end{pmatrix}T\begin{pmatrix}0&e^{S_{12}}\\\
e^{-S_{12}}&0\end{pmatrix}\widetilde{T}\begin{pmatrix}e^{-iS_{2\infty}}&0\\\
0&e^{iS_{2\infty}}\end{pmatrix}\begin{pmatrix}A^{+}_{\infty}\\\
A^{-}_{\infty}\end{pmatrix},$ (57)

where,

$S_{ij}=\int^{r_{j}}_{r_{i}}\left|W_{+}\right|^{\frac{1}{2}}dr,$ (58)

and we have used the fact that the factor $Q$ in (39) is the same for both the
$p^{+}$ and $p^{-}$ modes. Implementing the BCs in (56), the matrix equation
in (57) can be solved for the following condition,

$e^{2iS_{01}}+X=0,\qquad X=\frac{4+e^{-2S_{12}}}{4-e^{-2S_{12}}}.$ (59)

For real frequencies, the first phase term has unit modulus whereas the second
term is larger than 1. Hence, (59) is only satisfied for complex frequencies
$\omega_{\mathbb{C}}=\omega+i\Gamma$. Indeed, these modes are quasi-bound
states localised in the cavity $r\in[r_{0},r_{1}]$ which gradually leak out to
infinity by tunnelling through the barrier which exists in the region
$r\in[r_{1},r_{2}]$. These modes can either grow or decay in time as we shall
now see.

Let us assume that the characteristic growth/decay time is much longer than
the period of oscillation of the modes, i.e. $|\Gamma|\ll|\omega|$, allowing
us to perturb $S_{01}$ as,

$S_{01}(\omega_{\mathbb{C}})=S_{01}(\omega)+i\Gamma\partial_{\omega}S_{01}(\omega),$
(60)

which is a good approximation provided the barrier trapping the modes is large
(or $S_{12}\gg 1$). To leading order in $\Gamma$, (59) splits into two
separate conditions,

$\cos S_{01}(\omega)=0,\qquad\Gamma=-\frac{\log
X(\omega)}{2\partial_{\omega}S_{01}(\omega)}.$ (61)

The condition on the left in (61) is the Bohr-Sommerfeld quantisation
criterion which determines the real part of frequency. This condition
essentially states that a mode must have the correct wavelength to fit inside
the cavity, with an additional phase shift coming from each of the turning
points to account for the break down of WKB there. The expression on the right
of (61) tells us the imaginary component of the frequency once the real part
has been determined. The $\log X$ term is always positive since $X>1$. By
contrast, the sign of the term $\partial_{\omega}S_{01}$ coincides with the
sign of the norm of the modes trapped inside the cavity. In other words, if a
mode which solves (61) is on $\omega_{D}^{+}$ everywhere then it has
$\Gamma<0$ and the mode is stable, whilst if the mode tunnels from
$\omega_{D}^{+}$ outside the cavity to $\omega_{D}^{-}$ inside of it, then it
has $\Gamma>0$ and an instability occurs. From the expressions for
$\omega^{\pm}$ in (45), tunnelling between the positive and negative norm
branches of the dispersion relation only occurs for $0<m<2\ell$, which can
also be seen from Fig. 2. Hence any instabilities that exist must lie in this
range.

We comment in passing that there is an intuitive explanation for the form of
$\Gamma$ in (61). First, note that the time taken for the trapped mode to
complete a full orbit inside $r\in[r_{0},r_{1}]$ is given by $\tau=\pm
2\partial_{\omega}S_{01}$, where the $\pm$ coincides with the sign of the norm
of the mode in the cavity. Rearranging, we find that
$\exp(\Gamma\tau)=|\mathcal{R}|$, where $|\mathcal{R}|=X^{\mp 1}$ is the
modulus of the reflection coefficient associated with the barrier. Hence, the
amount of growth/decay that occurs over a complete orbit of a wave in the
cavity is determined by the size of the wave reflected back into it.

####  IV.8.3 Computing the instabilities

Figure 3: Examples of the vortex response in the time domain following
stimulation by a Gaussian pulse. In the upper panel, only the instability (the
growing mode) is present. The lower panel gives an example of a stable mode
which is present at early times and an unstable mode at late times. If the
magnitude of the growth/decay rate of these different modes is comparable in
magnitude, it is not possible to extract the imaginary component of the
frequency. Furthermore, if the growth of the mode is too slow, it is not
possible to measure over the window of observation.

We compute the unstable complex frequencies with $\omega>0$ by solving the
approximate resonance condition in (61). In Table 1, we display all unstable
modes that occur for $|m|\leq 4$ and $\ell\leq 4$. We find that the $\ell=1$
vortex is stable and that the $m=2$ mode is the most unstable for all the
other $\ell$ probed. This makes sense in light of the scattering picture in
Fig. 2 since the barrier is smallest for low $m$. We also see that, at fixed
$m$ and $\ell$, the modes with the largest frequency will be the fastest to
grow, since for small frequencies the barrier is larger which inhibits the
growth of the mode. As $\ell$ and $m$ increase, the size of the cavity
increases such that it becomes possible to fit more than one instability
inside, as in the case of $m=\ell=4$.

We cross validate our WKB results with a time domain simulation of the BdG
equation in (12). The equations are evolved forward in time via a trapezoidal
method with 5-point centred finite difference stencils for the spatial
derivatives (see Appendix C for more details). Using the behaviour of the
modes at the origin as discussed in Section IV.7, we apply a Dirichlet BC at
$r=0$ for $|m|\neq\ell$ and a Neumann BC for $|m|=\ell$. We apply a Dirichlet
BC at the outer boundary, which is placed sufficiently far away that
reflections from the boundary do not enter our window of analysis. For each
$\ell$ and $m$, we send a Gaussian pulse in towards the centre of the vortex
and after the initial pulse has passed, we are able to measure the response of
the vortex. Examples of the vortex response can be seen in Fig. 3.

The oscillation frequency is extracted by performing a time Fourier transform
of the signal and finding the frequency with maximum amplitude. The growth
rate is obtained by fitting the envelope of the signal with an exponential
function where possible. The results are displayed in Table 1 and display
excellent agreement in the real part of the frequency with those found using
(61). In particular, the agreement improves for larger $\omega$, as one would
expect of the WKB approximation. This confirms that our method can be used to
provide insight into the mechanism underpinning the instability. In all cases
except for the $m=2$ mode, we were unable to extract the value of $\Gamma$,
either because the growth rate was too slow to be resolved within our window
of observation or due to interference between multiple modes in the system;
namely between different instabilities, e.g. for $m=\ell=4$ or between the
unstable mode of interest and a decaying mode (not quoted in Table 1) which we
are not concerned with here.

$\ell$ | $m$ | $\omega_{\mathrm{WKB}}$ | $\omega_{\mathrm{sim}}$ | $\Gamma_{\mathrm{WKB}}$ | $\Gamma_{\mathrm{sim}}$
---|---|---|---|---|---
$2$ | $2$ | $0.3954$ | $0.4399$ | $2.543\times 10^{-3}$ | $2.431\times 10^{-3}$
$3$ | $2$ | $0.2475$ | $0.2871$ | $1.903\times 10^{-3}$ | $2.189\times 10^{-3}$
$3$ | $0.6613$ | $0.6648$ | $7.839\times 10^{-4}$ | n.a.
$4$ | $0.2161$ | $0.2462$ | $2.370\times 10^{-7}$ | n.a.
$4$ | $2$ | $0.1851$ | $0.2138$ | $1.593\times 10^{-3}$ | $1.832\times 10^{-3}$
$3$ | $0.4773$ | $0.4995$ | $8.667\times 10^{-4}$ | n.a.
$4$ | $0.1163$ | $0.1324$ | $3.406\times 10^{-8}$ | n.a.
| $0.7908$ | $0.7804$ | $1.711\times 10^{-4}$ | n.a.
Table 1: Unstable frequencies $\omega+i\Gamma$ of a vortex in an open system
are computed using the WKB result in (61) and by a time domain simulation of
(12). We display all the unstable frequencies for $|\ell|\leq 4$ and $|m|\leq
4$, and quote the values to four digit precision. At given $\ell$ and $m$,
modes with lower $\omega$ are deeper in the cavity and thus have a slower
growth rate, leading to difficulties in extracting $\Gamma_{\mathrm{sim}}$
from simulations. If the growth rate could not be reliably extracted, we
display n.a. in the last column.

###  IV.9 Closed systems

Truly open systems are problematic to probe both experimentally due to the
difficulty in imposing purely out-going BCs at the edge of the system. Indeed,
it is much more natural to impose a trapping potential $U$ which confines the
condensate to a spatially finite domain. The consequence of this is that any
out-going modes will be reflected back into the condensate, causing the
resonance condition to be modified. Instead of the quasi-bound states of the
previous section, the system will then exhibit true bound states which are the
normal modes of the condensate. However, it is still possible to have
instabilities in the spectrum as we will now show. A detailed study of the
real time dynamics of vortex decay as obtained by the GPE evolution can be
found in our companion paper [49].

####  IV.9.1 Boundary conditions

Let us approximate the potential as an infinite barrier located at $r_{B}$,
that it $U(r<r_{B})=0$ and $U(r_{B})=\infty$, with $r_{B}$ a large but finite
radius. Since the condensate wavefunction vanishes for $r>r_{B}$, we must have
$\Psi(r=r_{B})=0$. Taking $r_{B}\gg\ell$, the equation of motion for the
condensate wavefunction in (10) adopts the form,

$\partial_{r}^{2}Y+2Y(1-Y^{2})\simeq 0,$ (62)

near the boundary. The solution to this equation is,

$Y\approx\tanh(r_{B}-r),$ (63)

which drops quickly from the bulk value of $\approx 1$ to zero within a few
healing lengths (due to the quantum pressure term in (6) [51]). The question
of determining the correct BCs at $r_{B}$ then becomes equivalent to asking
how plane waves in a 1D condensate are reflected by a hard wall. In Fig. 4, we
verify that a pulse incident on a hard wall (comprised of wavelengths larger
than the healing length) reflects analogously to waves propagating through a
homogeneous condensate with a Neumann BC at $r_{B}$, that is
$A^{-}_{B}=A^{+}_{B}$. Since the instability is a low frequency phenomenon, we
will therefore assume in our WKB analysis that the density has the same form
of that of the free vortex all the way up to $r_{B}$ where we impose
$A^{-}_{B}=A_{B}^{+}$ on the propagating modes.

Figure 4: Initially left moving pulse, whose evolution is dictated by 1D BdG,
colliding with a hard wall at $x=0$. Three different times are shown, $t=0$
(top) $t=1250$ (middle) and $t=2500$ (bottom). The blue curve uses the true BC
at $x=0$ for $u_{+}$ (Dirichlet) with the background density profile
$\rho=\tanh^{2}x$. We find good agreement with the red dashed curved, which
uses a Neumann BC with $\rho=1$.

The remaining BCs are imposed by noting that the behaviour of the modes at the
origin is the same as in an open system, therefore at $r=0$ we have the same
BCs as in the case of the free vortex. Finally, since the $p^{t,b}$ mode will
also be completely reflected at $r_{B}$ and these modes do not scatter
anywhere else in the system, their amplitude must be zero everywhere since the
$p^{b}$ mode diverges at the origin.

In summary, the BCs for the excitations of a vortex in a bounded condensate
with a hard wall located at $r_{B}$ are,

$A^{t}_{B}=A^{b}_{0}=0,\qquad A^{-}_{0}=iA^{+}_{0},\qquad
A^{-}_{B}=A^{+}_{B}.$ (64)

####  IV.9.2 The resonance condition

The normal mode spectrum of the vortex can now be obtained from a resonance
condition similar to the one derived in Section IV.8.2. We focus our attention
on scenarios where there are, again, either two or three turning points using
the labelling convention for $r_{0}<r_{1}<r_{2}$ already established. Note
that in a closed system, there will also be normal modes in the frequency
bands containing either one or no turning points due to the box-like BCs,
however, these bands do not contain instabilities which are our main interest.
The matrix equation relating the amplitudes at $r_{0}$ and $r_{B}$ is the same
as in (57) (provided we interpret the amplitudes on the right of the
expression as the ones at $r_{B}$).

Applying the BCs in (64), the resulting condition can be cast into the
following form,

$4\cot(S_{01})\cot(S_{2B}+\pi/4)=e^{-2S_{12}}.$ (65)

We solve this condition for the case of $m=\ell=2$ for a variety of $r_{B}$
and display the spectrum of modes in Fig. 5.

Normal modes which are localised inside the cavity $r\in[r_{0},r_{1}]$ are
shown in red, whereas the black points are modes which are live outside in the
region $r\in[r_{2},r_{B}]$. The green points represent a resonant coupling
between these two different types of modes which leads to an instability (as
we will soon discuss). Also on Fig. 5, we display (in grey) the
eigenfrequencies as determined by a numerical resolution of the BdG equation
in (12). These are found as the eigenvalues of the BdG operator $\widehat{L}$,
where the radial derivatives in $\widehat{L}$ are expressed using 5-point
centred finite different stencils which encode the boundary conditions. The
WKB method agrees exceptionally well with the exact solution for modes
localised outside the cavity since the flow is slowly varying in this region.
The agreement for modes localised inside the cavity, although still
satisfactory, is reduced since the flow varies more rapidly in this region
(see e.g. the density profiles in Fig. 1). Crucially, however, the WKB
spectrum captures all the essential features of the exact spectrum (in
particular the crossing of positive and negative norm modes) meaning we can
use our WKB formula in (65) to interpret the instability mechanism.

The features of the plots in Fig. 5 can be understood purely from the form of
the condition in (65). If the zeros of the two cotangent functions are well
separated, the resonance condition is satisfied when one of these is
exponentially small (which we approximate as zero). Therefore, we either have
$\cot(S_{01})=0$, which gives the red data in Fig. 5, or
$\cot(S_{2B}+\pi/4)=0$, which gives the black dots. However, if the zeros of
the cotangents are close together, it is no longer correct to approximate
$\cot\simeq 0$ since only the product of the two needs to be exponentially
small. This leads to a deviation in behaviour away from the curves traced out
by the red and black data, which is precisely where the green dots appear.

When this coupling between the two types of modes occurs, the mode outside the
vortex transmits negative norm into the cavity when it scatters with the
barrier, leading to amplification of the mode outside the cavity. Similarly,
the mode inside the cavity transmits more positive norm to the outside of the
cavity as it reflects off the barrier, leading to amplification inside. This
continual amplification of the mode causes it to become unstable. Note that
during this process, norm conservation is never broken since positive and
negative norm are created in equal amounts. Indeed, an identity resulting from
the BdG equation (12) is,

$(\omega-\omega^{*})\int
d^{2}\mathbf{x}\left(|u_{+}|^{2}-|u_{-}|^{2}\right)=0,$ (66)

where the integral is the total norm from (16). This shows that for unstable
modes with $\mathrm{Im}[\omega]\neq 0$, the eigenmodes have zero total norm.

It has also been noted in the literature [31, 44] that in the limit that the
system becomes infinitely large, i.e. $r_{B}\to\infty$, the frequency of the
instability should asymptote to that in the open system. This is somewhat
counter intuitive since the eigenfrequencies of a system are highly dependent
on the BCs, and imposing open/closed BCs are completely different procedures.
Our condition in (65) offers some insight into this conundrum. Rearranging, we
find,

$\begin{split}e^{2iS_{01}}+X=&\ ie^{2iS_{2B}}\left(Xe^{2iS_{2B}}+1\right),\\\
\sim&\ \mathcal{O}\left(e^{-2\Gamma\tau_{B}}\right),\end{split}$ (67)

where $\tau_{\infty}=\partial_{\omega}S_{2B}$ is the time taken for the mode
to propagate from just outside the barrier to the outer boundary, and we have
assumed a slow growth rate with $|\Gamma|\ll|\omega|$. Hence, as
$r_{B}\to\infty$ we also have $\tau_{B}\to\infty$ and the resonance condition
for the closed system reduces to that of the open one in (59). Physically,
this is due to the fact that the instability originates inside the cavity, so
its amplitude at small $r$ will be exponentially larger than the mode
reflected by the boundary at large $r$. In the limit where the boundary is
infinitely far away, the amplitude of the reflected mode will be infinitely
smaller than the out-going mode and can therefore be discarded.

Finally, this resonant coupling picture also provides a useful perspective on
the vortex instability in the open system. In the vortex core, the normal
modes are those which roughly fit in the cavity whereas outside the vortex,
the spectrum of normal modes is continuous since the system in that case is
infinite. Therefore, the existence of the instability in the open system
depends only on the existence of a trapped mode in the cavity, since there is
guaranteed to be a mode outside which has the right frequency for a coupling
to occur.

Figure 5: The spectrum of the vortex for $m=\ell=2$ as a function of the
system size $r_{B}$. The upper panel gives the oscillation frequency. Modes
localised inside (outside) the cavity are shown in black (red). Note that the
total energy (norm multiplied by frequency) of the modes multiplied by the
frequency is positive for black and negative for red. When these two types of
modes cross each other, the resulting modes have zero total norm (green) [44]
and are unstable. The imaginary part of frequency of these unstable modes is
shown in the lower panel. This plot is a stability map, illustrating that the
instability occurs in bubbles of finite width. Eigenfrequencies computed from
a numerical solution of (12) are shown in grey for comparison.

##  V Black holes

In this section, we will study superradiance around BH spacetimes and
illustrate the connection with the vortex instabilities of the previous
section. Before looking at a gravitational (Kerr) BH, we first consider an
analogue BH in order to bridge the gap with the vortices we have already
looked at.

###  V.1 Draining vortex

Consider a BEC containing a vortex of the type (9), but now with a draining
component in the velocity field. It is well-known that the resulting system
shares many features with a rotating BH spacetime. In the region where the
density is approximately uniform, i.e. $\rho=1$, the continuity equation (7)
gives the following solution for the velocity field,

$\mathbf{v}=-\frac{d}{r}\,\vec{\mathbf{e}}_{r}+\frac{\ell}{r}\,\vec{\mathbf{e}}_{\theta},$
(68)

where we have adopted the units established in Section IV.1 and we take $d>0$
since the system is draining. To establish such a set up, one can imagine
pumping atoms out of the condensate at a rate $\dot{N}=-d$ in a small region
around the vortex axis $r=0$. Note that (68) is only an estimate and the full
fluid flow is more complicated in both classical [52, 53] and quantum [35]
fluids. However, this approximation is known to hold far from the vortex axis,
e.g. [21], and as such, can be treated as the leading order solution in a
large $r$ expansion.

Excitations of the flow satisfy the equations of motion in (18). Assuming that
$d\gg 1$, so that the characteristic draining length scale is much larger than
the healing length, the relevant frequencies to the problem will be smaller
than $\mathcal{O}(1)$ and the final term in (18b) (which arises from the
quantum pressure) can be neglected. The two equations in (18) now combine into
a single wave equation,

$D_{t}^{2}\phi-\nabla^{2}\phi=0.$ (69)

This wave equation describes long wavelength phonons propagating at speed
$c=\sqrt{\mu/M}$ (which is 1 in these units) independent of their frequency.
This approximation therefore neglects the effects of dispersion.

The wave equation in (69) is formally equivalent to the equation of motion for
a massless Klein-Gordon (KG) field propagating through a curved spacetime,

$\frac{1}{\sqrt{-g}}\partial_{\mu}\left(\sqrt{-g}g^{\mu\nu}\partial_{\nu}\phi\right)=0,$
(70)

where the effective metric entering this equation is given by,

$\begin{split}g^{\mu\nu}dx^{\mu}dx^{\nu}=&\
-\left(1-\frac{\ell^{2}+d^{2}}{r^{2}}\right)dt^{2}+\frac{2d}{r}dr\,dt\\\
&\qquad-2\ell\,d\theta\,dt+dr^{2}+r^{2}d\theta^{2},\end{split}$ (71)

with coordinates $(t,r,\theta)$. Here we take the metric signature as
$(-1,1,1,1)$ and $g$ is the determinant of the metric. The metric (71) can be
written in a form where the BH behaviour is more apparent by performing a
coordinate transformation $(t,\theta)\to(\tilde{t},\tilde{\theta})$ with,

$d\tilde{t}=dt+\frac{v_{r}dr}{1-v_{r}^{2}},\qquad
rd\tilde{\theta}=rd\theta+\frac{v_{\theta}v_{r}dr}{1-v_{r}^{2}},$ (72)

where we have written $v_{i}=\vec{\mathbf{e}}_{i}\cdot\mathbf{v}$. The metric
becomes,

$\begin{split}\tilde{g}_{\mu\nu}d\tilde{x}^{\mu}d\tilde{x}^{\nu}=&\
-\left(1-\frac{r_{e}^{2}}{r^{2}}\right)d\tilde{t}^{\,2}+\frac{dr}{1-r_{h}^{2}/r^{2}}\\\
&\qquad-2\ell\,d\tilde{\theta}\,d\tilde{t}+r^{2}d\tilde{\theta}^{2},\end{split}$
(73)

with $\tilde{x}=(\tilde{t},\tilde{\theta},r)$ and we have defined,

$r_{e}=\sqrt{\ell^{2}+d^{2}},\qquad r_{h}=d,$ (74)

which satisfy $\tilde{g}_{\tilde{t}\tilde{t}}(r_{e})=0$ and
$\tilde{g}_{rr}(r_{h})\to\infty$. These locations are ergosphere $r_{e}$, the
boundary of the region in which a phonon will be forced to co-rotate with the
vortex as compared to a static point at infinity (i.e. the ergoregion), and
the horizon $r_{h}$, the radius inside which all phonons are forced to move
inward (i.e. toward $r=0$). This is markedly different from the non-draining
vortex discussed earlier, which had an ergoregion but no horizon. We will
shortly use the WKB method to make the connection between $r_{e}$ (defined
using the metric above) and the ergoregion of the non-draining vortex (which
we defined as the region where excitations could have negative energy).

####  V.1.1 Conventional method

Here we recap the conventional analysis used to study BH superradiance, which
has been applied to this particular draining vortex set-up on multiple
occasions. More details can be found in e.g. [5, 6, 7, 8, 9, 9, 10, 11, 12,
13].

To begin, we split the field into the different $\omega$ and $m$ modes,

$\begin{split}\phi=&\
\frac{R(r)}{\sqrt{r}}e^{im\tilde{\theta}-i\omega\tilde{t}},\\\ =&\
\frac{R(r)}{\sqrt{r}}\exp\left(-i\int\frac{v_{r}\tilde{\omega}\,dr}{1-v_{r}^{2}}\right)e^{im\theta-i\omega
t},\end{split}$ (75)

with a generic wave consisting of a sum of these modes, as in (19), and we
have defined,

$\tilde{\omega}=\omega-\frac{m\ell}{r^{2}}.$ (76)

Next, one introduces a new coordinate (which in GR is called the tortoise
coordinate),

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

thereby mapping $r\in(r_{h},\infty)$ onto $r_{*}\in(-\infty,\infty)$. Plugging
(75) and (77) into (69), the problem of a wave scattering with a BH is
converted into 1D scattering problem,

$-\partial_{r_{*}}^{2}R+VR=0,$ (78)

in the presence of an effective potential barrier,

$\begin{split}V=&\
-\left(\omega-\frac{m\ell}{r^{2}}\right)^{2}+\left(1-\frac{d^{2}}{r^{2}}\right)\\\
&\qquad\qquad\qquad\qquad\times\left(\frac{m^{2}-\frac{1}{4}}{r^{2}}+\frac{5d^{2}}{4r^{4}}\right).\end{split}$
(79)

The scattering equation in (78) has the asymptotic solutions,

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

where $\tilde{\omega}_{h}=\omega-m\Omega_{h}$ is the quantity in (76)
evaluated on the horizon. Furthermore, an equation of the form (78) admits a
conserved quantity called the Wronskian,

$\partial_{r_{*}}\left(R^{*}\partial_{r_{*}}R-R\,\partial_{r_{*}}R^{*}\right)=0,$
(81)

which is equivalent to the conserved norm in (23). Finally, plugging the
solutions (80) in (81) one finds that the relation,

$|\mathcal{R}|^{2}+\frac{\tilde{\omega}_{h}}{\omega}|\mathcal{T}|^{2}=1,$ (82)

is satisfied, where the reflection $\mathcal{R}$ and transmission
$\mathcal{T}$ coeffcients are defined,

$\mathcal{R}=A^{+}_{\infty}/A^{-}_{\infty},\qquad\mathcal{T}=A_{h}/A^{-}_{\infty}.$
(83)

Thus, one deduces that amplification of the incident wave, i.e.
$|\mathcal{R}|>1$, occurs for $\tilde{\omega}_{h}/\omega<0$ or, for positive
frequency modes, when,

$\omega<m\Omega_{h},$ (84)

is satisfied. This is the phenomenon of rotational superradiance [22].

####  V.1.2 WKB method

Now, let us see how superradiance arises in the WKB approximation, since this
will allow us to make contact with our study of vortex stability in Section
IV. There are two possible approaches which turn out to be equivalent in the
limit of large $m$: either we apply the WKB approximation directly to (78) or
start from the full wave equation in (69). The second approach leads to the
conventional dispersion relation and as such, is the most natural to compare
to the vortex analysis in Section IV.5. However, we shall discuss both for
completeness.

Taking the second route, we write a WKB ansatz,

$\phi=\mathcal{A}(r)e^{i\int p(r)dr+im\theta-i\omega t},$ (85)

which we then plug into the wave equation in (69). Taking
$\mathcal{A}^{\prime}\ll p\mathcal{A}$, the leading order result is the
dispersion relation $\Omega^{2}=k^{2}$ with,

$\Omega=\omega-\frac{m\ell}{r^{2}}+\frac{pd}{r},\qquad
k=\sqrt{p^{2}+\frac{m^{2}}{r^{2}}}.$ (86)

The dispersion relation then separates into two branches,

$\omega_{D}^{\pm}=\frac{m\ell}{r^{2}}-\frac{pd}{r^{2}}\pm|k|,$ (87)

where modes with $\omega=\omega_{D}^{+}$ ($\omega=\omega_{D}^{-}$) have
positive (negative) norm, as can be verified by plugging (85) into (23) and
discarding the dispersive terms. At next to leading order, we obtain an
expression for the amplitude $A=\sqrt{r}\mathcal{A}$,

$A\propto|v_{g}^{r}\Omega|^{-\frac{1}{2}},$ (88)

where $v_{g}^{r}=\partial_{p}\omega_{D}^{\pm}$ is the radial component of the
group velocity. The dispersion relation in this case is quadratic in $p$ so
there are two different solutions,

$p^{\pm}=\frac{\tilde{\omega}d/r\pm\sqrt{-V_{|m|\gg 1}}}{1-d^{2}/r^{2}},$ (89)

where the term under the square root is the large $m$ limit of the potential
in (79), which has the effect of setting $m^{2}-1/4+5d^{2}/4r^{2}\to m^{2}$.
For the remainder of this section, we drop the $|m|\gg 1$ subscript and
implicitly assume that this limit has been applied to $V$. Following the
procedure outlined in Section IV.6, scattering occurs when $p^{\pm}$ change
from propagating modes into evanescent modes. This happens at the turning
points $r_{tp}$ which satisfy $V(r_{tp})=0$.

Let us briefly check what would happen had we applied the WKB approximation
directly to (78). Writing $R\sim A(r_{*})e^{i\int P(r_{*})dr_{*}}$ and taking
$A^{\prime}\ll PA$ and $P^{\prime}\ll P^{2}$, one finds,

$R\sim|V|^{-\frac{1}{4}}e^{i\int\sqrt{-V}dr_{*}}.$ (90)

Plugging this back into (75) and applying $|m|\gg 1$ yields precisely (85),
thereby confirming the equivalence of the different approaches for large
$|m|$. As a brief aside, using (90) we have $P^{2}+V=0$. By comparing this
with the $p^{\pm}$ modes in (35a) which obeyed $p^{2}+W_{+}=0$, we see that we
can think of $W_{+}$ as the effective potential for vortex excitations. We
will revisit this point when we look at Kerr BHs later on.

Returning to the task at hand, we now want to study scattering of the WKB
modes to see how superradiance arises in this picture. For this purpose, let
us work with (90), although identical results are obtained when using (85) as
shown in [24, 25]. Inspection of the form of $V$ shows that there are 0, 1 or
2 turning points for $r>r_{h}$. In the case of no turning points, in-going
modes are not scattered (at our level of approximation) and we have
$|\mathcal{R}|=0$. The case of 1 turning point is the marginal case where two
turning points merge into one. However, to apply the transfer matrices in (43)
we need to have propagating modes on one side of $r_{tp}$ and evanescent modes
on the other. In this case, one needs to build the transfer matrix using a
different asymptotic expansion than the one we have used here, see e.g. [54,
48]. Proceeding to the case of 2 turning points $r_{tp}\in[r_{1},r_{2}]$, we
have propagating modes for $r<r_{1}$ and $r>r_{2}$ and evanescent modes for
$r_{1}<r<r_{2}$ under the potential barrier $V$. To see where the zeros
$r_{1,2}$ of the $V$ occur for different frequencies, we may write,

$V=-(\omega-\omega^{+})(\omega-\omega^{-}),$ (91)

where these curves are $\omega^{\pm}=\omega_{D}^{\pm}(p_{tp})$. Note in
passing the similarity of (91) with the relation in (46), where in that case
the fact that two potentials are required to make the combination on the right
hand side is a consequence of dispersion. In the current case, the
$\omega^{\pm}$ are given by the following expressions,

$\omega^{\pm}=\frac{m\ell}{r^{2}}\pm\sqrt{\left(1-\frac{d^{2}}{r^{2}}\right)\frac{m^{2}}{r^{2}}}.$
(92)

These curves are precisely analogous to those defined in the vortex case in
(45). In particular, the intersections $\omega=\mathrm{const}$ and the
$\omega^{\pm}$ give the locations of the $r_{tp}$. When $\omega$ crosses
$\omega^{+}$ ($\omega^{-}$), there will be a turning point on the positive
(negative) norm branch of the dispersion relation. By looking at the example
in Fig. 6, we see that tunnelling to the negative norm branch occurs for
positive frequencies satisfying the superradiance condition in (84). Therefore
superradiance occurs when there is tunnelling between the positive and
negative norm branches of the dispersion relation. This was exactly the
mechanism that gave rise to the vortex instabilities in Section IV.8. Hence,
we find that superradiance and vortex instabilities are controlled by the same
physics.

In conventional BH physics, relativists talk about the existence of the
ergosphere as the mechanism underpinning superradiance. It is a simple matter
to see how this ties in with the WKB picture we are using, where the principal
tools are the functions $\omega^{\pm}$. From Fig. 6, it is clear that the
condition $\omega^{-}=0$ determines the radius inside which positive frequency
modes can be on the negative norm branch of the dispersion relation. Solving
this condition gives precisely the ergosphere $r_{e}$ defined in (74).
Therefore the existence of an ergosphere is equivalent to the existence of a
surface inside which there are $\omega>0$ modes are $\omega_{D}^{-}$, i.e.
modes with negative energy density. Note that the $m$-independence of the
ergosphere happens only in the non-dispersive regime. In the case of vortex
excitations, it also makes sense to talk about an ergosphere but as an
$m$-dependent notion. Indeed, the lack of ergosphere for vortex modes above a
critical $m$ is precisely the reason why those modes are stable.

Figure 6: The functions $\omega^{\pm}$ in (92) for the parameters $m=1$,
$\ell=d=5$. Superradiant amplification occurs for $0<\omega<m\Omega_{h}$ which
is the range in which positive frequency modes tunnel to the negative norm
branch of the dispersion relation. This identifies the mechanism underlying
superradiance as the same one causing the vortex instabilities.

As the final step, we can compute the reflection coefficient in the WKB
approximation to estimate the amount of amplification caused by BH
superradiance. The method follows the scattering computation in (57). However,
in this case, the calculation is simpler; we only need to evaluate the
scattering coefficients immediately outside the turning points since the WKB
approximation guarantees that the norm of each mode is fixed as it moves off
into the asymptotic regions. We write,

$\begin{pmatrix}A^{+}_{1}\\\
A^{-}_{1}\end{pmatrix}=\left|\frac{V(r_{2})}{V(r_{1})}\right|^{\frac{1}{4}}T\begin{pmatrix}0&e^{S_{12}}\\\
e^{-S_{12}}&0\end{pmatrix}\widetilde{T}\begin{pmatrix}A^{+}_{2}\\\
A^{-}_{2}\end{pmatrix},$ (93)

with,

$S_{ij}=\int^{{r_{*}}_{j}}_{{r_{*}}_{i}}|V|^{\frac{1}{2}}dr_{*}.$ (94)

Assuming $\omega>0$, the absorbing BC at the horizon fixes $A_{2}^{+}=0$ for
$\tilde{\omega}_{h}>0$ and $A_{2}^{-}=0$ for $\tilde{\omega}_{h}<0$. To solve
for the reflection coefficient, we then note that
$|V(r_{2})|^{\frac{1}{2}}A_{2}^{j}=|V(\infty)|^{\frac{1}{2}}A_{\infty}^{j}$
and we apply the definition in (83) to obtain,

$\mathcal{R}=e^{-\frac{i\pi}{2}}\left(\frac{4-e^{-2S_{12}}}{4+e^{-2S_{12}}}\right)^{\frac{\tilde{\omega}_{h}}{\omega}}.$
(95)

As expected, this is less than 1 for $\tilde{\omega}_{h}>0$ and greater than 1
for $\tilde{\omega}_{h}<0$.

In summary, our analysis demonstrates that superradiance occurs when in-coming
positive norm modes tunnel to the negative norm branch of the dispersion
relation. The comparison between Figs. 2 and 6 makes the connection to the
stability of quantised vortices apparent. With this connection, the vortex
instability can be interpreted as a superradiating bound state trapped in the
vortex core. Alternatively, with our understanding of vortex instabilities in
closed systems, superradiance can be interpreted as a coupling between
positive norm modes in the asymptotic region and negative norm modes near the
horizon, which act to lower the energy of the BH. Unlike the vortex where this
coupling only occurred for special system sizes, the spectrum of BH modes is
continuous due to the open boundary conditions and hence, the coupling always
occurs in the frequency band given by (84). Finally, although we have
neglected dispersion, note that the same picture emerges when taking
dispersion into account, as has recently been demonstrated in [25], although
the superradiance criterion in (84) picks up corrections in that case.

####  V.1.3 Stability

In the (analogue) BH case, incident modes can be amplified but no instability
occurs. This is due to the presence of absorbing BCs at the horizon and
infinity which means that perpetual amplification is prevented. In the vortex
case, there was a reflective BC at the origin which allowed for the amplified
mode to become trapped. This suggests that adding a generic drain component to
a vortex will stabilise it, since the drain can act like a horizon and
dissipate away the negative norm mode in the core. One might then reasonably
ask what happens if the drain is not perfectly absorbing and instead a small
amount of reflection occurs. In this section, we search for the conditions
under which a quantised vortex can be stabilised by a drain. This is also of
interest in view of recent reports that multiply-quantised vortices are stable
in polariton-fluid experiments when a drain is added to the flow [36].

Let us assume that we have an open system extending to infinity and that the
location $r_{h}$ is no longer a true horizon in the sense that modes are
reflected there with reflection coefficient,

$\mathcal{Z}=A_{h}^{-}/A_{h}^{+}.$ (96)

Here we have restricted ourselves to looking at modes with
$\tilde{\omega}_{h}<0$ (i.e. superradiant ones) so that $A_{h}^{+}$ is the in-
going mode and $A_{h}^{-}$ is the out-going one. We can now perform our
scattering computation, taking the matrix equation in (93) as our starting
point. Then applying $A^{-}_{2}=0$ (since we send no waves into the system)
and transporting the modes at $r_{1}$ to the horizon through
$A_{h}^{\pm}=|V(r_{1})/V(r_{h})|^{\frac{1}{4}}e^{\mp iS_{h1}}A^{\pm}_{1}$, we
obtain the following resonance condition,

$e^{2iS_{h1}}+\mathcal{R}\mathcal{Z}=0,$ (97)

where the $\mathcal{R}$ is the superradiant expression in (95). If the
growth/decay rate is small compared to the oscillation frequency, i.e.
$|\Gamma|\ll|\omega|$, this splits into,

$S_{h1}=\pi\left(\nu+\tfrac{1}{4}\right)+\chi,\qquad\Gamma=\frac{\log\left|\mathcal{R}\mathcal{Z}\right|}{2\left|\partial_{\omega}S_{h1}\right|},$
(98)

where $\nu=0,1,...$ indexes the different possible bound states and we have
written $\mathcal{Z}=|\mathcal{Z}|e^{2i\chi}$. Note that the formulae in (61)
are recovered for the case with $\mathcal{Z}=i$.

From the expressions in (98) we learn that for the system to be stable (i.e.
$\Gamma<0$), the product of the reflection coefficients needs to be less than
unity, i.e.

$|\mathcal{R}\mathcal{Z}|<1.$ (99)

Recalling we are in the superradiant case where $|\mathcal{R}|>1$, this gives
us the intuitive result that the system is stable if the boundary at $r_{h}$
absorbs the wave more than the barrier at $r_{1}$ amplifies it. Repeating the
same analysis for the non-superradiant case, one finds the same stability
condition as in (99), which is always satisfied if the boundary at $r_{h}$ is
an absorber and not an amplifier.

There are a couple of reasons why one can expect the condition in (99) to be
satisfied in a realistic experiment. Firstly, the reflection coefficient for
amplified modes in BH-type rotational superradiance is not usually vastly
greater than unity, see e.g. [11, 13, 24, 25, 22]. By contrast, an efficient
drain would be expected to highly damp out in-coming modes, e.g. in [54] the
maximum value of $\mathcal{Z}$ was $\sim 0.2$ for counter-rotating modes.
Lastly, dissipation inherent in any real experiment (due to e.g. viscosity or
conversion into other types of waves [21]) will further suppress the amount of
superradiant amplification, such that it becomes more difficult for
$|\mathcal{R}\mathcal{Z}|$ to exceed unity.

As a final note, we point out that here we have only considered the
possibility of unstable modes in the case of an infinite system, i.e.
ergoregion instabilities. If the system were finite, it would also be possible
to have instabilities outside the ergoregion, i.e. black hole bomb modes. In
the next section, we discuss this type of instability in the context of an
astrophysical black hole.

###  V.2 Kerr spacetime

Finally, we come to study superradiance around an astrophysical Kerr BH. In
this section, we will focus on a superradiant instability that occurs for
massive fields in the Kerr spacetime, comparing and contrasting this to the
instability around quantised vortices. We work in natural units such that
$G=\hbar=c=1$, where $G$ is the gravitational constant, $\hbar$ is the reduced
Planck’s constant and $c$ is the speed of light.

####  V.2.1 The metric

The Kerr metric describes a rotating BH spacetime with mass parameter $M$ and
angular momentum $J=aM$. The popular representation of the metric involves
Boyer-Lindquist coordinates $x^{\nu}=(t,r,\varphi,\theta)$ where $\nu$ is a
spacetime index and, in contrast to previous sections where only two spatial
dimensions were considered, $\theta\in[0,\pi]$ is the polar angle and
$\varphi\in[0,2\pi)$ is the azimuthal angle. Note in particular that in this
section $\varphi$ is a coordinate and not a field. The Kerr metric
$g_{\nu\sigma}$ in these coordinates is [55],

$g_{\nu\sigma}dx^{\nu}dx^{\sigma}=-\frac{\Delta_{r}}{\rho^{2}}\hat{\omega}_{1}^{2}+\frac{a^{2}\sin^{2}\theta}{\rho^{2}}\hat{\omega}_{2}^{2}+\frac{\rho^{2}}{\Delta_{r}}dr^{2}+\rho^{2}d\theta^{2},$
(100)

with the metric components,

$\hat{\omega}_{1}=dt-a\sin^{2}\theta
d\varphi,\qquad\hat{\omega}_{2}=dt-\frac{r^{2}+a^{2}}{a}d\varphi,$ (101)

The metric functions are,

$\Delta_{r}=r^{2}+a^{2}-2Mr,\qquad\rho^{2}=r^{2}+a^{2}\cos^{2}\theta,$ (102)

where $\rho$ should not be confused with the density in earlier sections. The
inner and outer horizons, $r_{-}$ and $r_{+}$ respectively are the locations
where $\Delta_{r}=0$,

$r_{\pm}=M\pm\sqrt{M^{2}-a^{2}}.$ (103)

Since we are only concerned with the outer horizon from here on, let
$r_{h}=r_{+}$ for labelling purposes. To avoid having a naked singularity, the
rotation parameter must satisfy $a<M$. The ergosphere is the region inside
which an observer is forced to co-rotate with the BH with respect to an
observer located at infinity [55]. It is located inside $r_{h}<r<r_{e}$ where,

$r_{e}=M+\sqrt{M^{2}-a^{2}\cos^{2}\theta}.$ (104)

From here on, we will work in units where $M=1$ so that $a<1$.

####  V.2.2 The Klein-Gordon field

The equation of motion for a massive scalar field $\phi$ with mass $\mu$ is
the Klein-Gordon (KG) equation,

$\frac{1}{\sqrt{-g}}\partial_{\nu}\left(\sqrt{-g}g^{\nu\sigma}\partial_{\sigma}\phi\right)-\mu^{2}\phi=0,$
(105)

where $g^{\nu\sigma}$ and $g$ are the inverse and determinant of the metric
respectively, and repeated indices are summed over in accordance with the
Einstein summation convention. Since $t$ and $\varphi$ are symmetries of the
metric (or in the language of GR, $\partial_{t}$ and $\partial_{\varphi}$ are
Killing vectors of the spacetime) the field $\phi$ is split into the different
$\omega$ and azimuthal $m$ components,

$\phi=\int d\omega\sum_{m}R_{\omega m}(r)S_{\omega
m}(\theta)e^{im\varphi-i\omega t},$ (106)

and KG equation splits into two separate equations. The angular part of the
field is governed by,

$\begin{split}&\sin\theta\partial_{\theta}\left(\sin\theta\partial_{\theta}S_{\omega
m}\right)-\big{[}\left(\omega a\sin^{2}\theta-m\right)^{2}\\\
&\qquad\quad+\sin^{2}\theta\left(\mu^{2}a^{2}\cos^{2}\theta-j_{\omega
m}\right)\big{]}S_{\omega m}=0,\end{split}$ (107)

and the radial part obeys,

$\Delta_{r}\partial_{r}\left(\Delta_{r}\partial_{r}R_{\omega
m}\right)+K^{2}R_{\omega m}=0,$ (108)

with,

$K^{2}=\left[\omega(r^{2}+a^{2})-am\right]^{2}-\Delta_{r}\left(j_{\omega
m}+\mu^{2}r^{2}\right).$ (109)

The separation constant $j_{\omega m}$ can be obtained as the eigenvalue of
(107) when imposing regular BCs at the poles $\theta=\pm\pi/2$. Indeed, as
well as solving numerically, one can apply the by now familiar WKB methods to
(107) to find approximate expressions for the $j_{\omega m}$ that give some
insight into it’s behaviour. Such an analysis can be found in Appendix E. In
the non-rotating $a=0$ case, we have $j_{\omega m}=l(l+1)$, where $l\geq|m|$
is a positive integer called the polar “quantum number”. In the rotating case,
$l$ can also be used to index the different $j_{\omega m}$. From here on, we
assume that (107) has been solved and the value of $j_{\omega m}$ determined.

Starting from (108), one can follow a similar analysis to that outlined in
Section V.1.1 to deduce the occurrence of superradiant amplification of
incident $\omega>0$ modes in the frequency band,

$\mu<\omega<m\Omega_{h},\qquad\Omega_{h}=\frac{a}{r_{h}^{2}+a^{2}},$ (110)

where $\Omega_{h}$ is the angular velocity of the spacetime on the horizon
[55]. Since we have already drawn the comparison between superradiance and the
vortex instability mechanism using the analogue BH, we will not reproduce the
full analysis here (the interested reader can find a full treatment in e.g.
[56, 57]). Instead, here we choose to focus on another aspect of the Kerr
spacetime particular to massive fields: the BH bomb instability.

###  V.3 Bound states of a Kerr BH

The possibility of having (quasi-)bound states around the BH arises due to the
non-zero rest mass $\mu\neq 0$ of the field, which acts as a mirror at large
$r$ reflecting the out-going part of the field back into centre. If this
happens to a superradiant mode, then the amplified mode is reflected back into
the BH where it is further amplified, leading to an instability. In the
literature, an instability of this type is known as a BH bomb [22]. We shall
shortly see that although, like the vortex instability, these instabilities
owe their existence to the coupling of positive norm modes to negative norm
ones, they are unlike the vortex one in the sense that they occupy the
positive norm region of the system. These unstable BH frequencies have been
computed using precise numerical algorithms in e.g. [43]. For the purposes of
this work, a WKB approximation of the modes will be sufficient since this
permits an easy comparison with the vortex case.

####  V.3.1 WKB method

To apply the WKB method to (108), we write,

$R_{\omega m}(r)=\mathcal{A}(r)\,e^{\pm i\int K(r)\,\Delta_{r}^{-1}dr},$ (111)

for $\Delta_{r}\partial_{r}K\ll K^{2}$ and
$\Delta_{r}\partial_{r}\mathcal{A}\ll K\mathcal{A}$. The leading order
equation is then exactly (109) and at next to leading order we find
$\mathcal{A}\sim|K|^{-\frac{1}{2}}$. It is useful at this point to introduce
the effective radial potential,

$\begin{split}V=&\ -K^{2}/(r^{2}+a^{2})^{2},\\\ =&\
-\left(\omega-\frac{ma}{r^{2}+a^{2}}\right)^{2}+\frac{\Delta_{r}}{(r^{2}+a^{2})^{2}}\\\
&\qquad\qquad\qquad\qquad\qquad\times\left(j_{\omega
m}+\mu^{2}r^{2}\right),\end{split}$ (112)

which is analogous to the large $|m|$ limit of the draining vortex potential
(79) in the sense that (112) is the WKB limit (in this case large $l$) of the
radial scattering potential arising when not working in the WKB approximation.
The turning points $r_{tp}$ are the locations where $V(r_{tp})=0$; for $V<0$
the modes are propagating whereas for $V>0$ they will be evanescent. In
particular, for $\omega^{2}<\mu^{2}$ the solution is evanescent at infinity.

####  V.3.2 Boundary conditions

The asymptotics of the solutions can be found by converting to the tortoise
coordinate $r_{*}$ defined by,

$dr_{*}=\frac{r^{2}+a^{2}}{\Delta_{r}}dr.$ (113)

Using this coordinate, the WKB solution in (111) may be written,

$R_{\omega m}\sim(r^{2}+a^{2})^{-\frac{1}{2}}|V|^{-\frac{1}{4}}e^{\pm
i\int\sqrt{-V}dr_{*}},$ (114)

which is notably similar to (90). The extra factor at the front is because in
this case we have $A=\sqrt{r^{2}+a^{2}}\mathcal{A}$ (rather than the relation
we had for the 2D vortex which was $A=\sqrt{r}\mathcal{A}$).

Plugging in the limiting form for $V$ in the two limits $r\to r_{h}$ and
$r\to\infty$, modes with $\omega^{2}<\mu^{2}$ become,

$R_{\omega
m}\sim\frac{1}{\sqrt{r^{2}+a^{2}}}\begin{cases}A_{h}e^{-i\tilde{\omega}_{h}r_{*}},\
\quad\qquad r\to r_{h}\\\
A^{\downarrow}_{\infty}e^{-\sqrt{\mu^{2}-\omega^{2}}r_{*}},\quad
r\to\infty\end{cases}$ (115)

where we have adopted the labelling convention established in Section IV.6. We
have also defined,

$\tilde{\omega}_{h}=\omega-m\Omega_{h},$ (116)

for $\Omega_{h}$ in (110). The BC at spatial infinity has been imposed by
discarding the non-physical growing mode and on the horizon we have imposed
the absorbing BC by allowing only the in-going mode, which in terms of the WKB
solution in (114) is the one with the $-$ in the phase for
$\tilde{\omega}_{h}>0$ and $+$ in the phase for $\tilde{\omega}_{h}<0$.
Therefore the BC for the (quasi-)bound states of the Kerr spacetime can be
summarised as,

$A^{\epsilon}_{h}=A^{\uparrow}_{\infty}=0$ (117)

with $\epsilon=\mathrm{sgn}(\tilde{\omega}_{h})$.

####  V.3.3 The resonance condition

The BCs for (115) are only satisfied only for a discrete set of frequencies.
To obtain a condition for these frequencies, consider the following argument.

Since the solution propagates on the horizon and decays toward spatial
infinity there must be at least one turning point $r_{tp}$ in the intermediate
region to convert propagating modes into evanescent modes. However, if there
is only a single turning point, then the in-going mode at $r_{h}$ will be
completely reflected as it is traced back to $r_{tp}$ and there will also be
an out-going mode on the horizon, which contradicts the BC. The next
possibility is that there are an extra two turning points located somewhere
between the horizon and first turning point. Let these three be labelled
$r_{1}<r_{2}<r_{3}$ such that the solution propagates for $r<r_{1}$ and
$r_{2}<r<r_{3}$, and is evanescent for $r_{1}<r<r_{2}$ and $r>r_{3}$. The
existence of three (real) turning points is a necessary (but not sufficient)
condition for a solution which satisfies (117). Note that the next possibility
up, i.e. 5 turning points, is not possible since $K^{2}$ is a quartic
polynomial and has only 4 roots (one of which is at negative $r$ and therefore
not in the coordinate range).

To find the resonance condition, we perform another scattering computation,

$\begin{pmatrix}A^{+}_{1}\\\
A^{-}_{1}\end{pmatrix}=\left|\frac{V(r_{3})}{V(r_{1})}\right|^{\frac{1}{4}}T\begin{pmatrix}0&e^{S_{12}}\\\
e^{-S_{12}}&0\end{pmatrix}\widetilde{T}\begin{pmatrix}e^{-iS_{23}}&0\\\
0&e^{S_{23}}\end{pmatrix}T\begin{pmatrix}A^{\downarrow}_{3}\\\
A^{\uparrow}_{3}\end{pmatrix},$ (118)

with,

$S_{ij}=\int^{r_{j}}_{r_{i}}|K|\,\frac{dr}{\Delta_{r}}.$ (119)

Since we are interested in the unstable modes (which can only exist when there
is superradiance) we take the superradiant BCs in (117) to get
$A^{-}_{1}=A^{\uparrow}_{3}=0$. Inserting these into (118) we solve to find
the resonance condition,

$e^{2iS_{23}}+|\mathcal{R}|^{-1}=0,$ (120)

where the reflection coefficient $\mathcal{R}$ is the same as the one defined
in (95) but using the definition of $\tilde{\omega}_{h}$ and $S_{ij}$ in (116)
and (119) respectively. This condition is solved by complex frequencies
$\omega_{\mathbb{C}}=\omega+i\Gamma$. Assuming $|\Gamma|\ll|\omega|$ splits
(120) in two in the usual fashion,

$S_{23}(\omega_{\nu})=\pi\left(\nu+\tfrac{1}{2}\right),\qquad\Gamma_{\nu}=\frac{\log|\mathcal{R}(\omega_{\nu})|}{2\partial_{\omega}S_{23}(\omega_{\nu})},$
(121)

where $\nu=0,1...$ indexes the different bound states. The expressions in
(121) describe modes which are in a cavity for $r\in[r_{2},r_{3}]$ located
outside the potential barrier at $r\in[r_{1},r_{2}]$. As verified in Appendix
D, these modes have positive norm inside the cavity and negative norm as they
transmit into the horizon. They also no longer radiate to infinity as in the
vortex case because the field mass acts as a perfect mirror and reflects them
back in. Since modes in (121) are in a superradiant frequency band, each time
they scatter at $r_{2}$ they will be superradiantly amplified, leading to
perpetual amplification. The negative norm mode which is transmitted through
the barrier then dissipates into the horizon, lowering the mass and angular
momentum of the BH as with normal superradiant scattering.

From (121), we can already draw comparisons with the vortex case. Comparing
with (61), and realising in that case $X=|\mathcal{R}|$, the most apparent
difference is the minus sign in front of the expression for $\Gamma$, which is
absent for the present case. This is because the term in the denominator of
$\Gamma$ has the same sign as the norm of the mode in the cavity, which is
positive in the Kerr case and negative for the vortex. Hence, the overall
effect is that the sign of $\Gamma$ coincides with that of
$\log|\mathcal{R}|$.

The instability rates predicted by (121) are plotted as a function of the mass
coupling for the $m=l=1$ mode for several values of $a$ in Fig. 7. In
particular, the most unstable mode (which is the one with the largest
$\Gamma$) occurs at a mass coupling just below $\mu\sim 0.5$ when the BH is
almost maximally rotating. These findings agree qualitatively with precise
numerical computations in the literature, see e.g. Fig. 6 of [43], with the
approximation used here having a tendency to over-estimate the size of
$\Gamma$. This is expected since the WKB method is known to overestimate the
reflection coefficient in similar calculations, e.g. [25].

Figure 7: Growth rates (black lines) for the $m=l=1$ mode. The black dots
indicate the value of $\mu$ above which the instability disappears. The red
line is the small $\mu$ approximation $\Gamma\sim\mu^{9}/24$ found in [58].

####  V.3.4 Hydrogenic spectrum

It is well known that the frequency spectrum of the Kerr bound states has a
similar form to the electronic spectrum of the Hydrogen atom. Here, we will
show how this analogy can be derived from (121).

The derivation begins by noting that as $\omega\to\mu$ from below, the turning
point $r_{3}$ is sent to infinity,

$r_{3}\overset{\omega\to\mu}{\sim}\frac{2\mu^{2}}{\mu^{2}-\omega^{2}}.$ (122)

Since the cavity becomes infinitely wide in this limit, we would expect to be
able to fit an infinite number of states inside of it. For $\omega\sim\mu$,
the function $K^{2}$ is well approximated by discarding the term which is
linear in $r$ and the constant term, and the zeros $r_{2,3}$ can be
approximated as the two zeros of the truncated $K^{2}$. Using $\Delta_{r}\sim
r^{2}$ for large $r$, the integral $S_{23}$ in (120) is,

$S_{23}(\omega)\approx\frac{\pi\mu^{2}}{\sqrt{\mu^{2}-\omega^{2}}}-\pi\left(l+\tfrac{1}{2}\right),$
(123)

where a detailed derivation of this result can be found in Appendix E.
Combining (123) with the condition (121), the spectrum of bound state
frequencies close to $\mu$ is approximately given by,

$\omega_{\nu}=\mu\left(1-\frac{\mu^{2}}{2\bar{\nu}^{2}}\right),$ (124)

where the “principle quantum number” $\bar{\nu}=\nu+l+1$ is assumed large,
since $\nu$ is large as $\omega\to\mu$ and the WKB approximation assumes large
$l$. This has precisely the form of the (semi-classical) electronic spectrum
of the Hydrogen atom [59]. Note that (124) can also be found by approximating
the KG equation at large $r$, which is then exactly in the form of the
Schrödinger equation for a charged particle moving in a Coulomb potential,
e.g. [60, 59, 43].

####  V.3.5 Comparing the Kerr and vortex potentials

Figure 8: Comparison of the radial potentials for the Kerr spacetime and the
quantised vortex. For the Kerr diagram, the parameters are $\omega=0.4075$,
$m=1$ and $\mu=0.42$ and $a=0.99$, whereas for the vortex, they are
$\omega=0.3954$, $m=2$ and $\ell=2$. These parameters correspond to an
instability in both cases. In these potentials, the modes have negative norm
to the left of the peaks of $V$ and $W_{+}$ and positive norm to the right.
The instability occupies the cavity (blue region). In both cases, the
oscillation frequency is determined by the phase integral over the blue
region, that is, (61) for the vortex and (121) for the BH. The growth rate is
given by barrier integrals $S_{12}$ appearing in those equations, which are
the related to the pink areas. The difference between the two systems is that
the Kerr instability lives in the positive norm region of the flow and tunnels
into the BH whereas the vortex instability lives in the negative norm region
and tunnels out to infinity. The insets indicate the corresponding 2D profiles
on a Cartesian grid.

Already from the resonance conditions in (61) and (121), we can see that the
features of the different instabilities are determined in a similar manner by
the properties of the two systems. To expand upon the comparison, we compare
the radial potential for the BH and the vortex to see explicitly where in the
system the instability becomes trapped.

In the BH case, the potential is the one in (112) which, as seen from (114),
obeys $p_{*}^{2}+V=0$ where $p_{*}$ is the component of the wavevector in the
tortoise coordinate $r_{*}$. In the vortex case, the curves $\omega^{\pm}$
appear as the more natural objects to determine the locations where scattering
occurs in the system. However, as seen from (35a), and also discussed around
(46), the function $W_{+}$ obeys $p^{2}+W_{+}=0$ and can therefore be thought
of as an effective potential. $W_{+}$ is therefore a suitable object to
compare with the Kerr potential $V$.

In Fig. 8, the two potentials are displayed for particular choices of
parameters corresponding to an instability. These plots paint a similar
picture of the BH and non-draining vortex instabilities. Namely, the unstable
mode occupies a cavity in the effective potential (blue region) where it grows
over time, then tunnels through the barrier (pink region). In both cases, the
area in the blue determines the oscillation frequency and the pink area gives
the growth rate (really it is area of the same region under the square root of
the black curve which gives the barrier integrals $S_{12}$). The plots also
help understand the key differences. In both cases, the modes have positive
norm to the right of the barrier and a negative norm to the left. Therefore,
the Kerr instability lives in the positive norm region of the spacetime
outside the ergosphere and tunnels into the BH. By contrast, the vortex
instability occupies the centre of the system where it has a negative norm.

Finally, this picture provides a natural interpretation of the different
dynamics resulting from the two instabilities. The Kerr instability saps
rotational energy from the system, spinning the BH down until it is no longer
rotating, leaving the scalar field to eventually dissipate through the
emission of gravitational waves [61]. This is made possible since the
instability lives outside the BH. By contrast, the vortex instability grows in
the centre of the system, forcing the $\ell$-charged vortex to fragment into
$\ell$ unit charged vortices which are pushed apart, e.g. [32, 33, 34] (see
also our companion paper where we study this splitting for $\ell=2$ [49]). In
this sense, the vortex is torn apart from the inside.

####  V.3.6 Analogy with $\alpha$-decay

In Section V.3.4, we saw that the spectrum of BH bound states in the limit
$\omega\to\mu$ have an analogue in atomic physics, namely, the electronic
spectrum of the Hydrogen atom. It is natural then to enquire as to whether the
vortex instabilities also have an atomic analogue. Having seen in Fig. 8 that
Kerr instabilities live outside the BH whereas vortex instabilities live
inside the core, one might expect that the vortex states should bear some
resemblance to certain nuclear excitations. We show that this expectation is
in a sense correct, using the example of $\alpha$-decay.

The phenomenon of $\alpha$-decay involves a radioactive isotope that fragments
into an $\alpha$-particle (i.e. two protons + two neutrons with charge $2e$,
where $e$ is elementary charge, and mass $M_{\alpha}$) and a daughter nucleus
(charge $Z_{n}e$, where $Z_{n}$ is the proton number, and mass $M_{n}$). It is
known that the lifetime of the radioactive species is inversely proportional
to the square root of the energy of the ejected particle, which is the so-
called Geiger-Nutall law [62, 63]. The following simple model was proposed by
Gamow [64] to explain this dependence. The Schrödinger equation for the two
particle system can be reduced to,

$i\hbar\partial_{t}\psi=-\frac{\hbar^{2}}{2\mu}\nabla^{2}\psi+V(r)\psi,$ (125)

where $\mu=(M_{\alpha}M_{n})/(M_{\alpha}+M_{n})$ is the reduced mass, $V(r)$
is a central potential and $r$ is the distance from the centre of mass.
Gamow’s model is based on a simple potential in which the Coulomb interaction
acts to repel the particles at large $r$ whereas the nuclear force leads to a
strong attraction at small $r$. In particular,

$V=-V_{0}\,\Theta(r_{1}-r)+\frac{q}{r}\,\Theta(r-r_{1}),$ (126)

where $\Theta(r)$ is the step function, $q=2Z_{n}e^{2}>0$ is the charge
coupling and $r_{1}$ is the radius inside which the nuclear attraction
dominates. The attractive nuclear force is assumed to have a fixed value of
$V_{0}=\mathrm{const}$. This has a potential well in the region $0<r<r_{1}$
inside which particles with $E<V(r_{1})$ can become trapped. If the particle
has $E>0$, it can escape to infinity with a probability determined by the
width of the barrier. Decomposing the wavefunction as
$\psi=r^{-1}R(r)Y^{m}_{l}(\theta,\varphi)e^{-iEt/\hbar}$ with $Y^{m}_{l}$ a
spherical harmonic, this becomes,

$\partial_{r}^{2}R(r)=\frac{l(l+1)}{r^{2}}R+(V-E)R,$ (127)

where we have scaled by the reference length scale
$r_{0}=\hbar/\sqrt{2V_{0}\mu}$, which amounts to setting $\hbar,\mu$ and
$V_{0}$ equal to 1. In the following, we consider only the $l=0$ mode for
simplicity.

The excited states in the nucleus and their associated lifetimes can be
estimated using the WKB approximation. In order to employ our method, we
require that $V$ be a smooth function at the turning points. We therefore make
the replacement $\Theta(r)\to\tanh(\kappa r)$ in (126), where $\kappa$ is a
large number controlling the steepness of the step. We show an example for $V$
in Fig. 9, taking $\kappa=10$. Writing $R\sim|p|^{-\frac{1}{2}}\exp(i\int
p\,dr)$ with $p^{2}+V=E$, the WKB matrix calculation is equivalent to the one
shown in (57), with $Q=p$ and $S_{ij}=\int^{r_{j}}_{r_{i}}|p|dr$. The four
turning points $r_{i}$ are $r_{0}=0$, $r_{2}\simeq q/E$ (provided $\kappa$ is
large enough), $r_{\infty}=\infty$ and $r_{1}$ is a free parameter. To obtain
the resonance condition, we apply the no-radiation BC ($A^{-}_{\infty}=0$) and
take a Neumann BC at the origin ($A^{-}_{0}=A^{+}_{0}$). In the limit
$e^{-S_{12}}\ll 1$, the complex eigenvalue $E+i\Gamma$ is given by,

$\cos\left(S_{01}(E)+\frac{\pi}{4}\right)=0,\qquad\Gamma=\frac{\log|\mathcal{R}(E)|}{2\partial_{E}S_{01}(E)}.$
(128)

where the reflection coefficient is $|\mathcal{R}|=X^{-1}$ with $X$ defined in
(59). Note that by assuming a large barrier with $r_{1}\ll r_{2}$, one finds
that the integral over the barrier scales as $S_{12}\sim E^{-\frac{1}{2}}$ to
leading order in $r_{1}/r_{2}$, which leads to the Geiger-Nutall rule
$\log|\Gamma|\sim E^{-\frac{1}{2}}+...$.

The likeness of this resonance to the unstable vortex mode become apparent,
firstly, by comparing Fig. 9 with the vortex potential on the right panel of
Fig. 8. In both cases, the resonant mode is trapped in a cavity centred on the
origin, then tunnels through a barrier out to infinity rendering the
eigenvalue complex. This analogous confinement suggests we can think of the
vortex instability as occupying the “nucleus” of the vortex. The similarities
are further apparent by comparing the resonance conditions for the
$\alpha$-particle in (128) to those for the vortex instability in (61). The
key difference is that in the case of the $\alpha$-particle there is no
amplification mechanism, i.e. $|\mathcal{R}|<1$. The consequence of this is
that $\partial_{E}S_{01}>0$ and therefore $\Gamma<0$, meaning the
$\alpha$-resonance is stable.

Figure 9: Example of the potential in (126) for $V_{0}=1$, $q=3$, $r_{1}=3.5$
and $\Theta(r)\to\tanh(10\,r)$. There is a clear similarity with the vortex
potential shown on the right panel of Fig. 8. In this sense, the confinement
of the vortex instability is analogous to the nuclear resonance in
$\alpha$-decay.

##  VI Conclusion

In this work, we have studied the instability of quantised vortices in BECs
from the perspective of the WKB approximation. As well as yielding good
agreement with more accurate numerical computations of the instabilities, this
method allowed us to explicitly demonstrate that the unstable modes originate
inside the ergoregion. Our analysis also revealed that the instabilities of
the open system are recovered when taking the infinite limit of the closed
system. We then studied an analogue black hole by taking our BEC vortex in the
non-dispersive approximation and adding a draining component to the flow.
There, we showed that the possibility of superradiant scattering is a result
of the same tunnelling of positive norm modes to the negative norm branch of
the dispersion relation, which was the mechanism at the heart of the
instability in the non-draining case. Based on these considerations, we also
presented a stability condition for a draining vortex in an infinite (open)
condensate when the drain is not a perfect absorber. This condition will be
modified in trapped condensates where the amplified modes are reflected back
into the vortex, making mode growth possible outside the ergoregion.

In the last part of this work, we considered the BH bomb instability of a
(gravitational) Kerr BH and compared this to this vortex instability. Our
method demonstrated that both processes result from the perpetual
amplification of superradiant modes, but the location of the cavity occupied
by these modes differs. In the case of a vortex, the cavity is inside the
ergoregion whereas for the black hole, it is outside. This leads to a natural
explanation of the differing long term dynamics in the two cases: the
ergoregion instability grows inside the vortex and ultimately causes it to
fragment, whereas the BH bomb grows outside the ergoregion where it eventually
saps all the angular momentum away from the BH. Inspired by the hydrogenic
approximation for the spectrum of BH bomb modes, we also argued that the
vortex instability has a similar origin to the nuclear resonances involved in
$\alpha$-decay.

By identifying the confinement mechanisms, our results provide deeper insight
into the origin of instabilities in rotating systems. The method we presented
could also be used to better understand the properties of draining superfluid
vortices, which have been argued to be stable [35, 36]. Our analysis also has
implications for the evolution of vortex states in symmetric trapping
geometries. In particular, in our companion paper [49] we study the
fragmentation of a doubly quantised vortex that ensues from the ergoregion
instability, finding that nonlinear modulations of the vortex separation can
occur due to a back-and-forth energy exchange between sound and vortices. We
also recently studied the propagation of sound around an $\ell=29$ vortex
after it decays into a cluster of $\ell=1$ vortices [37]. An analysis of the
instability spectrum of the initial vortex, of the kind we detailed here,
could be used to gain insight into the resulting vortex configurations.

As a final note, we would like to point out that interesting dynamics
resulting from ergoregion instabilities are not unique to vortices in BECs. In
fact, it is known that classical free surface vortices exhibit rotating-
polygon instabilities in containers with the correct dimensions. It has been
argued in [65] that this instability is also connected to superradiance, or
over-reflection as it is called in fluid dynamics. There, the shape of the
polygon is related to the $m$-mode which is unstable and the system size plays
the same role in it’s ability to stabilise the system as in the case of the
closed system studied here. The similarity between our stability map and those
in [65] is a strong indication that the same mechanism is at play in both
systems and the different dynamics (polygon vs. vortex splitting) can be
understood by the different behaviour of the free surface/density profiles. It
would be interesting to apply our analysis to the polygon-instability to
understand how the dispersive nature of surface water-waves can influence the
ensuing dynamics.

## References

  * Maxwell [1890] J. C. Maxwell,  _The Scientific Papers of James Clerk Maxwell_ , Vol. 2 (University Press, 1890).
  * Bekenstein [1973] J. D. Bekenstein, Phys. Rev. D 7, 2333 (1973).
  * Bekenstein [1974] J. D. Bekenstein, Phys. Rev. D 9, 3292 (1974).
  * Hawking [1974] S. W. Hawking, Nature 248, 30 (1974).
  * Basak and Majumdar [2003a] S. Basak and P. Majumdar, Classical and Quantum Gravity 20, 3907 (2003a).
  * Basak and Majumdar [2003b] S. Basak and P. Majumdar, Classical and Quantum Gravity 20, 2929 (2003b).
  * Berti  _et al._ [2004] E. Berti, V. Cardoso, and J. P. S. Lemos, Phys. Rev. D 70, 124006 (2004).
  * Basak [2005] S. Basak, arXiv preprint gr-qc/0501097 (2005).
  * Slatyer and Savage [2005] T. R. Slatyer and C. M. Savage, Classical and Quantum Gravity 22, 3833 (2005).
  * Federici  _et al._ [2006] F. Federici, C. Cherubini, S. Succi, and M. P. Tosi, Phys. Rev. A 73, 033604 (2006).
  * Richartz  _et al._ [2015] M. Richartz, A. Prain, S. Liberati, and S. Weinfurtner, Phys. Rev. D 91, 124018 (2015).
  * Churilov and Stepanyants [2019] S. Churilov and Y. Stepanyants, Phys. Rev. Fluids 4, 034704 (2019).
  * Demirkaya  _et al._ [2020] B. Demirkaya, T. Dereli, and K. Güven, Physica Scripta 95, 055001 (2020).
  * Barcelo  _et al._ [2011] C. Barcelo, S. Liberati, and M. Visser, Living Reviews in Relativity 14, 3 (2011).
  * Unruh [1981] W. G. Unruh, Phys. Rev. Lett. 46, 1351 (1981).
  * Weinfurtner  _et al._ [2011] S. Weinfurtner, E. W. Tedford, M. C. J. Penrice, W. G. Unruh, and G. A. Lawrence, Phys. Rev. Lett. 106, 021302 (2011).
  * Euvé  _et al._ [2016] L. P. Euvé, F. Michel, R. Parentani, T. G. Philbin, and G. Rousseaux, Phys. Rev. Lett. 117, 121301 (2016).
  * Steinhauer [2016] J. Steinhauer, Nat. Phys. 12, 959 (2016).
  * de Nova  _et al._ [2019] J. R. M. de Nova, K. Golubkov, V. I. Kolobov, and J. Steinhauer, Nature 569, 688 (2019).
  * Kolobov  _et al._ [2019] V. I. Kolobov, K. Golubkov, J. R. M. de Nova, and J. Steinhauer, arXiv preprint arXiv:1910.09363 (2019).
  * Torres  _et al._ [2017] T. Torres, S. Patrick, A. Coutant, M. Richartz, E. W. Tedford, and S. Weinfurtner, Nat. Phys. 13, 833 (2017).
  * Brito  _et al._ [2020] R. Brito, V. Cardoso, and P. Pani,  _Superradiance_ (Springer, 2020).
  * Dicke [1954] R. H. Dicke, Phys. Rev. 93, 99 (1954).
  * Patrick and Weinfurtner [2020] S. Patrick and S. Weinfurtner, Phys. Rev. D 102, 084041 (2020).
  * Patrick [2020a] S. Patrick, arXiv preprint arXiv:2011.04620 (2020a).
  * Fetter and Svidzinsky [2001] A. L. Fetter and A. A. Svidzinsky, Journal of Physics: Condensed Matter 13, R135 (2001).
  * Feynman [1955] R. Feynman (Elsevier, 1955) pp. 17–53.
  * Onsager [1949] L. Onsager, Il Nuovo Cimento (1943-1954) 6, 279 (1949).
  * Fetter [2009] A. L. Fetter, Rev. Mod. Phys. 81, 647 (2009).
  * Barenghi and Parker [2016] C. F. Barenghi and N. G. Parker,  _A primer on quantum fluids_ (Springer, 2016).
  * Takeuchi  _et al._ [2018] H. Takeuchi, M. Kobayashi, and K. Kasamatsu, Journal of the Physical Society of Japan 87, 023601 (2018).
  * Shin  _et al._ [2004] Y. Shin, M. Saba, M. Vengalattore, T. A. Pasquini, C. Sanner, A. Leanhardt, M. Prentiss, D. Pritchard, and W. Ketterle, Phys. Rev. Lett. 93, 160406 (2004).
  * Isoshima  _et al._ [2007] T. Isoshima, M. Okano, H. Yasuda, K. Kasa, J. Huhtamäki, M. Kumakura, and Y. Takahashi, Phys. Rev. Lett. 99, 200403 (2007).
  * Okano  _et al._ [2007] M. Okano, H. Yasuda, K. Kasa, M. Kumakura, and Y. Takahashi, J. Low Temp. Phys. 148, 447 (2007).
  * Zezyulin and Konotop [2014] D. A. Zezyulin and V. V. Konotop, arXiv preprint arXiv:1411.7948 (2014).
  * Alperin and Berloff [2021] S. N. Alperin and N. G. Berloff, Optica 8, 301 (2021).
  * Geelmuyden  _et al._ [2021] A. Geelmuyden, S. Erne, S. Patrick, C. Barenghi, and S. Weinfurtner, arXiv preprint arXiv:2105.11509 (2021).
  * Johnstone  _et al._ [2019] S. P. Johnstone, A. J. Groszek, P. T. Starkey, C. J. Billington, T. P. Simula, and K. Helmerson, Science 364, 1267 (2019).
  * Oliveira  _et al._ [2014] L. A. Oliveira, V. Cardoso, and L. C. B. Crispino, Phys. Rev. D 89, 124008 (2014).
  * Oliveira  _et al._ [2018] L. A. Oliveira, L. J. Garay, and L. C. B. Crispino, Phys. Rev. D 97, 124063 (2018).
  * Schutz and Friedman [1975] B. F. Schutz and J. L. Friedman, The Astrophysical Journal 199, L157 (1975).
  * Cardoso  _et al._ [2008] V. Cardoso, P. Pani, M. Cadoni, and M. Cavaglià, Phys. Rev. D 77, 124044 (2008).
  * Dolan [2007] S. R. Dolan, Phys. Rev. D 76, 084001 (2007).
  * Giacomelli and Carusotto [2020] L. Giacomelli and I. Carusotto, Phys. Rev. Research 2, 033139 (2020).
  * Berry and Mount [1972] M. V. Berry and K. E. Mount, Reports on Progress in Physics 35, 315 (1972).
  * Bühler [2014] O. Bühler,  _Waves and mean flows_ (Cambridge University Press, 2014).
  * Torres  _et al._ [2018] T. Torres, A. Coutant, S. Dolan, and S. Weinfurtner, Journal of Fluid Mechanics 857, 291 (2018).
  * Patrick [2020b] S. Patrick, arXiv preprint arXiv:2007.06671 (2020b).
  * Patrick  _et al._ [2021] S. Patrick, A. Geelmuyden, S. Erne, C. F. Barenghi, and S. Weinfurtner, arXiv preprint arXiv:2111.02567 (2021).
  * Bogolyubov [1947] N. N. Bogolyubov, J. Phys. (USSR) 11, 23 (1947).
  * Pethick and Smith [2008] C. J. Pethick and H. Smith,  _Bose–Einstein condensation in dilute gases_ (Cambridge University Press, 2008).
  * Andersen  _et al._ [2003] A. Andersen, T. Bohr, B. Stenum, J. J. Rasmussen, and B. Lautrup, Phys. Rev. Lett. 91, 104502 (2003).
  * Stepanyants and Yeoh [2008] Y. A. Stepanyants and G. H. Yeoh, Journal of Fluid Mechanics 604, 77 (2008).
  * Torres [2020] T. Torres, Philosophical Transactions of the Royal Society A 378, 20190236 (2020).
  * Visser [2007] M. Visser, arXiv preprint arXiv:0706.0622 (2007).
  * Unruh [1974] W. G. Unruh, Phys. Rev. D 10, 3194 (1974).
  * Ford [1975] L. H. Ford, Phys. Rev. D 12, 2963 (1975).
  * Zouros and Eardley [1979] T. J. Zouros and D. M. Eardley, Annals of Physics 118, 139 (1979).
  * Detweiler [1980] S. Detweiler, Phys. Rev. D 22, 2323 (1980).
  * Fox [2006] M. Fox,  _Quantum optics: an introduction_ , Vol. 15 (OUP Oxford, 2006).
  * East [2018] W. E. East, Phys. Rev. Lett. 121, 131104 (2018).
  * Geiger and Nuttall [1911] H. Geiger and J. Nuttall, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 22, 613 (1911).
  * Ren and Ren [2012] Y. Ren and Z. Ren, Phys. Rev. C 85, 044608 (2012).
  * Gamow [1928] G. Gamow, Zeitschrift für Physik 51, 204 (1928).
  * Mougel  _et al._ [2017] J. Mougel, D. Fabre, L. Lacaze, and T. Bohr, Journal of Fluid Mechanics 824, 230 (2017).

##  Appendix A Density of an $\ell$-charged vortex

To solve (10) for $Y(r)$, we use a shooting method. First we rewrite the
second order differential equation for $Z$ as two couple first order
equations,

$\begin{split}Y^{\prime}=&\ Z,\\\ Z^{\prime}=&\
-\frac{1}{r}Z-\left(2-\frac{\ell^{2}}{r^{2}}\right)Y+2Y^{3},\end{split}$ (129)

with ${}^{\prime}=\partial_{r}$. We evolve these equations from the starting
point $r_{1}$ using a standard differential equation solver (we used Matlab’s
inbuilt function ode45). Since setting $r_{1}=0$ causes a divergence in the
equations, we specify our initial conditions a small distance from the origin
$r_{1}=\varepsilon$. Using the asymptotics in (11), these are chosen to be
$Y(\varepsilon)=Y_{0}J_{\ell}(\varepsilon)$ and
$Y^{\prime}(\varepsilon)=Y_{0}J^{\prime}_{\ell}(\varepsilon)$.

The value of $Y_{0}$ is iteratively improved until the solution converges to
the desired value at $r_{N_{r}}$. Since we are ultimately interested in $Z(r)$
profiles in infinite systems, one would ideally set $r_{N_{r}}$ to some large
value and find the value of $Y_{0}$ such that $Z_{N}=0$. However, the shooting
problem for (10) is inherently unstable and if $r_{N}$ is chosen too large
then the value of $Y_{0}$ needs to be specified to an unworkable precision.
Instead, we look for the solution satisfying a Dirichlet BC, $Y_{N}=0$, which
corresponds to the wavefunction of a condensate inside a trapping potential of
the form $U(r<r_{N})=0$ and $U(r=r_{N})=\infty$. Once obtained, we then use
(10) to find the leading terms in the Taylor expansion of $Y$ in the limit
$r\to\infty$,

$\begin{split}Y(r)\overset{\infty}{\sim}&\
1-\frac{\ell^{2}}{4r^{2}}-\left(1+\frac{\ell^{2}}{8}\right)\frac{\ell^{2}}{4r^{4}}\\\
&\qquad-\left(8+2\ell^{2}+\frac{\ell^{4}}{16}\right)\frac{\ell^{2}}{8r^{6}}+\mathcal{O}\left(\frac{1}{r^{8}}\right).\end{split}$
(130)

This is stitched onto the numerical solution using a ramp function over the
range $r\in[a_{1},b_{1}]$. In order to guarantee the correct limiting
behaviour at the origin, we also stitch our numerical solution onto the
asymptotic solution $BJ_{\ell}(r)$ at the origin, using a ramp function over
the range $r\in[a_{2},b_{2}]$, with the parameter $B$ chosen to best fit the
numerical solution in the region $0<r<b_{2}$.

In our solver, we set the initial value of $Y_{0}$ to 1. The effects of
different initial points were checked using the values
$\varepsilon=10^{-4},10^{-3},10^{-2}$ and no noticeable difference was found.
We set the values $r_{N}=20$, $a_{1}=7$ and $b_{1}=15$ to minimise the time
taken to find a convergent solution in the numerics whilst still yielding a
good match in the stitching region. Near the origin, the stitching was
performed using $a_{2}=0.25$ and $b_{2}=0.5$.

##  Appendix B WKB amplitude

In this appendix, we derive the expression for the WKB amplitude in (37). The
WKB modes and their $y$ derivatives up to fourth order are,

$\begin{split}f&\ =\mathcal{A}e^{i\int qdy},\\\ f^{\prime}&\
=\left(iq\mathcal{A}+\mathcal{A}^{\prime}\right)e^{i\int qdy},\\\
f^{\prime\prime}&\
=\left(-q^{2}\mathcal{A}+iq^{\prime}\mathcal{A}+2iq\mathcal{A}^{\prime}\right)e^{i\int
qdy}+...\\\ f^{\prime\prime\prime}&\
=\left(-iq^{3}\mathcal{A}-3qq^{\prime}\mathcal{A}-3q^{2}\mathcal{A}^{\prime}\right)e^{i\int
qdy}+...\\\ f^{\prime\prime\prime\prime}&\
=\left(q^{4}\mathcal{A}-6iq^{2}q^{\prime}\mathcal{A}-4iq^{3}\mathcal{A}^{\prime}\right)e^{i\int
qdy}+...\end{split}$ (131)

where only terms involving up to one derivative on either $q$ or $\mathcal{A}$
have been retained, since these are the terms which contribute up to next-to-
leading order (NTLO) in the WKB expansion. Using just the first term from each
of (131) gives the leading order (LO) equation in (33). The NTLO terms give an
equation for the evolution of the amplitude. To find this equation, we can
combine the two equations in (29) into a single equation for $f$. The NTLO
contributions are,

$\begin{split}&\zeta
f^{\prime\prime\prime\prime}+2\zeta^{\prime}f^{\prime\prime\prime}-\left(2\widetilde{m}^{2}+4\rho
r^{2}\right)\zeta f^{\prime\prime}\\\
&\qquad\qquad\qquad\qquad-2(\zeta\widetilde{m}^{2})^{\prime}f^{\prime}+...=0,\end{split}$
(132)

where $\zeta=1/\Omega r^{2}$ has been defined for conciseness. Using the
expressions in (131) and extracting the NTLO terms, one finds,

$\begin{split}0=&\ 6\zeta q^{2}q^{\prime}\mathcal{A}+4\zeta
q^{3}\mathcal{A}^{\prime}+2\zeta^{\prime}q^{3}\mathcal{A}+2(\zeta\widetilde{m}^{2})^{\prime}q\mathcal{A}\\\
&\qquad\qquad\quad+\left(2\widetilde{m}^{2}+4\rho
r^{2}\right)\zeta\left(q^{\prime}\mathcal{A}+2q\mathcal{A}^{\prime}\right),\\\
=&\ \mathcal{A}^{-1}\big{[}\partial_{y}(2\zeta
q^{3}\mathcal{A}^{2})+\partial_{y}(2\zeta\widetilde{m}^{2}q\mathcal{A}^{2})\\\
&\qquad\qquad\quad+4\rho
r^{2}\zeta\partial_{y}(q\mathcal{A}^{2})\big{]},\end{split}$ (133)

where in going to from the first to the second line, we have multiplied
everything by $\mathcal{A}$ and collected terms under the $y$ derivatives. The
LO equation in (33) can be written in the form,

$\zeta^{2}\kappa^{4}+4\rho
r^{2}\zeta^{2}\kappa^{2}=4,\qquad\kappa^{2}=q^{2}+\widetilde{m}^{2}.$ (134)

Using this notation and multiplying the equation above by
$\kappa^{2}\zeta\mathcal{A}$ we find,

$\begin{split}0=&\
2\zeta\kappa^{2}\partial_{y}(\zeta\kappa^{2}q\mathcal{A}^{2})+4\rho
r^{2}\zeta^{2}\kappa^{2}\partial_{y}(q\mathcal{A}^{2}),\\\ =&\
2\zeta\kappa^{2}\partial_{y}(\zeta\kappa^{2}q\mathcal{A}^{2})-\zeta^{2}\kappa^{4}\partial_{y}(q\mathcal{A}^{2})+4\partial_{y}(q\mathcal{A}^{2}),\\\
=&\
\partial_{y}(\zeta^{2}\kappa^{4}q\mathcal{A}^{2})+4\partial_{y}(q\mathcal{A}^{2}),\end{split}$
(135)

Dividing by 4 and using $q=pr$ and $\kappa=kr$, the right hand side becomes,

$\begin{split}\partial_{y}\left[\left(1+\tfrac{1}{4}\zeta^{2}\kappa^{4}\right)q\mathcal{A}^{2}\right]=&\
\partial_{y}\left[\left(1+\frac{k^{2}}{4F}\right)rp\mathcal{A}^{2}\right],\\\
=&\
r\partial_{r}\left[\left(\frac{\rho+k^{2}/2}{\rho+k^{2}/4}\right)rp\mathcal{A}^{2}\right],\end{split}$
(136)

from which we obtain the result for the amplitude (37) in the main text.

##  Appendix C Simulations of the BdG equation

We evolve the BdG equation (12) using a Method of Lines approach. First, we
define a discrete $(r,t)$ grid with spacing $\Delta r=0.4$ and $\Delta
t=\Delta r^{2}/2$. We then write
$\partial_{t}|U\rangle=-i\widehat{L}|U\rangle$, where the derivative terms in
$\widehat{L}$ (i.e. $\partial_{r}$ and $\partial_{r}^{2}$) are expressed as
5-point (centred) finite difference stencils which encode the correct boundary
conditions at $r=0$ (as discussed in Sec. IV.7). We evolve in time using a
trapezium method, i.e.

$|U(t+\Delta t)\rangle=\left(\frac{1-i\widehat{L}\Delta
t/2}{1+i\widehat{L}\Delta t/2}\right)|U(t)\rangle,$ (137)

where the term on the denominator is understood as a multiplication by the
inverse matrix. For our initial condition, we take
$u_{+}(t=0)=e^{-(r-r_{p})^{2}/2\sigma}$, with $r_{p}=15$ and $\sigma=2$, and
$u_{-}(t=0)=0$. The outer boundary (where we impose Dirichlet BCs on both
$u_{+}$ and $u_{-}$) is located at $r_{D}=400$, which is sufficiently far away
that no reflections re-enter our window of analysis. To extract the frequency
content, we consider the signal $S(t)=u_{+}(r_{c},t)$, where $r_{c}$ is the
point inside the vortex core where the wave amplitude is largest. We compute
the Fourier transform $\tilde{S}(\omega)=\mathscr{F}[H(t)S(t)]$, where $H(t)$
is the Hamming window function which minimises the secondary lobes in
$\tilde{S}(\omega)$. The locations of the peaks of $\tilde{S}(\omega)$, which
are obtained via parabolic interpolation, give the real part of the unstable
frequency. To find the imaginary part, we fit a straight line to $\log|S(t)|$
in the region where the instability dominates the signal. It is not always
possible to extract the imaginary part, for example, when two instabilities
inference with each other or an instability interferes with a slowly decaying
mode.

##  Appendix D Klein-Gordon norm

The Klein-Gordon equation in (105) can be derived by minimising the following
action,

$\mathcal{S}_{\mathrm{KG}}=\int
d^{3}\mathbf{x}\sqrt{-g}\left(-\partial^{\nu}\phi\,\partial_{\nu}\phi^{*}-\mu^{2}\phi\,\phi^{*}\right).$
(138)

Since this invariant under phase rotations of $\phi$, it admits a conserved
norm,

$j^{\nu}=\sqrt{-g}g^{\nu\sigma}i\left(\phi\,\partial_{\sigma}\phi^{*}-\phi^{*}\partial_{\sigma}\phi\right),$
(139)

satisfying the conservation equation $\partial_{\nu}j^{\nu}=0$. Using the
decomposition of $\phi$ in (106) the main text with the angular functions
normalised as,

$2\pi\int^{\pi}_{0}\sin\theta\,S^{2}_{\omega m}(\theta)\,d\theta=1,$ (140)

the $j^{r}$ component is fixed with $r$, leading to,

$\partial_{r}\left(\Delta_{r}[R_{\omega m}\partial_{r}R^{*}_{\omega
m}-R^{*}_{\omega m}\partial_{r}R_{\omega m}]\right)=0,$ (141)

where the quantity in round parentheses is the conserved Wronskian associated
with (108).

The time component of (139) $j^{t}$ is the norm (density). In Section IV.5 it
was possible to find a concise expression for the norm of the vortex
excitations using the WKB approximation, which then allowed us to deduce the
sign of the norm in different regions of the flow. Things are not as simple
here, since the $\theta$ dependence in the $S_{\omega m}$ functions contained
in $\phi$ complicates the expressions. Instead, we consider the following
argument.

First, we plug the metric functions (100) and the $\phi$ decomposition (106)
into the expression for $j^{t}$ which results in,

$\begin{split}j^{t}=&\
\frac{2}{\Delta_{r}}\Big{[}\left\\{(r^{2}+a^{2})^{2}-\Delta_{r}a^{2}\sin^{2}\theta\right\\}\omega\\\
&\qquad\quad-
am\left(r^{2}+a^{2}-\Delta_{r}\right)\Big{]}\sin\theta\,|\phi|^{2}.\end{split}$
(142)