Superradiant phononic emission from the analog spin ergoregion in a two-component Bose-Einstein condensate

We make use of an analog gravity perspective to obtain a physical understanding of hydrodynamic instabilities stemming from the presence of quantized vortices in two-component atomic condensates and of their relation to ergoregion instabilities of rotating massive objects in gravitation. In addition to the localized instabilities related to vortex splitting, configurations displaying dynamically unstable modes that extend well outside the vortex core are found. In this case, the superradiant scattering process involves phonon emission into the much wider ergoregion of spin modes, so the physics most closely resembles the one of rotating massive objects. Our results confirm the potential of two-component condensates as analog models of rotating space-times in different regimes of gravitational interest.


Introduction
Analog gravity [1] is based on the use of table-top set-ups to simulate the physics of quantum fields on curved spacetimes. This approach is based on the formal equivalence of the equations of motion governing the dynamics of the two systems under suitable conditions. Since the original proposal by Unruh [2], curved space-time geometries have been realized in a large variety of systems, including surface waves on classical fluids [3], ultra-cold atoms [4], polariton fluids [5] and optical systems [6,7]. The phenomena considered in this context typically involve strong gravity environments that escape direct astrophysical investigations, as in the vicinity of a black hole or during the fast expansion of the early universe. Most remarkable examples are Hawking emission from black hole horizons [8][9][10][11][12][13], cosmological particle creation and dynamical Casimir effect [14][15][16][17][18][19], Penrose rotational superradiance [20][21][22][23].
In this context, the term rotational superradiance refers to the amplified scattering of waves from a rotating object at the expenses of the rotational energy of the object itself [24]. It is a fully classical phenomenon and is predicted to occur in a variety of spacetime geometries, in particular in the surroundings of massive rotating objects. In this case, dragging of the spacetime by the rotating object gives rise to a so-called ergoregion, namely a spatial region supporting negativeenergy light and matter modes. The process of superradiance can be microscopically explained as the amplified reflection of a positive-energy wave impinging onto the boundary of the ergoregion, accompanied by the transmission of a negative-energy wave inside the ergoregion. In the presence of some reflecting element, self-amplification of one of the outgoing modes may take place, leading to an exponentially growing perturbation. Depending on the location of the reflecting element on either the inner or the outer side of the ergosurface, the instability takes the name of ergoregion instability [25,26] or black hole-bomb instability [27].
Rotating Bose-Einstein condensates (BEC) are promising platforms to study superradiant phenomena [24,28]. A number of theoretical studies have investigated superradiance in configurations featuring a single vortex located at the center of a cylindrically symmetric trap as well as in planar geometries involving synthetic magnetic fields [29]. In particular, the analog gravity perspective has offered transparent physical explanations for the hydrodynamic instability of multiply-charged single vortex configurations [30][31][32][33] by connecting them to ergoregion instabilities of rotating massive objects in gravitation [34].
In the last decade, several authors have started investigating two-component BECs as a promising new platform for analog gravity [35][36][37][38]. Such systems are characterized by the existence of two independent branches of collective excitations, associated to perturbations of the total density and of the spin density, featuring very different values of the speed of sound and of the healing length [39]. As a consequence of the much faster speed of density-sound, configurations displaying horizons for spin-sound waves can be engineered while keeping the superflow sub-sonic with respect to density-sound and therefore robust against perturbations coupling to the density. Furthermore, the direct experimental observation of the system dynamics is facilitated by the significantly larger size of the spin-healing length and by the availability of interferometric techniques based on a coherent mixing of the two components which allow to image all quadratures of the quantum field. On the long run, all this suggests the possibility of investigating quantum entanglement features of the superradiant emission [40].
Within this general framework, the goal of this paper is to explore ergoregion instability phenomena in two-component BECs displaying a single quantized vortex in both components. Thanks to the very different value of the density-and spin-sound speeds, the ergosurface for spin waves is no longer bound to sit inside or in the close vicinity of the vortex core as it instead happens in single component BECs, but can be pushed far away from the vortex core into the external region where the density is approximately constant. This is of great interest for analog gravity as it leaves a wider space to investigate the dynamics of the quantum field in the ergoregion and disentangle the different effects at play.
Beyond analog gravity, an active interest for this physics is also coming from a purely quantum gas perspective: a recent numerical work [41] has pointed out the rich phenomenology of vortex splitting in harmonically trapped two-component BECs: singly-charged vortices are dynamically unstable for sufficiently strong repulsive inter-component interactions, whereas doubly-charged vortices dispose of several decay channels. Experimentally, the splitting of a singly-charged vortex into a pair of half-quantized vortices has been observed in spinor exciton-polariton superfluids [42] and antiferromagnetic spinor BECs [43].
Here, we take advantage of the analog gravity point of view to develop a transparent physical interpretation of vortex physics in two-component BECs in terms of superradiance effects. Use of this formalism allows us to highlight a number of features that are specific to two-components BECs and characterize them in view of forthcoming experimental studies. In addition to recovering the main predictions of previous works, e.g. on vortex splitting instabilities, we identify regimes where a clean analog of the ergoregion instability of massive rotating objects is visible.
The structure of the paper is the following: in Secs.2 and 3 we give a brief review of the theoretical tools that we are going to use in the following of the paper. In Sec. 4 we tackle the vortex stability problem within the local density approximation (LDA) that provides a clear qualitative understanding of the physical mechanisms. A more quantitative analysis using a Bogoliubov characterization of the collective modes and a Gross-Pitaevskii (GP) simulation of the timedependence is then provided in Secs.5 and 6, respectively. A summary of our conclusions and of further extensions of our work is finally given in Sec.7. Additional information on the numerical techniques adopted and on the independence of the results on the shape of the confinement potential are given in the Apppendices.

