Quasi-normal modes and fermionic vacuum decay around a Kerr black hole

We analyze the instability of the non-rotating fermion vacuum in Kerr spacetimes. We describe how the co-rotating Fermi sea is formed as a result of a spontaneous vacuum decay. Most significantly, and drawing upon intuition gained from analogous electrodynamic processes in supercritical fields, we show that this decay process is encoded entirely in a subset of quasi-normal fermion modes.


I. INTRODUCTION
Nothing can escape from the event horizon of a black hole, not even light, and for that reason, black holes are often thought of as perfect absorbers. Despite this, around rotating black holes, bosonic fields can extract angular momentum from the black hole and subsequently see their energy amplified. For this to occur, their angular frequency ω must satisfy ω < mΩ H , where m is their azimuthal number and Ω H is the angular velocity of the black hole. This process, known as 'superradiance', has been studied thoroughly since its discovery (see Ref. [1] for a recent review) and has even been tested experimentally in a water-wave black-hole analogue [2]. If the black hole is surrounded by a perfect mirror, this amplification process is turned into a dynamical instability: waves propagate back and forth from the mirror to the black hole, amplifying in energy each time and leading to an exponential blow up called the 'black hole bomb'. Nature "can also provide its own mirror" [3]: when the field is endowed with a non-zero mass, low-frequency modes are trapped, thereby triggering a black hole bomb instability [4].
On the other hand, fermionic fields around rotating black holes cannot be amplified [5,6] due to the Pauli exclusion principle. However, for a massless fermion, modes satisfying the condition ω < mΩ H display a vacuum instability: the ergoregion spontaneously emits a steady flux of particles in that frequency range. This is known as the Unruh-Starobinski radiation [7,8]. This dual behavior of superradiance is reminiscent of the Klein paradox [9], where amplification occurs only for bosons, but both bosons and fermions can experience a spontaneous emission in strong external fields. When considering fermions with a non-zero mass µ around a rotating black hole, it has been argued that modes in the range ω < mΩ H and ω < µ form a Dirac sea outside the black hole [10]. This Dirac sea carries a non-zero angular momentum vacuum expectation value, aligned with the one of the black hole.
In this work, we show that the Dirac sea is the result of a spontaneous vacuum decay, which constitutes the fermionic pendant of the black hole bomb instability. Indeed, shortly after a rotating black hole is formed, the fermion field is in a vacuum state resembling the Unruh vacuum and carrying zero angular momentum. Subsequently, the Unruh-Starobinski radiation sets in and fills a set of trapped modes outside the black hole until Pauli's exclusion principle stops the process. Through this mechanism, the initial, non-rotating vacuum decays into a new vacuum (the Dirac sea), which co-rotates with the black hole. As we will show, this decay process is entirely encoded in the set of quasi-normal modes (QNM) that satisfy both conditions ω < mΩ H and ω < µ. Notice that the structure of the vacuum states for massless fermions have been studied by various authors [11,12], and there is no vacuum decay in this case.
The paper is organized as follows. In section II, we quickly review the Dirac equation in the Kerr spacetime. We then use a WKB approximation to discuss the structure of eigenmode solutions and the set of QNMs in the range ω < mΩ H and ω < µ. In section III, we present a toy model for the vacuum decay instability, comprising a charged fermion in a strong electric field and tailored to resemble the Kerr problem. In section IV, we discuss the dynamics of the vacuum decay in detail, as well as the main properties of the initial and final states. Throughout the paper, we work in units such that c = = G = 1.

