# Quantum vortex stability in draining fluid flows

Sam Patrick<EMAIL_ADDRESS>Department of Physics, King’s College
London, University of London, Strand, London, WC2R 2LS, UK

###### Abstract

Quantum vortices with more than a single circulation quantum are usually
unstable and decay into clusters of smaller vortices. One way to prevent the
decay is to place the vortex at the centre of a convergent (draining) fluid
flow, which tends to force vortices together. It is found that whilst the
primary splitting instability is suppressed in this way (and completely
quenched for strong enough flows) a secondary instability can emerge in
circular trapping geometries. This behaviour is related to an instability of
rotating black holes when superradiantly amplified waves are confined inside a
reflective cavity. The end state of the secondary instability is dramatic,
manifesting as a shock wave that propagates round the circular wall and
nucleates many more vortices.

##  I Introduction

Quantum vortices are topological defects in the order parameter
$\Psi(\mathbf{x},t)$ describing a quantum fluid, e.g. Bose-Einstein
condensates (BECs), polariton fluids, superfluid 4He, photon superfluids etc.
Their dynamics is of widespread interest due to their role in various
important phenomena, including superfluid turbulence [1], topological phase
transitions [2] and superconductivity [3]. The circulation of a quantum vortex
is $\ell\kappa$, that is, $\ell$ units of circulation quantum $\kappa=h/M$,
with $h$ Planck’s constant and $M$ the mass of (bosonic) particles in the
fluid. The integer $\ell$ counts the number of times the phase
$\mathrm{arg}(\Psi)$ winds around the defect in the order parameter. The most
basic quantum vortex has $|\ell|=1$ and is called a singly quantised vortex
(SQV), whereas a vortex with $|\ell|>1$ is called a multiply-quantised vortex
(MQV).

Since the energy of a vortex scales with $\ell^{2}$, MQVs tend to split into
clusters of SQVs as a result of an energetic instability. Hence, any mechanism
that provides a reservoir to dissipate energy (e.g. coupling to a thermal
bath) will cause an MQV to decay [4]. Furthermore, even in ideal conservative
systems, vortices can couple to sound waves (phonons) in the fluid to produce
a dynamical instability [5, 6]. Because of this, many studies of quantum
vortices typically employ configurations of SQVs.

There are, however, certain conditions under which MQVs do not decay. For
example, energetic favorability of MQVs is found in rotating condensates with
anharmonic (quartic) traps [7] or pinning potentials [8]. Furthermore, the
dynamical splitting instability can be suppressed in finite size systems which
possess a discrete phonon spectrum [9]. Another stabilisation mechanism
involves filling the vortex core with a second condensate, such that there is
an energy barrier associated with splitting the MQV [10]. A well known
property of vortices (which follows from Kelvin’s circulation theorem [11]) is
that they tend to move under the influence of the local fluid velocity. Hence,
a natural stabilising mechanism involves placing an MQV at the focal point of
a convergent fluid flow, e.g. by draining fluid from the system in a small
region in the centre. If the drain is strong enough, one would expect the
convergence of the fluid to override the tendency of the MQV to split,
rendering the vortex stable against fragmentation.

Convergent fluid flows in superfluid 4He were realised in [12, 13, 14] and the
flow was reported to contain $\mathcal{O}(10^{4})$ vortex quanta. Motivated by
these experiments, [15] simulated the evolution of vortex filaments around the
drain at non-zero temperature, unveiling a formation mechanism for the vortex
bundle above the drain. 3D simulations based on a Gross-Pitaevskii model
containing a low number of vortices were performed in [16], which showed the
convergent flow tends to focus vortices over the drain hole, twisting them
into a bundle that funnels the flow into the drain below. More recently, a
novel method of surface wave spectroscopy was used to constrain the core size
of a macroscopic draining vortex containing $\mathcal{O}(10^{4})$ vortex
quanta [17].

MQV stability has also been demonstrated in exciton-polariton condensates,
systems with inherent particle losses due to the finite lifetime of the
polariton quasiparticles [18]. Simulations in [19] showed that a ring-shape
laser beam induces particle fluxes toward the centre, capable of stabilising
the MQV situated there. Furthermore, wavevector dependent losses can induce a
convergent flow centred on each vortex, leading to an attractive force between
that causes vortices to merge [20]. Stability of an $\ell=15$ vortex was
recently demonstrated in the experimental set-up of [21], although in that
case, the fluid velocity is forced by the pump laser and the convergent flow
is weak since the polariton decay is not restricted to a central region.

In atomic BECs, convergent flows can be realised using an electron beam to
ionise particles and deplete the condensate in a localised region. Numerical
simulations in [22] showed that an $\ell=2$ vortex can be stabilised by strong
localised dissipation, and that the vortex splits once dissipation is switched
off. Giant vortex aggregates were reported in a rotating condensate in which a
focused laser beam was used to remove atoms on the rotation axis [23]. These
vortex aggregates were attributed to the combined effect of the Coriolis force
and the radial influx of atoms.

In light of these results, one may expect that MQVs are unambiguously
stabilised under a strong enough convergent flow. The purpose of this work is
to demonstrate that this is not necessarily the case. It is shown that even
when the primary splitting instability is quenched, secondary instabilities
capable of driving the system far from equilibrium can arise. These secondary
instabilities are related to the black hole bomb instability of Press and
Teukolsky [24]. A rotating black hole is surrounded by an ergoregion, where
fluctuations in a certain frequency range possess negative energies. Incoming
radiation can scatter into these negative energy states, causing the escaping
fraction of radiation to become amplified. This process is called rotational
superradiance [25]. The black hole bomb proposal involves encircling the
system with a reflective mirror that scatters amplified radiation back into
the black hole where it is further amplified, leading to a runaway process. A
quantum vortex, which from the perspective of low frequency sound waves can be
viewed as an analogue rotating spacetime [26], also possesses an ergosphere
inside its core. Indeed, the MQV splitting instability can be viewed as a
negative energy bound state inside the ergosphere, which couples to sound
waves outside the vortex to release energy [9, 27, 28]. Dissipation inside the
vortex core leads to damping of this negative energy mode, akin the classical
formulation of black hole superradiance where the negative energy is absorbed
by the horizon [29]. We will show that, when the vortex is located at the
centred of a circularly symmetric trapping potential, superradiantly amplified
sound can trigger the black hole bomb instability.

The structure of this paper is as follows. In Sections II and III, we describe
the stationary draining vortex profiles of a quantum fluid. In Section IV, we
analyse the eigenmodes of the system for some illustrative parameters,
demonstrating the effect of the convergent flow on stability. Section V
provides an interpretation of the instabilities using a WKB approximation of
the resonance formula for unstable eigenmodes. In Section VI, numerical
simulations of the fully nonlinear equations of motion are described,
revealing the dramatic effect of the secondary instability on the late time
dynamics. Section VII concludes and implications of the results are discussed.

##  II Set-up

For concreteness, we consider a tightly confined atomic BEC with quasi-2D
dynamics in the $\mathbf{x}=(x,y)$–plane. In the mean field approximation, the
system is governed by the Gross-Pitaevskii equation (GPE),

$i\hbar\partial_{t}\Psi=\left[-\frac{\hbar^{2}\nabla^{2}}{2M}+V(\mathbf{x})-i\Gamma(\mathbf{x})+g|\Psi|^{2}-\mu\right]\Psi,$
(1)

where $\Psi$ is the order parameter, $M$ is mass of the particles (which are
bosons), $\mu$ is the chemical potential, $g$ is the 2D interaction parameter,
$V(\mathbf{x})$ is a trapping potential and $\Gamma(\mathbf{x})$ is a
dissipation term. In the context of the BEC, localised dissipation can be
achieved by illuminated the condensate with an ionising beam [30]. The GPE
applies equally to polariton condensates [21], photon superfluids [31] and can
even be considered as a phenomenological model of superfluid 4He on small
scales [32]. We consider traps with a hard circular wall located $r=r_{B}$.
Specifically, we use the function,

$V=\frac{V_{0}}{1+(V_{0}-1)e^{a(r_{B}-r)}},$ (2)

where $V_{0}$ is the height of the trap and $a$ is a smoothing parameter (we
mostly use $V_{0}=a=5$ and $r_{B}=25$ unless specified otherwise). In the
Madelung representation, $\Psi=\sqrt{n}e^{iM\Phi/\hbar}$ where $n(\mathbf{x})$
is the particle number density and $\Phi(\mathbf{x})$ a phase whose gradient
gives the velocity field $\mathbf{v}=\bm{\nabla}\Phi$. The GPE then separates
into a Bernoulli equation,

$\partial_{t}\Phi+\frac{1}{2}\mathbf{v}^{2}+\frac{gn+V(\mathbf{x})-\mu}{M}=\frac{\hbar^{2}}{2M^{2}}\frac{\nabla^{2}\sqrt{n}}{\sqrt{n}},$
(3)

where the term on the right is the quantum pressure modification and a
continuity equation,

$\partial_{t}n+\bm{\nabla}\cdot\left(n\mathbf{v}\right)=-2\Gamma(\mathbf{x})n.$
(4)

From the last equation, it is apparent that the function $\Gamma$ acts a sink
for the density. We consider a dissipation function localised at the centre of
the trap of the form,

$\Gamma=\Gamma_{0}\frac{1+\tanh[2(r_{0}-r)]}{2},$ (5)

where $\Gamma_{0}>0$ is the strength of dissipation and $r_{0}$ is the radius
of region where particles are lost. When $\Gamma_{0}=0$, the GPE conserves the
total number of particles in the condensate $N=\int d^{2}\mathbf{x}\,n$. For
finite $\Gamma_{0}$, $N$ decreases according to,