Two-component condensates and Bogoliubov theory of density and spin collective excitations
In this work, we are interested in the dynamics of a two-component condensate subject to an external potential V (r ) and to a resonant coherent coupling of strength Ω ≥ 0 between the two components. Let us focus on a spin-symmetric configuration of equal masses m 1 = m 2 ≡ m and intraspecies interactions g 1 = g 2 ≡ g > 0. Moreover, we assume to deal with a miscible mixture, |g 12 | < g , where g 12 defines the interspecies interaction strength.
The dynamics of such a system is ruled by two coupled GP equations [44]: where K = −ħ 2 ∇ 2 /2m is the kinetic energy operator, ψ j is the order parameter of the j -th component, and n j = |ψ j | 2 is the associated atomic density, normalized to the number of particles in the j -th state, N j . As long as N 1 = N 2 , the ground state of the system is everywhere unpolarized, with identical order parameters ψ o 1 (r ) = ψ o 2 (r ) and density profiles n 1 (r ) = n 2 (r ) for the two components of the mixture.
In addition to the Z 2 spin-symmetry related to the exchange of the two components, for a vanishing coherent coupling Ω = 0 the system has two additional U (1) symmetries associated to the independent conservation of the number of particles N 1,2 in the two states. The presence of the coherent coupling allows particles to be transferred from one component to the other and therefore breaks one of these two U (1) symmetries. In this case, we are left with a single U (1) symmetry related to the conservation of the total atom number N = N 1 + N 2 only.
A deep insight on the physical properties of the system is provided by the Bogoliubov theory of small perturbations of the collective excitations on top of the ground state. Thanks to the Z 2 spinsymmetry of the problem under ψ 1 ↔ ψ 2 , excitations can be classified as symmetric, density (d ) modes, and antisymmetric, spin (s) modes [45].

Uniform condensate
For a uniform condensate of total density n, the energy of collective excitations is given by the Bogoliubov dispersion relations as a function of the excitation wavevector k [45]: where µ d = (g + g 12 )n/2 and µ s = (g − g 12 )n/2 are the density and spin chemical potentials.
In the absence of Rabi coupling Ω = 0, both dispersions are linear at low momenta, with speed of sound equal to c d = µ d /m and c s = µ s /m in the density and spin channel, respectively; the presence of two gapless Goldstone modes is associated with the conservation of the total and relative numbers of particles, N 1 ± N 2 . For a finite Rabi coupling Ω = 0, instead, the breaking of the U (1) symmetry related to the relative phase of the order parameters of the two components is responsible for the opening of a gap of size ω p = Ω(2µ s /ħ + Ω) in the dispersion of spin modes, which acquire a finite effective mass [39].

