Hunting for fermionic instabilities in charged AdS black holes

Fermions scattering on a black hole background cannot develop an instability sourced by superradiance. However, in a global (or planar) AdS4-Reissner-Nordström background fermions can violate the AdS2 fermionic mass stability bound as measured by a near horizon observer at zero temperature. This suggests that AdS-Reissner-Nordström black holes might still be unstable to Dirac perturbations. Motivated by this observation we search for linear mode instabilities of Dirac fields in these backgrounds but find none. This is in contrast with the scalar field case, where a violation of the near-horizon Breitenlöhner-Freedman stability bound in the AdS-Reissner-Nordström background triggers the already known scalar condensation near-horizon linear instability (in the planar limit this is Gubser’s instability that initiated the holographic superconductor programme). We consider both the standard and alternative AdS/CFT quantizations (that preserve the conformal invariance of AdS). These are reflective boundary conditions that have vanishing energy flux at the asymptotic boundary.


Introduction
It is a well known fact that bosonic waves impinging charged or rotating black holes can be amplified via superradiant scattering (see e.g. the review [1] and references therein). It follows that black holes perturbed by bosonic fields in the presence of a gravitational potential well -provided by, for example, an asymptotic anti-de Sitter (AdS) potential, a physical cavity at finite radius or by the mass of the field -can develop a superradiant instability. However, this is not the only instability that can be present in such systems. Indeed, the family of superradiant black holes always has a configuration with zero temperature. Typically, such extreme black holes have a near-horizon geometry that JHEP04(2020)196 is the direct product (or a fibration) of a base space (e.g. a sphere) and an AdS 2 space [2]. These extreme (and near-extreme) black holes, when confined in a gravitational wall, are unstable if the 'effective mass of the perturbation (as seen by an AdS 2 observer) violates the 2-dimensional Breitenlöhner-Freedman (BF) bound for stability [3][4][5][6] (even though the asymptotic AdS 4 BF bound is obeyed). Recall that perturbations with a mass below this bound are normalizable (i.e. they have finite conserved energy) but their energy is negative and, consequently, they can trigger an instability. The superradiant and near-horizon instabilities have a different physical nature. But, in a non-extremal black hole they are usually entangled. However, if L is the typical dimension of the gravitational well (e.g. the radius of the cavity or the AdS radius), they disentangle for small dimensionless horizon radius r + /L since the near-horizon instability is suppressed for r + /L 1, while the superradiant instability is still present [7][8][9][10][11][12][13].
Not less well known is the fact that fermionic waves, unlike bosons, cannot suffer from superradiant scattering amplification [1,[14][15][16][17][18]. Perhaps less familiar is the fact that fermions, like bosons [7][8][9][10], can also violate the 2-dimensional stability bound of the near-horizon geometry of near-extremal black holes, if the fermion charge is high enough (for fermions the stability bound is lower than the BF bound). In particular, this can happen in a Reissner-Nordström (RN) black hole in an asymptotically planar AdS background [9,[19][20][21][22][23][24][25] (see in particular section V of [9] and sections 4.4 and 4.5 of [24]) or in an asymptotically global AdS geometry (section 4 below). This suggests that such a RN black hole might be unstable to condensation of a fermionic cloud around the horizon. Moreover, if the features from the bosonic field extend to the fermion case, then this might be a linear instability.
Motivated by these considerations, in the present manuscript we will search for linear mode instabilities of Dirac fields in a global AdS 4 RN black hole. We will not find any such instabilities. The absence of linear mode instabilities in the global AdS RN background is consistent with the fact that they are also not present in the planar AdS limit, r + /L → ∞, as found previously in [9,[19][20][21][22][23][24][25]. In view of these 'no-go' findings, in the conclusion remarks of section 5, we will argue that the absence of linear mode instabilities in the AdS RN background still leaves room (but not necessarily) for the following possibility: for a large number of fermions and in the semiclassical limit, the violation of the AdS 2 BF bound might signal a non-linear instability of the system.
In this manuscript, we also take the opportunity to improve the understanding of the near-horizon condensation instability of scalar fields. In particular, following a similar analysis for scalar fields in a Minkowsky cavity [13], we will explicitly show that the BF bound criterion for instability in AdS-RN is quantitatively sharp (section 2.4). This scalar field analysis will fit smoothly in our presentation since reviewing its details will also allow to make a direct comparison with the Dirac field case and pinpoint major differences between them. We further take the opportunity to show that the unstable modes belong to a family of near-extremal modes that are connected to the normal modes of AdS when the dimensionless horizon radius shrinks to zero. (This is an interesting observation because in de Sitter black holes the 'normal mode de Sitter family' is distinct from the 'near-extremal family of modes', where we are using the nomenclature of [26][27][28]).

JHEP04(2020)196
Necessarily, we will also clarify some misleading analyses and interpretations that were presented in previous literature. Namely, the authors of [29] 1 missed that the vanishing energy flux boundary conditions at the asymptotic boundary of AdS -that they propose to be "novel" -are nothing else but the AdS/CFT correspondence no-source boundary conditions -often denoted as 'standard ' or, if allowed, 'alternative' quantizations (discussed, for s = 1/2, originally in [25,[31][32][33][34] and specially in [3][4][5][6]35]; this discussion applies both to bosonic and fermionic fields). The vanishing energy flux boundary condition is fundamental to guarantee that the energy is conserved. If and only if this is the case, the Schrödinger operator that describes the (bosonic or fermionic) wave equation in AdS is Hermitian. 2 It follows that finiteness of energy then boils down to simply require that the wavefunctions of the system have finite norm in the usual quantum mechanical sense, i.e. that the solutions are normalizable (square integrable). In these conditions, the Schrödinger operator of the system is self-adjoint (the associated matrix is Hermitian) and we have a well-posed initial value problem (after imposing regularity at the inner boundary). That is, the dynamical evolution of the system is deterministic.
If the energy is positive, the evolution is stable (this happens if we are above the stability mass bound [6,35]); on the other hand, if the energy is negative we should have a dynamical evolution that develops an instability (this is the case if the mass of the perturbation is below the stability mass bound [6,35]). 3 The aforementioned homogeneous Dirichlet (standard) or Neumann (alternative) AdS/CFT boundary conditions are special in the sense that, by construction, they yield zero energy flux normalizable modes that preserve the conformal symmetry group of AdS (and thus do not deform the boundary conformal field theory). The zero-flux boundary conditions of [29,30] are nothing but these single-trace homogeneous boundary conditions [25,[31][32][33][34][35][36][37]. Besides these, the AdS/CFT correspondence literature identifies, for a certain range of the boson/fermion masses, other normalizable modes (finite conserved energy modes and thus with vanishing energy flux). This is the case of the inhomogeneous Dirichlet and Neumann boundary conditions but also of the often denoted mixed, Robin or multi-trace boundary conditions (see [6] for bosons and [35][36][37] for fermions). These zero-flux boundary conditions break the AdS conformal symmetry while still preserving its Poincaré symmetry subgroup.
The plan of this manuscript is the following. In sections 2.1-2.2, we will review the Dirac equation in a AdS-RN background and we will do the necessary field redefinitions in the (physical) Dirac field that allow it to separate and even decouple. Then, in section 2.3 we will analyse in detail the AdS/CFT standard and alternative quantizations of a Dirac field. In particular, we will check that the requirement that the source vanishes implies that the energy flux at the conformal boundary also vanishes. In this sense, these boundary conditions can also be denoted as reflective boundary conditions. For a massive Dirac field these no-source boundary conditions translate into homogeneous Dirichlet or Neumann boundary conditions in the auxiliary decoupled Dirac radial fields. However, for a massless JHEP04(2020)196 fermion (a Weyl field) these no-source boundary conditions -which are homogeneous Dirichlet or Neumann conditions on appropriate projections of the original physical Dirac field -translate into mixed (Robin) boundary conditions for the auxiliary decoupled Dirac radial fields. The misleading focus on the boundary condition for the auxiliary fields (and associated consequences) occurs recurrently. It is the case in [29] and it is similar to the one taken on the boundary conditions of the Regge-Wheeler-Zerilli master fields (aka Kodama-Ishibashi fields) in the case of gravitational perturbations of AdS black holes (as discussed in [38,39]). 4 In this context it is also important to clarify that in [41,42] massless Dirac quasinormal modes in Schwarschild-AdS were computed imposing Dirichlet boundary conditions in the auxiliary decoupled fields. These boundary conditions do not have vanishing energy flux at the asymptotic boundary, the energy of the system is thus not conserved, and it is not known what deformation they produce. Finally, in section 4.4 we will describe our strategy to search (unsuccessfully) for linear mode instabilities (eventually sourced by the 2-dimensional stability bound violation) of Dirac fields in global AdS 4 -RN black holes. We consider both the standard and alternative quantizations and we will highlight the differences between the scalar and fermion systems. For a reader interested in a future detailed analysis of the frequency spectra, we also compute (analytically) the normal modes of massive and massless Dirac field in global AdS.
2 Global AdS Reissner-Nordström black hole and the Dirac equation

AdS-RN black holes and an orthogonal vierbein
The gravitational g µν and Maxwell A µ fields of the AdS-RN BH are described by 5 A µ dx µ = A(r) dt, A(r) = − Q r + C; (2.1) where dΩ 2 2 is the line element of a unit radius 2-sphere, M and Q are the mass and charge parameters. We will find convenient to replace M and Q by the event horizon radius r + (where f (r + ) = 0) and chemical potential µ. The relation between these two pairs of parameters is In (2.1), C is an arbitrary integration constant and fixing it amounts to choosing a particular gauge. One common gauge choice is C = 0 where one has A| ∞ = 0 and the chemical potential is µ = −A| r + = Q/r + (we will typically use this one when presenting our results). Another gauge that is also commonly used is C = µ whereby A| r + = 0 and A| ∞ = µ. 4 For a discussion that AdS/CFT no-source boundary conditions for bosonic fields yield ('reflective') solutions with vanishing energy (and momentum) flux at the conformal AdS boundary see appendix A of [40]. 5 This is a solution of the Lagrangian L = √ −g R − 2Λ − 1 2 F 2 with F = dA. Note that if we rescale A → κA then the charge q of the perturbation field (to be discussed in later sections) rescales as q → κ −1 q so that qA, and thus the gauge covariant derivative, remain invariant.
The components of any tensor in the coordinate basis {dx µ } can be obtained from the components on the tetrad basis using the projectors e (a) and e (a) . For example, σ T (a)(b) . This example also illustrates that we can have mixed-index tensors with mixed components in the coordinate and tetrad bases.
The spin connection of non-coordinate based differential geometry can be introduced in terms of the affine connection Γ β µν of coordinate based differential geometry as with λ (a)(b)(c) = e α (a) ∂ β e (b)α − ∂ α e (b)β e β (c) , which allows to compute the spin connections without the use of the affine connection.
For a multi-index tensor with tetrad and coordinate indices the mixed-index covariant derivative is defined as Onwards, we take the affine connection Γ ν µα to be given by the Christoffel symbols that covariantly conserve the metric, ∇ α g µν = 0. It follows that the spin coefficients ω µ(a) (b) defined in (2.6) are such that the vierbein is also covariantly conserved, ∇ α e ν . Further note that the spin coefficient is anti-symmetric in the tetrad pair of indices, ω µ(a)(b) = −ω µ(b)(a) .
To discuss the Dirac equation one necessarily needs to introduce the (coordinate independent) Dirac gamma matrices γ (a) . Let σ i be the Pauli matrices and I n the n × n JHEP04(2020)196 identity matrix. We choose to work with the Weyl (chiral) spinor representation of the 4-dimensional Clifford algebra: (2.8) which indeed obeys the anti-commutation relations of the Clifford algebra: where {A, B} = AB + BA is the usual anti-commutator, as well as the relations (γ (0) ) 2 = −I 4 and (γ (i) ) 2 = I 4 for i = 1, 2, 3. Let us also introduce the pseudoscalar which obeys the relations (γ (5) ) 2 = 1 and {γ (5) , γ (a) } = 0. 6 The components of the (coordinate dependent) Dirac gamma matrices in the coordinate basis can be obtained from the tetrad basis components (2.8) using 11) and they obey the covariant Clifford algebra {γ µ , γ ν } = 2g µν I 4 .

Dirac equation in the AdS-RN background
One of our main goals will be to compare scalar (spin-0) and fermionic (spin 1/2) perturbations on the AdS-RN background. Therefore, below we briefly review the equations that these fields have to obey in a curved background, e.g. in the AdS-RN spacetime (2.1). For more detailed discussions see [43][44][45][46][47].
To start, consider spin-s fields in Minkowski spacetime. The spin of a field can be identified looking into how the field transforms under a Lorentz transformation, x µ →x µ = Λ µ ν x ν . Let M αβ = −M βα be the generators of Lorentz transformations (i.e. a basis of six 4 × 4 antisymmetric matrices obeying the Lorentz Lie algebra). A finite Lorentz transformation is described by Λ = exp 1 2 Ω αβ M αβ where Ω αβ = −Ω βα are six parameters describing the particular transformation Λ (boost, rotations) of the Lorentz group SO (3,1).

JHEP04(2020)196
In Minkowski spacetime, under Lorentz transformations the derivative of a spin-s field transforms as ∂ µ Ψ → ∂ µΨ (x) = Λ ν µ S[Λ]∂ ν Ψ(Λ −1 x). If we want to couple the spin-s field to a curved background, while preserving general covariance, one needs to promote the partial derivative ∂ µ to a covariant derivative D µ . This promotion is chosen such that any function of Ψ and D µ Ψ that is a scalar under Lorentz transformations in Minkowski spacetime remains a scalar -under general coordinate transformations and local changes in the vierbein -in the curved background. This is the case if, under an arbitrary Lorentz transformation, the covariant derivative D µ still transforms as a derivative of a spin-s field: It follows that the covariant derivative of a spin-s field that preserves Lorentz invariance in a curved background is where we took the opportunity to allow the spin-s field to have a charge q (that couples to the Maxwell background field A µ ), and Γ µ is a covariant spin connection which depends on the spin of the field it acts on since the generators S (a)(b) depend on whether we are looking into, for example, the scalar or spinor representation of the Lorentz group. In more detail, for a scalar (spin-0) field Ψ ≡ Φ the Lorentz group generator is simply S (a)(b) = 0. Therefore S[Λ] = 1 and the scalar covariant derivative (2.13) is simply D µ = ∂ µ − iqA µ . The action for a massive charged complex scalar field is given by where * stands for complex conjugation. The factor of √ −g is introduced to ensure that the Lagrangian L Φ is a scalar density and thus the action S Φ = M d 4 xL Φ is a scalar. Varying this action w.r.t. Φ * one gets the Klein-Gordon equation for the scalar field and similarly for Φ * . On the other hand, for a (spin-1 2 ) Dirac 4-spinor field 8 Ψ ≡ ψ, out of the gamma matrices γ (a) (2.8) that satisfy the covariant Clifford algebra (2.9) one can build the commutator

JHEP04(2020)196
The action that is Lorentz invariant and describes the coupling of a spin-1 2 fermion field ψ to a curved background M is where we have introduced the Dirac adjointψ = ψ † γ (0) with ψ † = (ψ * ) T being the Hermitian adjoint of the multi-component field ψ. One needs to work with the Dirac adjoint because the Fermi bilinearsψψ andψγ µ ψ transform covariantly (as a scalar and as a vector, respectively) under the Lorentz group (while the Hermitian partner objects do not).
Varying the action S D w.r.t.ψ and ψ, respectively, one gets the Dirac equations To find solutions of (2.18) it is advantageous to write the Dirac 4-spinor ψ in terms of the left-handed and right-handed 2-spinors Ψ − and Ψ + , respectively, as The chiral 2-spinors Ψ ± emerge naturally when we note that the pseudoscalar γ (5) defined in (2.10) obeys (γ (5) ) 2 = 1. Therefore we can introduce the Lorentz invariant projection operators P ± that project the Dirac 4-spinor ψ into the chiral spinors: and such that P 2 ± = P ± and P + P − = 0. 10 Moreover, the task of finding solutions of the Dirac equation in the AdS-RN black hole gets considerably simplified by the fact that under the separation anstaz [48,49]: the Dirac equations (2.18) reduce to a set of equations where the radial and angular functions of the fermion field are decoupled. This separation ansatz exploits the fact that ∂ t and ∂ φ are Killing vectors of the background AdS-RN solution. This allows to do a Fourier decomposition in these directions which introduces the frequency ω and azimuthal angular momentum m φ of the fermionic wavefunction. 11

JHEP04(2020)196
Concretely, the radial functions R 1 (r), R 2 (r) obey the coupled system of first order ODEs where λ is a separation constant, while the angular functions S 1 (θ), S 2 (θ) satisfy the coupled system of first order ODEs (2.23) Furthermore, the coupled pair of first order radial equations (2.22) can be decoupled in a pair of second order ODEs, one for R 1 (r) and the other for R 2 (r). For that we solve the first (second) equation in (2.22) w.r.t. R 2 (R 1 ) and replace it in the second (first) equation. We end up with two decoupled second order ODEs for R 1 and R 2 , where * denotes complex conjugation and we have defined Of course, we are only interested in solutions of (2.24) that also solve the original first order system (2.22). The requirement that (2.22) is solved imposes extra constraints on solutions of (2.24). This is best illustrated if we consider the Taylor expansion about the boundaries of the integration domain: the ODE pair (2.24) has four integration constants about each boundary but only two of them are independent when we further require that the solution solves the two first order ODEs (2.22); see discussion of (2.27) below. Similarly, the coupled pair of first order ODEs for S 1,2 (θ) can be written as a decoupled set of two second order ODEs for S 1 (θ) and S 2 (θ). They are hypergeometric equations and S 1,2 (θ) are the spin-1 2 weighted spherical harmonics. Regularity at θ = 0 and θ = π quantizes the angular separation constant as ( is a harmonic number related to the number of zeros of the wavefunction) with the azimuthal number being constrained as m φ ≤ .

JHEP04(2020)196
Unfortunately, the radial ODEs cannot be solved analytically. 12 We can however do a Frobenius analysis about the asymptotic boundary r → ∞ to find the asymptotic behaviours of R 1 (r) and R 2 (r). One finds that (for m = 0, 1 2 ) 13 where we used (−gf ) − 1 4 | r→∞ ∼ L 1/2 r −3/2 and, anticipating the AdS/CFT discussion below, we have introduced the conformal dimensions As expected for a coupled system (2.22) of two first order ODEs, there are two independent arbitrary constants (α 1 , β 1 ) in the asymptotic decay (2.27), that is to say, the decays of R 2 are fixed by the equations of motion as a function of (α 1 , β 1 ). The dots in (2.27) represent subleading terms that depend only on α 1 (in the ∆ − contribution) or β 1 (in the ∆ + terms). Before proceeding, one unavoidably needs to discuss the range of Dirac fermion masses that allow for normalizable solutions, i.e with conserved finite energy. We also have to distinguish the positive energy solutions (which are stable) from those negative energy states (which should trigger an instability). It was proven in section II/appendix B of [35] (see also [6,36,37]) that the fermionic bound for stability (in any dimension) is given by with the lower bound being the solution for which ∆ + = ∆ − in (2.28). 14 To understand this bound it is useful to rewrite the radial Dirac equation (2.24) as a Schrödinger equation [6,35]. Without further conditions, the associated Schrödinger operator is not self-adjoint (hermitian). It becomes self-adjoint if and only if we impose as boundary condition that the energy-momentum flux at the asymptotic AdS boundary vanishes. That 12 For global AdS, i.e. M = 0 = Q these ODEs are hypergeometric equations and can be solved analytically: see section 4.2. 13 For m = 1/2 one of the two independent solutions decays asymptotically as a power law in r and the other as a power law multiplied by a log r. For this reason (since a similar logarithmic solution appears in the scalar field case when m 2 = m 2 BF ), this case is often called the BF solution of the Dirac system. We do not discuss further this special case (see [19,35] for more details). It is however important to emphasize that for the scalar field, m 2 = m 2 BF corresponds to ∆+ = ∆− and is thus also the bound for stability while in the Dirac case, the mass stability bound is (2.29) not the BF mass m = 1/2.
14 So, for m 2 ≥ 0, ∆± are real; otherwise they are complex numbers. Note that for a scalar field the configuration ∆+ = ∆− corresponds to the BF bound where one of the independent solutions is logarithmic. However, for the Dirac field, the state ∆+ = ∆− is not the BF logarithmic solution (which occurs instead for m = 1/2). If follows that for the scalar field case the BF bound coincides with the bound for stability, m 2 ≥ m 2 BF , but not in the Dirac case. Moreover, in the scalar case, there is a 1-parameter family of boundary conditions that yield stable normalizable solutions for m 2 BF ≤ m 2 < m 2 BF + 1/L 2 and a unique boundary condition that generates stable normalizable solutions for m 2 ≥ m 2 BF + 1/L 2 [3][4][5][6]. However, in the Dirac case, normalizable stable states exist for: 1) a 1-parameter choice of boundary conditions for 0 ≤ m 2 < m 2 BF (with m 2 BF = 1/4), and 2) a unique boundary condition for m 2 ≥ m 2 BF [36]. Further note that, unlike in the scalar case, the Dirac stability mass bound is independent of the dimension of the spacetime.

JHEP04(2020)196
is to say, it becomes Hermitian if and only if the energy is conserved. In these conditions looking for (conserved) finite energy solutions boils down to look for normalizable states in the standard quantum mechanical sense. That is to say, normalizable solutions are those that are square integrable.
For m 2 < 0 there are normalizable solutions but they have negative energy. In a mathematical language, if m 2 < 0, the Schrödinger operator of the Dirac equation is unbounded below and thus it does not allow for a positive self-adjoint extension [6,35]. Alike in any other negative energy Schrödinger states, this signals the existence of an instability. We will explore further this in section 4.1.
On the other hand, if the mass is real, i.e. if it satisfies the bound (2.29), there are stable normalizable Dirac fermion solutions that are selected by a choice of boundary conditions. We will discuss in detail this issue of the boundary conditions in the next section. The upshot is that if mL ≥ 1/2 there is an unique complete set of normalizable modes (and the non-normalizable modes must be fixed by boundary conditions; e.g. nosource/homogeneous boundary conditions that eliminate them) [35]. On the other hand, for 0 ≤ mL < 1/2 there is a non-unique set of normalizable modes and thus a wider band of boundary conditions that yield normalizable solutions (e.g. the no-source/homogeneous Dirichlet or Neumann boundary conditions that we will use later but also more general multi-trace boundary conditions) [35]. Further note that if we take m → −m, we simply trade the role of the ∆ ± contributions while preserving condition (2.29). Therefore onwards we assume, without any loss of generality, m ≥ 0 in our discussion.
For our purposes, but without loss of generality, we will be particularly interested in the lower bound case of (2.29). For this m = 0 case and choosing the gauge A| ∞ = 0, a Frobenius analysis of the first order equations of motion about the asymptotic boundary yields 15 i.e. we can take the two independent integration constants associated to the coupled pair of first order ODEs to be α 1 and β 1 and the equations of motion then fix the decay of R 2 as a function of α 1 and β 1 .

Boundary conditions for the Dirac spinor in AdS-RN
To find the solution of the Dirac spinor field ψ and its Dirac adjointψ in the AdS-RN background we have to solve a system of two equations that are first order, namely (2.22), subject to boundary conditions imposed at the event horizon r = r + and at the asymptotic boundary r → ∞. Before imposing boundary conditions, such a system of two first order differential equations necessarily has two independent constants at the horizon boundary and another two independent constants at the asymptotic boundary (namely, α 1 and β 1 in (2.27)), which can be identified doing a Frobenius analysis at these two boundaries. To have a well posed formulation of the elliptic problem one should impose two boundary conditions that fix two of the independent constants and solve the equations of motion to find the other two. We certainly want the Dirac solutions to be regular at the event horizon: this boundary condition fixes one of the constants. 16 One should then fix one of the asymptotic constants α 1 or β 1 (or a relation between them) with an appropriate boundary condition [35][36][37]. But we certainly cannot fix both asymptotic independent constants: once the first is fixed, the second one must be found by solving the equations of motion in the bulk subject to the two aforementioned boundary conditions. This poses the question: how do we choose a boundary condition at the asymptotic boundary that is physically relevant? We should choose one that conserves the energy and thus yields a self-adoint Schrödinger operator for the system that ensures that we have a well-posed hyperbolic evolution if we let the perturbed system evolve in time. Next, we will review how two boundary conditions with these properties can be identified. They single out in the AdS-CFT context because they are single-trace (no-source) boundary conditions that preserve the conformal symmetry group of AdS (and thus do not deform the boundary conformal field theory) [25,[31][32][33][34][35][36][37] . Dirac spinor fields ψ are intrinsically quantum fields. The dynamics of such fields can be naturally described by a path integral formulation whereby one sums over all possible field configurations in configuration space to get the transition amplitude between two states. In particular, the partition function Z (i.e. the generating functional of correlation functions between operators) can also be naturally computed using the path integral formulation. Schematically one has, S cl [ψ,ψ] is the action evaluated on a solution of the classical equations of motion, that follow from the variation δS = 0 subject to the boundary conditions. As emphasised in [25,[31][32][33][34][35][36][37] this statement that the action must be stationary when evaluated on a classical solution severely constrains the type of boundary conditions we can impose on 16 Alternatively, since we have a ODE system, we could use two boundary conditions to fix the two asymptotic independent constants and solving the equations of motion would yield the behaviour of the Dirac fields at the event horizon. However, this is not a good strategy because in general these solutions would not be smooth at the event horizon.

JHEP04(2020)196
the field ψ. Indeed, if δS = 0 then it is not necessarily true that δ(S + B) = 0 where B is the boundary term describing the desired boundary conditions (i.e. a total derivative term that does not change the equations of motion). That is to say, the physical choice we make for the boundary conditions must be such that δB = 0. In particular, in the context of the AdS/CFT correspondence, this condition fixes the form of the boundary term that must be added to the standard Dirac action (2.17) to have stationary solutions. Vice-versa, this boundary term fixes the boundary field theory.
To determine the boundary term B, one first notes that the "radial" Dirac gamma matrix γ (1) defined in (2.8) satisfies (γ (1) ) 2 = I 4 and γ (1) = γ (1) † . It follows that we can decompose the Dirac spinor as where ψ ± (ψ ± ) are 4-eigenspinors of γ (1) with eigenvalue ±1 (∓1). 17 Using this property, including the associated properties listed in footnote 17, one finds that the terms in the Dirac action (2.17) that contain radial derivatives of the spinor are where we used the fact that D r = f 1/2 ∂ r . It follows that if we vary the Dirac action (2.17) w.r.t. ψ + and ψ − one gets, after integration by parts, where the bulk terms describe a contribution that vanishes when the equations of motion -which are equivalent to (2.18) -are satisfied and δS bdry is a boundary term resulting from integrating by parts the radial derivative terms (2.33) given by As discussed above, to have a well-posed boundary value problem, after requiring that the solution is regular at the event horizon we no longer have the freedom to fix both ψ + and ψ − at the asymptotic boundary (these are the two independent asymptotic constants of our pair of first order ODEs). Instead, we can either fix ψ + at the asymptotic boundary (in which case δψ + = 0) or fix the asymptotic value of ψ − (in which case δψ − = 0 at the boundary).
Suppose we want to fix ψ + at the asymptotic boundary (a similar analysis would apply if we wanted to fix ψ − ). In order to have a well-defined variational problem one should add a boundary term that cancels the contributionψ + δψ − in (2.35). Adding the boundary term [25,[31][32][33][34]

JHEP04(2020)196
produces the desired effect since the total on-shell action becomes which indeed vanishes when δψ + = 0 (and thus δψ + = 0). We can also compute the momentum conjugate to ψ + andψ + by varying S tot w.r.t. ψ + andψ + , respectively, yielding In terms of the functions R 1 (r), R 2 (r) and S 1 (θ), S 2 (θ) introduced in the separation ansatz (2.21), the 4-spinors ψ ± are given by From the asymptotic decays of R 1,2 in (2.27) (valid for m = 0, 1 2 ) or in (2.30) (valid for m = 0) one finds that ψ ± decay as (2.42) where α 1 , β 1 are the free constants introduced in (2.27) or (2.30) and the constants a(α 1 ), b(β 1 ),ã(α 1 ) andb(α 1 ) are fixed as functions of α 1 or β 1 (as described by their argument) by the equations of motion (details are irrelevant for our aim). The asymptotic decays of the Dirac adjointsψ ± follow straightforwardly from (2.40) with the exchange α 1 →ᾱ 1 , β 1 →β 1 , etc. For mL ≥ 1 2 the only normalizable mode (i.e. with finite energy) is ψ + [5,9,19,22,[35][36][37]. In the context of the AdS/CFT correspondence, the leading term of the asymptotic expansion lim r→∞ r ∆ − ψ + = 2α 1 is then identified with the source of the dual operatorŌ which has mass dimension ∆ + . We have a well-posed boundary value problem if we impose smoothness of ψ + at the event horizon and a Dirichlet boundary condition for α 1 at the asymptotic boundary. In particular, if we do not want to deform the boundary field theory we impose the no-source/homogeneous Dirichlet boundary condition: α 1 = 0. We have no freedom left to fix asymptotically ψ − i.e. β − . Instead, β − and thus ψ − | ∞ is determined by solving the Dirac equations subject to the above boundary conditions. The expectation value Ō of the dual operator is given by the conjugate momentum Π + defined in (2.38): Ō ∝ lim r→∞ r −∆ − Π + ∝β 1 .

JHEP04(2020)196
On the other hand for 0 ≤ mL < 1 2 both modes ψ ± are normalizable [5,9,19,22,[35][36][37]. Thus we can still impose the standard quantization where we identify the lim r→∞ r ∆ − ψ + ≡ ψ (0) + as the source of the dual operatorŌ. In particular, the no-source/homogeneous standard boundary condition for all possible masses: But, since for this range of masses both modes are normalisable, 18 we can also impose the so-called alternative quantization; where we identify the lim r→∞ r ∆ + ψ − ≡ ψ (0) − as the source of the dual operator O with mass dimension ∆ − . In particular, if we do not want to deform the boundary field theory we impose the no-source alternative boundary condition: The two quantizations (2.43) and (2.44) yield two distinct boundary conformal field theories [5,9,19,35,37]. For m = 0, note that the Dirichlet boundary condition on ψ + , ψ (0) + = 0, implies the Neumann condition in ψ − (i.e. the next-to-leading order term in the expansion for ψ − vanishes) and vice-versa. This follows straightforwardly from an inspection of (2.41).
We emphasize that the no-source standard and alternative boundary conditions (2.43)-(2.44) that do not deform the boundary theory imply that the energy flux and fermion particle flux vanish at the asymptotic boundary (this is also the case for more elaborated normalizable AdS/CFT boundary conditions [35,37]). In this sense we can regard these as 'reflective' boundary conditions. The Dirac action (2.17) (and (2.37)) is left invariant if we rotate the phase of the Dirac spinor, ψ → e −iα ψ. The Dirac current associated to this symmetry is j µ =ψγ µ ψ and one can check that it is conserved, ∇ µ j µ = 0 after using the first order equations of motion (2.18). This is an internal vector symmetry since ψ ± transform in the same way under this symmetry. This current gives the charge flux or particle number flux of fermions. The associated conserved charge is Q = V dx 3 √ γj µ ξ µ = V dx 3 √ γψ † ψ where V is the volume of a constant t hypersurface, γ ab is the associated induced metric, and ξ = ∂ t is the Killing vector describing time translations. In particular, j r | r→∞ gives the radial flux of particles at the asymptotic spacelike boundary Σ. One can also show that the energy flux across a spacelike boundary is proportional to the Dirac current. The energy flux across the asymptotic boundary Φ ∂t | ∞ is proportional to the particle flux j r | ∞ and is given by 19 (2.45)

JHEP04(2020)196
Inserting the asymptotic decays (2.27) for R 1,2 this yields That is, the energy flux at the asymptotic boundary vanishes if we impose the above discussed no-source Dirichlet boundary conditions α 1 = 0 or, for the alternative quantization, β 1 = 0 which do not deform the boundary conformal field theory. On the other hand, for m = 0, inserting the asymptotic decays (2.30) for Again, this flux vanishes if we impose the standard (2.43) or alternative (2.44) quantizations, β 1 = −i (±λ + ωL) α 1 (and thus β * 1 = i (±λ + ωL) α * 1 ). Here it is important to recall the clarification about AdS/CFT boundary conditions and vanishing flux conditions presented in the Introduction. The standard and alternative boundary conditions that we use have, by construction, zero energy flux at the asymptotic boundary, as reviewed above and originally discussed in [31][32][33][34][35][36][37]. Without noticing, these standard/alternative boundary conditions are also the boundary conditions used in [29,30] where the "generic physical principle of zero energy flux" was used to motivate the boundary conditions originally established in [31][32][33][34][35][36][37] (using precisely the same rationale). But there is a broader family of zero-flux boundary conditions. The AdS/CFT standard and alternative quantizations are a special class of zero-flux boundary conditions that, additionally, preserve the conformal symmetry group of AdS [31][32][33][34][35][36][37]. It is this property that singles them out among other zero-flux boundary conditions that break this conformal symmetry [35][36][37]. Further note that zero-flux boundary conditions are sometimes denoted as 'reflective' boundary conditions in some literature and both set of words encode the familiar idea that 'AdS behaves as a confining box' (under these boundary conditions). 20 It is also important to emphasize that in the AdS/CFT language the standard classification of Dirichlet/Neumann/Robin boundary conditions applies to the physical fields that obey the original differential equation (in the present s = 1/2 case, the Dirac equation). Often this classification does not then translate into the same type of boundary conditions on auxiliary (or even gauge invariant) fields that one might introduce. The classification should focus on the physical fields and not on auxiliary fields (we can fabricate many of these), unlike what is done for s = 1/2 in [29,30]. For example, for a massless Dirac fermion, no-source Dirichlet/Neumann boundary conditions ψ (0) ± = 0 translate into β 1 = −i (±λ + ωL) α 1 not α 1 = 0 or β 1 = 0. Facts like this are often missed: 21 zero-flux boundary conditions that preserve conformal symmetry require ψ (0) ± to vanish not R 1,2 | ∞ . 22 20 Note however that in AdS/CFT there are other sets of boundary conditions that yield a well-defined boundary value problem but do not correspond to zero-flux boundary conditions (e.g. mass deformations describe sourced solutions with important physical interpretations where gauge field(s) have a non-vanishing asymptotic flux). 21 This is e.g. the case in [41,42] where massless Dirac quasinormal modes of Schwarschild-AdS are computed with the Dirichlet boundary condition α1 = 0. This choice of boundary condition is not one of the AdS/CFT zero flux boundary conditions for a massless Dirac field. 22 Further note that there are other boundary conditions (e.g. β1 = −i (±i λ + ωL) α1) that make the flux (2.47) vanish. These should correspond to multi-traced (i.e. mixed or Robin) AdS/CFT boundary conditions [35,37] which deform the boundary theory in a way that might be interesting for other studies.

Near-horizon geometry of the extreme AdS-RN black hole
The near-horizon geometry of the extremal AdS-RN black hole will play an important role in our discussions in sections 3 and 4. Therefore, we review it here. The limiting procedure described below was first presented in [2]. The extremal AdS-RN black hole is given by (2.1) with µ = µ ext given by (2.4). To obtain the near-horizon geometry, it is convenient to work in the gauge A| r + = 0 (C = µ; otherwise we can do a gauge transformation in the end). One first zooms in around the horizon region by making the coordinate transformations: where L AdS 2 is the AdS 2 radius (to be defined below). Now the near-horizon geometry is obtained by taking ε → 0 which yields This geometry is the direct product of AdS 2 ×S 2 and has a Maxwell potential that is linear in the radial direction. Remarkably, in spite of the limiting procedure, it is still a solution of the 4-dimensional Einstein-Maxwell-AdS theory. On the other hand, the AdS 2 metric solves the 2-dimensional Einstein-AdS equations, R µν = −L −2 AdS 2 g µν , if L AdS 2 is identified as a function of the AdS 4 radius L and the horizon radius r + as indicated in the first line of (2.49).

Scalar fields in a AdS-RN background and their instabilities
Scalar fields confined inside the gravitational potential (like the AdS potential or a box in an asymptotically Minkowski background) of a black hole can condense creating nearhorizon linear instabilities [7][8][9][10][11][12] (for planar AdS, this instability triggered the holographic superconductor programme [7][8][9]). Essentially this happens because we can have scalar fields that obey the asymptotically AdS 4 UV Breitenlöhner-Freedman (BF) bound but violate the 2-dimensional BF stability bound associated to the AdS 2 × S 2 near-horizon geometry of the extremal black hole of the system. As we shall discuss in section 4.1, a similar violation of the 2-dimensional stability bound can occur for Dirac fields. In spite of this, as we will find in section 4.4, it turns out that Dirac fields are not linearly unstable to the near-horizon condensation mechanism. Therefore, before we discuss the fermionic case, it is important to revisit the scalar field case. This will allow to: 1) motivate the search of linear instabilities due to Dirac fields done in this manuscript, 2) eventually identify differences between the two spins that could help in understanding the opposite outcomes. We also take the opportunity to demonstrate: i) how remarkably sharp the near-horizon instability bound (3.7) is by comparing it with the numerical solutions of the

JHEP04(2020)196
Klein-Gordon equation, and ii) that the unstable modes are both peaked near the horizon but also connected to the AdS normal modes (that is to say, in the language of [26][27][28] the AdS and near-extremal families of modes coincide and describe the unstable modes).
Using the fact that the AdS-RN background (2.1) is static and spherically symmetric we can consider a separation ansatz for the scalar field (with mass m and charge q) with the Fourier decomposition where Y (θ) are the familiar (spin-0) spherical harmonics -which are regular when the separation constant of the system is quantized as λ = ( + 1), = 0, 1, 2, · · · -and |m φ | ≤ is the azimuthal quantum number. The Klein-Gordon equation yields the following equation for the radial function φ(r): A Taylor expansion around the asymptotic boundary yields the two independent solutions being the conformal dimensions of the field. Such a scalar field in AdS 4 is normalizable as long as its mass obeys the AdS 4 Breitenlöhner and Freedman (BF) bound, m 2 ≥ m 2 BF ≡ − 9 4 1 L 2 [3,4]. Such a scalar field that is stable in the UV region can however be unstable in the IR region. This is best understood if we take the near-horizon limit of (3.2). Concretely, applying the near-horizon coordinate transformation (2.48) together with the near-horizon frequency transformationω → ω ε/L 2 AdS 2 (so that e −iωt → e −i ωτ ) followed by the near-horizon limit ε → 0 yields the radial Klein-Gordon equation in the near-horizon geometry (2.49): This is nothing else but the Klein-Gordon equation for a scalar field around AdS 2 with an electromagnetic potential A τ = α ρ. A Frobenius analysis of (3.4) yields which determines the effective mass of the scalar field from the perspective of a near-horizon observer, Now, a scalar field with mass (3.6)

JHEP04(2020)196
Note that scalar fields can also induce instabilities due to another mechanism that is known as superradiance. Unlike the near-horizon instability -which is suppressed in the limit r + /L → 0; indeed (3.7) goes as q 2 L 2 ≥ L 2 8r 2 + +O(1) -the superradiant instability is present for small r + /L 1 black holes. For example, for m = 0, from the perturbative results of [12] one finds that the superradiant instability in extremal AdS-RN 4 is present for scalar charges 23 Next, we solve the Klein-Gordon equation numerically to confirm that the near-horizon and superradiant instabilities are indeed present and to find how sharp the instability bounds (3.7) and (3.8) are. We present results for scalar masses above the unitarity bound m 2 BF + 1 = −5/4 so asymptotically we impose the Dirichlet boundary condition a = 0; see (3.3). 24 On the other hand at the horizon we require that the solution is regular in the future horizon which discards outgoing modes. To present the results, note that our system has a scaling symmetry [12] which means that the physical dimensionless quantities that are relevant for the problem are (this effectively sets L ≡ 1) r + L , µ; mL, qL, ωL, . First, we are interested in finding the onset of the instabilities namely, the scalar field charge q onset (µ, r + /L; mL, ) above which the system is unstable. This onset occurs when the frequency satisfiesω = ω−qµ = 0. The Klein-Gordon equation (3.2) is then solved as an eigenvalue problem for q = q onset . For concreteness, we fix = 1 (we need m φ ≥ 1 to have an instability). In the left plot of figure 1 we set m = 0 and we plot the dimensionless onset charge q onset L as a function of the dimensionless horizon radius r + /L for different values of the chemical potential µ = µ ext (1 − 10 −x ) that increasingly approaches the extremal value. From top to bottom, the green numerical curves describe chemical potentials with x = 2, 3, 6, 15. We see that as we get closer to extremality these onset curves increasingly approach (for values of r + /L larger than ∼ 0.25) the red dashed curve which describes the near-horizon bound (3.7). This strongly suggests that the instability, for large values of the horizon radius and near extremality, can be understood as due to the violation of the AdS 2 BF and that the associated near-horizon bound (3.7) is sharp (i.e. it is attained at extremality). On the other hand, as pointed out before, the near-horizon red dashed curve diverges as r + /L → 0. However, figure 1 shows that q onset L is finite for small r + /L. Actually, in this regime the numerical onset curves are well described by the superradiant bound (3.8) (blue dashed curve with negative slope). This suggests that for small horizon radius and near extremality the instability has a superradiant nature and the superradiant bound (3.8) is sharp. For finite values of r + /L, i.e. away from the r + /L → 0 and r + /L → ∞ 23 This bound can be obtained from the expression for the frequency obtained in section III.D of [12].
Namely, the onset charge (3.8) of the superradiant instability is obtained by settingω = 0 and µ = µext in equation (55) of [12] and solving for the charge q. For further discussions between the entanglement of the superradiant and near-horizon instabilities and their different nature we ask the reader to see [12] and [13]. 24 For m 2 BF < m 2 < m 2 BF + 1 both modes are normalizable and thus we could also impose the Neumann boundary condition b = 0 (the so-called alternative quantization in the context of AdS/CFT) [3,4]. regions, the superradiant and near-horizon instabilities are entangled. These features are not unique to the massless case. For example, the onset charge plot for a scalar mass of mL = 2 is shown in the right panel of figure 1. Again, as extremality is approached the numerical green curves increasingly approach the near-horizon onset bound (3.7) (in this case we do not show the curve corresponding to the perturbative superradiant curve because it was not computed in [12] but we see that the behaviour of the onset curves for r + /L 1 is similar to the massless case). To compare with what happens in the Dirac field case, it is enlightening to do the following exercise whose results are summarized in figure 2. We pick a particular AdS-RN background with chemical potential µ = 0.99µ ext and horizon radius r + /L = 0.5. We also fix the scalar mass to be mL = 0 and the scalar field harmonic number = 1. Then we solve the Klein-Gordon equation to find the imaginary and real parts of the frequency ωL as a function of the dimensionless scalar field charge qL: these are shown in the left and right panels, respectively, of figure 2. From the left panel we see that, in accordance with the conclusions of figure 1, for small qL the system is stable (since Im(ωL) < 0) but there is a critical charge q L ∼ 1.863 (vertical dotted line) above which the system becomes unstable. Precisely at this critical onset charge one has Re(ω) − qµ = 0 and this quantity is negative (positive) for q < q (q > q ). The inset plot of figure 1 zooms-in the region q < q . In section 4.4 we will find that the partner plot for Dirac fields is substantially different.

JHEP04(2020)196
We take also the opportunity to understand better the frequency spectrum of scalar fields in AdS-RN. For global AdS RN black holes there are two quasinormal mode families [51][52][53]: one whose imaginary part grows negative without bound as the horizon radius r + /L decreases, and another whose imaginary part vanishes as r + /L → 0 and whose real JHEP04(2020)196 Figure 2. Scalar field frequency as a function of the dimensionless scalar charge qL for a AdS-RN black hole with µ = 0.99µ ext and r + /L = 0.5 (mL = 0 and = 1). Left panel : Imaginary part of the dimensionless frequency, Im(ωL). The system becomes unstable for q > q where q L ∼ 1.863. Right panel : Real part of the dimensionless frequency, Re(ωL), measured with respect to qµL. This quantity changes sign at q = q , i.e. when Im(ωL) changes sign. In both plots, the black dashed curves describe the analytic prediction of the asymptotic matching expansion (A.10). We find very good agreement with the numerical results (blue curves) for qL < 1.1 (say). This is a further justification of the crude assumption that we should match with 0 in the overlapping region.
However we find that these modes connect with the AdS normal modes as r + → 0.
part approaches the normal modes of AdS. The unstable modes are found in this second family. This could well be the complete story. However, in de Sitter black holes there is a third family of quasinormal modes -called the near-extremal family -whose wavefunctions are spatially peaked near the horizon and that is distinct from the de Sitter family (as the name suggests, the latter is connected to the normal modes of de Sitter when the black hole shrinks). This naturally raises the question: could it be that in AdS one also has a near-extremal family of quasinormal modes that is not connected to the AdS family? If so, do the near-horizon unstable modes with bound (3.7) fit in this near-extremal family? We find a negative answer to these questions: the unstable modes belong to the AdS family of modes and the near-extremal family coincides with the AdS family. To arrive to this conclusion we first use a matching asymptotic expansion similar to the one used in de Sitter [26][27][28]54] to find the frequency spectrum of the near-extremal family of quasinormal modes. This is done in appendix A and here we just quote the final result: near-extremality and for small scalar field charge one finds that near-extremal modes have the frequency (for the lowest radial overtone p = 0)

JHEP04(2020)196
where R + = r + /L, e = qL and σ = r + −r − r + measures the distance away from extremality with r − (r + , µ, L) being the inner (Cauchy) horizon for which f (r − ) = 0. In figure 2, this analytical near-extremal frequency (with = 1) is described by the dashed black curve. We find that it matches quite well the numerical result for small scalar charge. This indicates that the unstable modes fit into the near-extremal family of modes. But they also fit into the AdS family of normal modes. That is to say, unlike in the de Sitter case, in AdS the near-extremal and AdS family of modes coincide. To see this is indeed the case we pick two solutions in figure 2 that have qL = 2.5 (orange diamond) and qL = 2 and (keeping µ = 0.99µ ext fixed) we follow this family of unstable modes as r + /L decreases to zero. 25 This is done in figure 3 for qL = 2.5 and figure 4 for qL = 2. In both cases we find that, as r + /L → 0, Im(ωL) → 0 and Re(ωL) → 3, which is indeed the normal mode frequency of AdS with = 1 (and lowest radial overtone).
In figure 3 and figure 4 the magenta dashed lines departing from the normal mode of AdS describe the frequency that one obtains when we consider a perturbative expansion in r + /L and near-extremality about global AdS (and = 1, m = 0). This result is taken from [12] (we already mentioned it to get the bound (3.8)): where e = qL. So we see that not only the unstable modes approach the normal modes of AdS but they also do it at the expected rate in an expansion in r + /L. The matching of our numerical results with the perturbative results (3.10) and (3.11) represents a non-trivial check of our results and illustrates the regime of validity of the perturbative results. Now that we have highlighted the key features of the near-horizon (and superradiant) instabilities due to scalar perturbations in AdS-RN, we can proceed to the study of perturbations of Dirac fields in AdS-RN.

Searching for an instability of Dirac fields in the AdS-RN background
In section 3 we have seen that scalar fluctuations in the AdS-RN background give rise to the near-horizon scalar condensation instability. Moreover, we have seen that this instability is closely associated to the violation of the AdS 2 scalar BF stability bound. So much that the associated stability bound (3.7) for the onset of the instability is sharp. This naturally invites the questions: in the fermionic case can we also have a range of parameters where the AdS 2 fermionic stability bound is violated? If so what is the equivalent bound to (3.7) for the onset of the instability? 25 Note that qL = 2.5 is well above the onset curves of figure 1 for any r+/L while the qL = 2 line is above the onset curves only above a certain horizon radius. So, for the latter charge, the system is unstable only above a critical value of r+/L, as shown in figure 4. Figure 3. Scalar field frequency as a function of the dimensionless horizon radius r + /L for a AdS-RN black hole with µ = 0.99µ ext and qL = 2.5 that is always above the near-horizon bound (3.7) (mL = 0 and = 1). Left panel : Imaginary part of the dimensionless frequency, Im(ωL). Right panel : Real part of the dimensionless frequency, Re(ωL). In both plots, the magenta dashed curves describe the analytic prediction of the perturbative expansion in r + /L about AdS. The unstable modes are thus connected to the AdS normal modes when r + → 0 (brown disk). For reference, the blue disk with r + /L = 0.5 is shown (which makes contact with figure 2).

JHEP04(2020)196
In this section we will address these questions. We will find that a near-horizon analysis of the Dirac equation indeed indicates that the AdS 2 fermionic stability bound can be violated near-extremality if the charge of the fermion is above a critical value (subsection 4.1). Encouraged by this result we will do a numerical analysis that will search for unstable modes in the region of parameters of interest (subsection 4.3). However, we will find no trace of instabilities, unlike in the scalar field case.

Argument for a near-horizon instability of Dirac fields
The Dirac equation in the near-horizon geometry (2.49) of the extreme AdS-RN black hole can be obtained taking the near-horizon limit of section 2.4 directly on the Dirac equation (2.24) for the extreme AdS-RN black hole. Concretely, applying the near-horizon coordinate transformation (2.48) together with the near-horizon frequency transformatioñ ω → ω ε/L 2 AdS 2 (so that e −iωt → e −i ωτ ) followed by the near-horizon limit ε → 0 yields the Dirac equation in the near-horizon geometry (2.49): 26 where the AdS 2 radius L AdS 2 and the Maxwell near-horizon parameter α are defined in (2.49) and λ is the angular eigenvalue quantized as in (2.26). Also, recall that m and q are the mass and charge of the fermionic field. Asymptotically, as ρ → ∞, a Frobenius analysis of (4.1) finds that the solution R 1 (ρ) decays as where α 1 , β 1 are two arbitrary constants and we have introduced the AdS 2 conformal dimensions 3) The s = 1/2 stability bound is independent of the spacetime dimension and still given by (2.29), m 2 ≥ 0 [35,37]. Thus, the 2-dimensional fermionic stability bound is obeyed if m 2 eff ≥ 0 in (4.3). It follows that we can have situations where the Dirac field obeys the 4-dimensional fermionic stability bound (2.29), m 2 ≥ 0, but violates the 2-dimensional stability bound. When this happens, i.e. when m 2 eff < 0, one might expect an instability. This condition can be rewritten: the 2-dimensional stability bound is violated if the charge of the fermion is larger than The equality applies strictly to the extremal case; as we move away from extremality, by continuity the instability should still be present but a higher fermion charge is needed to trigger it. Figure 5 illustrates the regions where (4.4) predicts instability/stability. Figure 5. Predicted Dirac field charge for the onset of an instability as a function of the horizon radius for an extremal AdS-RN black hole (µ = µ ext ). In both plots the red dashed curve is the lower bound of (4.4). The plot in the left (right) panel is for fermion mass mL = 0 (mL = 4) and harmonic number = 1/2. The near-horizon analysis of the 2-dimensional stability bound violation leading to (4.4) predicts that region B should be unstable while region A should be stable (at least with respect to the stability mass bound mechanism). Note that for small r + /L the system is not unstable because there is no superradiance for fermions and the predicted near-horizon instability is also suppressed. These Dirac figures can be (qualitatively) compared with figure 1 for the scalar field.