$\partial_{t}N=-2\int d^{2}\mathbf{x}\,\Gamma n.$ (6)

Particles then flow inward to replenish the depleted condensate inside the
dissipation region, establishing a convergent fluid flow. This is the quantum
analogue of the classical draining bathtub flow that forms over a physical
fluid outlet [33, 34] and the dissipation region can be viewed at the analogue
of the drain hole. In order to establish a stationary state in the simulations
of the following sections, particles will be resupplied near the edge of the
trap to hold $N$ fixed.

##  III Vortices

We search for stationary states containing vortices at the centre of the trap.
The condensate phase and velocity field are,

$\Phi=\frac{\hbar}{M}\left[\ell\theta+\Theta(r)\right],\quad\mathbf{v}=\frac{\hbar}{M}\left[\frac{\ell}{r}\hat{\mathbf{e}}_{\theta}+\partial_{r}\Theta\hat{\mathbf{e}}_{r}\right],$
(7)

where $\mathbf{v}$ has the profile of a draining vortex. The radial part of
the phase $\Theta(r)$ is an unknown function that must be solved for. Note
that our assumption of a circular trap with the draining vortex located at the
centred implies the stationary background state described by $(n,\mathbf{v})$
is a function of $r$ only.

We start by considering the non-draining case when $\Gamma_{0}=0$ and
$v_{r}\equiv\mathbf{e}_{r}\cdot\mathbf{v}=0$. In that case, the stationary
solutions are the typical quantum vortices, whose velocity field is fully
determined by the constraint of an irrotational fluid flow and the density $n$
is the only function one has to solve for. The size of the vortex core is set
by the healing length,

$\xi_{0}=\hbar/\sqrt{M\mu_{0}},$ (8)

where $\mu_{0}=\mu(\Gamma_{0}=0)$. We can then simplify the calculation by
defining the following dimensionless quantities,

$\begin{split}\frac{\mathbf{x}}{\xi_{0}}\to&\
\mathbf{x},\quad\frac{\mu_{0}t}{\hbar}\to t,\quad\frac{\mu}{\mu_{0}}\to\mu,\\\
&\frac{gn}{\mu_{0}}\to n,\quad\frac{M\Phi}{\hbar}\to\Phi.\end{split}$ (9)

In these units, the chemical potential for the non-draining vortex is unity by
definition. In equilibrium, the density is determined by the equation,

$\frac{1}{2}\nabla_{r}^{2}z+\left[1-\frac{\ell^{2}}{2r^{2}}-V(r)\right]z=z^{3},$
(10)

where $z=\sqrt{n}$ and $\nabla^{2}_{r}=\frac{1}{r}\partial_{r}r\partial_{r}$.

The presence of a radial component of $\mathbf{v}$ complicates the flow
pattern. To simplify the expressions, we work with the dimensionless
quantities defined in (9). In equilibrium, (4) implies
$\partial_{r}(rnv_{r})=-2r\Gamma n$ which is solved by,

$\begin{split}v_{r}(r)=&\
-\frac{2\int_{0}^{r}dr^{\prime}r^{\prime}\Gamma(r^{\prime})n(r^{\prime})}{n(r)r},\\\
=&\ \begin{cases}-D/r,\qquad\qquad\,r\ll r_{0}\\\ -\Gamma_{0}r/(\ell+1),\quad
r\gg r_{0},\end{cases}\end{split}$ (11)

where $D=2\int dr\,r\Gamma n/\mu$ is the drain constant and the integral is
performed over the region where $\Gamma(r)$ has support. Hence, approaching
the centre from the boundary, the radial velocity magnitude increases as the
fluid converges on the dissipation region. Inside the dissipation region,
particles are lost and mass conservation causes a decrease in the radial
velocity magnitude, with $v_{r}$ vanishing at the very centre. Using these
properties of $v_{r}$, (3) implies the following asymptotics for the density,

$\begin{split}n(r)=&\ \begin{cases}\mu-(\ell^{2}+D^{2})/2r^{2},\qquad\quad
r\gg r_{0}\\\ r^{2\ell}\left[A+Br^{2}+\mathcal{O}(r^{4})\right],\quad r\ll
r_{0}\end{cases}\\\ B=&\ A[A\delta(\ell)-\mu]/(\ell+1),\end{split}$ (12)

where the value of $A$ can only be determined from the full solution. For
$\ell=0$, $A>0$ since it is the value of the density at $r=0$. It will turn
out that $A<\mu$ so that $B<0$, i.e. $n$ has negative curvature at $r=0$. This
implies that the density is a non-monotonic function of $r$. When $\ell\neq
0$, $A$ is also positive since it determines the rate at which the density
rises from zero at the vortex centre.

To search for equilibrium states when $\Gamma_{0}\neq 0$, we apply the
following prescription. Starting from the non-draining solution, we switch on
a finite value of $\Gamma_{0}$ so that the density in the region $r<r_{0}$
begins to deplete. We then resupply the depleted particles near the edge of
the trap. The over density near $r_{B}$ will cause the condensate to flow
inward toward the centre, establishing a convergent radial flow. The fixed
particle number $N$ is accommodated by a shift in the chemical potential. The
details of this procedure are described in Appendix A.1.

Figure 1: Stationary draining flow profiles for $r_{0}=5$. In the top panels,
we display the density $n$ and radial velocity magnitude $|v_{r}|$ for three
illustrative values of $\Gamma_{0}$. On the lower panels, we show the
variation of the chemical potential $\mu$, the drain rate $D$ and the horizon
radii $r_{h}$ with $\Gamma_{0}$.

In Fig. 1, we display some examples of the density and radial velocity
profiles when there is no vortex present. The density exhibits a depression
moving toward the centre followed by a sudden rise at the very centre. This is
the behaviour anticipated in (12). The magnitude of the radial velocity
displays the properties anticipated in (11), although the behaviour in the
intermediate region can be intricate due to the non-monotonic nature of the
density. As expected, the chemical potential adjusts as we increase
$\Gamma_{0}$ but keep $N$ fixed. The general increase in $\mu$ results from
$v_{r}$ lowering the density in the dissipation region, thereby pushing atoms
into the bulk and raising the asymptotic density level. However, the variation
of $D(\Gamma_{0})$ reveals a curious feature that beyond a certain
$\Gamma_{0}$, stronger dissipation can weaken radial flow, i.e. the more one
tries to drain the fluid, the less one manages to achieve it. This is related
to the macroscopic Zeno effect studied in [22]. We also display the locations
$r_{h}$ satisfying $(n-v_{r}^{2})|_{r=r_{h}}=0$. In the region where $n$ is
close to its asymptotic value, this represents the location where radial flow
exceeds the sound speed $\sqrt{n}$, i.e. the acoustic horizon [35]. Note that,
since $v_{r}$ decreases toward the origin, the horizon necessarily comes as a
black hole-white hole horizon pair.

In Fig. 2, we illustrate the effect of placing a vortex (in this case one with
$\ell=2$) in the centre of the dissipation region. Since the density drops to
zero on the vortex axis, there are fewer particles present in the dissipation
region to be damped, hence, the radial velocity is smaller when compared with
the $\ell=0$ case. The black lines represent $r_{0}=5$ (same as Fig. 1)
whereas the red lines are for $r_{0}=2$, illustrating that the draining flow
is weaker when the dissipation region is smaller as one would expect.

Figure 2: In the top two panels, we show the density and radial velocity
profiles when there is an $\ell=2$ located at $r=0$. The legend gives the
dissipation function parameters in the format $(r_{0},\Gamma_{0})$. The radial
velocities are much weaker than the $\ell=0$ case since the vortex reduces the
density in the dissipation region. In the lower panels, we show the dependence
of $\mu$, $D$ and $r_{h}$ on the dissipation strength for two different values
of $r_{0}$. Decreasing the size of the dissipation region weakens the draining
flow.

##  IV Eigenmodes

To find out whether the stationary profiles obtained above are stable, we
analyse their perturbations. The perturbed order parameter may be written,

$\Psi=e^{i\Phi}\left[\sqrt{n}+\sum_{m\geq
0}\left(ue^{im\theta}+v^{*}e^{-im\theta}\right)\right],$ (13)

where $u=u(r,t;m)$ and $v=v(r,t;m)$ are the radial eigenfunctions. Under this
decomposition, the linearised equations may be written as,

$\begin{split}&i\partial_{t}\begin{pmatrix}u\\\
v\end{pmatrix}=\begin{bmatrix}\mathcal{D}_{+}-i\widetilde{\Gamma}&n\\\
-n&-\mathcal{D}_{-}-i\widetilde{\Gamma}\end{bmatrix}\begin{pmatrix}u\\\
v\end{pmatrix},\\\
&\mathcal{D}_{\pm}=-\frac{1}{2}\nabla^{2}_{r}+\frac{(\ell\pm
m)^{2}}{2r^{2}}+\frac{v_{r}^{2}}{2}\mp iv_{r}\partial_{r}+2n+V-\mu,\\\
&\widetilde{\Gamma}=\Gamma+\bm{\nabla}\cdot\mathbf{v}/2.\end{split}$ (14)

This implies,

$\partial_{t}\mathcal{N}=-2\int d^{2}\mathbf{x}\,\Gamma\rho_{n},$ (15)