Spatially inhomogeneous condensate
For spin-symmetric, but spatially inhomogeneous systems, the Bogoliubov spectrum can only be obtained numerically (details are given in the Appendix): it is first neceassary to compute the density profile of a stationary state, that is, of an eigensolution of the stationary GP equation, where Φ(r) is the order parameter of both components, ψ o 1 (r) = ψ o 2 (r) = Φ(r) and µ is the oscillation frequency of the matter field. Of course µ takes the meaning of chemical potential when we consider the lowest-order stationary state of (4) corresponding to the ground state of the condensate. In practice, the stationary state is obtained through imaginary-time evolution.
The frequencies ω i and spatial shape [u i (r), v i (r)] of the collective excitation modes on top of the stationary state are then computed as the eigenvalues and eigenvectors of the Bogoliubov operators in the density d and spin s sectors: where the operators on the diagonal read: and the interaction energies acquire a spatial dependence due to the inhomogeneity of the density profile, µ d ,s (r) = (g ± g 12 ) n(r)/2. For each collective mode, the time-evolving perturbation of the order parameter of each species then has the form: where the two components u i and v i of the eigenvector are often referred to as particle and antiparticle components. Depending on which of the two dominates over the other, excitation can have positive, negative or zero Bogoliubov norm: In this work we exploit the analysis of the Bogoliubov spectrum [46] to study the stability of stationary states, in particular vortex configurations. In stable configurations, all frequencies are real-valued and, for each excitation mode, have the same sign as the norm. Modes with positive frequency and negative norm (or vice versa) have instead a negative energy and their existence is a signature of an energetic (or thermodynamic) instability, whereas the configuration is dynamically unstable if some modes have a complex frequency. Thanks to the pseudo-hermiticity of the Bogoliubov matrix, these zero-norm modes always come in pairs, whose frequencies share the same real part and have opposite imaginary parts, ω = η ± i Γ. The imaginary part Γ represents the rate of exponential growth of the dynamically unstable mode.

Analog gravity and vortices in two-component BECs
Analog gravity is based on the formal equivalence between the Klein-Gordon equation for the dynamics of a scalar field on a curved spacetime and the linearized GP equations in the hydrodynamic regime in which all excitations propagate with the same speed of sound. In the particular case of excitations on top of a two-component BEC, the effective metric tensors of density and spin excitations are different, but share the general form [1]: where v(r) is the (irrotational) velocity field associated to the background condensate flow: Let us now focus on the specific case of two-dimensional systems. Even though condensation is strictly speaking not possible in 2D, effectively two-dimensional atomic superfluids can be realized experimentally by tightly confining the condensate along one direction, say z, imposing a much weaker confinement in the orthogonal plane (x, y), and working at a sufficiently low temperature for the coherence length to exceed the finite size of the system. The main effect of the vertical confinement is a geometrical renormalization of the interaction strength.
In the following we focus our attention on a radially symmetric two-dimensional BEC confined in the (x, y) plane. Thus r = (x, y) and polar coordinates are defined as θ = arctan(y/x) and r = x 2 + y 2 . For such configurations, the analog spacetime defined by (10) may display a density and a spin ergosurface, defined by |v| = c d ,s , while analog event horizons are found if the radial fluid velocity equals the speed of sound, |v · u r | = c d ,s and then exceeds it. Of course, a radial flow requires the presence of some particle non-conservation mechanism. A study of this physics goes beyond this paper, which is restricted to purely tangential flows with no radial component.
The ground state of a non-rotating condensate is well approximated by the Thomas-Fermi (TF) solution: In the presence of rotation, a spin-symmetric quantized vortex is a solution of the stationary Gross-Pitaevski equation (4) of the form: where the radial function f (r ) defines the vortex profile and L is the integer-valued vortex charge.
The density profile of a vortex differs from the TF one only in the region of the vortex core, while n(r ) → n TF (r ) at larger distances. For a charge-L vortex, the core has a size on the order of L ξ d , where ξ d = ħ/mc d is the density healing length. The spin healing length is defined in an analogous way as ξ s = ħ/mc s .
The plot is obtained using the numerically calculated density profile of a vortex of charge L = 2; interactions are set to g 12 As usual, the velocity flow of the vortex is tangential and inversely proportional to the distance r from the center: v(r ) = ħL Assuming the two speeds of sound c d ,s are slowly varying, the analog ergosurfaces for density and spin modes, defined by |v| = c d ,s , are approximately located at r ∼ Lξ d ,s . An example of promising configuration is shown in Fig.1, where the local values of c d ,s have been computed using the exact density profile of an L = 2 vortex in a very wide trap and the interactions constants have been chosen such that the large-r asymptotic values c s c d and ξ s ξ d . This condition guarantees that the the spin ergoregion is located well outside the vortex core, where the density profile is approximately flat. The absence of a radial flow imposed by particle condensation prevents the system from having an analog event horizon, that would act as an absorbing boundary: the vortex is thus the ideal background to study analogs of the ergoregion instabilities of rotating space-times.
Taking advantage of the cylindrical symmetry, the Bogoliubov excitations of frequency ω and angular momentum M on top of a charge-L vortex can be written in the factorized form and the stability of the vortex configuration can be assessed by solving the effectively 1D radial Bogoliubov problem for u i (r ) and v i (r ). The Bogoliubov matrices (5) involve the kinetic energy operator. In polar coordinates this reads: where the + (−) sign must be chosen when applying the operator to the particle (antiparticle) component u (v).