JHEP04(2020)196
At this level, we see that the near-horizon analysis of the possible violation of the AdS 2 stability bound for a Dirac field parallels very much the partner analysis done for a scalar field in section 3, with the minimum value for the charge (3.7) for the scalar case just replaced by the fermionic minimum value (4.4). In the scalar field case, we found (through a numerical study of linear perturbations in the AdS-RN background) that the violation of the 2-dimensional stability bound translates into the existence of a linear scalar condensation instability. Moreover, the near-horizon scalar bound (3.7) turns out to be very sharp, as best illustrated in figure 1. This scalar condensation linear instability indicates that non-linearly the AdS-RN black hole, when perturbed by a scalar field evolves towards a new configuration -a hairy black hole (with a scalar condensate floating above the horizon) -with the same UV asymptotics (since the 4-dimensional stability bound is satisfied) but with a different near-horizon geometry where the 2-dimensional stability bound is no longer violated [11,12,55,56].
These considerations motivate the study done in this manuscript for a Dirac field. In this case the AdS 2 stability bound can also be violated: at extremality this occurs for a fermion charge that saturates (4.4). From the lessons learned in the scalar field case one might well expect that the AdS-RN black hole, when perturbed by a Dirac field, is linearly unstable. To confirm whether this is the case, in the rest of this section we will solve numerically the Dirac equation in the AdS-RN background to hunt for a signature of the near-horizon linear instability. However, unlike the scalar field case, we will not find any evidence of a linear instability.