where $\mathcal{N}=\int d^{2}\mathbf{x}\,\rho_{n}$ is the norm and
$\rho_{n}=|u|^{2}-|v|^{2}$ is the norm density, i.e. the norm isn’t conserved
in the presence of a dissipative mechanism. Unlike the previous section where
we resupplied $N$ near $r_{B}$, we do not replenish the dissipated mode norm
since our reinjection strategy in Section VI will assume all lost particles
are resupplied to the axisymmetric background solution. Notice that the norm
density (and hence the norm) can be either positive or negative.

Figure 3: Left panels: Eigenvalue trajectories through the complex plane as
$\Gamma_{0}$ is varied. The trajectories start at the solid dots
($\Gamma_{0}=0$) and end on crosses ($\Gamma_{0}=0.5$). We show results for
$m=1,2$ and $r_{0}=2,5$. In the top panel, the inset shows in detail the
trajectories of the modes inside the grey rectangle. Note the nonuniform
spacing on the vertical axis. Right panels: Unstable modes for $r_{0}=2$ (top)
and $r_{0}=5$ bottom. The conventions for the line style and color correspond
to those on the left panels. For most modes in the spectrum (including the MQV
instability) the localised dissipation tends to make the mode decay. However,
certain modes attain a positive imaginary part when $\Gamma_{0}>0$,
representing a new class of unstable modes. Fewer secondary instabilities are
found for larger $r_{0}$.

In the stationary state, the operator on the right hand side of (14) is
$t$-independent and the equation is separable. Writing
$(u,v)^{\mathrm{T}}=a(t)(\tilde{u},\tilde{v})^{\mathrm{T}}$, one finds,

$\omega\begin{pmatrix}\tilde{u}\\\
\tilde{v}\end{pmatrix}=\begin{bmatrix}\mathcal{D}_{+}-i\widetilde{\Gamma}&n\\\
-n&-\mathcal{D}_{-}-i\widetilde{\Gamma}\end{bmatrix}\begin{pmatrix}\tilde{u}\\\
\tilde{v}\end{pmatrix},$ (16)

and $i\partial_{t}a=\omega a$, which is solved by $a(t)=a(0)e^{-i\omega t}$
where $\omega=\omega^{r}+i\omega^{i}\in\mathbb{C}$. The equation for the norm
in (15) implies the modulus of $a$ behaves like,

$|a(t)|=|a(0)|e^{\omega^{i}t},\qquad\omega^{i}=-\frac{1}{\mathcal{N}}\int
d^{2}\mathbf{x}\,\Gamma\rho_{n},$ (17)

where $\omega^{i}$ is a spatial averaging of the dissipation function weighted
by the norm density. For uniform dissipation, we find $\omega^{i}=-\Gamma$ and
all modes are damped irrespective of the sign of their norm. When $\Gamma$
varies in space, the sign of $\omega^{i}$ depends on the sign of the norm
density in the dissipation region. In that case, both growing (i.e. unstable)
and decaying modes can exist.

Note that the naive expectation that the dissipation function should damp
positive energy modes and enhance negative energy ones (where
$\mathcal{H}=\mathrm{Re}[\omega]\mathcal{N}$ is the mode energy) is not
satisfied here. Unlike damping found in stochastic formulations of the GPE
[36], the dissipative term in (1) decreases the particle number rather than
the energy, which is the statement of (6). The important thing for modes is
whether $\rho_{n}$ in the dissipation region has the opposite sign to
$\mathcal{N}$, as expressed by (17).

To solve (16) for the eigenvalues $\omega$ and radial eigenfunctions
$(\tilde{u},\tilde{v})^{\mathrm{T}}$, we diagonalise the matrix in square
parentheses using standard numerical algorithms (see Appendix A.2). Fig. 3
shows the eigenvalues $\omega$ obtained in this way for the $\ell=2$ vortex.
The left panels show the trajectories of the eigenvalues through the complex
plane as $\Gamma_{0}$ is varied, for $m=1$ (top panel) and $m=2$ (bottom
panel). We also show the effect of varying the size of dissipation region by
displaying results for $r_{0}=2$ (red lines) and $r_{0}=5$ (black lines). As
$\Gamma_{0}$ increases, most modes in the spectrum tend to migrate down into
the lower half plane where $\mathrm{Im}[\omega]<0$. However, certain modes can
acquire a small positive imaginary part when the dissipation is turned on,
with more unstable modes arising for smaller dissipation regions. This is
further exemplified by the panels on the right, which displays the variation
of the eigenvalues $\mathrm{Im}[\omega]>0$ with $\Gamma_{0}$. The mode with
$\mathrm{Im}[\omega]>0$ at $\Gamma_{0}=0$ is the MQV splitting instability,
which stabilises as the dissipation is increased. However, other modes in the
spectrum can be unstable for $\Gamma_{0}\neq 0$, with the $m=1$ modes
displaying the largest tendency for instability.

##  V WKB method

There are three key features of the spectra in Fig. 3 that require
explanation: (1) the quenching of the MQV splitting instability and (2) the
existence of secondary instabilities which are (3) fewer in number for larger
$r_{0}$. To explain this behaviour, we apply the WKB method of [28], the
application of which to our problem is detailed in Appendix B. The basis of
the approximation is to write fluctuations as a superposition of in- and out-
going plane waves with a radial wavenumber $p$ which varies with $r$ due to
the inhomogeneous flow field. Schematically these two modes of wave motion
oscillate as,

$e^{i\int(p^{\pm}+i\Gamma\partial_{\omega}p^{\pm})dr+im\theta-i\omega t},$
(18)

where $p^{\pm}$ are the two solutions to local dispersion relation, which
asymptotically looks like the standard Bogoliubov one, i.e.
$(\omega-\mathbf{v}\cdot\mathbf{k})^{2}=nk^{2}+k^{4}/4$ with $k=|\mathbf{k}|$.
Note that this gets modified inside the vortex core (see Appendix B for
details). Throughout the analysis, $\omega$ is treated as a real quantity and
the imaginary part obtained at the end of the procedure as a next-to-leading
order correction once the boundary conditions have been imposed.

The term in the exponent of (18) proportional to $\Gamma$ describes the decay
of the wave amplitude in the direction of propagation. The $e^{i\int
p^{\pm}dr}$ has different behaviour in three separate regions: two propagating
regions where $\Omega=\omega-\mathbf{v}\cdot\mathbf{k}$ is positive and
negative respectively and a zone separating the two where the wave tunnels.
These regions are depicted on Fig. 4. Since $\rho_{n}\propto\Omega$ in the WKB
approximation, a mode with a given $\omega,m$ has a positive (negative) norm
density when it is in the green (pink) area. Since the energy density is
proportional to $\mathrm{Re}[\omega]\rho_{n}$, a frequency that crosses from
the green to the pink region undergoes superradiant amplification (see [27,
28] for details on this point).

We can derive a resonance condition for modes which satisfy the boundary
conditions. For those $\omega$ which cross the tunnelling region in Fig. 4, it
is given by,

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

where we have defined,

$\begin{split}2S_{ab}=&\
\int^{r_{b}}_{r_{a}}\left(1+i\Gamma\partial_{\omega}\right)(p^{+}-p^{-})dr,\\\
I_{12}=&\ \int^{r_{2}}_{r_{1}}|\mathrm{Im}[p^{+}]|dr.\end{split}$ (20)

The derivation of this formula is given in Appendix B. The $S$ integrals can
be split into real and imaginary parts $S_{ab}=S^{r}_{ab}+iS^{i}_{ab}$, where
these two terms are respectively the phase integral from $r_{b}$ to $r_{a}$
then back again, and the total damping accumulated along the trajectory. When
the frequency needed to solve this condition has a small imaginary part, we
make the replacement $S^{r}_{ab}\to
S^{r}_{ab}+i\omega^{i}\partial_{\omega}S^{r}_{ab}$, whilst $S^{i}_{ab}$ is
evaluated for $\omega^{r}$ since $\Gamma_{0}$ is assumed small. The $I$
integral is evaluated in a region where $p^{\pm}$ are complex conjugate pairs,
and the effect of damping has been ignored in computing it since its
contribution to (19) is already exponentially small.

Figure 4: The WKB modes in (18) can exhibit three different behaviours:
propagating with positive $\rho_{n}$ (green region), propagating with negative
$\rho_{n}$ (red region) or tunnelling (grey region). These regions are
delineated by the turning points $r_{0,1,2}$. The curves shown are for the
parameters $\Gamma_{0}=0.05$, $r_{0}=5$, $\ell=2$ and $m=1,2$. Note that for
$m=\ell$, the inner turning point is $r_{0}=0$ by convention.

Before delving into the novel features of the spectrum when $\Gamma_{0}>0$, we
first we show how (19) encodes the physics of the basic MQV splitting
instability in the case $\Gamma_{0}=0$. This condition describes two types of
eigenmode in the system. The phonon frequencies (i.e. modes which live outside
the vortex) are exponentially close to the zeros of $\cot(S_{2B}+\pi/4)$.
Similarly, excitations of the vortex core have frequencies close to the zeros
of $\cot S_{01}$. However, when the zeros of the two cotangent functions are
close, the exponential on the right of (19) leads to a coupling between the
two types of eigenmode. Let $\cot S_{01}(\omega_{v})=0$ and
$\cot[S_{2B}(\omega_{p})+\pi/4]=0$. In that case, we expand (19) around the
zeros,

$(\omega-\omega_{p})(\omega-\omega_{v})+\epsilon|\mathcal{T_{\bar{\omega}}}|^{2}/T_{p}T_{v}\simeq
0,$ (21)