Local density approximation
In general, the Bogoliubov problem on top a spatially inhomogeneous vortex configuration can hardly be solved analytically. Before proceeding with the numerical diagonalization of the  (19).
Bogoliubov matrices which will be the subject of the next Section, we start here with an intuitive discussion of the physics on the basis of a local density approximation (LDA). Let us assume that the condensate is enclosed in a large box trap, so that the density can be considered approximately uniform far from the vortex core, n(r Lξ d ) n ∞ . Within a small spatial region around r, excitations of angular momentum M can be approximated as plane waves with wavevector k = ku r + (M /r )u θ , so that the Bogoliubov wavefunctions u(r), v(r) ∝ e i kr e i M θ . By exploiting the LDA, the k-dependent Bogoliubov dispersion relation in the radial direction is obtained by taking the uniform result (3) and replacing k 2 → k 2 + M 2 /r 2 . The effect of the finite M is thus to open a gap in the dispersion relation and is equivalent to having an additional space-dependent coherent coupling Ω M (r ) = ħM 2 /2mr 2 . On top of this, one needs to include the overall M -dependent Doppler shift δ(r ) ≡ k · v = ħLM /mr 2 due to the tangential velocity flow. Due to the absence of in-going flow, this Doppler term does not tilt the dispersion relation as a function of the radial k, but only rigidly shifts it vertically along the frequency axis. This shift is crucial to bring negative-norm (positive-norm) modes up (down) to positive (negative) frequencies.
Combining all these elements together, we obtain the final formula where the sign ± refers to the positive and negative branches, respectively, and the effective chemical potential and coupling strength are µ ∞ s = (g − g 12 )n ∞ /2 andΩ(r ) = Ω + Ω M (r ). From this formula, the space-dependent centrifugal gap turns out to be ∆(r ) ≡ Ω (r ) Ω (r ) + 2µ ∞ s /ħ . Because of the Doppler shift and the centrifugal gap, for a generic positive frequency ω, there can be positive norm, negative norm or no modes available depending on the radial position r . Three specific examples for different radial positions r are shown in Fig.2(a-c) for the most relevant Ω = 0 case; including a finite Ω would simply result in a larger centrifugal gap. A summary of the results for generic (r, ω > 0) is displayed in panel (d): here, the red and gray shaded areas indicate the (r, ω) regions where positive and negative norm modes are available, respectively. The boundaries of these regions identify two frequency-dependent positions in real space, r ± (ω), indicated by the red and black solid lines in Fig.2(d). In particular, r − (ω) can be thought as a sort of effective frequency-dependent ergosurface position.
An ergoregion instability occurs when a negative norm mode in the inner part of the system is coupled to a positive norm mode living in the outer region, so that the two combine into a pair of zero-norm modes with complex frequencies ω±i Γ [34]. In other words, the dynamical instability can be thought as resulting from a tunnelling process between the red-and gray-shaded regions in Fig.2(d), with a tunneling-mediated instability rate determined by the real-space width of the forbidden white region separating them. According to this interpretation, modes with smaller frequency, that is, hydrodynamic excitations, have a more extended antiparticle component, but also a smaller instability rate since they have to cross a much wider forbidden region.
In general, a necessary (but not sufficient) requirement for the occurrence of an ergoregion instability is the presence of negative norm modes in the inner region, which translates in the condition δ(r ) ≥ ∆(r ) between the Doppler shift and the centrifugal gap. The limiting radius r E for which such condition is satisfied can be thought as the boundary of the ergoregion, that is the ergosurface. In the absence of Rabi coupling (Ω = 0), the effective ergosurface for angular momentum M modes is located at a radial position This value is slightly reduced with respect to the hydrodynamic prediction Lξ s , which is recovered in the M → 0 limit, with the correction being due to the non-linear behaviour of the dispersion relation. For larger M values, the ergoregion shrinks in space and eventually disappears for large angular momentum modes M ≥ 2L. As mentioned above, the presence of a coherent coupling Ω = 0 results in a larger centrifugal gap [white region in Fig.2(d)], so the matching between the external positive-norm modes (redshading) and the negative-norm ones in the ergoregion (gray-shading) modes is harder to obtain and the unstable modes, if present, are more localized and display a weaker instability rate.
While this discussion based on the LDA provides an intuitive understanding of the microscopic process underlying the ergoregion instability, it completely misses all those features that stem from the quantization of the radial wavevector k in a finite-size system 1 . Since the radial wavevector is associated to the radial kinetic energy contribution ħ 2 k 2 /2m, we expect that for tight enough ergoregions, the minimum value of the kinetic energy set by the quantization of the negative norm mode may increase the size of the effective gap so much that it prevents the development of the instability. Assessing all these features in a quantitative way requires going beyond the LDA and solving the full Bogoliubov problem with numerical tools. This will be the subject of the next Section.