Dirac normal modes of global AdS
Before looking for potential instabilities (or frequency spectrum of damped oscillations) of Dirac modes in the global AdS-RN black hole it is convenient to first compute the normal mode spectrum of Dirac fields in global AdS. Indeed, some families of AdS-RN perturbations must reduce to these in the limit where the horizon shrinks to zero. Massive (section 4.2.1) and massless (section 4.2.2) Dirac fields require a distinct analysis.

Massive normal modes
For massive fermions in global AdS, it is not easy to solve directly the Dirac equations to get the radial functions R 1,2 . There is however an appropriate combination of R 1,2 that yields equations of motion that are explicit hypergeometric equations. The linear combination for R 1,2 that we use below is motivated by a similar analysis done to compute the massive normal modes of fermions for de Sitter in [57].
For m = 0 and in global AdS, we introduce the new radial variable y = −ir/L and make the following field redefinitions where f 1,2 (y) are functions to be determined. In these conditions, the coupled system of Dirac equations (2.22) yields Adding and subtracting these two ODEs yields This pair of coupled first order ODEs can be straightforwardly rewritten as a decoupled pair of second order ODEs for f 1 and f 2 . Moreover, if we introduce the new radial coordinate z = y 2 and the field redefinitions each of the ODEs becomes a hypergeometric ODE with the standard form with parameters a i , b i and c i given by The most general solutions of (4.8) are [58] is the Gaussian (ordinary) hypergeometric function and A 1,2 , B 1,2 are arbitrary amplitudes. We can now plug (4.11) into (4.7) and then into (4.5) to get the most general solution for R 1 (r) and R 2 (r). Finally, we can insert this most general solution for R 1,2 (r) into (2.39) to get the most general solution for the Dirac fields ψ ± (r). These are the physical fields that have to be regular everywhere and this constrains some of the amplitudes A 1,2 and B 1,2 and the frequencies. Namely, at the origin, r = 0, one finds that both ψ ± have two divergent terms of the form B 2 /r +3/2 and (2A 2 −B 2 )/r +1/2 . Regularity at the origin thus requires that one sets A 2 = 0 and B 2 = 0 and the other two amplitudes A 1 and B 1 are left arbitrary. It follows that the regular normal eigenmodes are We have not yet imposed the asymptotic boundary condition. A Frobenius analysis of (4.12) near the conformal boundary together with the use of (2.39) finds that ψ ± behaves JHEP04(2020)196 as (2.40) or (2.42) with For m > 0 (m = 1/2) the no-source standard boundary condition (2.43) requires α 1 = 0. Using Γ[−p] = ∞ for p = 0, 1, 2, · · · this quantizes the frequency as ωL = + 2 + mL + 2p or ωL = −( + 1 + mL + 2p) , (standard quantization) (4.14) For 0 < mL < 1 2 we can also impose the alternative quantization (2.44), i.e. β 1 = 0. This quantizes the frequency spectrum as (also with radial overtone p = 0, 1, 2, · · · )