where $|\mathcal{T}_{\bar{\omega}}|\simeq e^{-I_{12}(\bar{\omega})}$ is the
modulus of the transmission coefficient across the tunnelling zone evaluated
at the central frequency $\bar{\omega}=(\omega_{p}+\omega_{v})/2$,
$\epsilon=\mathrm{sgn}(S^{\prime}_{01}S^{\prime}_{2B})$ with the prime
denoting an $\omega$ derivative, and $T_{p}=2|S^{\prime}_{2B}|$ and
$T_{v}=2|S^{\prime}_{01}|$ are the crossing times for phonons and vortex
modes. In our case, $\epsilon=+1$ for $\omega>0$ (which superradiate) and
$\epsilon=-1$ for $\omega<0$ (which do not superradiate). The eigenfrequencies
are then approximately given by,

$\omega\simeq\bar{\omega}\pm\sqrt{\Delta\omega^{2}-\epsilon\gamma^{2}}.$ (22)

Hence for non-superradiant states, there is an avoided crossing whereas for
superradiant states, we obtain a complex conjugate pair (one of which is
unstable) when the mode coupling
$\gamma=|\mathcal{T}_{\bar{\omega}}|/\sqrt{T_{p}T_{v}}$ is larger than the
frequency splitting $\Delta\omega=(\omega_{p}-\omega_{v})/2$. This coupling of
the negative energy vortex mode to a positive energy phonon is precisely the
reason for the dynamical splitting instability, as the MQV radiates energy
into sound when it splits (see [27, 28] for a comparison between the exact
spectrum and the predictions of (19)).

For the discussion of $\Gamma_{0}\neq 0$, let us rewrite (19) in the two
following equivalent forms,

$\displaystyle X+e^{2iS_{01}}+ie^{2iS_{2B}}\left(1+Xe^{2iS_{01}}\right)=$
$\displaystyle\ 0,$ (23) $\displaystyle
1+iXe^{2iS_{2B}}+e^{-2iS_{01}}\left(X+ie^{2iS_{2B}}\right)=$ $\displaystyle\
0,$ (24)

where we have defined,

$X=\frac{1+e^{-2I_{12}}/4}{1-e^{-2I_{12}}/4},$ (25)

which is the modulus of the reflection coefficient for superradiant modes and
the inverse of the same for non-superadiant modes.

Condition (23) is amenable for studying the MQV splitting instability when the
propagation time between the vortex and $r_{B}$ is sufficiently large that
$|e^{2iS_{2B}}|\ll 1$ (e.g. for low frequencies the relevant condition is
$r_{B}\gg 1/\omega^{i}$). Note that for this exponential to be small, we must
have an unstable mode with $\omega^{i}>0$, which only occurs for
$\omega^{r}>0$. If we assume $\omega^{r}\gg\omega^{i}$, the oscillation
frequency is determined to leading order by $\cos S^{r}_{01}(\omega^{r})=0$,
describing a mode in the vortex core, and the growth rate by,

$\omega^{i}=\frac{\log X-2|S^{i}_{01}|}{T_{v}}.$ (26)

The first term shows that these modes are unstable due to superradiance, since
$X>1$ and therefore the logarithm contributes positively. The second term
results from dissipation and acts to reduce the amplitude of the instability,
which is consistent with Fig. 3. Physically this makes sense since the
dissipation acts in the region where the mode lives, and from (15) we expect
dissipation to deplete the norm density (or by proxy the amplitude) rather
than the energy. Therefore, in this case, superradiance tries to make the mode
grow whereas dissipation tries to stop it growing. This explains the first
feature of Fig. 3 identified at the beginning of this section.

To understand the origin of the secondary instabilities, we study the second
condition (24). Consider $\omega^{r}>0$ modes in the limit where dissipation
is sufficiently strong that $|e^{-2iS_{01}}|\ll 1$. To satisfy this condition,
the vortex core should be large enough that the mode propagating inside it has
enough time to be damped significantly, which is achieved in practice for
large $\ell$. For a small growth rate, the leading order condition for the
oscillation frequency is $\cos(S^{r}_{2B}(\omega^{r})+\pi/4)=0$ describing
phonons trapped outside the vortex. The growth rate is given at leading order
by,

$\omega^{i}=\frac{\log X-2|S^{i}_{2B}|}{T_{p}},$ (27)

where the logarithmic term again contributes positively due to superradiance.
When dissipation is confined to the region where the mode propagates in the
vortex core, $S^{i}_{2B}=0$ and phonons become unstable, explaining the second
finding of Fig. 3. Physically, the negative energy part of the mode
responsible for superradiant amplification gets damped out inside the
dissipation region and (to leading order) the part of the wave which reflects
back can be neglected. This is the same mechanism as the black hole bomb
instability, where superradiant modes outside the black hole are confined by
an external mirror and the horizon provides a perfectly absorbing inner
boundary. When there is dissipation in the green region of Fig. 4,
$S^{i}_{2B}$ becomes non-zero and acts to weaken the instability, completely
quenching it for large enough values. This explains the third finding of Fig.
3 that secondary instabilities are suppressed for large dissipation regions.
Note that although we have considered here two specific limits (first large
$r_{B}$ then large $\ell$) we expect the same qualitative features to emerge
in the spectrum for general $r_{B}$ and $\ell$.

Figure 5: Intensity plots of $n(x,y)$ from three different simulations, which
differ by the dissipation strength and the initial vortex locations. We used
$r_{0}=5$ and $r_{B}=25$ and the $x$ and $y$ limits for each panels are
$[-r_{B},r_{B}]$. Panels (a1-4) show a low $\Gamma_{0}$ simulation, and the
insets ($r\in[0,r_{B}]$ on the horizontal and a.u. on the vertical axis) show
that the $m=2$ waveforms during vortex recombination and splitting are in good
agreement with the corresponding eigenmodes in the spectrum. Panels (b1-2)
show that for higher dissipation, the vortices do not separate at late times.
Panels (c1-4) show that for an asymmetric initial condition, the inspiral can
trigger a secondary instability in the $m=1$ mode. The inset compares the
measured waveform with corresponding $m=1$ eigenmode. Panels (d1-3) show the
shape of the unstable mode near the boundary in the angular direction. At late
times, the wave becomes sharply peaked, triggering a short wavelength
instability at the leading edge. Panels (b3-4) show the total energy and the
most highly occupied $m$-components near the centre of the trap as functions
of $t$.

##  VI Numerical simulations

In this section, we explore the predictions of Fig. 3 by simulating the full
nonlinear equation in (1), demonstrating the effect of the convergent flow on
a pair of initially separate vortices. We take as our initial condition a
stationary flow with $\Gamma_{0}\neq 0$ but no vortex. We then imprint a pair
of SQVs (both rotating counter-clockwise) near the edge of the trapping
potential. Particles lost due to dissipation are replenished at the edge of
the trap such that $N$ is conserved. Details of the procedure are described in
Appendix A.3. Fig. 5 shows the results of three such simulations.

Figure 6: Panels (a1-4) depict a continuation of panel (c4) in Fig. 5 with the
colour axis increased to clearly show the overdense region. Panels (b1-4) show
a vector map of the velocity field at the same times, with regions where
$v_{r}>0.06$ ($v_{\theta}>0.45$) highlighted in magenta (yellow). On panels
(c1-4), we see a shock wave being generated at the leading edge of the $m=1$
wave eventually propagating all the way round the trap edge. Once the shock
wave reaches large enough amplitudes, corotating vortices nucleate near the
boundary (panel (a3)) and move into the centre to form a giant vortex mass
(panel (a4)).

In panels (a1-4) of Fig. 5, the amount of dissipation is small and we start
with the vortices located symmetrically about the centre at
$(x,y)=(12.7,12.7)$ and the antipodal point. The vortices quickly spiral into
the middle and, once they are separated by around three healing lengths, the
motion is well described by the conjugate mode to the MQV splitting
instability, which describes vortex recombination. Between panels (a2) and
(a3), we Fourier transform $n_{m=2}$ over the interval $t=[317.2,410]$ and
plot the $\omega$ component with the largest amplitude (black curve), which
agrees well with the recombination mode predicted from (16) (dashed red line).
The frequency is computed by performing a quadratic interpolation around the
peak value of the Fourier transform at $r=1.4$ (where the amplitude of the
spatial waveform exhibits a maximum). The growth rate is found by performing a
linear regression on $\log|n_{m=2}(t,r=1.4)|$ over the same time interval. We
obtain a measured eigenfrequency $0.40(3)-0.0206(4)i$ 111For the measured
eigenfrequencies, the error in the real part is half the frequency resolution,
whilst the error in the imaginary part is quoted from the $95\%$ confidence
bounds. The predicted values are quoted to 3 s.f., which we compare against
$0.439-0.0293i$ from Fig. 3. The agreement is acceptable given that at this
stage, the system is in a highly non-equilibrium situation. Indeed, since $N$
is fixed, the value of $\mu$ will shift to accommodate the non-zero occupation
of various excited states, leading to differences with the spectrum of Fig. 2
which was computed in equilibrium. After a long time, the vortices separate by
a small amount since the dissipation is not sufficient to quench the MQV
splitting instability. The measured eigenfrequency of splitting mode (obtained
using the same method as above over the interval $t\in[1000,1400]$) is
$0.44(1)+0.00150(3)i$ whereas the predicted value is $0.440+0.00365i$. The
poor agreement in the imaginary part can be attributed to the effects of $N$
conservation and nonlinearities for high mode amplitudes. We checked this by
performing a simulation starting from equilibrium with a small initial
amplitude in the splitting mode, finding that, once the mode grew large
enough, the growth rate dropped to around the measured value above. The black
curves on panels (b3-4) are the total energy $H$, which is a slowly increasing
function at late times, and the $m=2$ mode amplitude which, after an initial
transient, decreases then increases exponentially as expected from the
discussion above.