II. FERMIONS IN KERR
We consider a rotating black hole, whose spacetime is described by the Kerr metric. Using the Boyer-Lindquist coordinate system, the line element reads 2M r ρ 2 dt 2 − 4aM r sin 2 θ ρ 2 dt dφ + ρ 2 ∆ dr 2 + ρ 2 dθ 2 + r 2 + a 2 + 4a 2 M r sin 2 θ ρ 2 sin 2 θ dφ 2 , where ρ 2 = r 2 + a 2 cos 2 θ and ∆ = r 2 − 2M r + a 2 . This spacetime possesses many interesting features and has been studied extensively in the literature [13][14][15]. In particular, there is an event horizon at r = r + = M + √ M 2 − a 2 (and a Cauchy horizon at r = r − = M − √ M 2 − a 2 , although the latter won't play any role in what follows). In addition, the asymptotically time-like Killing field ∂ t becomes space-like close to the black hole, when 2M r > ρ 2 . This defines the ergoregion, where any time-like trajectory must co-rotate with the black hole. To characterize the rotation of the black hole, we define its angular velocity As we shall see, Ω H plays an important role in characterizing the vacuum instability that we will analyze. In addition, it is quite convenient to introduce the tortoise coordinate r * , adapted to scattering problems from the horizon to infinity. It is defined via We see that r * goes from −∞ to +∞ when r runs from r + to +∞.

A. The Dirac equation in Kerr
We now consider non-interacting fermions of mass µ in this spacetime. These are described by a Dirac field ψ, which obeys the equation where D ν is the spinor covariant derivative. It turns out that, in the Kerr metric, both the massive and massless Dirac equations admit a complete separation of variables. This was first shown by Unruh [16] for the massless case and later by Chandrasekhar for the massive case [13,17]. In this paper, we shall use the same conventions as in Ref. [18]; namely, we shall work in the Weyl/Chiral representation, where in which I 2 = diag(1, 1) is the two-dimensional unit matrix and the σ j (j = 1, 2, 3) are the Pauli matrices. We refer the reader to that work [18] for a detailed derivation of the Dirac equation and its separation, or our Appendix A, where we summarize the various conventions.
Due to the symmetries of the metric, the eigenmodes of the Dirac equation (4) can be written in the form where the complex variable z = r + ia cos θ. In addition, ω is the frequency, and m is a half integer giving the azimuthal number. The above ansatz leads to an ordinary differential equation for both the radial and angular parts. First, the angular parts obey the coupled equations j,m,P + aµ cos θ S 2 , (7a) where λ (ω,µ) j,m,P is the separation constant, which is closely related to the total angular momentum j ∈ N + 1/2 (see Ref. [19] for a detailed discussion of the angular equation). In the non-rotating limit, i.e. when a → 0, λ ω reduces to P(j + 1/2), where P = ±1 is the parity, which we recognize in terms of the eigenvalues of the spin-orbit coupling. However, in this work, we are instead interested in sufficiently rapidly rotating black holes, and λ (ω,µ) j,m,P is, in general, a complicated function of j and P, as well as aω and aµ. As we shall see, the modes responsible for the vacuum decay essentially satisfy 0 < ω µ. It turns out that the angular equations for ω = µ can be solved exactly [19], and one can use these solutions to obtain leading-order expressions. In this regime, the separation constant is given by To simplify our notation, we hereafter drop the indices on λ. Once λ is determined, the radial part obeys the coupled equations with The aim now is to solve equation (9) with appropriate boundary conditions.

B. Complete mode basis
Once the mode equations (7) and (9) have been solved, general solutions of the Dirac equation (4) are obtained by linear superpositions. In particular, the field operator decomposes into a sum of eigenmodes, where the associated operator-valued amplitudes are the familiar creation and annihilation operators. Formally, the decomposition has the form Once (j, m, P) is chosen, the modes are determined by solutions of equation (9) with specific boundary conditions. (In the following, we drop the superscript (j, m, P), as we have done for the eigenvalue λ.) There are now two families of modes. The first one is defined at infinity (r → ∞) as incoming towards the black hole. These modes are asymptotically onshell and hence exist only above the mass gap ω > µ. Following standard conventions, we refer to them as in-modes. Moreover, we have The second family is defined near the horizon (r → r + ) as emanating from it, and we refer to these as up-modes. The positive and negative continua are related by Therefore, the field operator has the general eigenmode decomposition whereâ ω andb ω satisfy the canonical anti-commutation relations. These anti-commutation relations are guaranteed when using normalized modes with respect to the (conserved) scalar product, given by Using the ansatz (6), the norm of a mode separates into Notice that the radial part takes the canonical form R dr * (|R 1 | 2 + |R 2 | 2 ) using the tortoise coordinate (3). Since the stationary modes form a dense set, we must use a Dirac normalization, that is, we require with a vanishing overlap for differing values of (j, m, P). As usual, when dealing with a dense set of modes, it is simpler to normalize the conserved current. Hence for the radial part, we simply require As we show in detail in appendix B, the conservation law [18] with dΩ ρ 2 J t = R dr * (|R 1 | 2 + |R 2 | 2 ) guarantees that the current normalization (18) is equivalent to the norm (17b).

C. Eigenmodes and QNMs in the WKB approximation
To understand the structure of the eigenmode solutions of equation (9), we shall employ a WKB method. The WKB approximation is particularly good in the regime M µ 1 [18]. Since, from now on, we shall mainly focus on the radial part, we use the condensed notation To proceed with the analysis of equation (9), we assume and that the amplitudes R 0 j vary much more slowly than the phase dr k(r). This leads to the 2 × 2 linear system This system admits non trivial solutions if the determinant of the coefficient matrix vanishes. This gives the dispersion relation or, equivalently, the Hamilton-Jacobi equation Moreover, equation (22) gives a linear relation between R 0 1 and R 0 2 . The additional requirement of unit current (18) leads to the full expression of the amplitudes for a given k; namely, As we see in the Hamilton-Jacobi equation (23), depending on the sign of the right-hand side, the modes are either propagating or evanescent. A convenient way to discuss its sign is to recast equation (23) in the form [20] Around a certain radius r, the modes propagate if ω > ω + or if ω < ω − . Hence, ω ± conveniently delimitate the positive-and negative-frequency continua. The explicit expressions for ω ± follow from the Hamilton-Jacobi equation (23), and we have In Fig. 1, we have represented ω ± as functions of r. We first notice the behavior of ω ± on both asymptotic regions (i.e. near the horizon and at infinity): A remarkable fact of rotation is that frequencies satisfying 0 < ω < mΩ H can tunnel from the lower continuum close to the horizon to the upper continuum outside the potential barrier. If such a frequency is also above the mass gap ω > µ, it tunnels out and propagates to infinity. This gives rise to the Unruh-Starobinski radiation. In this work, we are interested in the frequencies below the mass gap that are also able to tunnel to another propagating region (such that ω > ω + ). This happens if ω is above the local minimum of ω + (see Fig. 1), which we call ω min . These frequencies therefore satisfy ω min < ω < µ and ω < mΩ H . In this case, there are three turning points r j (j = 1, 2, 3) such that ω ± (r j ) = ω, and the modes are propagating for r < r 1 and r 2 < r < r 3 . On the horizon side, for r < r 1 , the modes are superpositions of right movers and left movers with a scattering phase shift e 2iδω [21], i.e.

R(r)
Notice that, at this point and onward, k denotes only the positive root of the Hamilton-Jacobi equation (23). For r 2 < r < r 3 , the modes are again a superposition of left and right movers, with an overall amplitude A ω , which depends on the tunneling probability. Explicitly, we have (29) Using connection formula on the turning points [22,23] (see also appendix C), we obtain both the scattering phase shift and the amplitude in the trapped region: In these equations, I ω and S ω are the classical actions in the forbidden and allowed regions, that is The expression (30b) for A ω shows a set of resonant frequencies near cos(S ω ) = 0. These resonances are poles of A ω in the complex plane. Since they are also poles of the phase shift e 2iδω , we see that they are associated with in-going boundary conditions at the horizon and are therefore quasi-normal frequencies 1 . Since they are trapped due to the condition ω < µ, they are usually called quasi-bound states (QBS) [18] and are non-radiative, contrary to usual QNMs. It is convenient to write the quasi-normal frequencies w as a sum of real and imaginary parts, viz. w = ω R + iω I . Using the WKB analysis outlined above, and in the limit of small tunneling |ω I | |ω R |, we obtain an approximate expression for the quasinormal frequencies. The real part is obtained by solving a Bohr-Sommerfeld condition, and the imaginary part follows from the tunneling amplitude: where n is an integer and is the time needed for a semiclassical wave packet to propagate at a local group velocity v g (r) = 1/∂ ω k(r) from r 2 to r 3 . As we shall see, this particular set of QNMs encodes a vacuum instability.  (8) with (j, P) = (1, 1). The line of constant ω (dashed line) crosses ω ± three times: ω + twice at r 2 and r 3 , which forms the trapping region, and ω − once at r 1 . Right side: magnified view of the trapping region.

III. ELECTRIC TOY MODEL
To illustrate the mechanism of the vacuum decay in Kerr, we shall first consider a simple toy model in an electric field. Analogous toy models have been used in the past to study the black hole bomb of scalar fields, see, e.g., the appendix of Ref. [25]. Our model comprises a Dirac fermion in 1 + 1 dimensions, which is subject to an external electric potential V (x) and has a local mass µ(x). The two functions will be chosen to mimic the Kerr problem, while yielding much simpler calculations. Moreover, this electric model allows us to make connection with the literature on fermions in overcritical fields [26,27], which bare many similarities with black hole physics [28]. The Dirac equation for the mode of energy ω is Note that the coordinate x here is the analogue of r * in the BH case. As in previous sections, we work in the Weyl/Chiral representation. This means, in particular, that where I 2 = diag(1, 1) is the two-dimensional unit matrix and σ 3 = diag(1, − 1) is the third Pauli matrix. We assume that V (x) and µ(x) are piece-wise constant. Alternatively, one could work with slowly varying, smooth V (x) and µ(x), and use a WKB approximation, as employed for the Kerr problem in Sec. II. As we shall see, the piece-wise constant problem is already quite similar to the Kerr problem, without generating extra technicalities. In each region where both V and µ are constant, the mode solutions of equation (34) have the form where k is the wavenumber. The dispersion relation can be found by acting on equation (34) from the left with the conjugate Dirac operator. This gives We now consider a problem with two main regions in which the field is massless. For x < 0, the potential takes the constant value V 0 , and it vanishes (V = 0) for x > 0. Moreover, we assume that there is a perfectly reflecting boundary at x = L. V 0 is the electric analogue of mΩ H in the Kerr case, and the reflecting condition at x = L mimics modes in Kerr below the mass gap (ω < µ), which are trapped in a finite region near the black hole. In addition, we add a delta-barrier potential at the transition, i.e. near x = 0. It is, however, delicate to define a delta-barrier potential for fermions, due to the first-order character of the Dirac equation [29]. To obtain such a barrier, we assume that there is a small region −D < x < D where the field has a large mass µ 0 . We then take the width of the barrier to zero and the height to infinity such that the product remains constant. The parameter ξ allows us to control the 'strength' of the barrier: in the limit ξ → 0, there is complete transmission and, in the limit ξ → ∞, the barrier is impermeable. Note also that the perfectly reflecting boundary condition at x = L is obtained with the same method by taking ξ → ∞. The effective dispersion relation of the problem is represented in Fig. 2. The Dirac mode equation (34) has four linearly independent solutions: comprising both spin projections of the positive-and negative-frequency branches. Applying the thin-barrier limit described above, we obtain the boundary condition at x = 0. This boundary condition relates the value of the mode on each side. Explicitly, we have Taking the limit ξ → ∞ of this expression also yields the reflecting boundary condition at Notice that the boundary conditions mix the two chiralities, which reside in the upper and lower pairs of elements in the four-component spinor. Recalling that the spin operator is Σ = diag(σ 3 , σ 3 ), we see that the energy eigenstates are also spin eigenstates, with the two spin polarizations evolving independently. We now consider the setup where the central barrier is initially infinite, i.e. ξ = ∞ for t < 0, and is suddenly lowered to a finite (but large) value for t > 0. Initially, the eigenmodes compose a continuous set on the left side, reflected on the impermeable barrier, and a discrete set between 0 < x < L. We denote by (ϕ in ω ) ω>V 0 (resp. (ψ in ω ) ω>V 0 ) the continuous modes of positive (resp. negative) energy, and (ϕ in n ) n∈N 0 (resp. (ψ in n ) n∈N 0 ) the discrete ones of positive energy (resp. negative). Since the barrier is infinite, these two sets are uncoupled. When the barrier is lowered, the two regions become coupled. In particular, discrete modes of positive energy with ω n < V 0 are coupled to the continuous set of negative energy. As a consequence, the initial vacuum will spontaneously populate these discrete states with particles by emitting anti-particles to the left. This process is the decay of the initial neutral vacuum into the true ground state: a charged vacuum. After the barrier is lowered, for t > 0, the eigenmodes compose only a continuous set extending over x ∈]−∞; L] and denoted (ϕ in ω ) ω>V 0 for positive energies and (ψ in ω ) ω>V 0 for negative ones. As we will now show, the decay process described before is encoded in the overlap between this continuous set and the discrete set of initial eigenmodes.
In Appendix D, we present the full set of eigenmode solutions of equation (34). Since we are interested only in the overlap between the initial discrete modes and the outgoing continuum modes, we can restrict our attention to the initial discrete modes ϕ in n of the cavity [0, L] and their charge conjugates ψ in n , as well as the outgoing continuum modes on the interval [0, L]. Moreover, since the two spin eigenstates decouple, it is sufficient to consider only one spin polarization. The discrete set of in-modes is given by where The continuous set of outgoing modes has a structure very similar to the one in Kerr [see equations (29) and (28)]. On the left side (x < 0), it is a superposition of left and right movers with a scattering phase shift. However, we are mostly interested in these modes in the trapped region (0 < x < L), where they read with A ω = L π e ξ sin(ωL) + i cos(ωL)e 2ξ .
We can now decompose the field operator for both time regions aŝ Pre-empting the analysis of the Kerr problem that follows in Sec. IV, we draw attention to the similarity of the ϕ in and ϕ out modes on the interval [0, L]. Most importantly, when ξ 1, the amplitude A ω is maximal for ω ∼ ω n , such that the main distinction between the ϕ in and ϕ out is in their normalization.
We require the field operator to be continuous at t = 0. This gives a linear relationship between the two representations of the field operator, otherwise known as a Bogoliubov transformation. We have, for instance, that b out where, for brevity, we have included only the contributions from the nearby modes with which ψ out * ω (x) has significant overlap. The coefficients β n (ω) and α ωω describe the Bogoliubov transformation between the two representations. The mode basis is orthonormal with respect to the conserved scalar product of the one-particle Hilbert space, i.e.
Moreover, as we discuss above, the most relevant overlap (see also appendix D) is thanks to the identity We now have all the ingredients that we need to understand the decay of the vacuum state. To describe this decay, we first define an instantaneous occupation number for the discrete modes. This will allow us to follow their evolution from initially empty to occupied, at the end of the decay process. The natural definition of the instantaneous mean occupation number of the state |ϕ in n is [26,27] where we have defined the instantaneous annihilation operator of the (one-particle) state |ϕ in For t < 0, |ϕ in n is a stationary state, andâ n (t) evolves with the phase e −iωnt . For t > 0, |ϕ in n is no longer a stationary state, andâ n (t) is instead associated with a resonance. As a result, the vacuum state for t < 0 is no longer stationary for t > 0, and it decays to the true ground state. If one reinstates the barrier (ξ → ∞) after some time T ,â n (t) corresponds again to a stationary state, and N n (T ) will give the mean number of particles created in the state n between t = 0 and t = T . This justifies our definition of the instantaneous occupation number.
We anticipate that the state |ϕ in n will initially be empty for 0 < ω n < V 0 , subsequently decaying to an occupied state and meaning that N n (t) = 1 − (decaying terms). Hence, it is more convenient to use the anti-commutation relations to write We now use equation (51) to obtain substituting for b out ω from equation (46). At leading order in 1/ξ, the overlap ϕ in n |ϕ out ω is negligible, because ω n < V 0 < ω, and we find such that From equations (48) and (49), we see that β n (ω) is proportional to sinc[(ω −ω n )L/2]. Hence, at leading order, β n (ω) * β n (ω) ∝ δ nn , and we can drop the sum in the above equation. Using the change of variable ω → 2V 0 − ω and equation (49), it follows that All that remains is to compute the integral over ω by the residue theorem. We see from equation (48) that all the singularities of the overlap are those of F (ω) = |A ω | 2 . Once extended in the complex plane, the poles and corresponding residues of F (ω) are, to leading order in e −ξ , given by Res n (F ) = − 1 2iπ .
For t > 0, we pick up the poles in the lower-half plane. Moreover, the factor sinc[(ω−ω n )L/2] in the overlap selects only the pole w n whose real part ω R n is close to ω n . We therefore obtain and everything combines to give the simple result which holds only for t ≥ 0. We see that the modes with 0 < ω n < V 0 are initially unoccupied. On lowering the barrier, the vacuum state decays to the true ground state, and these modes become occupied. In the next section, we will return to the Kerr problem, where we will find an analogous result.