Massless normal modes
In this section we find the normal modes in global AdS for a massless fermionic field. These have been previously discussed in [29,59] but these references have not identified the full spectra of frequencies.
For m = 0 in global AdS, introducing the change of coordinates and field redefinition the radial equation (2.24) can be rewritten as a hypergeometric ODE in the standard form Its most general solution is [58]  Introducing this into (4.16) one gets R 1,2 (r) (note that R 2 = R * 1 as discussed in the next section). Plugging this into (2.39) one finds the most general solution for the Dirac fields ψ ± (r). At the origin, r = 0, these ψ ± have a divergent term proportional to C 2 r − −3/2 . Regularity at the origin thus requires that we set C 2 = 0. Now we need to impose the JHEP04(2020)196 asymptotic boundary condition. One finds that asymptotically ψ ± decays as (2.41) with As explained previously, for m = 0 we can impose either the standard or alternative boundary conditions. The no-source standard boundary condition (2.43),  .20) and (4.21) were computed in [29] using vanishing flux boundary conditions that, as explained in the end of section 2.3, are exactly the AdS/CFT standard and alternative boundary conditions. However, [29] missed the existence of half of the normal mode spectrum, namely the half part that has negative frequencies. The relevance of the full spectrum (and associated relations between standard/alternative quantizations) is further analysed in the discussion of figure 10. Further note that in RN, the four families of modes that reduce to (4.20)-(4.21) in the AdS limit become completely independent (i.e. they are not related by complex conjugation and the "degeneracy" is broken). This is further discussed in the next subsection.