Bogoliubov spectrum
The numerical solution of the Bogoliubov problem is performed as follows (more details can be found in the Appendix). For the chosen value of the vortex charge L, we first find the radial function f (r ) that describes the density profile of the stationary vortex and the associated oscillation frequency µ. These are then used to build the Bogoliubov matrices for the density and spin excitations, whose diagonalization gives the eigenstates and the corresponding eigenvalues, that is, the Bogoliubov spectrum on top of the stationary state. In what follows we present results for harmonically trapped systems, but no qualitative difference is observed for box traps (see Appendix). In both cases we choose the parameters so that the radial size of the condensate R is much larger than the vortex size. Fig.3 and Fig.4 show the Bogoliubov spectra in the spin channel for different angular momenta M as a function of the interaction strength ratio g 12 /g for vortices of charge L = 1, 2 in the Ω = 0 case: in correspondence of the crossing points between a positive and a negative energy mode, we observe the appearance of dynamically unstable modes (red dots) with the associated imaginary frequency bubble. The insets of Fig.3(panels c, f ) and Fig.4(panels c,d,g) show examples of the spatial profiles of the particle and antiparticle components of the dynamically unstable modes: the common feature consists in the localization of the antiparticle component (red solid line) in the inner portion of the system, while the particle component (black solid line) is spread throughout the whole volume occupied by the fluid.
In the case of singly-charged L = 1 vortices, in the density channel a negative energy density mode is always present, but its frequency remains below that of positive norm modes (not shown); hence a dynamical instability never develops, in agreement with previous works [34]. On the other hand, the spin channel shows energetic stability for attractive interspecies interactions g 12 < 0 (not shown) and alternated intervals of dynamical stability and instability for repulsive interactions g 12 > 0 [see Fig.3(a-b)]. Similarly to what happens for density fluctuations of multiplecharged vortices in single-component BECs, this originates from the crossing of a spectrally isolated negative norm mode with the band of positive norm ones; the discreteness of the band, due to the finite size of the cloud, explains the existance of instability bubbles separated by intervals of dynamical stability [34]. Moreover, since the frequency of the isolated mode can be arbitrarily high, the associated instabilities are typically well localized in the vortex core [see Fig.2(d)]. modes is shown in panels (c) and (d) for g 12 = 0.75g and g 12 = 0.95, respectively. The main difference between the two modes is in the localization of the antiparticle component, which dominates up to r − ∼ 3ξ d for the former, and up to r − ∼ 14ξ d for the latter; indeed, according to our LDA treatment, lower frequency modes are more extended. Remarkably, the low frequency mode is only present for g 12 ∼ g . The M = 3 channel shows a single instability: an example of the real-space profile of the mode is shown in panel (g) for g 12 = 0.81g .
As it is well known [30], when multiply-charged L > 1 vortices are considered, dynamical instabilities of the same nature appear in the density channel too: since the energy of a vortex with charge L is larger than that of L singly-charged vortices, multiply charged vortices may be unstable against splitting into several lower-charge vortices. For instance, for L = 2 vortices, we find, as a function of the interparticle interaction, alternate intervals of dynamical stability and instability in the M = 2 channel of the density Bogoliubov problem (not shown), again in agreement with previous work [34].
An additional class of dynamical instabilities occur in the spin channels: besides the highfrequency isolated modes crossing the positive energy band discussed above, which are present for all angular momenta 1 ≤ M < 2L [see Fig.3(d), Fig.4(a,e)] and which are related to the deformation of the vortex core and/or its splitting, low-frequency and spatially extended instabilities may also appear. As visible in Fig.4(a), such an instability may originate from the crossing of the two bands of closely spaced positive and negative norm modes that extend throughout the ergoregion. As such it belongs to the class of superradiant ergoregion instabilities closer to the hydrodynamic regime. The possibility of observing instabilities of this superradiant kind is one of the main interests of replacing single-component systems with spin mixtures with a much wider ergoregion.
Due to their different origin, these additional modes appear only around zero frequency for g 12 g close to the demixing point: most importantly, their antiparticle component extends well outside the vortex core [see Fig.4(d) as compared with Fig.4(c)]. Physically, this means that long-wavelength spin waves are involved in the instability mechanism. From a quantum gas perspective, the instability can be understood as a demixing instability whose appearance is facilitated by the presence of a vortex which lowers its threshold. From the point of view of the gravitational analogy, this instability is of superradiant nature: as it happens for ergoregion instabilities around massive rotating objects, the positive energy of the wave in the outer region is compensated by the negative energy of the wave inside the ergoregion.
Analogously, in the presence of a coherent coupling, low-frequencies instabilities appear close to the ferromagnetic phase transition point g 12 g + 2Ω/n, where the massive gap in the spin dispersion closes and the spin healing length diverges. This will be the subject of future work.
In general, depending on the parameters' values, multiple dynamically unstable modes can be found for the same configuration. When observing the time evolution of the system starting from an unperturbed L-charged vortex, after a short transient in which all the unstable modes compete, the mode with the largest imaginary frequency will eventually win over the others and dominate the intermediate-time evolution. In most cases, the dominating mode is a highfrequency localized one. However, it is possible to play with the system size R and with the interaction constants ratio g 12 /g to make all localized modes dynamically stable, while keeping a spatially extended mode unstable via superradiant mechanisms. In this case the dynamical instability is associated to the generation of long-wavelength phonon within the ergoregion, as expected from the gravitational analogy with ergoregion instabilities of massive rotating objects [24].
At very long times, nonlinear effects set in, resulting in complex mode mixing phenomena and additional instabilities. A study of this physics will be the subject of the next Section.