In panels (b1-2) of Fig. 5, we show results of a simulation prepared with the
vortices at the same locations but a larger value of $\Gamma_{0}$. The two
SQVs combine much quicker in this case and residual fluctuations dissipate
away, in keeping with the fact that the spectrum contains no unstable modes
for $m=1,2$. The energy (blue curve on panel (b3)) decreases very slowly at
late times whilst $|n_{m=2}|$ decreases exponentially, exhibiting a beating
effect due to the superposition of two decaying co- and counter-rotating
modes. These results support the idea that a dense quantum vortex cluster can
be prevented for splitting by a sufficiently strong dissipation mechanism.

In the third row of Fig. 5, we pick a configuration where the strongest
unstable mode in the spectrum is one of the secondary instabilities with
$m=1$. Two SQVs are placed asymmetrically around the centre at
$(x,y)=(12.7,12.7)$ and $(11.5,-10.2)$ to promote the excitation of this mode
when they spiral into the middle and combine. Between the second and third
panels, we display the $m=1$ component of the density (which is the maximum of
the time Fourier transform over the interval $t\in[500,800]$) and find good
agreement with the predicted waveform. The measured eigenfrequency of this
mode is $0.083(11)+0.00068(1)i$ whilst the predicted value is
$0.0767+0.00236i$. The real part of the frequency agrees within the estimated
error. The discrepancy in the imaginary part is again attributed to the large
amplitude during the measurement window (see Appendix A.4 for a discussion of
this point). Measurement of $\mathrm{Im}[\omega]$ at earlier times where the
amplitude is smaller is complicated by the presence of a decaying counter-
rotating mode, which causes the beating effect in the (red) $m=1$ curve on
panel (b4). Panel (b3) shows that there is a drastic increase in the total
energy, which is due to the damping of the negative energy part of the wave in
the dissipation region. Panels (d1-3) illustrate the shape of the growing
$m=1$ mode in the $\theta$ direction near the edge of the trap. At late times,
the wavefront becomes sharply peaked, which can be attributed to the nonlinear
excitation of the $m=2$ wave (also shown on panel (b4)) which has roughly
twice the frequency of the $m=1$ component. In the last panel, a short
wavelength instability (a shock wave) appears at the leading edge of the
overdense region and propagates around the trap wall. The observed self-
steepening dynamics followed by the appearance of a shock wave is similar to
that studied in the 1D set-up of [38].

Further development of this state is depicted in detail in Fig. 6. We show the
density (panels (a)), the velocity field (panels (b)) and the angular profile
of the density near $r_{B}$ (panels (c)) at four different instants of time.
Eventually, the angular shock wave becomes large enough that additional
vortices nucleate at the trap wall, before being dragged into the centre by
the draining flow, which is most apparent from panels (a). The vector maps in
the (b) panels represent the velocity field (with the arrows scaled with
respect to the largest one in the plot). We highlight regions where
$v_{r}>0.06$ ($v_{\theta}>0.45$) in magenta (yellow). The overlapping regions
appear as red. We find that the radial velocity is directed outward in front
of the $m=1$ wave. This supplies additional particles to the leading edge of
the overdense region resulting in a sharp peak that continues to grow. The
angular velocity of the condensate is largest in the overdense region and
decreases suddenly at the leading edge, in a manner consistent with known
shock wave formation mechanisms [38]. By the end of the simulation, the mass
of vortices in the centre was precessing in the counter-clockwise direction
whilst more SQVs continued to nucleate from the trap edge.

For the presented data in Fig. 5 panels (c,d) and Fig. 6, the trap height was
increased to $V_{0}=10^{3}$ to ensure the shock wave near the trap edge was
not a consequence of the mirror image condensates arising from the periodic
boundary conditions. Qualitatively similar results were observed to those
performed with $V_{0}=5$. As a further check that the shock wave was not a
consequence of the numerical resolution, we checked the effect of (firstly)
doubling the space outside the $2r_{B}\times 2r_{B}$ condensate zone whilst
keeping the same number of grid points, and (secondly) changing the number of
grid points to $N_{x}=256$ (i.e. half the spatial resolution) and the time
step to $\Delta t=10^{-2}$. The shock wave was found to form in both cases,
with minor difference expected due to the change in the precise seeding
conditions at different resolutions.

##  VII Discussion

Multiply quantised vortices (MQVs) are generally expected to decay into
clusters of singly quantised vortices (SQVs). A natural mechanism purported to
enhance the stability of MQVs involves a dissipation mechanism which induces a
convergent (draining) flow in the system [22, 19, 20, 17]. We have
demonstrated that, in circularly symmetric trapping geometries (realisable
experimentally [39, 40]) stability is not guaranteed, challenging the
perspective that a convergent fluid flow has a stabilising influence on MQVs.

We considered a particle dissipation term used to model losses in Bose gases
illuminated with a laser/electron beam [23, 30, 22], polariton condensates
[18, 19, 21] and nonlinear optical condensates (the so called photon
superfluids) [41, 31]. The model also captures phenomenological aspects of
superfluid 4He flows in the thin film regime, with the density acting as a
proxy for the height field and the dissipation representing loss of fluid from
the plane, e.g. due to a physical outlet.

The observed tendencies of our system are a consequence of the dissipation
reducing the particle number rather than the total energy which, when combined
with an amplification mechanism, can result in secondary instabilities. Our
WKB method revealed that the relevant factor is where the confined mode lives
relative to the localised dissipation region. Since the MQV splitting
instability is confined to the vortex core, where dissipation acts, it can be
quenched. When superradiant phonon modes are confined outside the dissipation
region, the negative energy transmitted into the vortex core is absorbed
whilst the mode occupation outside the vortex increases. However, if the
dissipation region extends outside the vortex core, the occupation of phonon
modes can similarly be reduced. We then demonstrated that secondary
instabilities can have drastic implications for the late time dynamics of the
systems, growing so large as to generate a shock wave at the edge of the
condensate which seeds the nucleation of many vortices near the trap wall.

Observation of these features is within current experimental capabilities of
existing platforms implementing atomic BECs. For example, in [42], a 2D planar
geometry containing vortices was realised in a 87Rb condensate with a lifetime
of $28\mathrm{s}$. In this system, the healing length is on the order of 1µm
and the natural timescale is $\tau=\xi/c\sim 1\mathrm{ms}$, where $c$ is the
sound speed. If dissipative losses are induced with an electron beam (e.g.
[22, 30]) our results suggest that the secondary instability would onset on a
timescale of $10^{2}-10^{3}\tau$, well within the condensate lifetime.

For polariton condensates, stability of an MQV with 15 circulation quanta was
found in [21]. In that case, dissipation was uniform through the system,
hence, superradiant phonons are dissipated as well as the MQV splitting mode.
The system size and phonon decay length were $\sim 70$ and $\sim 8$ healing
lengths respectively, explaining the absence of secondary instabilities. In
our dimensionless units, their dissipative parameter is $\Gamma\sim 0.1$, for
which the WKB treatment of Section V seems reasonable. Hence, our methods
could be applied to provide insight into properties of superradiance in these
systems.

Finally, since Gross-Pitaevskii based models are a useful phenomenological
tool for studying the microscopic dynamics of quantum liquids [1], our results
may help elucidate the behaviour of superfluid vortex experiments. In
superfluid 4He, the core of a macroscopic vortex cluster is either a
depression in the liquid’s surface over the drain hole (which provides the
effective dissipation) or a fully formed throat which plunges through it [14,
17]. Studies of the problem in 3D reveal the formation of a complex tangle of
vortex lines over the drain [15, 16]. A detailed understanding of the simpler
2D problem may help identify the key mechanisms at play, and may be of further
use in studying the effectively 2D dynamics of waves on the superfluid
interface.

Acknowledgements.— I would like to thank Maxime Jacquet and Ruth Gregory for
helpful discussions, and Perimeter Institute for Theoretical Physics for
hospitality whilst this research was conducted. This work was supported by the
Science and Technology Facilities Council through the UKRI Quantum
Technologies for Fundamental Physics Programme [Grant ST/T005858/1]. Research
at Perimeter Institute is supported by the Government of Canada through the
Department of Innovation, Science and Economic Development Canada and by the
Province of Ontario through the Ministry of Research, Innovation and Science