Setup of the numerical problem
In this section we solve numerically the Dirac equation and search for linear instabilities of the Dirac solution in the AdS-RN background. Before proceeding it is important to note that: 1) the Dirac radial equation (2.24) for R 2 (r) is just the complex conjugate of the radial equation for R 1 (r) so if R 1 (r) is a solution one automatically has R 2 (r) = R 1 (r) * , and 2) the Dirac angular equations for S 1,2 are related by the symmetry θ → π − θ so if S 1 (θ) is a solution then S 2 (θ) = S 1 (π − θ). Therefore, we just need to find the solutions R 1 (r) (S 1 (θ) are just the spin-weighted s = 1/2 spherical harmonics with quantum number ).
Further note that if R 1 has charge q then R 2 = R * 1 has charge −q, and complex conjugation maps quasinormal modes to quasinormal modes. It follows that if ω = ω r + iω i is a linear mode frequency of R 1 then −ω * = −ω r + iω i is a linear mode frequency of R 2 = R * 1 . Thus, if we compute the frequency spectrum of R 1 , we have the spectrum of R 2 too. It also follows that there is no loss of generality in assuming that qQ > 0 in our analysis: results for qQ < 0 are obtained simply by reversing the sign of the real part JHEP04(2020)196 of the frequencies. Finally note that when we compute ω we have to allow both positive and negative values of ω r , i.e. if ω r = Re(ω) is a frequency of R 1 there is no symmetry in the system that requires −ω r to be also a frequency of R 1 . The only exception is if µ=0 (or e = 0) and m = 0 i.e. a massless Dirac field in Schwarzschild-AdS. In this case For concreteness, we will set the mass of the fermion to zero, i.e. we solve the Dirac equation (2.24) for R 1 with m = 0 (and in the gauge A| ∞ = 0 where the frequency of the fermionic wave is ω; see footnote 15) subject to the physically relevant boundary conditions. The asymptotic decay of R 1 (r) is given in (2.30). For reasons discussed previously, we impose the asymptotic boundary condition (2.43) (standard quantization, ψ . At the horizon, for a non-extreme black hole, a Frobenius analysis finds that the two pairs of independent solutions are (4.22) Rewriting this in ingoing Eddington-Finkelstein coordinates (v, r, θ, φ), with v = t + f −1 dr, which are smooth across the future event horizon H + , we find that regularity of R 1 (r) at H + requires that we impose the boundary condition B out = 0. 27 For the numerical solution it is convenient to redefine 23) and to work with the compact radial coordinate such that the horizon is now at z = 0 and the asymptotic infinity at z = 1. This has the advantage that analytical solutions that are smooth at the horizon simply have to obey the horizon Neumann boundary condition q (z = 0) = 0. The asymptotic boundary condition for q(z) follows straightforwardly from (2.43)-(2.44) and (4.23)-(4.24).
The numerical methods that we use are very well tested [10,11,39,40,[62][63][64][65][66][67][68] and reviewed in [69]. To discretize the field equations we use a pseudospectral collocation grid on Gauss-Chebyshev-Lobbato points. The eigenfrequencies and associated eigenvectors are found using Mathematica's built-in routine Eigensystem. For a given , this method has the advantage of finding several modes (i.e., from distinct families and with distinct radial overtones) simultaneously. However, to increase the accuracy of our results at a much lower computational cost we use a powerful numerical procedure which uses the Newton-Raphson root-finding algorithm discussed in detail in section III.C of the review [69]. All our results have the exponential convergence on the number of gridpoints, as expected for a code that JHEP04(2020)196 uses pseudospectral collocation. In particular, all the results that we present are accurate at least up to the 10th decimal digit.
The Dirac equation in AdS-RN also has the scaling symmetry that determines that the physical dimensionless quantities are those listed in (3.9).

Main results
As discussed in section 4.3, for fermion mass m = 0, we can have two independent homogeneous boundary conditions that yield normalizable modes: the standard (ψ (0) + = 0) and alternative (ψ (0) − = 0) boundary conditions. Moreover, for each of these boundary conditions, the eigenvector R 1 can have negative or positive real part of the frequency. It follows that, for a given harmonic , m φ (and m = 0) we have a total of two frequency spectra to discuss for each one of the two possible boundary conditions. Before proceeding to the actual physical analysis of the frequency spectrum and instabilities of the system, in appendix C we first test our numerical code by comparing the associated numerical results with some analytical perturbative expansions that are derived in appendix B. This confirms that our numerical code is generating physical data and we can now proceed and discuss our main physical findings.
Our aim is not to present the full spectrum of frequencies of a Dirac field in AdS-RN black hole. Instead, we are motivated to search for unstable modes, i.e. on eventually finding modes that, in some range of parameters, have Im(ωL) > 0. There is a wide window of parameters to explore although the instability, if it exists, should appear near extremality for fermion charges q above a critical value. Thus one needs a good strategy to hunt efficiently for unstable modes. We proceed as follows. From the near-horizon bound (4.4) arguing for the existence of an instability, we see that this bound is lower if we set m = 0 and = 1/2. So, first, we either: 1) fixed m = 0, = 1/2 and µ close to µ ext , and varied {r + /L, qL}, or 2) fixed m = 0, = 1/2 and qL and varied {r + /L, µ}. In both cases, as described in the end of section 4.3, we solved our system as an eigenvalue problem for ωL. This finds "all" the solutions of the system (as long as the hierarchies do not grow large, e.g. |ωL| 1, which makes the numerical problem hard). This allows to eventually identify unstable modes with Im(ωL) > 0 or, in the worst case, to identify modes with Im(ωL) < 0 that are closest to the marginal case for instability (Im ω = 0). Once these interesting modes are identified we then used a Newton-Raphson root-finding algorithm to follow efficiently the modes to other values of the parameter space.
In spite of our efforts, we have found no sign of an instability. Recall that for mL = 0 both the standard (2.43) and the alternative (2.44) boundary conditions yield normalizable modes. In general, we do find that the stable modes with smallest |Im(ωL)| are those that reduce to the alternative normal modes of AdS (4.21) or to the standard AdS normal modes (4.20), when the horizon radius shrinks to zero. Among these, we further find that the modes with smallest |Im(ωL)| are, for both quantizations, the ones that reduce to the positive normal mode frequencies when r + /L → 0, i.e. ωL = 3/2 (alternative boundary condition) and ωL = 5/2 (standard quantization). Therefore, to avoid distraction from the main point, in the rest of this manuscript we only discuss these two families of modes. Probably the plots that best illustrate the main conclusions of our Dirac study are those of figure 6 (for alternative quantization) and of figure 9 (for standard quantization). Recall that in the best case scenario the expectation is that, close to extremality, modes should become unstable above a fermion charge q that should be higher than the nearhorizon bound (4.4). Thus, in these figures we fix the black hole horizon to be r + /L = 0.5 and choose a chemical potential close to extremality, µ = 0.99µ ext . Starting from qL = 0, where Im(ωL) < 0, we then increase this charge to see if there is a critical value above which Im(ωL) becomes positive. (That is to say, we adopt a similar strategy as the one followed in the scalar field case to get figure 2).