Long-time dynamics
It is well-known that the linearized Bogoliubov theory holds as long as the the excited modes are small perturbations on top of the stationary vortex state. When dealing with a dynamical instability, this approximation is valid for a limited amount of time, roughly given by the inverse of the growth rate Γ of the unstable mode. After that, the cylindrical symmetry of the problem is likely to be broken by complex nonlinear processes, and mixing phenomena start occurring.
In order to access the long-time dynamics of the spin mixture, it is necessary to simulate the full GP equations (1) in two dimensions. The numerical protocol is the following (more details can be found in the Appendix): the stationary state corresponding to a vortex of charge L sitting at r = 0 is first obtained via a conjugate gradient algorithm [47]. In order to trigger the instability, we perturb the stationary state with some weak random noise, independently applied on the total and relative density. The temporal evolution of the system under the GPE (1) is then simulated using a split-step method.
Examples of the simulated time-evolution are presented in Figs.5-8 for vortices of charge L = 1, 2. As already discussed, L = 1 vortices have a single M = 1 potentially unstable mode in the spin channel. On the other hand, L = 2 vortices can feature up to five unstable modes, four of which in the spin channel (one for M = 1, two for M = 2, one for M = 3). For the sake of simplicity, the simulation parameters are chosen in a way to have a single dynamical instability in the spin channel: in practice, this is done by keeping the interaction ratio g 12 /g fixed and tuning the system size R to select the desired mode.
In all these figures, the short-and intermediate-time dynamics is consistent with our expectations based on the Bogoliubov theory: the instability develops initially as a well visible perturbation characterized by 2M lobes in the azymuthal direction, where M is the angular momentum of the unstable mode. For high-frequency instabilities, this excitation has a larger-amplitude part localized around r = 0 and a smaller-amplitude one visible as interference fringes at larger Figure 5. Dynamics of the vortex splitting instability. Real time evolution of a L = 1 vortex with R = 21ξ d and g 12 = 0.8g . For these parameters, the growth rate of the M = 1 unstable mode, whose real-space profile is shown in panel (c), is Γ/µ 0.014. Panels (a1-a6): total density n = n 1 + n 2 , measured in units of the peak Thomas-Fermi density 2n TF (r = 0) = 2µ/(g + g 12 ). Panels (b1-b6): polarization, Z = (n 1 − n 2 )/n. Each column is computed at the time indicated by one of red dots in panel (d), which shows the temporal evolution of the standard deviation of the spin density δn = n 1 − n 2 (black solid line), compared with the analytical exponential growth exp(Γt ) (blue dashed line). from its ground state vanishing value, defined as std(δn) = 1 where Q × Q is the number of points in the 2D numerical grid. As expected, this quantity is exponentially growing at a rate equal to the imaginary part Γ of the frequency of the dynamically unstable mode, whereas the total density remains almost unperturbed. At longer times, t Γ 1, for the intermediate values of g 12 /g considered in Figs.5-7, the high-frequency dynamically unstable modes lead to the splitting of the original vortex into smaller vortices and/or pairs of half-quantized vortices. Two examples are shown in Figs.5-6 for L = 1 and L = 2 vortices, respectively, and closely resemble the results of Ref. [41]. Interestingly, the splitting is often followed by a recombination process; similar dynamics was predicted for multiply-charged vortices in a single-component BEC in [48]. We attribute this phenomenon to energy conservation and to the finite size of the system: once the vortex is split, the excess of energy is released in the form of spin excitations, whose interference, after bouncing on the trap walls, reverses the splitting process. This typically occurs multiple times, before, in the absence of a dissipation mechanism, the evolution becomes turbulent, with additional vortices being nucleated from the boundaries and/or deformation of the (otherwise circular) cloud [see panels (a6-b6) of Fig.6]. All these complex nonlinear features go beyond this work and will be the subject of future studies.
The physics is more intriguing in the g 12 g regime close to the demixing point, where vortex splitting is much harder to observe and, if it does so, only occurs at much longer times scales: in addition to the reduced value of the instability rate Γ, the softness of spin modes prevents in fact the cloud from absorbing the excess energy that would derive from vortex splitting. This result is confirmed by looking at the Bogoliubov spectra in Figs.3-4: for g 12 /g 1, the instability bubbles disappear or become tiny, and instability rates drop by roughly an order of magnitude. As a consequence, when approaching the demixing point, vortices in condensate mixtures remain stable for a very long time, while their spin dynamics is ruled by linear Bogoliubov theory. Two examples of the real-time dynamics in this regime are shown in Figs.7-8 for L = 1, 2 vortices, respectively, and g 12 /g = 0.97. In the former case, the unstable Bogoliubov mode is spatially localized around the vortex and corresponds to a displacement of the cores in the two components, with the consequent appearance of a net polarization. This process, however, does not lead to a full splitting of the vortex. We rather observe a sequence of intervals of suppression and revival of the instability, as witnessed by the oscillation in time of the standard deviation (STD) of the relative density. The underlying mechanism is rooted in the finite-size of the system and is analogous to the one leading to the sequence of vortex splitting and recombination events seen above.
The temporal evolution associated to a superradiant instability is finally visible in Fig.8. At intermediate times [panel (b2)], a spatially extended M = 2 modulation is clearly visible in association to the vortex core deformation. Its characteristic size corresponds to the spatial profile of the unstable mode shown in panel (c) and, in particular, its antiparticle component matches the extension of the ergoregion. This is the smoking gun of the superradiant nature of the instability.
At later times, the dynamics develops instead a more complicated behaviour that involves strong nonlinear mode-mixing effects: the perturbation of the spin density, due to the development of the instability, is large enough to excite the M = 2 waves with higher radial momentum and several radial nodes that are visible in panel (b3). These different excitation patterns are in competition, eventually leading to a marked deformation of the vortex and a significant modulation also of the total density [see panel (a6)]. Once again, the accurate analysis of these complex mode-mixing processes goes beyond the scope of the present work.