## References

  * [1] Barenghi C F, Donnelly R J and Vinen W F 2001 Quantized vortex dynamics and superfluid turbulence vol 571 (Springer Science & Business Media)
  * [2] Kosterlitz J M 2017 Rev. Mod. Phys. 89 040501
  * [3] Abrikosov A A 2004 Rev. Mod. Phys. 76 975
  * [4] Barenghi C F and Parker N G 2016 A primer on quantum fluids (Springer)
  * [5] Shin Y, Saba M, Vengalattore M, Pasquini T A, Sanner C, Leanhardt A E, Prentiss M, Pritchard D E and Ketterle W 2004 Phys. Rev. Lett. 93 160406
  * [6] Isoshima T, Okano M, Yasuda H, Kasa K, Huhtamäki J A, Kumakura M and Takahashi Y 2007 Phys. Rev. Lett. 99 200403
  * [7] Lundh E 2002 Phys. Rev. A 65 043604
  * [8] Simula T P, Virtanen S M M and Salomaa M M 2002 Phys. Rev. A 65 033614
  * [9] Giacomelli L and Carusotto I 2020 Phys. Rev. Research 2 033139
  * [10] Patrick S, Gupta A, Gregory R and Barenghi C F 2023 Phys. Rev. Res. 5(3) 033201 URL https://link.aps.org/doi/10.1103/PhysRevResearch.5.033201
  * [11] Lautrup B 2011 Physics of continuous matter: exotic and everyday phenomena in the macroscopic world (CRC press)
  * [12] Yano H, Ohyama K, Obara K and Ishikawa O 2018 Observation of the spiral flow and vortex induced by a suction pump in superfluid 4he J. Phys. Conf. Ser. vol 969 (IOP Publishing) p 012002
  * [13] Matsumura I, Ohyama K, Sato K, Obara K, Yano H and Ishikawa O 2019 J. Low Temp. Phys. 196 204–210
  * [14] Obara K, Matsumura I, Tajima N, Ohyama K, Yano H and Ishikawa O 2021 Phys. Rev. Fluids 6 064802
  * [15] Inui S, Nakagawa T and Tsubota M 2020 Phys. Rev. B 102 224511
  * [16] Ruffenach W, Galantucci L and Barenghi C F 2023 J. Low Temp. Phys. 1–14
  * [17] Švančara P, Smaniotto P, Solidoro L, MacDonald J F, Patrick S, Gregory R, Barenghi C F and Weinfurtner S 2023 arXiv preprint arXiv:2308.10773
  * [18] Carusotto I and Ciuti C 2013 Rev. Mod. Phys. 85 299
  * [19] Alperin S N and Berloff N G 2021 Optica 8 301–307
  * [20] Solnyshkov D, Septembre I and Malpuech G 2023 arXiv preprint arXiv:2309.06269
  * [21] Delhom A, Guerrero K, Calizaya P, Falque K, Brady A J, Agullo I and Jacquet M J 2023 arXiv preprint arXiv:2310.16031
  * [22] Zezyulin D A and Konotop V V 2014 arXiv preprint arXiv:1411.7948
  * [23] Engels P, Coddington I, Haljan P C, Schweikhard V and Cornell E A 2003 Phys. Rev. Lett. 90 170405
  * [24] Press W H and Teukolsky S A 1972 Nature 238 211–212
  * [25] Brito R, Cardoso V and Pani P 2020 Superradiance vol 10 (Springer)
  * [26] Slatyer T R and Savage C M 2005 Class. Quant. Grav. 22 3833
  * [27] Patrick S, Geelmuyden A, Erne S, Barenghi C F and Weinfurtner S 2022 Phys. Rev. Research 4 043104
  * [28] Patrick S, Geelmuyden A, Erne S, Barenghi C F and Weinfurtner S 2022 Phys. Rev. Research 4 033117
  * [29] Cardoso T A S and Richartz M 2022 Phys. Rev. A 106 063310
  * [30] Barontini G, Labouvie R, Stubenrauch F, Vogler A, Guarrera V and Ott H 2013 Phys. Rev. Lett. 110 035302
  * [31] Braidotti M C, Prizia R, Maitland C, Marino F, Prain A, Starshynov I, Westerberg N, Wright E M and Faccio D 2022 Phys. Rev. Lett. 128 013901
  * [32] Volovik G E 2003 The universe in a helium droplet vol 117 (OUP Oxford)
  * [33] Andersen A, Bohr T, Stenum B, Rasmussen J J and Lautrup B 2003 Phys. Rev. Lett. 91 104502
  * [34] Stepanyants Y A and Yeoh G H 2008 J. Fluid Mech. 604 77–98
  * [35] Basak S and Majumdar P 2003 Class. Quant. Grav. 20 3907
  * [36] Cockburn S P and Proukakis N P 2009 Laser Physics 19 558–570
  * [37] For the measured eigenfrequencies, the error in the real part is half the frequency resolution, whilst the error in the imaginary part is quoted from the $95\%$ confidence bounds. The predicted values are quoted to 3 s.f.
  * [38] Damski B 2004 Phys. Rev. A 69 043610
  * [39] Johnstone S P, Groszek A J, Starkey P T, Billington C J, Simula T P and Helmerson K 2019 Science 364 1267–1271
  * [40] Navon N, Smith R P and Hadzibabic Z 2021 Nat. Phys. 17 1334–1341
  * [41] Zezyulin D A, Konotop V V, Barontini G and Ott H 2012 Phys. Rev. Lett. 109 020405
  * [42] Gauthier G, Reeves M T, Yu X, Bradley A S, Baker M A, Bell T A, R D H, Davis M J and Neely T W 2019 Science 364 1264–1267

##  Appendix A Numerical methods

###  A.1 Stationary profiles

First, we describe our numerical technique to find the stationary draining
vortex profiles. We exploit the symmetry of the desired solution by factoring
out the vortex winding $e^{i\ell\theta}$ so that the equation to solve is a
function of $r$ only. Using the units of (9), we write (1) in the form,

$i\partial_{t}z=\left(-\frac{\partial_{r}^{2}}{2}-\frac{\partial_{r}}{2r}+\frac{\ell^{2}}{2r^{2}}+V-i\Gamma+|z|^{2}\right)z,$
(A.1)

where $z(r)=\sqrt{n(r)}e^{i\Theta(r)-i\mu t}$. The radial grid is
$r_{j}\in[\Delta r,r_{B}+d]$, and $j=1...N_{r}$ with $N_{r}=256$ and $d=5$.
For the radial derivatives, we use 5-point finite different stencils accurate
to $\mathcal{O}(\Delta r^{4})$, where $\Delta r$ is the radial grid step size.
The behaviour in (12) implies that correct boundary conditions at $r=0$ are
Neumann for $\ell=0$ and Dirichlet for $\ell\neq 0$. The boundary condition at
the edge of the grid is arbitrary since it is applied deep in the region where
$z$ goes to zero under the potential barrier. We choose it to be a Dirichlet
boundary condition.

Due to the dissipation in the centre, any excitations of the ground state
decay away. An implicit assumption here is that there are no unstable modes
with $m=0$, which turns out to be true (at least for the scenarios we are
interested in). Therefore, starting from an arbitrary initial state, we should
eventually arrive at a $z$ corresponding to the ground state provided we
evolve (A.1) for long enough. The evolution is enacted using a fourth order
Runge-Kutta algorithm. $\Delta t=2\times 10^{-3}$ is taken for the time step,
satisfying $\Delta t<\Delta r^{2}/2$ which ensured stability. This meant that
going to higher spatial resolution was challenging since decreasing the
spatial grid step size by a factor of 2 would increase the simulation time for
the stationary states by a factor of 8, whilst the 2D simulations described in
Appendix A.3 would take 16 times longer. After each time step, we compute the
number of particles $\Delta N$ lost due to dissipation. These particles are
then resupplied at the edge of the trap before the next time step according to
the following prescription,

$\begin{split}&z\sqrt{1+\Delta N\beta(r)/|z|^{2}}\to z,\\\
&\beta=\beta_{0}\begin{cases}r-r_{1},\quad\,r_{1}<r\leq r_{2}\\\ 1,\
\qquad\quad r_{2}<r\leq r_{3}\\\ r_{B}-r,\quad r_{3}<r\leq r_{B}\\\
\end{cases}\end{split}$ (A.2)

where $r_{1}=r_{B}-4$, $r_{2}=r_{B}-3$ and $r_{3}=r_{B}-1$ and $\beta_{0}$ is
a normalisation such that $2\pi\int dr\,r\beta=1$. This ensures that the total
number of particles is conserved throughout the simulation.

We estimate the chemical potential at each time step $t_{k}$ using
$\mu_{k}=-\mathrm{Im}[\log(z_{k}/z_{k-1})]/\Delta t$ which gives a function
over the $r$ grid. The evolution is continued until the maximum value of
$\mu_{k}-\mu_{k-1}$ in the region $r<r_{B}$ is below $5\times 10^{-6}$. At
this point there will still be residual fluctuations in $z$. Since these are
decaying, we evolve for a further 2000 ($10^{4}$) time steps for
$\Gamma_{0}>0.05$ ($\Gamma_{0}\leq 0.05$) and take the value of $|z|$ as its
average during this time. The phase is measured at each point $r_{j}$ and fit
with the function $\Theta_{j}-\mu_{j}t$. The first term in this fit gives the
phase function in (7) from which obtain $v_{r}$ by enacting a $\partial_{r}$
stencil incorporating the Neumann boundary condition at the origin. The second
term gives the chemical potential as a function of $r$. We observe that this
is constant to a high precision and take its average over the region
$r<r_{B}$.

To generate the solutions in Figs. 1 and 2, we start by calculating the non-
draining solution with $\ell\neq 0$ using an imaginary time prescription (see
e.g. Appendix A of [27]). This is the solution for $\Gamma_{0}=0$. We then use
this as our initial condition for the largest value of $\Gamma_{0}$ simulated
(which was $0.5$). This was done since the strong dissipation here will
quickly damp the solution to the desired equilibrium state. We then generate
the solutions for all $\Gamma_{0}$ in descending order, each time taking the
result from the previous computation as the initial condition.

###  A.2 Eigenmodes

To solve (16) for the eigenmodes, we use the same radial grid $r_{j}$ and
finite difference stencils defined in the previous section. We choose a
Dirichlet boundary condition for the edge of $r$ grid, although the results
are insensitive to this choice in the frequency range probed. The boundary
condition at $r=0$ is more delicate. By expanding $(u,v)$ around $r=0$ and
inserting into (16), we deduce that the leading order behaviour is $u\sim
r^{|m+\ell|}$ and $v\sim r^{|m-\ell|}$. Hence, we apply a Dirichlet boundary
condition whenever the exponent is different from zero. Since the next term in
the series is multiplied by an extra factor of $r^{2}$, we apply a Neumann
boundary condition for a zero exponent. (16) is then solved using a standard
matrix diagonalisation routine (we apply Matlab’s `eig` function).