JHEP04(2020)196
For the alternative quantization, the left panel of figure 6 shows that, starting from q = 0, as qL grows, Im(ωL) < 0 increases and approaches Im(ωL) = 0 very closely. However, no matter how large qL is we never reach a situation where Im(ωL) ≥ 0. Interestingly, there is a critical value of q, namely qL = q max L ∼ 0.9390 (vertical brown dashed line) where Im(ωL) reaches a maximum value of Im(ωL) ∼ −0.000548 (see the inset plot which zooms-in around this maximum). But increasing qL further, Im(ωL) becomes again increasingly more negative (instead of becoming positive). The Dirac field system behaves therefore substantially distinctly from the scalar field case of figure 2 (left panel) where there was a critical qL above which Im(ωL) becomes positive. To complete the analysis, in the right panel of figure 6, we plot Re(ωL) − qµL. We find that for small qL this quantity is positive but becomes negative above q = q ∼ 0.9344/L. Interestingly, this occurs at a charge that is smaller than q max where the maximum of Im(ωL) is reached ( dashed line): this is better seen in the inset plot which zooms-in the relevant region. Again we note the difference to the scalar field case displayed in the right panel of figure 2 where Re(ωL) − qµL changes sign precisely at the critical value of qL where Im(ωL) = 0. Further note that these plots also show that for a Dirac field we do not have a value of qL for which we simultaneously have Re(ωL) − qµL = 0 and Im(ωL) = 0. Therefore, we cannot set ωL = qµL in the equations of motion and solve these as an eigenvalue problem for the instability onset charge. That is to say, unlike the scalar field case, we do not have an onset charge that would produce the partner plots of the scalar field onset plots of figure 1. The predictions of figure 5 do not hold (at least at the linear mode level).
We have done similar experiments as those of figure 6 for other black hole parameter values µ and r + /L. Keeping µ fixed, black holes with distinct r + /L have plots similar to figure 6 with the feature that larger values of r + /L reach the maximum of Im(ωL) (but remaining negative) at smaller critical values of q = q max . On the other hand, keeping r + /L fixed, black holes with distinct µ also have similar plots to figure 6 with the property that larger values of µ reach the maximum of Im(ωL) (but still negative) at smaller critical values of q = q max and this maximum of Im(ωL) is increasingly closer to zero as µ approaches the extremal value µ ext .
To have a complementary perspective of the system's properties, in figure 7 and in figure 8 we illustrate other attempts we have made to find an instability. In figure 7, we keep the alternative quantization and fix the chemical potential at µ = 0.999µ ext , and plot Im(ωL) as a function of r + /L for five different values of qL, namely, qL = 0.7, 0.8, 0.9 (from bottom to top in the left panel) and qL = 1, 1.1 (right panel). (The two plots are needed for the presentation of the results because the relevant qL = 1 case in the right panel reaches a maximum that is approximately two orders of magnitude higher than the JHEP04(2020)196 first three cases on the left panel). The main feature in these plots is the typical presence of a local minimum and local maximum (bump). As we increase the fermion charge from zero to a value slightly above 1, the relative minimum and relative maximum of Im(ωL) raise and shift to lower values of r + /L. But the local maximum always has Im(ωL) < 0, i.e. there is no instability. However, for charges qL above a value that is in between 1 and 1.1, the local minimum and maximum are no longer present and Im(ωL) decreases monotonically with r + /L (see e.g. qL = 1.1 displayed as the yellow curve in the right panel; higher values, qL ≥ 1.1, have a similar monotonic behaviour). As yet another illustration of experiments we made, in figure 8 we fix the fermion charge to be qL = 1 (which was already analysed in figure 7 for µ = 0.999µ ext ) and we study the effect that changing the chemical potential has by considering a total of 5 curves with 5 different values of µ. Namely, in the left plot we consider the cases µ = 0.9µ ext and µ = 0.95µ ext . These cases have no bump (no local maximum) and illustrate that it only appears close to extremality. In the right panel we show three more cases where we fix µ = 0.99µ ext , µ = 0.995µ ext and µ = 0.999µ ext (from bottom to top). The bump is now present and the local maximum increases as one approaches extremality but never becomes positive. For the case µ = 0.999µ ext this local maximum is at Im(ωL) ∼ −5.93 × 10 −6 .
So far we have focused our discussion of the results for the alternative quantization case because, typically, for the same values of black hole parameters this is the case where Im(ωL) approaches Im(ωL) = 0 the most. Nevertheless, we have also tried hard to find an instability in the standard boundary condition (2.43) case. Again without success. To illustrate briefly this conclusion, in figure 9 we give the partner plot of figure 6 but this time for the standard quantization. Although the features of figure 9 are clearly more JHEP04(2020)196 Figure 9. Dirac field frequency with standard quantization (2.43) as a function of the dimensionless scalar charge qL for a AdS-RN black hole with µ = 0.99µ ext and r + /L = 0.5 (also, mL = 0 and = 1/2). Left panel : Imaginary part of the dimensionless frequency, Im(ωL) which attains a maximum of Im(ωL) ∼ −0.0037491 for q = q max ∼ 3.3555/L (vertical brown dashed line). The main inset plot zooms-in around this maximum and shows that Im(ωL) < 0 for any qL. The secondary inset plot shows the detail of the curve around qL ∼ 1.45 to show that the apparent cusp in the main plot is smooth. Right panel : Real part of the dimensionless frequency, Re(ωL), measured with respect to qµL. This quantity changes sign at q = q ∼ 3.2873/L with q < q max , i.e. for a smaller qL than the one where Im(ωL) attains its maximum value (vertical brown dashed line): this is better seen in the inset plot which zooms-in the relevant region. elaborated than those of figure 6 (e.g. there are several local maxima and minima), the main conclusions are still the same: i ) one always has Im(ωL) < 0; ii ) there is a q = q max where the solution approaches Im(ωL) = 0 the most (vertical brown dashed line); iii ) Re(ωL)−qLµ changes sign at q = q < q max . It follows that we find no sign of an instability and the standard boundary condition case, much like the alternative quantization case, also gives results that are very different from the scalar field case of figure 2.

Discussion and conclusions
A scalar field in an asymptotically AdS 4 Reissner-Nordström black hole can satisfy the asymptotically AdS 4 UV Breitenlöhner-Freedman (BF) stability bound but violate the infrared 2-dimensional BF stability bound associated to the AdS 2 ×S 2 near-horizon geometry of the extremal black hole of the system, as reviewed in section 3. When this is the case,