Conclusions and perspectives
In this work, we have theoretically analysed the stability of quantized vortices in symmetric twocomponent atomic Bose-Einstein condensates (BEC) and we have interpreted the results within an analog gravity context in terms of ergoregion instabilities. In addition to instabilities related to high-frequency Bogoliubov modes localized in the vortex core and associated to the distortion or the splitting of the vortex core like in single-component BECs [34], suitable regimes are found where the physics is determined by low-frequency and long-wavelength spin excitations. In this case, the superradiant scattering process underlying the ergoregion instability involves the excitation of a Doppler-shifted, negative-energy spin-sound wave spread over the ergoregion and the simultaneous emission of positive-energy spin waves into the outer part of the BEC, in close analogy with the ergoregion instability of space-time around a rotating massive object.
On the theoretical side, the most challenging open questions concern the late time dynamics of unstable configurations. First studies of this physics focused on the simpler case of black hole laser instabilities [49] and first hints of the remarkable complexity of the superradiant case are visible in the simulations reported in this work. In this long-term adventure, a special attention has to be paid to the role of nonlinear processes in saturating the instability, in analogy with related phenomena predicted in the gravitational context such as the growth of the so-called black hole "hair" [50][51][52][53].
On the experimental side, two-component condensates can be obtained in fluids of atoms and of light. In the atomic case, vortices can be generated by means of suitably chosen stirring potential [31,54] and the different components can be chosen within the atomic hyperfine structure in a way to obtain slow spin-sound waves [39]. In the optical case, the polarization degree of freedom can be used to obtain two-component condensates and arbitrary velocity patterns can be imprinted onto the fluid by suitably designing the phase profile of the pump beam so to generate rotating configurations [55]. In both cases, the control of the coherent coupling strength via either an external electromagnetic dressing of the atoms or the optical birefringence of the cavity material may be used to tune the mass of the quantum field.
The remarkable experimental advances in both fields make us confident that it will be soon possible to validate our conclusions using state-of-the-art set-ups. As a more ambitious challenge, the analysis of correlations in the spirit of [11,56,57], will allow to investigate superradiant processes at the quantum level, so to prove superradiant amplification is intrinsically connected to the spontaneous creation of correlated pairs of Bogoliubov modes with opposite energy at the ergosurface.