###  A.3 Nonlinear dynamics in 2D

For the simulations in Section VI of the full equation (1), we take as our
initial condition a stationary solution with $\Gamma_{0}\neq 0$ and no vortex,
say $\Psi_{0}$, and interpolate this onto a Cartesian $(x,y)$ grid. We take
$x_{j}\in[-(r_{B}+d),r_{B}+d]$ (and similarly for $y_{j}$) with $j=1...N_{x}$
and $N_{x}=2N_{r}=512$, so that the resolution is comparable to our radial
code. A pair of SQVs is then imprinted. For this process, we first calculate
the non-draining SQV density profile which asymptotes to unity, say
$n_{1}(r)$, then extend the solution to a high value of $r$ using the expected
asymptotic behaviour. This solution is recentred to the vortex coordinates
$(x_{\nu},y_{\nu})$ and the initial condition becomes,

$\begin{split}&\Psi_{\mathrm{init}}=\Psi_{0}(x,y)\\\
&\times\prod_{\nu}n_{1}(x-x_{\nu},y-y_{\nu})\frac{(x-x_{\nu})+i\ell_{\nu}(y-y_{\nu})}{|\mathbf{x}-\mathbf{x}_{\nu}|},\end{split}$
(A.3)

where the last factor ensures the defects have the correct phase winding. We
take $\ell_{\nu}=1$ with $\nu=1,2$ and $\mathbf{x}_{\nu}$ given in the main
text. Since the vortices merge in the centre during the simulation, we
multiply $\Psi_{\mathrm{init}}$ by a factor $(N/\int
d^{2}\mathbf{x}|\Psi_{\mathrm{init}}|^{2})^{1/2}$ so that the system has the
same total number of particles as our equilibrium solutions in Fig. 2 (which,
we recall, have the same $N$ as a non-draining MQV with unit chemical
potential by construction). The time evolution is then performed using a
split-step Fourier spectral method. That is, given the solution at time
$t_{k}$, we find $\Psi$ at $t_{k+1}$ according to,

$\begin{split}\widetilde{\Psi}_{1}=&\ e^{-i\Delta
t(V-i\Gamma+|\Psi_{k}|^{2})/2}\Psi_{k},\\\ \widetilde{\Psi}_{2}=&\
\mathcal{F}^{-1}\big{[}e^{-i\Delta
t\mathbf{k}^{2}/2}\mathcal{F}\big{[}\widetilde{\Psi}_{1}\big{]}\big{]},\\\
\Psi_{k+1}=&\ e^{-i\Delta
t(V-i\Gamma+|\widetilde{\Psi}_{2}|^{2})/2}\widetilde{\Psi}_{2},\end{split}$
(A.4)

where $\mathcal{F}$ and $\mathcal{F}^{-1}$ denote the 2D Fourier transform and
its inverse on the $(x,y)$ grid, with $\mathbf{k}$ coordinates in Fourier
space, and we set $\Delta t=4\times 10^{-3}$. At the end of each time step, we
resupply particles lost through dissipation at the edge of the trap with the
same prescription used to find our stationary solutions in the radial
simulations. In this way, the total particle number is conserved throughout
the evolution. We can think of this as representing a scenario where a fixed
number of particles remain in circulation.

Figure 7: Growth of the unstable $m=1$ mode (black line) with the same
parameters used for Fig. 5. The best fit to the absolute value over the window
$t\leq 500$ is also shown (red line). The measured eigenfrequency (displayed
in figure) is in excellent agreement with the predicted value of
$0.0767+0.00236i$.

###  A.4 Mode growth

In the simulations of Fig. 5, we found discrepancies between the measured
eigenfrequencies and the predicted ones, particularly in the imaginary part.
We attribute this to nonlinear effects at large mode amplitudes as well as our
resupply strategy, which keeps $N$ fixed at the expense of depleting the
background for large mode occupations. We illustrate an example of this on
Fig. 7 for the unstable $m=1$ mode. The same parameters were used as for the
(c) panels of Fig. 5 except, this time, the initial condition was the
stationary $l=2$ vortex state seeded with the (normalised) $m=1$ unstable mode
multiplied by an amplitude $0.5$. We measure the eigenfrequency for $t\leq
500$ and display the result on the figure. The growth rate and oscillation
frequency are in excellent agreement with the predicted values at early times
when the amplitude is small. At intermediate times, the growth rate starts to
deviate from the initial value, with the amplitude saturating at late times.
When the eigenfrequencies are measured in Fig. 5, the wave amplitudes are of
the order $0.1$ which is where the discrepancy starts in Fig. 7. This suggests
that deviations in the growth rate are consistent with a nonlinear mechanism.

##  Appendix B WKB method

In this appendix, we provide the details of the method leading to the
resonance condition (19) in the main text. First, it is helpful to define new
eigenfunctions $g,f$ according to,

$\sum_{m\geq 0}e^{im\theta}\begin{pmatrix}u\\\
v\end{pmatrix}=\begin{bmatrix}g(\mathbf{x},t)+if(\mathbf{x},t)/2\\\
g(\mathbf{x},t)-if(\mathbf{x},t)/2\end{bmatrix},$ (B.1)

which satisfy,

$\begin{split}(D_{t}-\widetilde{\Gamma})f+ng+\frac{1}{4}\mathcal{D}_{r}g=&\
0,\\\ (D_{t}-\widetilde{\Gamma})g-\mathcal{D}_{r}f=&\ 0,\\\
\mathcal{D}_{r}=-\nabla^{2}+\mathbf{v}^{2}+2(n+V-&\ \mu),\end{split}$ (B.2)

and $D_{t}=\partial_{t}+\mathbf{v}\cdot\bm{\nabla}$ is the material
derivative. In the WKB approximation, we write,

$\begin{bmatrix}f\\\
g\end{bmatrix}=\begin{bmatrix}\mathcal{A}(\mathbf{x},t)\\\
\mathcal{B}(\mathbf{x},t)\end{bmatrix}e^{i\mathcal{S}(\mathbf{x},t)}.$ (B.3)

The amplitudes are assumed to vary gradually compared to the phase. That is,
$|\partial\mathcal{A}|\ll|\mathcal{A}\partial\mathcal{S}|$ and
$|\partial^{2}\mathcal{S}|\ll(\partial\mathcal{S})^{2}$ (and similarly for
$\mathcal{B}$) for $\mathbf{x}$ and $t$ derivatives. Using the symmetry of the
system, we can write,

$\mathcal{S}=\int p(r)dr+m\theta-\omega
t\quad\begin{cases}\bm{\nabla}\mathcal{S}\equiv\mathbf{k}=(p,m/r),\\\
-\partial_{t}\mathcal{S}\equiv\omega\end{cases}$ (B.4)

where the decomposition into $\omega,m$ is exact since the background is
independent of $t,\theta$. Under these assumptions, (B.2) becomes,

$\begin{split}\Omega^{2}=&\ n\widetilde{k}^{2}+\widetilde{k}^{4}/4,\\\
\Omega=&\
\omega-\mathbf{v}\cdot\mathbf{k},\qquad\widetilde{k}^{2}=p^{2}+\widetilde{m}^{2}/r^{2},\\\
\widetilde{m}^{2}=&\ m^{2}+r^{2}\mathbf{v}^{2}+2r^{2}(n+V-\mu),\end{split}$
(B.5)

at leading order and,

$\begin{split}&\partial_{t}\left(\frac{\Omega\mathcal{A}^{2}}{F}\right)+\bm{\nabla}\cdot\left(\frac{\mathbf{v}_{g}\Omega\mathcal{A}^{2}}{F}\right)=-\frac{2\Gamma\Omega\mathcal{A}^{2}}{F},\\\
&F=n+\widetilde{k}^{2}/4,\end{split}$ (B.6)

at next to leading order. A similar equation for the second amplitude can be
derived using $\mathcal{B}=i\Omega\mathcal{A}/F$. In the derivation of (B.5)
and (B.6), we have assumed that dissipation is a weak effect that enters at
second order. Hence, the method applies only for low $\Gamma_{0}$.

The first line of (B.5) is the dispersion relation whose solutions
$\mathbf{k}=\mathbf{k}(\omega,m;r)$ represent the local value of the
wavevector for a given $\omega$ and $m$. The group velocity of the wave is
then determined by $\mathbf{v}_{g}=(\partial_{\omega}\mathbf{k})^{-1}$; in
particular, its radial component is
$\hat{\mathbf{e}}_{r}\cdot\mathbf{v}_{g}=(\partial_{\omega}p)^{-1}$. The
equation in (B.6) describes how amplitude is transported, with the term on the
right acting as a damping term. Evaluating the norm density for a WKB mode
gives $\rho_{n}=\Omega|\mathcal{A}|^{2}/F$ (where this expression is evaluated
on the real part of $\omega$), so that (B.6) describes how the norm is
transported across the system. In particular, the sign of the norm density
coincides with that of $\Omega$ and is transported in the direction of
$\mathbf{v}_{g}$. That is, a fluctuation with $\Omega<0$ has negative norm
density. The mode energy is defined as,