JHEP04(2020)196
the AdS-RN black hole is unstable to scalar condensation and the system evolves to a new configuration in the phase diagram of solutions that preserves both the UV BF bound and the near-horizon 2-dimensional stability bound. Such a solution is a hairy black hole with a charged scalar field floating above the horizon [8,9,11,55,56,70,71]. Coulomb repulsion balances the gravitational force and the system is static. There is no doubt that the violation of the AdS 2 stability bound is the physical mechanism responsible for the near-horizon scalar condensation instability since the associated minimum bound (3.7) on the scalar field charge that triggers the instability is sharp (at extremality) as best demonstrated by figure 1.
Given these considerations, the study done in this manuscript for Dirac field perturbations in the global AdS 4 RN black hole was motivated by the following observation. Dirac fields in AdS-RN can also preserve the UV fermionic stability bound (2.29) [35,37] but violate the near-horizon infrared fermionic stability bound, as seen in section 4.1. From the scalar field case lessons, this suggests that the system might be unstable to fermion condensation. However, in spite of our efforts to scan the relevant parameter space near extremality, we found no sign of a linear mode instability. The sharp distinction between the scalar and Dirac field cases is best illustrated comparing the scalar figure 2 with the Dirac figure 6 (for alternative quantization) or figure 9 (for standard quantization). Of course our numerical study does not prove linear stability but we did such a detailed scan that we are very confident that no linear instability is present. Our stability results are also consistent with the stability study of fermions in planar AdS, where no instability was found [9,19] (see also [20][21][22][23][24][25]). 28 Indeed, the planar AdS case is the r + /L → ∞ limit of the global AdS system. So, the planar AdS studies [9,[19][20][21][22][23][24][25] and our present study in global AdS establish that the violation of the 2-dimensional stability bound of a Dirac field in AdS-RN does not lead to a linear mode instability. However, such solutions correspond to negative energy Schrödinger states: without a positive self-adjoint extension for the Schrödinger operator the dynamical evolution of the system should develop an instability. . . In particular, the system might indeed still be unstable if non-linear effects play a role in the discussion. That is to say, if we perturb a AdS-RN black hole with a Dirac field in a region of parameter space where the infrared stability bound is violated, it could still be the case that the system evolves non-linearly to a new configuration that has a charged Dirac field floating above the horizon and that preserves both the UV and IR stability bounds. How difficult would it be to prove whether this scenario is correct?
One must proceed with caution. To begin with one needs to first formulate more precisely the setup of the problem. It is certainly much harder to find, if they exist, the proposed Dirac hairy black holes than it was to construct the scalar hairy black holes [8,9,11,55,56,70,71]. There is a fundamental difference between fermionic and bosonic fields. The fermion has no classical limit: Planck's constant is present in the stress tensor and associated equations of motion for a fermion. As discussed in detail in section 14.3

JHEP04(2020)196
of Wald's textbook [72], the absence of a classical limit means that in Einstein's equation we have to promote the differential Einstein and energy-momentum operators G and T to quantum operators and the quantum version of Einstein's that gives the back-reaction of a fermion on the gravitational field is G µν = 8π( T Max µν + T µν ), where T Max µν and T µν stands for the Maxwell and Dirac stress tensor contributions and · · · stands for the expectation value of the corresponding operator.
Thus, to find the backreaction that fermions induce on the gravitoelectromagnetic background one needs to first compute the expectation value of the fermion energy momentum tensor T µν . This is a highly non-trivial task. Even worse, once we consider the quantum backreaction of fermions one also needs to consider the quantum backreaction of gravitons and photons, i.e. one also needs to compute G µν and T Max µν [72]. In a best case scenario, where we have a large number N of Dirac fields, one might be able to assume that, roughly speaking, the effects of N Dirac fields are N times as relevant as that of the gravitons and photons [72]. For a 'fermionic hairy black hole', the fermionic condensate should be made of a large number of fermions. In these conditions, for large N , one might be able to neglect the quantum backreaction of gravitons and photons and work in the semi-classical limit whereby the backreaction of the Dirac field on the gravitoelectromagnetic background is simply governed by G µν = 8πT Max µν + 8πN T µν . This semi-classical system should be viewed as the leading term of a 1/N expansion of the full theory [72]. But this semi-classical computation still requires that one computes T µν . And this is still a remarkable task. An overview on the physical and technical tools required to accomplish this task can be found in [24,73,74] (and references there-in) where asymptotically planar AdS quantum electron stars are discussed as semi-classical solutions of Einstein-Maxwell theory.
Finally note that in the present manuscript we focused our attention on modes that could eventually become unstable. We have not studied in detail the full spectrum of quasinormal mode frequencies of the Dirac field in AdS-RN. Moreover, we focused on the case of a massless fermion because, as explained previously, this was enough for our purposes. However, the equations of motion and relevant boundary conditions for any fermion mass and any sector of perturbations are given in section 2. We have also computed the normal modes of massive fermions in AdS (previously only the massless spectrum was computed). It might be useful to have a more complete frequency spectra study for future studies/applications. It might also be interesting to look for perturbations of spin 3/2 Rarita-Schwinger fields about AdS-RN black holes. In this case, there are also normalizable solutions that become unstable for negative square masses [36]. Probably there will be no linear near-horizon instabilities when the effective 2-dimensional mass violates the AdS 2 stability bound but, as far as we know, this was never checked.

A Near-extremal modes in AdS
For a Reissner-Nordström de Sitter (Λ > 0) background, in [26][27][28] it was found that there is a family of quasinormal modes -denoted as the 'near-extremal' family of modesthat is distinct from the 'de Sitter' family quasinormal mode, where the latter connects to the normal modes of de Sitter when the horizon radius shrinks to zero size, r + → 0.
This naturally raises the question of whether there is also such a 'near-extremal' family of modes in AdS and, if so, wether they do or not coincide with the AdS family of modes. In this appendix, we address this question in the simplest case, namely in the case of a (charged) scalar field that obeys the Klein-Gordon equation (3.2). More concretely, we arrive to the near-extremal frequency (3.10) which is used in the main text (see section 3 and the discussion there-in of the dashed curves of figure 2) to show that in AdS the 'near-extremal' and AdS families of modes coincide (unlike in the de Sitter case).
The 'near-extremal' modes we seek obey (3.2) in the background (2.1) (we will work with the gauge choice C = 0) and, at least in the near extremal limit, are expect to be highly peaked near the horizon. So we want to simultaneously zoom into the horizon and approach extremality. For that we first introduce the dimensionless quantities For x 1 one is close to the outer horizon and for σ 1 the inner and outer horizon are very close, i.e. one is close to extremality. Next, we take the limit σ → 0 whilst keeping z = x σ fixed. From [26][27][28] the 'near-extremal' modes are expected to saturate the superradiant bound ω = qµ at extremality so onwards we measure the frequency difference δω with respect to this bound via the redefinition Using the condition f (r − ) = 0 for the location of the inner horizon one can find r − = r − (r + , Q, L) which is then inserted into (A.1) to express Q as a function of (r + , σ, L).
In these near-extremality conditions, we are ready to find the near-horizon solution of the Klein-Gordon equation. Concretely, introducing the above redefinitions into the Klein-Gordon equation (3.2), to leading order in σ, we obtain:
With the field redefinition 3) is rewriten as a standard hypergeometric ODE The regular (i.e. ingoing) solution at the future event horizon is given by where 1 F 2 (a, b, c; z) is the standard hypergeometric function and a − , a + are defined by: In the context of a matched asymptotic expansion, the near-region (near-horizon) solution (A.7) must now be matched with the far-region solution of (3.2) (in near-extremality conditions). As explained above we expect the 'near-extremal' modes we are looking into to have wavefunctions that die-off very quickly away from the black hole horizon (at least near-extremality). Therefore, as a first rude approximation we take the far-region to be described by a vanishing wavefunction. That is to say, in the overlapping region, we match the near-region solution (A.7) with φ = 0. In the end of the day, this approximation turns out to be quite good because the analytical approximation for the 'near-extremal' frequency that we obtain -see (A.11) -matches remarkably well the numerical solution of (3.2). This is best seen comparing the black dashed analytical curve of our expansion in figure 2 with the numerical blue dot results. For this reason, we do not try to improve further our matching asymptotic approximation.
Proceeding in these conditions, the leading order behaviour of the large R = r/L series expansion (z → −∞) of φ, namely φ ≈ (−z) ± √ 1+4η , needs to be matched with the farregion solution φ = 0. Before we can do it, we still need to distinguish the cases 1 + 4η ≥ 0 and 1 + 4η < 0. For our proposes (comparing with the numerical results of section 3), we want to consider the small scalar field charge case for which one finds that 1 + 4η ≥ 0 holds as long as e ≤ e c where e 2 c =

JHEP04(2020)196
(the reader also interested in the case 1 + 4η < 0 can follow the steps detailed in [28]). In these conditions, from the matching condition one finds that where p = 0, 1, 2, · · · is the radial overtone of the mode. Replacing this into (A.2) one finally finds that 'near-extremal' modes have a frequency given by: This is (3.10) in the main text when we set the radial overtone p = 0.

B Perturbative results
Although we solve the Dirac equation numerically in the main text, it is very good practice to testify the numerical results against analytical predictions that can be obtained within perturbation theory in some region of the parameter space. Therefore, in this appendix we find some useful analytical perturbative approximations for the Dirac frequencies. More concretely, this appendix is divided in two parts. In appendix B.1 we use a matching asymptotic expansion approach to find the frequency approximations (C.1)-(C.2) that, in figure 11 of appendix C, are compared against our numerical results. Then, in appendix B.2 we use a systematic perturbative expansion in the dimensionless horizon radius r + /L 1 (with no further approximations) to find the analytical frequency approximations (C.3)-(C.4) which, in figure 12 of appendix C, are also compared with our numerical results. In both cases, there is agreement between the analytical approximation predictions and the numerical results in the regime of parameter space where the former are valid.
In this appendix, as in the main text, we solve the Dirac equation (2.24) in a AdS-RN background (2.1) (gauge choice C = 0) with regular (ingoing) boundary conditions at the future event horizon and the standard (2.43) or alternative (2.44) boundary conditions at the conformal boundary. We will work exclusively with vanishing fermion mass, m = 0.

B.1 Matched asymptotic expansion
In this section we derive an analytical expression for the imaginary part of the Dirac frequency using the method of matched asymptotic expansion introduced in [75,76] (see also e.g. [77][78][79]). We assume r + L and split our spacetime into two regions; an asymptotic globally AdS far region where the effects of the black hole can be neglected and a near region about the black hole outer horizon where the effects of the cosmological constant can be neglected. In each region the associated perturbation equation can be solved analytically, then matching the near and far region solutions in their overlapping region will fix the integration constants as well as the imaginary part of the perturbation frequency ω. More concretely, the near region is defined by r − r + 1/ω and the far region is defined JHEP04(2020)196 by r − r + r + . It follows that the overlapping region exists for ωr + 1. A further assumption we must make is that the Coulomb interaction is weak, Qq 1, where q (Q) is the fermion (black hole) charge.

B.1.1 Far region solution
Since in the far region the effect of the black hole (BH) is assumed to be negligible, we effectively have a fermion field in the global AdS background. Thus, the general solution for the massless fermionic field R 1 is given by (4.18) that we reproduce here: where C 1,2 are two arbitrary amplitudes to be determined below. Asymptotically this solution decays as (2.30), namely This solution has to satisfy the asymptotic boundary condition (for a massless fermion). For the standard quantization this is (2.43) while for the alternative quantization the boundary condition is (2.44) which fix β 1 as a function of α 1 or, equivalently, C 1 as a function of C 2 . This yields where the upper sign refers to the standard quantisation (2.43) and the lower sign to the alternate quantisation (2.44). Note that we do not impose any boundary condition at a inner boundary since this far region solution does not extend till there.