Acknowledgements
AB aknowledges financial support from Q@TN, the joint lab between University of Trento, FBK-Fondazione Bruno Kessler, INFN-National Institute for Nuclear Physics and CNR-National Research Council. IC acknowledges financial support from the European Union H2020-FETFLAG-2018-2020 project "PhoQuS" (n.820392), from the Provincia Autonoma di Trento, and from the Q@TN initiative.

Numerical details
The numerical calculation of the Bogoliubov spectrum is performed as follows. Given a value of the vortex charge L, we first find the radial function f (r ) that describes the density profile of the stationary vortex and the associated oscillation frequency µ. This is achieved by means of an imaginary-time evolution . Once the vortex profile f (r ) is known, for each value of angular momentum 1 ≤ M < 2L, we build the density and spin Bogoliubov matrices and diagonalize them: the eigenvalues form the Bogoliubov spectrum of the vortex, while the eigenvectors give the real-space profiles of the modes.
The free parameters of the simulation are the interaction ratio g 12 /g , the Rabi frequency Ω, the radial size of the system R and the expected chemical potential in the TF regime (i.e. without vortices) µ T F . The remaining parameters, i.e. the number of particles N and the trapping frequency ω 0 , are chosen accordingly; in the case of a uniform mixture in a box of radius R, we set: N πR 2 = 2µ T F + ħΩ g + g 12 (21) while, for a system in an harmonic trap of TF radius R, we choose: While we can take advantage of the axial symmetry of the system to compute the Bogoliubov spectra, the same strategy is not applicable to the problem of determining the long-time dynamics of the mixture: we expect indeed the cylindrical symmetry to be broken by the development of the instability.
We thus solve the full GP equations in 2D, according to the following numerical protocol: we first identify, thanks to the Bogoliubov spectra, a set of parameters leading to a dynamical instability; we calculate the 2D profile of the stationary vortex via a conjugate gradient algorithm [47], where the azymuthal profile of the wavefunction phase is enforced at every step. The stationary state is perturbed with some weak random noise, applied independently on the total and relative density, to trigger the instability. We then simulate the time evolution of the system given by the GPE (1) using a split-step method.

Dependence on the size and on the external potential
In order to assess the dependence of dynamical instabilities on the size of the system and on the external potential applied to the atoms, we performed a calculation of the Bogoliubov spectra as a function of R for two experimentally relevant cases: the harmonic trap and the box. The results are shown in Fig.9 for L = M = 1: in agreement with Ref. [34], for both cases we find alternate intervals of dynamical stability and instability; as shown in panels (b) and (e), the instability bubbles acquire a larger width but smaller amplitude as R increases. This is the reason why in the main text, if possible, we simulate the real-time dynamics of systems with relatively small size (R ∼ 20ξ d ): in addition to the advantage of allowing for a better resolution in space, such configurations also feature a faster development of the instabilities. As a further check, panels (c) and (f) of Fig.9 show that there is no qualitative difference in the real-space profiles of the unstable modes given in the two examined configurations. We verified that this holds true also for other values of vortex charge L and angular momentum M (not shown).

Figure 9.
Bogoliubov spectra for M = 1 spin excitations over a vortex of charge L = 1 with g 12 = 0.93g . Panels (a-c) refer to an harmonically trapped system of TF radius R, while panels (d-f) refer to a uniform mixture in a box of size R. Panels (c) and (f ) show the realspace profile of the particle (black solid line) and antiparticle (red solid line) components of the dynamically unstable mode for R = 23ξ d ; the dashed line represents the rescaled vortex profile f (r ).