IV. VACUUM DECAY IN KERR
We now return to the original Kerr problem, set up in Sec. II. For the electric analogue of Sec. III, the boundary conditions were changed from those of the initial to the final states by instantaneously lowering the infinite barrier at t = 0. For the Kerr case, the emergence of the effective potential in Fig. 1 will depend on the way in which the Kerr black hole is formed. One might envisage, e.g., realistic scenarios where the initial condition corresponds to some asymmetrically collapsing mass distribution that forms the Kerr black hole or, alternatively, one might consider an initial Schwarzschild black hole that is spun up by accreting matter from an orbiting body. For our purposes, however, we can be far more pragmatic, since we do not expect the state to depend strongly on the initial conditions at sufficiently late times. What we do expect is that states associated with the upper frequency continuum (see Fig. 1) are initially empty. This initial state would resemble the 'Unruh vacuum' of massless fermions [11,30]. Note that this state normally radiates a thermal flux at the Hawking temperature T H , but because we mainly consider modes ω µ, and usually T H µ, we shall neglect the Hawking temperature. It would, however, be interesting to investigate this in more detail and, in particular, whether thermal effects could somehow accelerate the decay via stimulated emission.
In order to construct the initial state, we would like to build normalizable modes associated with the relevant QNM frequencies. This is a delicate and ambiguous procedure, since it is a well-known issue of QNMs that they are spatially growing and hence not normalizable [31][32][33]. Here, we circumvent the problem by constructing approximate eigen-states associated with the real part ω R n , and of finite support, in close analogy to the preceding section. An alternative way of dealing with these modes and constructing a pseudo-orthonormal basis would be to define them as approximations to the exact eigenstates by means of complex distributions [34,35]. Importantly, we assume that the in modes constructed this way form a complete basis, at least in the relevant range of frequencies ω min < ω < min(µ, mΩ H ), exactly as in equation (45) for the electric case.
To explicitly construct the in basis, first notice that the modes described in Sec. II C are a set of trapped modes in the region r 2 < r < r 3 that can tunnel to the continuum set living in r < r 1 (see Fig. 1). By analogy with the electric case, the tunneling amplitude plays a role similar to the barrier strength λ. There, we can construct approximate in modes ϕ in n by assuming that they have support only in the region r 2 < r < r 3 . We select the discrete set satisfying the Bohr-Sommerfeld quantization condition (32a) and normalize them using the integral of the scalar product (16) between the turning points r 2 and r 3 . In this way, the radial part of the initial and final states are given in the region r 2 < r < r 3 by and (61) In the above definition of in modes (and of k n ), the frequencies are the real parts ω R n of the QNM frequencies given in equation (32). T cl is defined in equation (33).
We now show how the initial vacuum, where the states |ϕ in n are unoccupied, decays spontaneously by filling them in. Proceeding in exact analogy to the electric toy model in Sec. III, we can define the instantaneous occupation number of the state |ϕ in n as As before, this occupation number is governed by the analytic structure of the overlap ϕ in n |ϕ out ω . Since we work in the WKB regime, we shall estimate this overlap by exploiting the fact that the relevant exponentials oscillate rapidly. By means of the separation in Eq. (16), we have where we have neglected the highly-oscillatory and therefore subdominant contributions proportional to e +i dr [k(r )+kn(r )] and e −i dr [k(r )+kn(r )] . The remaining integrals are also highly oscillatory except when ω − ω R n ∼ 0, such that the overlap has dominant support for ω ∼ ω R n . We can therefore expand the phase to leading order in ω − ω R n . Doing so, we have where is the time lapse (in the Boyer-Lindquist coordinate system) necessary to propagate from a reference point (here r 2 ) to r. In this way, we obtain By estimating the amplitudes R 0 1,2 at leading order in ω − ω R n (which is justified in the WKB approximation), we recognize that r 2 +a 2 ∆ |R 0 n,1 | 2 + |R 0 n,2 | 2 is nothing else than 1/v g (r). This is a direct consequence of the normalization of the in modes ϕ in n , and the correspondence between the norm and the conserved current (see appendix B and specifically equation (88)). This leads to Now, to evaluate the integral giving the evolution of the occupation number, i.e. equation (62), we need to deform the contour in the lower-half complex ω-plane (t > 0). Doing so, we pick up all the poles of |A ω | 2 , in perfect analogy with the electric case. These poles are nothing else than the QNM frequencies, and, as in Sec. III, only the QNM corresponding to the same index n contributes significantly. For the Kerr case, the residue of |A ω | 2 at its pole is given by and hence, after making use of Eq. (67), we find By means of equation (62), we then arrive at the simple result, analogous to equation (59), that the occupation number at the time t is given by This shows what was announced earlier: the initial vacuum slowly decays by filling in the discrete set of energy levels corresponding to the QNM satisfying ω R < min(µ, mΩ H ). Since these modes are trapped, this process is not accompanied by any radiation emitted to infinity. The decay process we describe does not arise on top of the Unruh-Starobinski radiation, but it replaces it. In fact, we can decompose the process in three steps. First, following initial transients, the field is in the non-rotating vacuum. Then, as long as the exponentials e −iwnt are approximately linear in t, there is a steady flux of radiation emitted from the ergoregion: the Unruh-Starobinski radiation. This radiation fills in the discrete set of trapped modes, until every level is filled and the process stops. One is left with a Dirac sea, as described in Ref. [10], which now possesses a non-zero expectation value of the angular momentum. To see this, we can evaluate the total angular momentum inside the trapped region. Assuming that the in modes used earlier are complete (at least approximately), we have Notice that the above sum means that we sum over the (discrete) set of QNMs satisfying the condition ω R < min(µ, mΩ H ). In particular, we sum over the quantum numbers j and P, thereby taking into account the contributions of all values of the angular momentum and the spin. Note also that the stress-energy tensor is renormalized at early times, which amounts to normal order with respect to theâ in 's. We then conclude that the true vacuum possesses a mean angular momentum given by Of course, if one computes the total angular momentum on the whole space, one finds that it is conserved and vanishes. This means that the spontaneous acquisition of an angular momentum vacuum expectation value is accompanied by a net flux in the horizon of the black hole, thereby slowing it down, as in bosonic superradiant processes. Although the total angular momentum (73) can be a significant fraction of the black hole angular momentum (see section 3.1.2 in Ref. [10]) depending on the parameter M µ, after a finite time, only a fraction of states have decayed. The extraction of angular from the black hole is a mainly a collective effect, since each state carries a small angular momentum. Hence, a proper estimate of the angular momentum transfer after a time t requires to work with equation (72b) rather than (73).
To finish this section, we point out that we have carried out this analysis in the regime M µ 1, wherein the time necessary for the decay is exponentially long (see equation (32b)). However, as in the electric case, the same effect exists for any value of the parameter M µ. The modes that contribute to the effect are those satisfying |ω I | |ω R | (which justifies the construction of the initial vacuum) and the condition ω R < min(µ, mΩ H ). In the opposite regime, when M µ 1, the life-time is also large, increasing as a power law in M µ (see the works of Refs. [36,37]). Based on the results of Ref. [19], the effect is expected to be the quickest in the regime M µ ∼ 1, in close similarity to the bosonic black hole bomb instability.