B.1.2 Near region solution
In the near region we can approximate ∆(r) = r 2 f (r) by: where r − ≈ Q 2 2r + . This follows from the assumption that r + L and therefore in the near region we have r ∼ O(r + ) L, so we can neglect the r 2 /L 2 term in f (r). Further applying the near region assumptions to the Dirac equation (2.24) we can neglect terms of order ωr + or higher powers. Other terms appear which are dominated by a 1/∆ term in the small black hole (BH) approximation near the horizon; therefore we can evaluate the numerators of these terms at r ≈ r + .
With these approximations, and the coordinate transformation (the horizon r = r + is at z = 0) the Dirac equation is approximately given by the near region equation, Making the field redefinition  Using the property of the hypergeometric function 2 F 1 (a, b, c, 0) = 1 we find that the z → 0 behaviour of the near region solution is R near 1 (z) ≈ z 1 4 (αz −iσ + βz iσ ). Requiring regularity (only ingoing modes) at the horizon implies that we must set β = 0.

B.1.3 Matching
To find the large r (z → 1) behaviour of the near region solution (B.11) with β = 0, we use the z → 1 − z transformation law of the hypergeometric function [58]  In the overlapping region, one must have R near 1 large r = R far 1 small r . That is to say, we must match independently the r 1 2 + and r − 1 2 − terms of (B.13) with those of (B.14). This matching yields: where we have used (B.9) and (B.7) forσ to restore the explicit dependence on the frequency ω. The ratio C 2 C 1 follows straightforwardly from (B.3). We want to solve the transcendental equation (B.15) to get an analytical expression for the frequency. There is no closed form solution, unless we do some educated approximations that we now discuss. Since we are working with very small and weakly charged RN-AdS black hole, one expects that the mode frequencies in such a background are close to the massless Dirac normal modes frequencies of AdS 4 , already computed in (4.20) (standard quantization) or (4.21) (alternative quantization). Denote this normal mode frequency by ω AdS 4 . However, since the background now has a horizon, the system becomes dissipative and with respect to the normal modes of AdS 4 , the frequency of the system should acquire a small imaginary contribution. Denote it by i δ. In (B.15), it is thus a good approximation to replace the frequency ω by ω = ω AdS 4 + i δ with |δ| ω AdS 4 . Our target now is to solve (B.15) at leading order for δ ω AdS 4 ∼ O(1). A posteriori, we compare the prediction of our analytical computation with the numerical result to confirm that this approximation is valid. Equation (B.15) has the additional challenges that: 1) the frequency appears in the argument of the Gamma functions and 2) Γ [−1 − 2 ] in the numerator diverges for the allowed values (2.26) of the harmonic number (recall that Γ[−p] = ∞ for non-negative p). To deal with these obstacles, we use the Gamma function property Γ[z + 1] = zΓ [z] JHEP04(2020)196 and the assumptions of our problem, ωr + 1 and δ ω AdS 4 ∼ O(1). This allows to expand the Gamma functions whose argument depends on ω (and thus on δ) to extract δ out of the argument of the Gamma functions. In particular, this permits to find that the divergence of Γ [−1 − 2 ] in the numerator is cancelled by the Gamma function in the denominator that depends on δ.
In these conditions, one finds that the leading order solution for δ is whereK is a positive real number for each , n. For the alternative boundary condition (2.44) it is given by (for overtone p = 0) These are the analytical predictions we use in (C.1)-(C.2) of appendix C to compare against the numerical results: see figure 11. We expect these analytical predictions to be valid only for small horizon radius and away from extremality and for Qq 1, which we are able to confirm numerically in section 4.4. An added bonus for this method is that we have an expression for general . The method we present in the next appendix below has to be done for each individually, but it is more systematic than this one, since it only requires an expansion in r + /L 1.

B.2 Perturbative expansion in R +
In this section we find an analytical prediction for the frequency using a systematic perturbative expansion in r + /L. Unlike in the previous subsection, the only approximation that will be made is that the expansion parameter of this expansion is small, r + /L 1. We will do this expansion up to the order that finds the first correction (in the real part JHEP04(2020)196 of the frequency) to the global AdS normal mode frequency. Should we wish, we could go one order higher in the analysis and find also the correction to the imaginary part of the frequency (although this is computationally more demanding). For our purposes of comparing with the numerical results, it is enough to have the correction to the real part of the frequency (the results of appendix B.1 already allow us to test independently the imaginary part).
The systematic perturbative expansion in r + /L 1 used in this appendix was first introduced in [55] and further explored in [11-13, 80, 81] where the reader can find full details of the method (we will be very succinct in our exposition). In short, we split our spacetime into a near and far regions. We expand the frequency Ω = ωL and the field R 1 in each region in a power series in R + = r + /L: We now series expand the Dirac equation (2.24) in small R + . The leading, zeroth order equation is simply the Dirac equation in global AdS 4 for a massless fermion (B.1) (that we already studied in sections 4.2.2 and B.1.1). Not surprisingly, the small R = r/L expansion of this leading order solution breaks down at order R + /R. This motivates splitting our spacetime into a far region, R R + , and a near region, R + ≤ R 1. In the far region we work with the radial coordinate R but in the near region we work instead with the radial coordinate y = R/R + (since the far region small R expansion breaks down at order R + /R).
In the far region, at each order in R + , we impose the standard boundary condition (2.43) or the alternative boundary condition (2.44). In the near region we impose boundary conditions that only allow for ingoing waves at the horizon. We then perform a matching procedure at each order in R + , in the region where the far and near region overlap, to determine the frequency coefficients Ω (k) , as well as amplitudes that were not fixed by the two boundary conditions. At leading (zeroth) order, we fix Ω (0) to be the normal mode frequency for a massless fermion already obtained in (4.21) or (4.20) for the boundary conditions (2.44) or (2.43), respectively. We will do this for the mode with harmonic number = 1/2, and radial overtone to be p = 0. Our aim is then to find the first frequency correction Ω (1) due to the presence of the black hole. For concreteness, in most of our discussion below we only explicitly present details of the case where we impose the alternative boundary condition (2.44). We then present the final result also for the standard boundary condition (2.43).
In the far region the leading order R 0 + solution is (B.1) and imposing the asymptotic boundary conditions amounts to repeat mutatis mutandis the analysis done in (B.1)-(B.3). With our choice of = 1/2 and p = 0 this fixes the frequency at order zero to be Ω (0) = 3 2 ; see (4.21). To fix the normalization, we set the amplitude of the Dirac field at infinity to be 1 at all orders in R + : R far 1 | R→∞ = 1 + O(1/R). Introducing the near region radial coordinate y = R/R + , still at leading order R 0 + , the near region Dirac equation (for = 1/2) reads 1 2 (y − 1) 2y − µ 2 ∂ 2 y ψ near (0) + y −

JHEP04(2020)196
The solution which is regular at the horizon (y = 1) is ψ near (0) = α (0) cosh 2 log 2 y − 1 + 4y − 2µ 2 + iβ (0) sinh 2 log 2 y − 1 + 4y − 2µ 2 , (B.23) with We must now match the far and near regions solutions at order R 0 + in their overlapping region R + R 1. This procedure, typically fixes all other constants of the problem that were not fixed by the boundary conditions. The large R expansion of ψ near (0) is ψ near (0) | large R = α (0) + β (0) R + · · · whereas the small R expansion of ψ far (0) is ψ far (0) | small R = (−1) 3/4 R + · · · . Therefore, matching ψ near (0) | large R = ψ far (0) | small R requires that we set α (0) = 0 and β (0) = (−1) 3/4 . Collecting the results at order 0 for the alternative quantization (2.44) we have: ; For the standard quantisation one has a similar result with Ω (0) = 5/2. We can now proceed to the first order R 1 + contribution that enables us to find Ω (1) . We repeat the procedure outlined above for the far and near regions, imposing the alternative boundary condition in the far region and ingoing boundary condition at the horizon in the near region. The near region equation for ψ near (1) is the same as at order R 0 + (i.e. there is no source). After imposing the horizon boundary condition, we find that the large R expansion of the near region solution gives R near large R ≈ (−1) 3/4 R + β (1) R + + · · · . (B.26) The far region equation at order R 1 + can also be solved analytically for ψ far (1) . As usual, we find a solution with two integration constants, say C 1 and C 2 . These constants will in general depend on µ, e and Ω (1) . Firstly by imposing the alternative quantisation we find an expression for C 2 in terms of µ, e and Ω (1) . Then we again impose that our field has amplitude 1 at infinity, yielding an expression for C 1 in terms of µ, e and Ω (1) . After doing this we find the small R behaviour of the far region solution to be In the overlapping region R + R 1, (B.27) must match (B.26). A straightforward matching of the terms R + R 0 fixes the constant β (1) in (B.26). On the other hand, there is no R + /R term in the large R series expansion of the near region solution. It follows that Ω (1) must be such that it eliminates the corresponding R + /R term in the small R expansion of the far region solution (B.27). This fixes Ω (1)  The frequencies (B.29) and (B.30) are the analytical frequency approximations (C.3) and (C.4) that we reproduce in the main text and that we compare against the numerical data in figure 12 of section C.

C Comparing numerical results with analytical expansions
To test our numerical code we first compute the quasinormal mode frequencies for a massless Dirac field in global AdS and in Schwarzschild-AdS: see figure 10 for the standard (2.43) and alternative (2.44) quantizations. When r + = 0 our frequencies reduce to the normal mode frequencies of global AdS (4.20) and (4.21) for the standard and alternative quantizations, respectively. Recall that is the spin-weighted harmonic number and n is the radial overtone (related to the number of nodes of the wavefunction along the radial direction). So the smallest |ω| is obtained for = 1/2 and n = 0. For the alternative boundary condition these are ωL = 3/2 and ωL = −5/2, while for the standard quantization these are ωL = 5/2 and ωL = −3/2). All the numerical results we will present describe solutions with = 1/2 and n = 0. If there is an instability it should already be present in this sector of perturbations (see the argument in section 4.1).
On the other hand, for finite r + /L our numerical curves reproduce the values first computed in [29]. As explained in the end of section 2.3, this is because the AdS/CFT standard (2.43) and alternative (2.44) quantizations have vanishing energy flux and correspond precisely to the boundary conditions imposed in [29]. To complete the spectrum of Schwarzschild-AdS (and global AdS) note that frequencies that are the negative of the complex conjugate of the frequencies (−ω * ) plotted in figure 10 are also eigenvalues of the system (which was missed in [29]).
Next we test our numerical code for AdS-RN. As a first test, in appendix B.1 we use a matching asymptotic expansion method to find that for ωr + 1 and qQ 1 (and m = 0, = 1/2, n = 0) the imaginary part of the frequency is approximately given by Im(ωL) ≈ − 1 4π Alternative quantization; (C.1)

Standard quantization. (C.2)
In figure 11 we plot Im(ωL) as a function of r + /L, for µ = 1 2 µ ext and a fermion charge qL = 1 (blue dots). We compare these numerical results with the analytical result (C.1).
Both agree for small r + /L, i.e. in the regime were the matching expansion (C.1) is valid. The plots for the real part of the frequency are displayed in figure 12. The analytical perturbative results are (C.3) and (C.4). An important feature of this result is that for a fixed chemical potential, the slope of the real part of the frequency (the term proportional to R + ) changes sign for a certain electric charge. This coincides with our numerical findings.

JHEP04(2020)196
In the top panels we show the real part of the frequency for the boundary condition (2.44) with qL = 0.7, 0.8 and µ = µ ext (1 − 10 −3 ). In the bottom panels of figure 12 we show the real part of the frequency for the boundary condition (2.43) as a function of the horizon radius for fixed chemical potential µ = µext 2 and fermion charges qL = 1.9, 2 from left to right. The dashed red lines are the analytical results (C.4) and (C.3) respectively. Note that for these parameters the change of sign occurs for charges qL ∼ 1.95 and qL ∼ 0.71 respectively.
As a general comment, it is perhaps worthy to comment that we find that matching asymptotic and perturbative analysis like the ones provided in appendices B.1 B.2 are valid for smaller windows of r + /L in the Dirac field case when compared with similar analysis done for the scalar field.
We have confirmed that our numerical code is generating physical data. Our main physical results are reported in section 4.4 of the main text.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.