$\begin{split}\mathcal{H}=&\ \int
d^{2}\mathbf{x}\bigg{[}\frac{1}{2}|\partial_{r}u|^{2}+\frac{1}{2}|\partial_{r}v|^{2}+\frac{(\ell+m)^{2}}{2r^{2}}|u|^{2}\\\
&\
-\frac{iv_{r}}{2}\left(u^{*}\partial_{r}u-u\partial_{r}u^{*}-v^{*}\partial_{r}v+v\partial_{r}v^{*}\right),\\\
&\ +\frac{(\ell-m)^{2}}{2r^{2}}|v|^{2}+(2n+V-\mu)(|u|^{2}+|v|^{2})\\\ &\
+n(uv^{*}+u^{*}v)\bigg{]}.\end{split}$ (B.7)

which we evaluate in the stationary regime with the aid of (14) to obtain,

$\mathcal{H}=\omega\mathcal{N}+i\int d^{2}\mathbf{x}\,\Gamma\rho_{n}.$ (B.8)

From (17), it then follows that $\mathcal{H}=\mathrm{Re}[\omega]\mathcal{N}$.
That is, a $\mathrm{Re}[\omega]>0$ wave with $\Omega<0$ (and therefore
$\rho_{n}<0$) has negative energy density.

The dispersion relation (B.5) has four distinct solutions representing the
different allowed $\mathbf{k}$. When $\Gamma_{0}$ is small, only two of these
(which we call $\mathbf{k}^{\pm}$ such that
$\mathrm{Re}[p^{+}]>\mathrm{Re}[p^{-}]$) correspond to propagating solutions.
The other two are evanescent everywhere and enforcing the boundary conditions
sets their amplitudes to zero. Hence, the general solution can be written as,

$\begin{split}f=&\ \sum_{\omega,m}e^{im\theta-i\omega
t}\sum_{j=\pm}\frac{R^{j}(r)}{\sqrt{r}},\\\ R^{j}=&\
\alpha^{j}\sqrt{\frac{F^{j}\partial_{\omega}p^{j}}{\Omega^{j}}}e^{i\int(p^{j}+i\Gamma\partial_{\omega}p^{j})dr},\end{split}$
(B.9)

where $\alpha^{j}$ is an adiabatically conserved integration constant that
comes from solving (B.6), and superscript $j$ implies that the quantity is
evaluated for one of the $p^{\pm}$ solutions to the dispersion relation. The
integrals in the exponent are evaluated relative to a reference point where
the phase is known. Near $r_{B}$, $p^{+}$ corresponds to the out-going mode
with $\partial_{\omega}p^{+}>0$ whilst $p^{-}$ is in-going with
$\partial_{\omega}p^{-}<0$. Notice that the damping term in the exponent
implies that the amplitude decays in the direction of the group velocity
$(\partial_{\omega}p^{j})^{-1}$ as one would expect.

The locations where the $p^{j}$ change from real (propagating) solutions to
complex (evanescent) solutions are the turning points $r_{tp}$. At these
points, $p^{+}_{tp}=p^{-}_{tp}$ and the group velocity in the radial direction
vanishes, $(\partial_{\omega}p^{\pm}_{tp})^{-1}=0$. Intuitively, the wave
propagating in one direction comes to rest instantaneously before propagating
back in the other direction. Hence, the turning points are the locations where
waves scatter off the inhomogeneous background. On the other side of $r_{tp}$,
$p^{j}\in\mathbb{C}$ and the waves tunnel.

Looking at (B.9), we can see that the amplitude diverges as a consequence of
the radial group velocity vanishing at $r_{tp}$, signalling the breakdown of
the WKB approximation. This is not a problem but rather a necessity. The WKB
approximation encodes only the adiabatic variation of a given solution of the
dispersion relation, i.e. the various solutions do not exchange energy with
each other. To capture scattering (the transfer of energy between waves which
travel in opposite directions) the WKB approximation must breakdown at certain
locations.

There is a standard procedure to fix the breakdown of (B.9) which we now
describe. The full equations of motion are expanded in the vicinity of
$r_{tp}$ and an exact solution for the scattering waves is obtained. The
asymptotic form of this solution is matched onto the WKB modes either side of
$r_{tp}$ to provide a relation between them. When waves propagate for
$r<r_{tp}$ and are evanescent for $r>r_{tp}$, the relation is

$\begin{pmatrix}R^{+}_{tp}\\\
R^{-}_{tp}\end{pmatrix}=e^{\frac{i\pi}{4}}\begin{pmatrix}1&-\frac{i}{2}\\\
-i&\frac{1}{2}\end{pmatrix}\begin{pmatrix}R^{\downarrow}_{tp}\\\
R^{\uparrow}_{tp}\end{pmatrix},$ (B.10)

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

$\begin{pmatrix}R^{\uparrow}_{tp}\\\
R^{\downarrow}_{tp}\end{pmatrix}=e^{\frac{i\pi}{4}}\begin{pmatrix}\frac{1}{2}&-\frac{i}{2}\\\
-i&1\end{pmatrix}\begin{pmatrix}R^{+}_{tp}\\\ R^{-}_{tp}\end{pmatrix}.$ (B.11)

Here, a superscript $\uparrow$ ($\downarrow$) denotes a wave which grows
(decays) in the direction of increasing $r$.

In Fig. 4 of the main text, we depict the dependence of the various turning
points $r_{0,1,2}$ on $\omega$ for $m=1,2$. Between turning points, there is a
region where the waves are evanescent and tunnel across the flow. This
tunnelling zone separates the region where $\Omega>0$ ($\rho_{n}>0$) from that
with $\Omega<0$ ($\rho_{n}<0$). We then see that $\mathrm{Re}[\omega]>0$ modes
can have negative energy densities ($\mathrm{Re}[\omega]\rho_{n}<0$) inside
the vortex core, provided they are below a certain threshold (for $m=\ell$
this threshold is $\omega=\mu$). Hence, when such a wave impinges on the
vortex, it excites a state with negative energy in the core and (by energy
conservation) gets amplified upon reflection. This is the phenomenon of
rotational superradiant scattering familiar from black hole physics [25].

We can obtain a formula that describes the scattering of modes that cross this
tunnelling zone in the following way. Firstly, using (B.9), we define a number
that shifts one of the $R^{j}$ modes at $r_{b}$ to a smaller radius $r_{a}$,
such that the WKB amplitude does not diverge for $r_{a}<r<r_{b}$. That is,

$\begin{split}\mathcal{F}^{j}_{ab}=&\
\sqrt{\frac{\Omega^{j}_{b}F^{j}_{a}\partial_{\omega}p^{j}_{a}}{\Omega^{j}_{a}F^{j}_{b}\partial_{\omega}p^{j}_{b}}}\exp\left[-i\int^{r_{b}}_{r_{a}}(1+i\Gamma\partial_{\omega})p^{j}dr\right],\\\
R^{j}_{a}=&\ \mathcal{F}^{j}_{ab}R^{j}_{b},\end{split}$ (B.12)

and a single subscript indicates the location where a quantity is evaluated
e.g. $R^{j}_{a}=R^{j}(r_{a})$. We then combine (B.10), (B.11) and (B.12) to
define a relation between propagating amplitudes either side of the tunnelling
zone $r_{1}<r<r_{2}$,

$\begin{split}\begin{pmatrix}R^{+}_{1}\\\ R^{-}_{1}\end{pmatrix}=&\
\mathcal{N}_{12}\begin{pmatrix}R^{+}_{2}\\\ R^{-}_{2}\end{pmatrix},\\\
\mathcal{N}_{12}=&\
\mathcal{F}^{\downarrow}_{12}\begin{bmatrix}1+f_{12}^{2}/4&i(1-f_{12}^{2}/4)\\\
-i(1-f_{12}^{2}/4)&1+f_{12}^{2}/4\end{bmatrix},\end{split}$ (B.13)

where,

$f_{12}=\exp\left\\{-\int^{r_{2}}_{r_{1}}\mathrm{Im}[(1+i\Gamma\partial_{\omega})p^{\downarrow}]dr\right\\}.$
(B.14)

For a wide tunnelling zone, $f_{12}$ is assumed small and can be evaluated to
leading order for $\Gamma=0$. These relations can now be used to find a
resonance condition for modes that cross the tunnelling zone, by transporting
the solution from $r_{0}$ in the vortex core (see Fig. 4) to the outer
boundary at $r_{B}$,

$\begin{split}R^{+}_{0}=&\
\mathcal{F}^{+}_{01}\mathcal{F}^{\downarrow}_{12}\left(\sigma_{+}\mathcal{F}^{+}_{2B}R^{+}_{B}+i\sigma_{-}\mathcal{F}^{-}_{2B}R^{-}_{B}\right),\\\
R^{-}_{0}=&\
\mathcal{F}^{-}_{01}\mathcal{F}^{\downarrow}_{12}\left(-i\sigma_{-}\mathcal{F}^{+}_{2B}R^{+}_{B}+\sigma_{+}\mathcal{F}^{-}_{2B}R^{-}_{B}\right),\\\
\sigma_{\pm}=&\ 1\pm f_{12}^{2}/4.\end{split}$ (B.15)

The final step is to impose the boundary conditions. In the vortex core,
regularity of the solution requires us to have $R^{-}_{0}=iR^{+}_{0}$. At the
outer boundary, we can assume that the density goes to a constant value and
that the velocities are vanishingly small provided we impose a Neumann
boundary condition there; that is $R^{+}_{B}=R^{-}_{B}$. Both of this steps
are explained in detail in [28]. Inserting these relations into (B.15) and
rearranging yields the resonance condition (19) in the main text. Note, the
assumption of vanishingly small $v_{r}(r_{B})$ means that square root factors
in (B.12) cancel on the boundary when computing
$\mathcal{F}^{+}_{2B}/\mathcal{F}^{-}_{2B}$.