V. CONCLUSION
In this work, we have shown that massive fermions around a rotating black hole are subject to a vacuum instability. After the black hole is formed, the initial state is similar to the Unruh vacuum and has a zero angular momentum expectation value. As time progresses, this state decays into the true vacuum, which corresponds to the Kerr-Fermi sea identified by Hartman et al. [10]. We show that this decay is entirely governed by the set of QNMs satisfying the superradiance condition 0 < ω < mΩ H (although there is no amplification) and lying below the mass gap ω < µ. These modes correspond to a discrete set of bound states, which orbit well outside the ergoregion. The states are progressively filled with particles, with an exponential decay probability characterized by the imaginary part of the QNM frequency. Once every state is filled with one particle of each spin, the process stops. As usual, the vacuum decay is accompanied by a spontaneous symmetry breaking. Here, the true vacuum acquires a non-zero expectation value of angular momentum.
When the effect is applied to particles of the Standard Model, the time scale of the decay is extremely large. Using a conservative estimate for the mass of the light neutrino, and a solar mass black hole, one has M µ ∼ 10 9 , and therefore, the imaginary parts of the QNM frequencies are suppressed by a factor e −10 9 , making the decay time larger than the age of the Universe. This suggests that the state of neutrino fields around astrophysical black holes is likely to be the metastable, non-rotating one. However, this statement deserves a more precise analysis. Indeed, since the density of states is also very large in the limit M µ 1, by statistical effects, several states will be filled after a reasonable amount of time. It would be interesting to study the structure of the state as a function of time, when only a fraction of the bound states has been filled. Moreover, since the decay is accompanied by a decrease of the black hole angular momentum, this effect could be used to constrain the possible range of mass for light fermions, in the same spirit as what was done for light Bosons using the black hole bomb instability [25,[38][39][40]. One can also imagine various alternative scenarios where the quantity M µ becomes of order 1, and the vacuum instability becomes relevant. This could happen around a primordial black hole with a sufficiently small mass or in the presence of ultra-light fermions [25]. In this context, it is interesting to remark that the mass of the lightest SM neutrino does not have a lower experimental or observational bound.
In this work, we have concentrated on Kerr black holes. A similar effect is anticipated to be present also for charged black holes, as described by the Reissner-Nordström spacetime. In this case, it is anticipated that the true vacuum will acquire a non-zero charge instead of angular momentum, as in the electric model of Sec. III. This vacuum instability could also be of interest for holography, where similar effects have been studied in the context of holographic superconductivity [41,42], or for the Kerr-CFT correspondence [43].

A. DIRAC EQUATION IN CURVED BACKGROUND IN THE WEYL REPRESENTATION
For completeness, we provide in this appendix all the conventions for the Dirac equation in Kerr, as used in this work. We recall that we have chosen the same conventions as in Ref. [18] and we refer to the literature on the subject for more details [13]. To define the Dirac equation (4) fully, we need to build the spinor covariant derivative D ν and the curved spacetime Dirac matrices γ µ . To do so, we use the tetrad formalism. We first write the metric as g µν e µ a e ν b = η ab , where η ab is the Minkowski metric (see equation (11) in Ref. [18] for the explicit expressions in Kerr). We use lower-case Greek characters to label coordinates in the curved space and lower-case Roman characters to label those in the tangent space.
Using the tetrad e µ a , we construct the curved spacetime Dirac matrices γ µ from their flat spacetime counterpartsγ µ : Notice that γ µ and e µ a depend on the point in the spacetime manifold, whereas theγ a do not. In the Weyl/Chiral representation, the flat Dirac matrices are given by where σ j (j = 1, 2, 3) are the Pauli matrices The spinor covariant derivative can be written where is the spin connection in which λ abc = e µ a (∂ ν e bµ − ∂ µ e bν )e ν c .

B. EIGENMODE NORMALIZATION AND CONSERVED CURRENT
In this appendix, we show that the normalization of the conserved norm (17) is equivalent to that of the current (18). We shall present a very general argument.
Let's consider a field φ : R 2 → C n in 1+1 dimensions (t, x). Associated with the dynamics of the field is a conservation law of the form where the norm density D and the current J are quadratic forms. In the core of the paper, C n is the space of the 4 components of a Dirac spinor, and the radial coordinate r plays the role of x, once one has integrated over the angles. The conservation law above corresponds to equation (19). One can already draw two important consequences of equation (79): 2. If φ(t, x) is a stationary mode, meaning φ(t, x) = ϕ ω (x)e −iωt for some real frequency ω, then the current J[ϕ ω (x)] is conserved along x.
Now, we assume that any solution of the equation of motion that is localized in space can be written as a dense sum of eigen-frequency modes ϕ ω . For simplicity, we assume that the background is homogeneous, and hence eigenmodes are plane waves: where C ω is a vector in the internal space C n . However, the same result applies directly to WKB modes, simply by computing the various integrals using appropriate stationary phase approximations. We want to normalize the eigenbasis in the Dirac sense, that is 2 with N ω = ±1. The normalization conditions generally fix the amplitude C ω in equation (81). Our general claim is that this normalization is equivalent to that of the conserved current. More specifically, we will show that which is what was used in the main text. For this, we fix a frequency ω and consider a tight wave packet whose Fourier spectrum is peaked around ω, i.e.
where the function f is highly peaked about 0. The usual trick to describe tight wave packets is to consider only the vicinity of ω. For this we define ω = ω + , and expand phases to first order in : wheref is the Fourier transform of f . We recover here the standard result: φ essentially oscillates as e −iωt+ikωx , and since f is sharply peaked,f is a long envelope (i.e. containing many wavelengths defined by k ω ), moving at the group velocity v g = (∂ ω k) −1 . Notice that if one works with WKB modes, the argument of f is replaced by t − dx/v g (x) since the group velocity is defined locally. In the case of a field obeying a wave equation in a curved spacetime, this exactly means that tight wave packets follow null geodesics.
We will now use this expression in the conservation law (79). Noticing that it follows that Hence, we obtain the following identity, valid for any ω: This is the natural translation of the conservation law for a definite frequency: the current is simply the norm density times the moving velocity of wave packets v g . All that is left to do now is to use it in the normalization condition (82): Using the identify (88), we have leading directly to the result in equation (83).

C. CONNECTION FORMULA FOR THE DIRAC EQUATION
In this appendix, we briefly review the reasoning that leads to the boundary conditions used at the turning points to build the modes in Kerr (see equations (29) and (28)). The idea is rather standard from WKB methods [22,23] and amounts to determining how to connect two oscillating solutions on one side of the turning point to the growing and decaying exponentials on the other side. To see this, we consider a general solution of the left of the outer turning point r 3 . For r 2 < r < r 3 , it is a sum of oscillating exponentials obtained from equations (21) and (24): We see that this combination is the same as the usual case of a scalar function. The usual trick is to identify it with the asymptotic behaviour of an Airy function, which furnishes us with a regular solution across the turning point. Using the following asymptotics of the Airy function: and In order to obtain the Kerr modes of equations (29) and (28), one needs to apply the connection formula three times: at r 1 , r 2 and r 3 . Notice that we have to face the standard problem of the irreversibility of the connection formula [22]. However, as we saw, everything happens exactly as in the scalar case, exposed in detail in most treatise, and one can circumvent the irreversibility problem in the same way. As a final remark, we point out that what we used in this work is a combination of (three times) the single turning point formula. It is possible to obtain more general three-turning-point connection formulae [22]. Doing so is significantly more involved and barely improves the accuracy of our results.

D. FURTHER DETAILS OF THE ELECTRIC TOY MODEL
In this appendix, we include further details for the electric toy model in Sec. III, including the complete mode functions, their overlaps and the decomposition of the field operator.
We first compute the set of four linearly independent solutions. For x < −D, we have V (x) = V 0 , and µ(x) = 0. This gives four solutions Similarly, for D < x < L, we have V (x) = 0 and Note that the subscripts L and R, indicate the solutions to the left and right of the barrier and not chiralities. Lastly, inside the barrier, we have V (x) = 0 and µ(x) = µ 0 . We assume that |ω| < µ 0 . (We will ultimately take the limit µ 0 → ∞.) In this case, the modes are either growing or decaying but not oscillating. Defining κ = µ 2 0 − ω 2 > 0, we get Using these and assuming continuity of the field for the potential and mass of Fig. 2, we can extract from them the boundary conditions (39) and (40). The outgoing modes are given by a single continuum, composed of For completeness, we also give the various mode functions overlaps on the interval [0, L]: The dominant support of the overlaps occurs at the zeros of the arguments of the sinc functions, setting ω ∼ ω n .