Eigenvalue repulsions and quasinormal mode spectra of Kerr-Newman: an extended study

The frequency spectra of the gravito-electromagnetic perturbations of the Kerr-Newman (KN) black hole with the slowest decay rate have been computed recently. It has been found that KN has two families $-$ the photon sphere and the near-horizon families $-$ of quasinormal modes (QNMs), which display the interesting phenomenon of eigenvalue repulsion. The perturbation equations, in spite of being a coupled system of two PDEs, are amenable to an analytic solution using the method of separation of variables in a near-horizon expansion around the extremal KN black hole. This leads to an analytical formula for the QNM frequencies that provides an excellent approximation to the numerical data near-extremality. In the present manuscript we provide an extended study of these properties that were not detailed in the original studies. This includes: 1) a full derivation of a gauge invariant system of two coupled PDEs that describes the perturbation equations \cite{Dias:2015wqa}, 2) a derivation of the eikonal frequency approximation \cite{Zimmerman:2015trm,Dias:2021yju} and its comparison with the numerical QNM data, 3) a derivation of the near-horizon frequency approximation \cite{Dias:2021yju} and its comparison with the numerical QNMs, and 4) more details on the phenomenon of eigenvalue repulsion (also known as \emph{level repulsion}, \emph{avoided crossing} or \emph{Wigner-Teller effect}) and a first principles understanding of it that was missing in the previous studies. Moreover, we provide the frequency spectra of other KN QNM families of interest to demonstrate that they are more damped than the ones we discuss in full detail.


Introduction
When a black hole is (moderately) perturbed, it typically relaxes back to equilibrium by emitting gravitational waves with damped characteristic frequencies − the quasinormal mode (QNM) frequencies − that depend on the conserved charges of the black hole. It follows that these QNM frequencies may be used to determine the mass and angular momentum of a black hole. In fact, this is one way of measuring the mass and angular momentum of the final black hole [4] that emerges from the black hole binary coalescences observed in gravitational wave detector experiments [5][6][7][8][9][10].
Astrophysical black holes are expected to be described by Einstein gravity; more specifically, by its Kerr solution parametrized by the mass M and angular momentum J ≡ M a (where a is the rotation parameter) [11]. Therefore, all LIGO-Virgo [6,7] observations of events compatible with black hole binaries [8][9][10] have been described so far mainly under the working assumption that the coalescing objects can be modelled by the Kerr solution or parametrically small deviations thereof [4]. However, to discuss the physical interpretation of the observed data, we might also want to consider black hole solutions of the Einstein-Maxwell theory that have an electric charge Q, in addition to M and J. 1 In this case, the uniqueness theorems [14,15] guarantee that the Kerr-Newman black hole (KN BH) [16,17] parametrized by M , J and Q is the unique, most general, analytic, stationary asymptotically flat electro-vacuum black hole of Einstein-Maxwell theory. The Kerr [5], Reissner-Nordström (RN) [18,19] and Schwarzschild [20] black holes are then viewed as limiting cases of KN with Q = 0, a = 0 and Q = a = 0, respectively.
Although astrophysical black holes are expected to quickly lose any electric charge that they may have [21,22], one should nevertheless study the properties of KN black holes and compute their quasinormal mode frequencies. With this theoretical information at hand, we will be better equipped to analyse and interpret observational data to unequivocally establish that the observed system has no charge (or even to compute its charge in the lucky but unlikely event of observing a system during the short timescale where the discharge has not yet occurred). Furthermore, the QNM spectra of KN might be of interest for other interpretations of observational data and for applications in both ground and spacebased gravitational wave detectors [6,7,[23][24][25][26]. For example, it can be used to model gravitational wave emission [27], and it might even be useful for constraining some dark matter models [28] and modified gravity models [29]. For these reasons, in this manuscript we conclude a series of papers, started in [1,3], that compute the main families of QNMs of the KN BH and identify their key properties.
The QNM spectra of Schwarzschild, RN and Kerr black holes were determined many decades ago [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46] (see review [47]). This was possible at a relatively small computational cost because for these black holes the QNM spectrum turns out (remarkably) to be encoded in a single separable equation that effectively yields a pair of angular and radial ODEs that one can solve as an eigenvalue problem. For Schwarzschild and RN black holes this is known as the (odd mode) Regge-Wheeler and (even mode) Zerilli equations [30][31][32], while for the Kerr black hole this is known as the Teukolsky equation [38]. The existence of such a simplification allows one to find the QNM spectra, and in doing so, to establish evidence in favour of the linear mode stability of these solutions and to ultimately motivate a formal proof of the linear mode stability of the Kerr solution [48]. 2 The state of affairs is very different in the Kerr-Newman case. Generic gravitoelectromagnetic perturbations of KN are no longer described by a single separable equation. Thus, initial hints about the QNM spectra of KN were obtained only within perturbation theory about the RN or Kerr black holes: perturbative results in the small rotation parameter a about RN were discussed in [58,59], and perturbative results in the small charge parameter Q around Kerr were computed in [60].
To make further progress and compute the KN QNM spectrum for generic Q and J, one must solve the perturbed Einstein-Maxwell equation which is a coupled partial differential equation (PDE) system. Naïvely, one expects to find a system of nine coupled PDEs. However, working in the so-called phantom gauge, Chandrasekhar reduced the problem to the study of 'just' two coupled PDEs [40] (see also [60]). Despite this significant progress, finding the QNM spectrum and addressing the problem of the linear mode sta-bility of the KN BH has remained an open problem for several decades. Further progress was made in [1] where it was shown that generic gravito-electromagnetic perturbations of KN (except for those that change the mass and angular momentum of the solution) are described by a coupled system of two PDEs for two gauge invariant Newman-Penrose (NP) fields. Upon gauge fixing, these reduce to the coupled PDE system originally found by Chandrasekhar [40,60]. Moreover, in [1] a numerical search of KN modes was finally performed in regions of the KN parameter space that could be more prone to developing an instability, finding none and thus providing evidence for the linear mode stability of KN (further supported by the non-linear time evolution study of [61]). More recently, in [3], the numerical code of [1] was made computationally more efficient and extended to compute the frequency spectra, across the full KN 2-parameter space, of the most dominant (i.e. with slowest decay rate) gravito-electromagnetic QNM family. These are the modes that reduce -in Chandrasekhar's notation [40] -to the Z 2 (i.e. gravitational), = m = 2, n = 0 modes in the Schwarzschild limit (a = Q = 0), where the harmonic number gives the number of zeros of the eigenfunction along the polar direction and n is the radial overtone. In the process, [3] found that KN has not one but two main families of Z 2 = m = 2 QNMs which were coined the photon sphere (PS), and the near-horizon (NH) families, although the sharp distinction between the PS and NH modes is unambiguous only for small rotation a, i.e., when the KN black hole is close to the Reissner-Nordström family. Quite remarkably, [3] further found that as we evolve along the KN parameter space, the imaginary part of the frequency of these two PS and NH families intersect each other (however, the real part of the frequency is very similar for the PS and NH modes and, typically, does not display crossings). Sometimes this intersection of the imaginary part of the frequencies is a simple crossover where the modes simply trade dominance but, other times this interaction is much more intricate and displays a behaviour that suggests repulsions between the PS and NH modes. These "eigenvalue repulsions" were unexpected since they are not observed in the QNM spectra of neither Kerr nor Reissner-Nordström. 3 As a result of these repulsions, well away from the RN limit of the KN solution, the PS and NH families lose their individual identities and instead combine to yield what is more appropriately described as a PS−NH family of QNMs and its radial overtones.
In the current manuscript we complement and complete the studies of [1,3] in five main ways: 1. We use the Newman-Penrose (NP) formalism to derive the aforementioned coupled system of two PDEs for two gauge invariant NP variables, first presented in [1], that describes the most general gravito-electromagnetic perturbations of KN (except for those that change the mass and angular momentum of the solution) and that reduces, upon gauge fixing, to the Chandrasekhar PDE system [40,60]. This derivation was only very briefly sketched in [1] but we now give a detailed derivation of it in Sec-tion 2. We also take the opportunity to revisit a simple proof of isospectrality of the Schwarzschild and RN QNM spectra [1].
2. We can envisage solving the perturbation equations for the two gauge invariant NP fields in a WKB analysis at large |m| = 1. Similar to the Schwarzschild and Kerr cases, the leading order contribution of this analysis, known as the eikonal or geometric optics limit |m| = → ∞, is expected to be closely connected to the properties of unstable null circular orbits revolving around the KN black hole. In Section 3.1 we will compare this eikonal result with the numerical data for photon sphere modes to conclude that the eikonal frequency indeed provides a relatively good approximation to the PS frequencies that gets better as m grows.
3. There is a second class of QNMs that have eigenfunctions that, near-extremality, are very localized around the event horizon and quickly decay to zero away from the horizon. These are the near-horizon modes or the PS−NH modes that were already mentioned above. This suggests doing a 'poor-man's' matched asymptotic expansion (MAE) whereby we take the near-horizon limit of the perturbed equations to find the near-region solution and match with a vanishing far-region wavefunction in the overlapping region where both solutions are valid. Remarkably, this can be done because the perturbation equations, in spite of being a coupled system of two PDEs, can be solved analytically in the near-horizon region around the extremal (zero temperature) KN black hole using the method of separation of variables. Ultimately, this is possible because the near-horizon limit of the extremal KN BH is a warped circle fibred over AdS 2 (Anti-de Sitter in 1+1 dimensions) and thus its perturbations can be decomposed as a sum of known radial AdS 2 harmonics. The system of 2 coupled PDEs for the gauge invariant NP fields in the near-horizon region of the near-extremal KN geometry separates into a system of 2 decoupled radial ordinary differential equations (ODEs) and a coupled system of 2 angular ODEs. We can solve this near-horizon system, match it with the trivial far-region, and obtain an analytical expression for the NH and PS−NH frequencies. The final expression was presented in [3] but not the long derivation that leads to it. We will present this detailed derivation in Section 3.2 and show that it provides an excellent approximation to the numerical frequencies when we are close to extremality. 4. In the Reissner-Nordström background, there are exactly two distinct sectors of QNMs: the aforementioned PS and NH families (and their radial overtones). However, as we move away from this limit in the KN parameter space we find that this clear distinction between the two families is lost and the two families and their overtones combine in an intricate way to form what is more appropriately described as PS−NH modes and their radial overtones. This happens because the phenomenon of eigenvalue repulsion occurs. These eigenvalue repulsions were already reported in [3] but in Section 4 we will give a detailed description of these eigenvalue repulsions in the KN QNM spectra, and we will see how the frequency gaps between different QNM families develop and evolve. No less important, we will provide a first princi-ples understanding of this phenomenon that was not discussed in [3]. For that we will start by pointing out that eigenvalue repulsion is common in some eigenvalue problems of quantum mechanical systems where it is also known as level repulsion, avoided crossing or Wigner-Teller effect [64,65]. In Section 4.1 we will start by reviewing (following §79 of the Landau-Lifshitz textbook [64]) the simplest quantum mechanical two-level system with a self-adjoint Hamiltonian that exhibits avoided crossing. We will then extend the discussion of avoided crossing to the case where the perturbed Hamiltonian of the system is not self-adjoint, as is the case with the KN QNM system. Having understood that level repulsions should be present in the QNM spectra of KN, in Section 4.2 we will give a detailed description of eigenvalue repulsions in the frequency spectra of KN. The analysis of Section 4 together with the one of Section 3.2 will allow us to conclude that the complex frequencies ω of KN have level crossing (i.e. both the real and imaginary parts of the PS and NH modes cross each other) exactly at one, and only one, point in the 2-dimensional KN parameter space (we collect strong evidence to claim that this is the point at extremality where the PS modes reach Im(ω) = 0, which will be represented by a in Fig. 12). In all other KN black holes we either have no crossovers of the imaginary and real parts of the frequency or the imaginary part of the PS and NH frequencies cross, but not the real part of the frequencies. These features are in agreement with the predicted properties of the eigenvalue spectra of a 2-dimensional parameter space system with avoided crossing, as explained in Section 4.1. This analysis will also explain why avoided crossing is not observed in the 1-parameter family of Kerr solutions. Ultimately, the intricate QNM spectra of KN emerges from the fact that level crossing occurs only at one point but the system reacts to avoid crossings at other points. This leads to the observed elaborate features/repulsions when one is approaching the level crossing point of the system.

5.
After revisiting in Section 5 the properties of the Z 2 = m = 2 KN QNMs (first presented in [3]) that are expected to be the least damped ones, in Section 6 we will present the frequencies of some other relevant gravito-electromagnetic modes of KN. This will give solid, explicit, evidence that the Z 2 = m = 2 QNM is indeed the mode with the slowest decay rate in KN (as with the Schwarzschild, RN and Kerr black holes).

Derivation of the gauge invariant perturbation equations for KN
In subsection 2.1 we briefly review the Kerr-Newman black hole solution. Then, in subsection 2.2, we detail how the Newman-Penrose (NP) formalism can be used to derive a coupled system of two PDEs for two gauge invariant NP variables [1] that describes the most general gravito-electromagnetic perturbations of KN (except for those that change the mass and angular momentum of the solution). Finally, in subsection 2.3 we discuss the boundary conditions that allow one to solve the final eigenvalue problem to find the QNM frequencies of KN.

KN black hole: an algebraically special Petrov type D solution
The KN BH solution with mass M , angular momentum J ≡ M a and charge Q is most commonly expressed in standard Boyer-Lindquist coordinates {t, r, θ, φ} (time, radial, polar, azimuthal coordinates) [16,17], in which the metric takes the form with ∆ = r 2 − 2M r + a 2 + Q 2 and Σ = r 2 + a 2 cos 2 θ. Roots of the function ∆, namely correspond to the inner and outer event horizons, respectively. Physically, one is most interested in the outer event horizon (r = r + ), which is a Killing horizon generated by the Killing vector with angular velocity Ω H and temperature T H given by where we have used (2.2) to express M as a function of r + , a and Q. If r − = r + , i.e. a = a ext , the KN BH has a regular extremal ("ext") configuration with T ext H = 0, and maximum angular velocity Ω ext (2.5) Here, we are interested in linear gravito-electromagnetic perturbations about the KN background. Following Teukolsky [38,66], we work within the Newman-Penrose (NP) formalism [36]. We will not review the NP formalism here, but instead refer the reader to Chapter 7 of [67] for a comprehensive review. Suffice it to say that the NP formalism starts with a complex null frame or tetrad 4 and uses this tetrad to transform all quantities of interest (connection coefficients, Ricci, Weyl and Maxwell field strength components) into complex scalars. In such a manner, the Weyl tensor, for example is transformed into a set of five complex scalars: Ψ a (a = 0, 1 · · · , 4) or the Maxwell field strength into a set of three complex scalars: Φ a (a = 0, 1, 2) [40,67]. Furthermore, the existence of a NP frame in which a certain combination of the Weyl scalars vanishes determines the Petrov type of the background solution.
Teukolsky [38,66] showed that on an algebraically special vacuum background, which is defined to be one in which there exists a null frame so that Ψ 0 = Ψ 1 = 0, the linear perturbations of the background may be expressed in terms of a decoupled equation There is a spinor version of the NP formalism. However, here, we deal only with the Lorentzian version.
where O is some linear second-order differential operator and Ψ (1) 0 is the gauge-invariant perturbed value of Ψ 0 .
Fortunately, the Kerr BH, which was of principal interest for Teukolsky, is algebraically special. In fact it is Petrov type D (i.e. a NP frame exists such that the only non-vanishing Weyl scalar is Ψ 2 ). Furthermore, given global and hidden [68] symmetries of the Kerr BH, the coordinate dependence of the perturbations separate leading to a single ODE. Thus, the combined simplification of decoupling and separability on the Kerr BH allows one to study its linearized mode perturbations [38,47,66] and prove its linear mode stability [42] (see also footnote 2).
Like its vacuum cousin the Kerr BH, the KN BH is also Petrov type D. In particular, in a NP null frame {e (a) } = { , n, m,m} with (a = 1, 2, 3, 4) adapted to the principal null directions, given by wherer = r + ia cos θ, 5 the only non-zero Weyl scalar is 6 Moreover, the only non-zero Maxwell scalar is (2.11) However, importantly, the decoupling result of Teukolsky does not apply to the KN BH, since it is a non-vacuum solution. In fact, such a decoupling result does not seem possible for the KN BH (see e.g. [40]). The best that can be done is to derive a gauge-invariant coupled PDE system [1], which we now derive and which reduces to the Chandrasekhar system [40] under a particular gauge choice. 5 The standard notation for the complex conjugation in the NP formalism is to use a bar. We will stick to this notation as far as NP quantities are concerned. However, this should not be confused withr defined here, whose complex conjugation (r * ) we shall denote with a star. 6 Recall that the 5 complex Weyl scalars Ψa in the NP formalism encode the information in the 10 independent components Cµνρσ of the Weyl tensor, 8) and the 3 complex NP scalars Φa encode the information in the 6 independent components of the antisymmetric Maxwell field strength, F = dA,

Derivation of the gauge invariant perturbation equations
To discuss generic perturbations of the Kerr-Newman black hole one needs to find the perturbed Einstein-Maxwell equation which, a priori is a system of nine coupled PDEs. Although, a decoupling result cannot be obtained for the perturbations on the KN BH background, one can still reduce this perturbation system to a simpler set of two gauge invariant coupled PDEs [1] that, after gauge fixing, reduces to the Chandrasekhar coupled system of two PDEs [40,60]. In this section we give the details of the derivation of this system of PDEs. The linearised perturbations on any background satisfy the linearised Einstein equation on that background. Therefore, any perturbation equation that we derive must ultimately come from some operator acting on this linearised Einstein equation [69].
There are two common ways of deriving the Teukolsky equations: the original method relies on a particular manipulation of the NP Bianchi equations, which now comprise the non-trivial content of the Einstein equations [38]. Another method is a more straightforward contraction of the Penrose wave equation into the NP null frame (R µνρσ is the Riemann tensor) [70]. While the second method is more prescriptive and does not require much guesswork as to which equations to look at and how to manipulate them, the former method requires less calculation, once a strategy has been determined. Therefore, we shall derive the coupled equations using the Bianchi equations, which, of course, coincide with those derived from the Penrose wave equation. We derive the equations as follows. First, let us settle the notation. In this section all equations labelled as (7.xx) refer to equation (xx) in chapter 7 of [67]. Since these are long equations we do not reproduce them here and simply refer the reader to that reference. Further recall that the fundamental quantities in the NP formalism needed to study perturbations are the directional derivative operators [40,67] 13) and the 12 complex spin coefficients defined from linear combinations of the 24 background Ricci rotation connection coefficients γ cab = e µ (c) e ν (b) ∇ ν e (a) µ [40,67], (2.14) Their complex conjugates (denoted by a bar) correspond to the replacement 3 ↔ 4 in γ cab . 8 Consider the expressionδ (7.32d) −∆(7.32c) (2.15) as a first order perturbative equation. Let us consider the left hand side of this expression, 9 which involves second order in derivative quantities: 10 Now, the operator acting on Ψ 3 is a commutation operator. Therefore, we can use the equation (7.6c) to rewrite it in terms of first order in derivative quantities, which can themselves be turned into zeroth order in derivative quantities using equations (7.32c) and (7.32d). However, the third term involving Φ 2 cannot be similarly simplified. At best we can use the commutation relations to rewriteδ∆ in terms of∆δ. Therefore, it is clear already at this stage that a decoupled equation is not going to be possible on the KN background and at best we can only hope to derive a coupled equation involving Ψ 4 and Φ 2 . Further simplifying the first order in derivative terms using equations (7.32f ), (7.32g), (7.32j) and (7.32k) and using various NP equations (7.21a)-(7.21r), as well as the Maxwell equations (7.22)-(7.25), gives At this stage we find that perturbed quantities λ and ν are obstructions to a coupled equation involving Ψ 4 and Φ 2 . However, inspecting the Bianchi equations closely, we find that λ appears in equation (7.32c) with coefficient 3Ψ 2 + 2Φ 1Φ1 , while ν appears in (7.32d) with coefficient 3Ψ 2 − 2Φ 1Φ1 . Thus, we can solve for λ and ν in terms of differential operators on perturbed quantities Ψ 3 , Ψ 4 and Φ 2 . Significantly, the operator on λ in equation (2.17) and the form of (7.32c) means that Ψ 3 will end up with a second order derivative∆δ, the same operator that acts on Φ 2 in equation (2.17). Thus, we have the possibility of defining a perturbed quantity involving a particular combination of Φ 2 and Ψ 3 such that this quantity couples with Ψ 4 . Studying the form of the equations, it is not too difficult to conclude that such a quantity can be defined and is of the form (2.18) 8 KN is Petrov type D so, from the Goldberg-Sachs, one must have κ = σ = ν = λ = 0. Moreover, one has = 0 because we have chosen to be tangent to an affinely parametrized null geodesic µ ∇µ ν = 0. 9 Note that in equations (7.xx), the derivative terms are written on the left hand side, while the rest of the terms are placed on the right. 10 Generally, we will use the notation that NP scalars with superscript (0) refer to scalars in the KN background and the superscript (1) to first order perturbations of the scalar. However, in the equations below, for brevity, we suppress the superscripts. From the expressions above for the background NP scalars it should be clear what is a background quantity and what is a first order perturbed quantity.
The resulting equation is of the form The second coupled equation is derived in a similar manner, except that it is now easier, because we know that the perturbed quantity that couples to Ψ 4 is ϕ −1 as defined in (2.18). Thus we begin by considering (2.20) using the fact that an equation for (∆D −δδ)Ψ 3 may be obtained from and an equation for (∆D −δδ)Φ 2 may be obtained from ∆(7.23) −δ(7.25).
The strategy used to simplify the resulting equation is very similar to that used to derive equation (2.19). Therefore, without going through the details, we give the resulting coupled equation: In summary, we have derived from the NP equations two coupled PDEs, (2.19) and (2.21), satisfied by ϕ −2 = Ψ 4 and ϕ −1 = 2Φ 1 Ψ 3 − 3Ψ 2 Φ 2 , which are invariant under infinitesimal diffeomorphisms and tetrad rotations [40].
These NP scalars ϕ −2 and ϕ −1 are the ones relevant for the study of perturbations that are outgoing at future null infinity and regular at the future horizon. Note that we could have equally derived a set of coupled equations satisfied by Ψ 0 and 2Φ 1 Ψ 1 − 3Ψ 2 Φ 0 ; the positive spin counterparts of ϕ −2 and ϕ −1 . Such equations are simply obtained via Geroch-Held-Penrose (GHP) transformations [37] of equations (2.19) and (2.21) and are relevant for perturbations outgoing at past null infinity.
We can now substitute the background values of the NP quantities into equations (2.19) and (2.21) (recall footnote 10). Since ∂ t , ∂ φ are Killing vector fields of the KN background, its gravito-electromagnetic perturbations can be Fourier decomposed as e −iωt e imφ , where ω and m are the frequency and azimuthal quantum number of the mode. Moreover, we rescale the perturbed quantities, {ϕ −2 , ϕ −1 } → {ψ −2 , ψ −1 } as [1]: Having done this, we obtain the following coupled system of two PDEs (first presented in [1]): 11 where the second order differential operators {F, G, H} are given by with α ± ≡ 3(r 2 M −rQ 2 ) ± Q 2r * −1 , and the radial and angular Chandrasekhar operators [40] are defined The complex conjugate of these operators, namely D † j and L † j , can be obtained from D j and L j via the replacement K r → −K r and K θ → −K θ , respectively. 11 There is a set of two coupled PDEs -related to (2.23) by a Geroch-Held-Penrose [37] transformation -for the quantities ψ2 and ψ1 that are the positive spin counterparts of (2.23); however these would be relevant if we were interested in perturbations that were outgoing at the past null infinity.
Note that fixing a gauge in which Φ (1) (2.23) reduces to the Chandrasekhar coupled PDE system [40] (see also the derivation in [60]). Finally, note that in the limit Q → 0 equations (2.23) decouple yielding the familiar Teukolsky equation for Kerr [38].
Before finishing this section, we take the opportunity to discuss a property of QNMs of Schwarzschild and RN black holes that has raised a lot of attention in the literature. This is the fact that the spectra of the Regge-Wheeler (aka odd or axial) [30], and Zerilli (aka even or polar) [31] QNMs is isospectral (i.e. these two QNM families have exactly the same frequency) [40]. In the Schwarzschild limit, (2.23) decouples and we can look independently at the gravitational ψ −2 or electromagnetic ψ −1 perturbations. In particular, the decoupled equation for ψ −2 (ψ −1 ) corresponds to the original Teukolsky master equation [38] for gravitational (electromagnetic) perturbations in Kerr with a = 0. Thus, we can use the Teukolsky master equation to study the QNMs of the Schwarzschild black hole, instead of the Regge-Wheeler−Zerilli (RWZ) formalism. The two must give the same spectrum. The Teukolsky formulation has a single gauge invariant variable ψ −2 that must translate into two gauge invariant variables in the RWZ formulation, namely Regge-Wheeler's Φ v and Zerilli's Φ s eigenfunctions. Refs. [71][72][73] give the unique differential map that allows one to derive Φ v and Φ s from ψ −2 : see e.g., equation (4.16) of [73] (which holds for any cosmological constant). Isospectrality is the statement that Φ v and Φ s have the same QNM spectrum. Since Φ v and Φ s are constructed from the same Teukolsky NP gauge invariant variable ψ −2 , it follows that the eigenfrequencies of the Regge-Wheeler and Zerilli QNMs must necessarily be the same. This proves the isospectral property of QNMs in Schwarzschild and RN black holes and shows that this property is only non-trivial when viewed from the perspective of the Regge-Wheeler−Zerilli formalism. 12

Boundary conditions of the problem
To have a well-posed boundary value problem we must supplement the coupled PDE system (2.23) with appropriate (physical) boundary conditions. At spatial infinity, we require only outgoing waves, and at the future event horizon, we keep only regular modes in ingoing Eddington-Finkelstein coordinates. Moreover, we must require regularity at the north (south) pole θ = π (−π). In this subsection, we state what the conditions that these boundary conditions impose on the fields {ψ −1 , ψ −2 } are. 13 Recall that ω and m are the frequency and azimuthal quantum number m of the linear mode perturbations, respectively. The t − φ symmetry of the KN BH allows one to consider only modes with Re(ω) ≥ 0, as long as we study both signs of m. Then, to solve the coupled PDEs (2.22), we need to impose physical boundary conditions. At spatial infinity, a Frobenius analysis of (2.23) yields two independent solutions that at leading order behave as C ± e ±iωr . Imposing the boundary condition C − = 0, i.e. allowing only 12 The proof given in [1] and revisited here is for ψ−2 in the Schwarzschild black hole but it extends trivially to ψ−1 modes and the RN background. 13 The reader interested on a more detailed discussion of boundary conditions in perturbation problems about asymptotically flat backgrounds can see e.g. [74][75][76].
outgoing waves yields the decay: where s = −2, −1, and β s (θ) is a function of α s (θ) and its derivative fixed by expanding (2.23) at spatial infinity. At the horizon, the boundary condition must be such that only ingoing modes are allowed. A Frobenius analysis at this boundary gives two independent solutions, 26) where A in , A out are arbitrary amplitudes and Ω H , T H are the angular velocity and temperature defined in (2.4). To impose the correct boundary condition, we introduce the ingoing Eddington-Finkelstein coordinates {v, r, x, φ}, which extend the solution through the horizon. These are defined via The boundary condition is determined by the requirement that the metric and Maxwell field perturbations are regular in these ingoing Eddington-Finkelstein coordinates. This happens if and only if ψ s (r) behaves as where ψ EF s (r) is a smooth function. Thus, we must set A out = 0 in (2.26). To conclude, at the horizon, a Frobenius analysis whereby we require only regular modes in ingoing Eddington-Finkelstein coordinates, yields the expansion where b s (θ) is a function of a s (θ) and its derivative. At the north (south) pole x ≡ cos θ = 1 (−1), regularity dictates that the fields must behave as (ε = 1 for |m| ≥ 2, while ε = −1 for |m| = 0, 1 modes) where B + s (r)(B − s (r)) is a function of A + s (r)(A − s (r)) and its derivatives along r, whose exact form is fixed by expanding (2.23) around the North (South) pole.
The PDE system (2.23) subject to the above boundary conditions that describe the gravito-electromagnetic QNMs of the KN black hole with parameters {M, a, Q} has a useful scaling symmetry. When we scale the metric and Maxwell field strength as g µν → Λ 2 g µν and F µν → ΛF µν , for an arbitrary constant Λ, the equations of motion are left invariant. This means we can scale out one of the 3 parameters of the solution. Therefore, we can work with the adimensional parameters {ã,Q} ≡ {a/M, Q/M } (or {a/r + , Q/r + }) andω ≡ ωM . To find the frequency spectra of KN BHs we thus 'just' need to scan a 2-dimensional space.
To solve the PDE problem numerically, we use a pseudospectral method that searches directly for specific QNMs using a Newton-Raphson root-finding algorithm. We refer the reader to the review [77] and [74][75][76][78][79][80][81][82][83][84] for details. The exponential convergence of the method, and the use of quadruple precision, guarantee that the results are accurate up to, at least, the eighth decimal place.
3 Two families of QNMs: photon sphere and near-horizon modes The frequency spectra of gravito-electromagnetic perturbations of KN has two main families of QNMs: 1) the photon sphere (PS), and 2) the near-horizon (NH) families. Each of these families can dominate the frequency spectra (i.e. have lower |Imω|) depending on the region of the parameter space we look at. These two families are the natural extension to the rotating case (a = 0) of the PS and NH families of QNMs that are present in the Reissner-Nordström case, although this sharp distinction between the two families is unambiguous only for small rotation a parameter. There are particular regimes of parameter space where the frequency of each of these two families can be captured by perturbative expansions (WKB expansion and/or a near-horizon matched asymptotic expansion). This allows us to identify these two families of QNMs (thus providing the basis for their nomenclature), while providing also analytical formulae that give good approximations to the actual frequencies. Therefore, in this section we discuss in detail two useful perturbative analyses. In subsection 3.1 we consider a large m WKB expansion that identifies the PS QNMs, while, in subsection 3.2 we describe a simple but efficient matched asymptotic expansion that captures the NH modes.
Before considering the perturbative analyses, it is enlightening to identify the two families of QNMs in the simplest black hole where they co-exist. This is the Reissner-Nordström (RN) black hole. This identification will be our guide once we delve into the parameter phase space of KN away from the RN limit. It will also allow us to speculate about expectations for the KN QNM spectra that will then be discussed in the next subsections 3.1-3.2 and in section 4.
In Fig. 1 we plot the frequency spectra for = m = 2, n = 0 gravitational (Z 2 ) QNMs in the RN BH, 14 which are the QNMs with slowest decay rate in RN and KN as we will demonstrate in section 6. We see that there are two clearly distinct families of QNMs: 1) the PS family (orange diamonds) which dominates for a wide range of charge, namely for Q <Q RN c ∼ 0.9991342 (and reduces to the Schwarzschild QNM when Q = 0 [40,41]), and 2) the NH family (green circles) which becomes the slowest decaying mode forQ RN c <Q ≤ 1 and approaches Imω = 0 at extremality as best seen in the inset plot (in RN these NH modes have Reω = 0 for anyQ).
By continuity, once rotation is turned on but with smallã = a/M , we expect the KN spectra to be similar to Fig. 1 (NH modes should still approach extremality with Imω = 0 but this time, as we confirm later, with Reω = mΩ ext H = 0). Moreover, it also seems reasonable to expect the existence of a line − let us denote it asQ =Q c (ã) − that describes the intersection of the PS and NH surfaces and that eventually extends from Q c (ã = 0) =Q RN c (identified in Fig. 1) all the way up towards extremality. However, and ❬ ❭ ( where the PS mode is comfortably the dominant one. Right panel: Real part of the frequency as a function of charge. The NH mode has Reω = 0 in the RN limit (and only in this limit) and is not shown. Thus, the real part of the PS and NH frequencies have no crossing (unlike the imaginary part). interestingly, our full numerical results will prove that our expectations are only partially correct. Indeed, the KN frequency spectra for fixed but smallã is similar to Fig. 1. In particular, for 0 ≤Q <Q c (ã), the PS family has the lowest |Imω| and forQ c (ã) <Q ≤ Q ext it is the NH QNM that has slowest decay rate. Moreover, keeping a/a ext small, these two families trade dominance along their intersection lineQ =Q c (ã) with a simple crossover in the imaginary part of the frequency like the one observed in the inset plot of the left panel of Fig. 1 (the real part of the frequencies display no crossing as is clear in the right panel of Fig. 1 for a = 0). However, as we keep increasing the rotation a/a ext , we find that we enter a region of the parameter space (a window ofQ) where an unexpected change occurs: instead of having simple crossovers in Im(ω) where the PS and NH families should intersect, one starts observing intricate eigenvalue repulsions in Im(ω) that will be discussed in section 4 and associated Figs. 13-14, and the sharp distinction between PS and NH modes is lost (in this region PS and NH modes have similar Re(ω) with no crossings).
So much so that the QNM families will now be a combination (to be made precise later) of the two old modes in what will be more properly denoted as PS-NH families and their radial overtones.
After this brief summary of the findings to come, let us discuss the eikonal and nearhorizon analytical descriptions of the KN modes.
3.1 Photon sphere modes in the eikonal limit: analytical formula for the frequencies In the eikonal or geometric optics limit ∼ |m| 1, where a WKB approximation holds, there are QNM frequencies − known as "photon sphere"(PS) QNMs − that are closely related to the properties of the unstable circular photon orbits in the equatorial plane of the KN black hole. Namely, the real part of the PS frequency is proportional to the Keplerian frequency Ω c of the circular null orbit and the imaginary part of the PS frequency scales with the Lyapunov exponent λ L of the orbit [85][86][87][88][89][90][91][92][93][94]. The latter describes how quickly a null geodesic congruence on the unstable circular orbit increases its cross section under infinitesimal radial deformations.
The PS modes with an eikonal limit that we will consider are those with = m or = −m. This includes the = m = 2, n = 0 modes that have the slowest decay rate and that we typically display as orange diamond curves/surfaces (e.g. in Fig. 1 and Figs. 13-14, among others). And these PS modes of the KN BH are those that reduce to the wellknown QNM frequencies of Schwarzschild BH in the limit Q → 0 and a → 0 (typically identified as a dark-red point in our figures) first studied by Chandrasekhar (see Table  V, page 262 of [40]). Therefore, in this subsection we use geometric optics to compute an analytical approximation (to be denoted as ω eikn PS ) for the frequency of these PS modes in the KN background. A similar analysis was originally done in [2,95]. Although, the final analytical formula for the PS QNM frequencies is strictly valid in the WKB limit ∼ |m| → ∞, in practice we find that it matches reasonably well the PS frequencies even for values as small as = |m| = 2. Therefore, the eikonal limit allows us to identify the nature of this QNM family and, furthermore, it provides a check on our numerics.
The geodesic equation, describing the motion of pointlike particles around a KN BH, leads to a set of quadratures. A priori this is perhaps an unexpected result since KN only possesses two Killing fields, K = ∂/∂ t and ξ = ∂/∂ φ . We seem to be one Killing field short of an integrable system. However, there is another conserved quantity − the Carter constant − associated to a Killing tensor K ab , which saves the day [40].
The most direct way to identify this integrable structure is to consider the Hamilton-Jacobi equation [40]: where S is known as the principal function. One can obtain the motion of null particles by noting that, according to Hamilton-Jacobi theory, the principal function and the particle momenta are related by ∂S ∂x µ ≡ p µ and p µ = with τ denoting an affine parameter of the null geodesic. We can then take a separation ansatz of the form (x = cos θ, where θ is the polar angle) where the constants e and j are the conserved charges associated with the Killing fields K and ξ via 15 where the dot (˙) denotes a derivative with respect to the affine parameter τ . Substituting the ansatz (3.3) into the Hamilton-Jacobi equation (3.1) for null geodesics yields a coupled system of ordinary differential equations for R(r) and X(x) (the prime denotes a derivative w.r.t. the argument, r or x, respectively) where Q is a separation constant known as the Carter constant.
We are interested in matching the behaviour of null geodesics with that of QNMs with large values of = |m|, so we can restrict attention to the equatorial plane where x = 0.
where the potential is We are now interested in finding the photon sphere (region where null particles are trapped on circular unstable orbits), i.e. the values of r = r s and b = b s , such that V (r s , b s ) = 0 and ∂ r V (r, b)| r=rs,b=bs = 0. (3.10) From the first equation we get which we insert in the second equation of (3.10) to get a fourth order polynomial equation for r s : 12) where ∆(r) is defined below (2.1) and we are interested in solutions with r s > r + . Alternatively, we can solve (3.10) to get the black hole parameters M and Q that have circular orbits with radius r s and impact parameter b s , namely There are two real roots r s higher than r + which are in correspondence with two PS modes: the co-rotating one (with m = ) that maps to the eikonal orbit with radius r s = r − s and b s > 0 (and that has the lowest |Imω|) and the counter-rotating mode with m = − which is in correspondence with the orbit with radius r s = r + s and b s < 0, with r + s ≥ r − s ≥ r + . The two real roots r ± s higher than r + are displayed in Fig. 2. We can finally compute the orbital angular velocity (also known as Kepler frequency) of the null circular photon orbit, that is simply given by 14) Figure 2. The radii r ± s (with r + s ≥ r − s ≥ r + ) of the two unstable circular orbits in the equatorial plane of the KN black hole that ultimately yield the co-rotating m = (in the r − s case) and the counter-rotating m = − (in the r + s case) PS QNM frequencies in the eikonal limit. For a = 0, one has r + s = r − s , and at (Q,ã) = (0, 1) one has r − s = r + .
where we used (3.6) evaluated at r = r s and b = b s . We can also compute the largest Lyapunov exponent λ L , measured in units of t, associated with infinitesimal fluctuations around photon orbits with r(τ ) = r s . This can be done by perturbing the geodesic equation (3.8) with the potential (3.9) evaluated on an orbit with impact parameter b = b s and setting r(τ ) = r s + δr(τ ). One finds that small deviations decay exponentially in time as δr ∼ e −λ L t with Lyapunov exponent given by One finally obtains the approximate spectrum of the photon sphere family of QNMs in the WKB limit = |m| 1 using [85][86][87][88][89][90][91][92][93] where n = 0, 1, 2, . . . is the radial overtone. This is the eikonal approximation for the PS modes we were looking for. Note that this expression is blind to the spin of the perturbation, i.e. it is the same for scalar and gravito-electromagnetic perturbations. The eikonal analysis, although only based on a geodesic analysis, gives the same result as a (light blue surface) with the actual numerical frequencies (orange points) for co-rotating PS modes with m = = 6, n = 0. The former is given by (3.11)-(3.16) with r s = r − s of Fig. 2 and b s > 0. The brown curve at extremality has Imω = 0 and Reω = mΩ ext H . So, it turns out that m = 6 seems to be already within the WKB validity |m| 1. The dark-red point at (Q,ã) = (0, 0) coincides with the Schwarzschild QNM, ω 1.21200982 − 0.09526585 i, first computed in [40]. leading order |m| = → ∞ WKB analysis of the wave perturbation equations. Although the eikonal frequency is independent on the spin of the perturbation, the higher order frequency corrections in the 1/m WKB expansion should certainly depend on the spin of the perturbation. Recall that (3.16) is strictly valid in the geometric optics limit, |m| 1, with cor-  rections to Imω and Reω of order O (1/|m|) and O (1), respectively. However, Fig. 3 compares (3.16) (light blue surface) with the actual numerical frequency (red dots) of the co-rotating PS modes with m = = 6, n = 0 and already finds an excellent agreement (of course this agreement will improve for m > 6). Moreover, as Fig. 4 demonstrates, the eikonal approximation (3.16) (light blue surface) still provides a reasonably good qualitative approximation to the numerical co-rotating PS modes (orange dots) even in the m = = 2, n = 0 case. Taken together, this identifies the PS QNM family and validates our numerics.
As shown in Fig. 5, an important feature of the eikonal frequency (3.16) (the solid dark-blue line) is that it is in good agreement with Imω = 0 and Reω = mΩ H (dotted brown line) forã ext >ã eikn , but not so forã ext <ã eikn withã eikn = 1 2 (recall thatã ext = 0 andã ext = 1 in the RN and Kerr limits, respectively). This transition point is indeed observed in the numerical data of Fig. 3 (for m = 6) and Fig. 4 (for m = 2), where the values Imω = 0 and Reω = mΩ H are represented by the solid brown lines. As expected, the eikonal quantitative value ofã eikn = 1 2 is not yet a good approximation for m = 2, where numerically we findã 0.360, but it becomes a better approximation as m increases. For example, for m = 6 one hasã 0.463 and for m = 10 one hasã 0.480. We come back to this issue in the discussion of Fig. 11.
For completeness, in Fig. 6 we turn our attention to the counter-rotating PS modes, and compare the eikonal prediction (3.16) (grey surface) with the numerical data (blue/green points) for the counter-rotating PS modes Z 2 m = − = −2, n = 0. Although m = −2 is certainly outside the regime of validity of the geometric optics approximation, |m| 1, it turns out that it gives a relatively good approximation for the qualitative shape of the PS QNM surface (although less than for the m > 0 case). As expected, the quantitative eikonal prediction improves considerably as m grows more negative in the same way as for the m > 0 case. This is illustrated for the PS modes Z 2 m = − = −6, n = 0 modes in Fig. 7.

Near-horizon family of QNMs: analytical formula for the frequencies
Near-extremality, there is a family of KN wavefunctions that are very localized near the horizon and quickly decay to zero away from it. This suggests doing a 'poor-man's' matched asymptotic expansion (MAE) whereby we take the near-horizon limit of the perturbed equations (2.23), which can be solved analytically, to find the near-region solution and then match it with a vanishing far-region wavefunction in the overlapping region where both solutions are valid. 16 In fact, motivated by the result that the near-horizon limit of the extremal KN BH corresponds to a warped circle fibred over AdS 2 (Anti-de Sitter) [96], the perturbations of which can be decomposed as a sum of known radial AdS 2 harmonics, we can attempt to use separation of variables. It turns out that this can indeed be done, and the system of 2 coupled PDEs for {ψ −2 , ψ −1 } separates into a system of two decoupled radial ODEs and a coupled system of two angular ODEs. This is a non-trivial, remarkable property of perturbations on KN.
At extremality, the modes with slowest decay rate (independently of belonging to the NH or PS families or, as we will introduce and discuss later, to the PS−NH family) always approach Imω = 0 and Reω = mΩ ext H and the near-horizon matched asymptotic expansion analysis that we perform below produces a prediction for the frequencies of these modes that will prove to be an excellent approximation near-extremality.
After this preliminary outline, we are ready to formulate and perform in detail our matched asymptotic expansion to find the NH family of QNMs. The near-region is defined as the region r r + −1 1 and its wavefunction must be regular, in particular, at the horizon. The far-region is the zone r r + − 1 σ, with σ = 1 − r − r + being an off-extremality parameter, and its wavefunction must obey the outgoing boundary condition at r → +∞. The two wavefunctions must be simultaneously valid − and thus the free parameters of the two regions must be matched − in the matching region σ r r + − 1 1. We can guarantee that the latter overlap region exists if we take σ − which is our expansion parameter − to be small, i.e. if σ 1 and we are thus near-extremality r − r + . Under these conditions, in the near-region r r + − 1 1 we want to simultaneously zoom in around the horizon and approach extremality. For that we first introduce the dimensionless quantities where recall that r = r − and r = r + are the Cauchy and event horizon locations, respectively, that satisfy ∆ = 0 with ∆ = r 2 − 2M r + a 2 + Q 2 as defined in (2.1). Equivalently, we can also write ∆ = (r − r − )(r − r + ). Equating these two expressions and their derivatives we can express M and Q as a function of (r − , r + , a): From (3.17), one sees that for y 1 one is close to the event horizon and for σ 1 the Cauchy and event horizons are very close, i.e. one is close to extremality. Next, we take the limit σ → 0. From previous works on QNMs of RN, Kerr, KN [1-3, 46, 93, 97-101] and even de Sitter black holes [84,102,103], when we Fourier decompose the modes as e −iωt e imφ , the near-horizon modes are expected to saturate the superradiant bound ω = mΩ H at extremality (this will be further confirmed by our numerical results). Therefore, we expand the frequency about this bound via the redefinition Our task is to find δω. In (3.19) and hereafter, a and Ω H in expressions always refer to their extremal values, a ext and Ω ext H , although we drop the super/subscripts 'ext' for brevity. In these near-extremality conditions, we are ready to find the near-horizon solution of the KN gravito-electromagnetic perturbation equations (2.23). We substitute Next we attempt to separate variables. In the KN system − described by a coupled pair of PDEs − this might seem bound to fail. So it is enlightening to make a small diversion from our exposition to explain the motivation for even considering this possibility. An extremal KN BH has Q = r 2 + − a 2 . Similar to the Kerr case [96], the near-horizon limit of the extremal KN black hole (NHEKN) can be obtained by performing the coordinate in the KN solution (2.1) and taking the limit ε → 0 (recall that x = cos θ). This yields the near-horizon geometry of the extremal KN solution (which also solves the original Einstein-Maxwell equation): Surfaces of constant x are warped AdS 3 geometries; that is they correspond to a circle fibred over AdS 2 (parametrized by T and Z) with warping parameter The isometry group is SL(2, R) × U (1). Consequently, perturbations in NHEKN can be expanded in terms of the AdS 2 harmonics and thus they separate into a radial and angular part. This observation is relevant for our purposes because, returning to the full KN geometry, it suggests that near-extremality and near the horizon the two coupled PDEs for {Σ −2 , Σ −1 } might be amenable to a solution by separation of variables.
What is the meaning of (3.30) and (3.31)? Recall that in the case of the Teukolsky equation describing perturbations in the Kerr black hole [38], it is well known that the so-called Starobinsky-Teukolsky relations relate perturbations with spin s to those with spin −s [40] (see also Appendix of [83]). Thus, one interprets relations (3.30)-(3.31) as being effectively a kind of Starobinsky-Teukolsky relations for the KN perturbations. In this case they relate the wavefunction of spin s = −2 with that of spin s = −1 because the perturbations for these two spins are coupled.
After this long tour, we should recap what we have learned so far. The gravitoelectromagnetic perturbations of the KN black are described by a coupled system (2.23) of two PDEs for {ψ −2 , ψ −1 }. However, if we take its near-horizon limit near extremality, as described in (3.17)-(3.20), we get two near-horizon coupled PDEs for {Σ −2 , Σ −1 } that can be solved assuming the separation of variables (3.23). After using the Starobinsky-Teukolsky−like relations (3.30)-(3.31), we verify that the system indeed separates. We get two decoupled ODEs (3.28) for the radial wavefunctions Y 1 (y) and Y 2 (y). (This decoupling reflects the fact that in NHEKN the radial perturbations are exactly described by the AdS 2 harmonics as explained below (3.22)). Once we know the separation constant λ 2 , and thus λ 1 via (3.31), these two ODEs (3.28) can be solved independently as a quadratic eigenvalue problem for δω (for a given m). On the other hand, the angular equations for X 1 and X 2 -given by the curly brackets of (3.24) after using (3.30)-(3.31) -do not decouple. Thus we have to solve this coupled system of two ODEs (that are independent of δω) to find the eigenvalue λ 2 (and thus λ 1 given in (3.31a)). This can be done numerically as we discuss later. But we can also solve this coupled ODE system analytically in a large m WKB expansion. This is what we do next.
Substituting (3.25)-(3.27) and (3.30)-(3.31) into the curly brackets expressions of (3.24), we find that (3.24a) is a second order ODE for X 1 (hereafter we denote this as the 'first' angular equation) that also depends on X 2 and X 2 but not on X 2 . Similarly, (3.24b) is a second order ODE for X 2 (henceforth denoted as the 'second' angular equation) that also depends on X 1 and X 1 but not on X 1 . If we redefine where K 12 was first introduced in (3.30a), we can solve the equation for χ 1 to express χ 2 = χ 2 (χ 2 , χ 1 , χ 1 , χ 1 ). We substitute this relation and its derivative into the second angular equation to get a differential equation that can be solved to express χ 2 = χ 2 (χ 1 , χ 1 , χ 1 , χ 1 ). Substituting this back in the first angular equation we end up with a fourth order differential equation for χ 1 that no longer depends on χ 2 . This is a non-polynomial eigenvalue problem for χ 1 and λ 2 ; recall (3.31a). Perhaps remarkably, this fourth order differential equation can be solved analytically in a |m| 1 WKB expansion to find χ 1 and λ 2 . We substitute the WKB ansatz into the fourth order ODE and solve it order by order in a standard large m expansion, requiring that the solution is regular at x = ±1. The leading order WKB wavefunction is ) and the separation constant λ 2 is given by (3.33) with WKB coeficients Of course, now that we have the eigenpair λ 2 and χ 1 (x) (in the WKB approximation) we can straightforwardly obtain λ 1 and χ 2 (x) using (3.31a) and the aforementioned relation χ 2 = χ 2 (χ 1 , χ 1 , χ 1 , χ 1 ). This terminates our WKB analysis of the angular equations. We can also solve numerically the coupled pair of angular ODEs for X 1 and X 2 to check that the WKB result is indeed a good approximation, even for m = 2. For m ≥ 2, regularity at x = ±1 requires that we keep the X 1 , X 2 solution that behaves as (1 − x) 1 2 (s+m) at x = 1 and as (1+x) 1 2 (−s+m) at x = −1 where s = −2, −1 for X 1 , X 2 , respectively. We can impose these boundary conditions straightforwardly if we introduce the field redefinition and solve the two coupled second order ODEs for smooth Q 1 and Q 2 and the eigenvalue λ 2 , after using (3.31a). As explained above, we use a Newton-Raphson algorithm with pseudospectral discretization [77]. In Fig. 8 we compare the WKB result (3.33)-(3.35) with the numerical λ 2 . We see that, as expected, for large m, m = 10 (left panel), there is perfect agreement between the WKB result (continuous dark-blue curve) and the numerical result (blue circles). However, as the right panel demonstrates, the WKB approximation (continuous dark-green line) proves to be a good approximation to the exact result (green circles) already for m = 2. Also note that asâ increases fromâ = 0, λ 2 changes sign from positive to negative. This fact will be important later.
Having solved the angular equations we can now focus our attention on the radial ODEs (3.28)- (3.29). Recall that we can solve one of these, e.g. (3.28b) for Y 2 , and the solution for Y 1 is then straightforwardly given by the Starobinsky-Teukolsky differential map (3.30a). Further recall that (3.28b) is a quadratic eigenvalue problem in δω. This ODE turns out to be a standard hypergeometric equation with most general solution given by where c 1 , c 2 are arbitrary integration constants. At the event horizon, y = 0, this solution behaves as Y 2 | y=0 ∼ c 1 y i(1+â 2 )δω + c 2 y 1−i(1+â 2 )δω . Regularity in ingoing Eddington-Finkelstein coordinates at the future event horizon requires that we set c 1 = 0 to eliminate the outgoing modes. On the other hand, far away from the horizon, i.e. at large y, the regular solution at the horizon behaves as Assume for now that λ 2 > 0. From Fig. 8, this happens whenâ ext = 1 −Q 2 is small, which occurs for largeQ. For λ 2 > 0, at large y, the solution y diverges. 17 17 Note that the metric components that must be a regular 2-tensor behave as y In the context of a matched asymptotic expansion, the large behaviour of the nearregion (near-horizon) solution (3.38) must now be matched with the far-region solution (near extremality). As explained at the beginning of this section, we expect the nearhorizon modes we are looking into to have wavefunctions that die-off very quickly away from the black hole horizon (near extremality). This will be confirmed by our numerical analysis. Therefore, as a first approximation − that we henceforth call a 'poor man's matched asymptotic expansion (MAE) − 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 (3.38) with Y 2 | far 0. 18 It is important to emphasize that this is an ansatz or educated guess that we cannot argue for in a formal mathematical away that goes deeper than the above heuristics. It is ultimately only validated a posteriori by the fact that the final quantization agrees with the numerical results for the frequency spectra (indeed, Y 2 | far is never exactly zero and thus a small component of the divergent term in (3.38) should be used in a proper matching). This ansatz requires that we kill the solution y in (3.38) that diverges for large y. Since Γ(−n) → ∞ for n ∈ N 0 , this is the case if we quantize the frequency correction to be such that the argument of the gamma function in the denominator of the divergent term is a non-positive integer n: δω mâ Inserting this frequency correction into the frequency expansion (3.19) one gets the final expression for the frequency in units of r + :ω = mΩ H + σ δω. We can now convert this into units of M by multiplying this expression by M/r + (sinceωM/r + =ω = ωM ) and expanding it in terms of σ while keeping terms only up to O(σ) (since all our analysis is valid only up to this order). This yields the frequency quantization for the near-horizon (NH) QNMs which can be written as: 1 + 2n 1 +ã 2 − −λ 2 (m,ã) 4(1 +ã 2 ) 2 + O σ 2 , n = 0, 1, 2, 3, · · · (3.40) whereã in this expression must be evaluated at extremality, i.e.ã =ã ext , the off-extremal parameter σ is defined in (3.17), and we have defined √ z to be such that Re( √ z) > 0 (Im( √ z) > 0) for positive (negative) values of z. How good an approximation is (3.40)? It is in excellent agreement with the numerical NH frequencies near extremality, as will be discussed in the analysis of Figs. 13-14. This is further illustrated in the left panel of Fig. 9 where we take a KN BH family with Q/r + = 0.99 and compare the numerical results (green circles) with the red curveω MAE given by (3.40). It turns out that for very largeQ the agreement is excellent not only nearextremality but also far away from it down to smallâ. So much so that we can basically use (3.40) for any astrophysical application that requires the knowledge of the dominant NH 18 Ideally, we would also solve the far-region equations to obtain the sub-leading far-region solution but in the KN background we cannot do this analytically. frequencies for 0.99 < Q/r + < 1, say. Accordingly, the reader will later find that we have not felt the need to collect numerical data in the window 0.99 < Q/r + < 1 in our plots: see e.g. the gap between the green surface and extremal brown curve in Fig. 16 and the similar gaps in Figs. 18−20. Naturally, as we decreaseQ the approximation (3.40) becomes increasingly less accurate when we move away from extremality. This is demonstrated in the right panel of Fig. 9 where we do the comparison between (3.40) (red line) and the numerical data for Q/r + = 0.95.    In the final steps leading to (3.40), we assumed that λ 2 > 0. From Fig. 8, this happens whenã ext = 1 −Q 2 is small, which occurs for largeQ, as is the case in Fig. 9. This also includes the extremal RN limit, (Q,ã) = (1, 0) in which case (3.40) reduces to the expression first found in [2]. However, nothing impedes us from extending the application of (3.40) also to the case where λ 2 < 0. From Fig. 8, this happens for largeã ext = 1 −Q 2 , and thus for smallQ. In particular, this includes the extremal Kerr limit, (Q,ã) = (1, 0). Interestingly, when λ 2 < 0 (unlike for λ 2 > 0), one is effectively in a region of the parameter space where the PS family terminates at Imω = 0 and Reω = mΩ ext H and, quite importantly, it dominates over the NH family. Therefore, by construction (3.40) should be able to capture (also, or in this case) the frequency of the dominant PS modes near extremality. And indeed it does so, as illustrated in Fig. 10 where we compare (3.40) (black curve) against the numerical PS frequency (orange diamonds) for the KN families with Q/r + = 0.5 (left panel) and Q = 0 (right panel). The latter case is the Kerr solution, where (3.40) reduces to the expression first found in [46,100]. Thus,ω MAE in (3.40) (also) provides an analytical approximation for PS modes when they approach Imω = 0 at extremality that complements, and is independent of, the eikonal analytical approximation ω eikn PS given in (3.16). It has the added value of being very accurate near extremality already for m = 2 (i.e. well outside the |m| 1 eikonal regime of validity). Interestingly, the approximation (3.40) for the PS modes improves asQ decreases, as can be inferred from the two cases presented in Fig. 10.
Altogether, and to summarize, we find that (3.40) is an excellent approximation for the dominant modes (which always approach Imω = 0 and Reω = mΩ ext H at extremality) when we are close to extremality, i.e. when a/ã ext 1, independently of the QNM family that dominates, as best illustrated in Figs. 9-10. For largeQ the dominant modes are the NH modes and (3.40) describes them. For smallQ the dominant modes are instead the PS modes and (3.40) also describes them. This might sound a bit puzzling: how can it be that the near-horizon MAE analysis captures sometimes the PS modes? This is because, away from the RN limit, the distinction between the PS and NH families becomes less clean and actually the dominant QNM family is better described by a combination of the PS and NH modes (that we will denote as a PS−NH family) due to the phenomenon of eigenvalue repulsion. This statement will be clarified and made accurate when discussing the results of Figs. 13−14 so we postpone further discussion till then.
In the analysis of the eikonal expression (3.16) and associated Fig. 5, we have already pointed out that when we are at extremality, e.g., if we place ourselves on the extremal brown curve of Fig. 5 (or of Fig. 16) and move along it fromã ext = 1 down toã ext = 0 (or, equivalently, fromQ ext = 0 toQ ext = 1), there is a critical rotationã ext =ã (or, equivalently, a critical chargeQ = 1 −ã 2 ). Forã <ã ext ≤ 1 (i.e. 0 ≤Q ext <Q ) the PS family terminates at Imω = 0 and Reω = mΩ ext H at extremality (e.g., in the Kerr limit whereã ext = 1), but it fails to do so otherwise (e.g., in the RN limit whereã ext = 0 and Q ext = 1). Interestingly, we find that this transition point turns out to be very well approximated (if not exactly given) by the point where the separation constant λ 2 (m,ã ext ) in (3.40) vanishes: λ 2 (m,ã NH ) = 0. Forã ext <ã NH (or equivalently,Q ext >Q NH ) one has λ 2 > 0 and forã ext >ã NH we have λ 2 < 0. To get the accurate values forã NH − which are displayed as black 's in Fig. 11 − we use the numerical solution for λ 2 (m,ã ext ) as displayed in Fig. 8. Alternatively, since λ 2 has the WKB expansion (3.33b) and (3.35), we can use it to find a NH | WKB orQ NH | WKB , yielding Using our numerical data for λ 2 (Fig. 8) (3.41) is expected to be accurate only as m → ∞. To confirm this, we compute these critical rotations for m = 2 to m = 10 and Fig. 11 shows thatã NH | WKB as given by (3.41) (the solid blue line) indeed approaches increasingly the value ofã NH (the black 's) as m grows, with excellent agreement already for m = 10 (or even m = 6). Further note that as m increases,ã NH andã NH | WKB approach from below the eikonal value ã eikn = 1 2 (or from above the eikonalQ eikn = 1 − (ã eikn ) 2 = √ 3/2 0.866025) discussed in Fig. 5 (see its red ). So our numerical results indicate that the critical rotation/charge {ã ,Q } seem to be given to very good accuracy by the values {ã NH ,Q NH } discussed above and displayed in Fig. 11. This is further demonstrated in Fig. 12. In these plots we show the imaginary and real part of the PS frequency as a function of the rotation at extremality,ã ext , for Z 2 = m = 2, n = 0 modes. The black atã ext =ã NH 0.360 (i.e.Q ext =Q NH 0.932) is the point already displayed in Fig. 11. The set of black squares displayed only forã ext ∈ [0, 0.24] describe data we obtained by solving the gravito-electromagnetic PDEs directly at extremality (numerically, it is very hard to extend the computation for higherã ext ; recall that at extremality we have a degenerate horizon and thus the boundary conditions differ from the non-extremal case). On the other hand, the auxiliary grey line that joins these black squares and connects to the black point atã ext ∼ 0.360, is an interpolation curve built from the black square and points. Finally, the PS modes closest to extremality that we found using our non-extremal code are identified with orange diamonds (withã ext 0.2 since it is hard to obtain data whenã ext → 0). For 0 <ã ext <ã 0.360 they are just below the interpolation grey line. Altogether this indicates that PS modes indeed terminate at the grey interpolation line for 0 ≤ã ext ≤ã . On the other hand, forã <ã ext ≤ 1, the grey horizontal line displayed in Fig. 12 has Imω = 0 and Reω = mΩ ext H . The orange diamonds in this region are the closest PS modes we obtained using the non-extremal PS numerical code; these points are at 99% of extremality. Again we see that they indeed approach the grey horizontal line. To find even further approach we need to extend our collection of data closer to extremality, say up to 99.9% of extremality. We did this for a few constant charge families (not shown) to confirm it is indeed the case (and these are very accurately described byω MAE in (3.40) as discussed previously). 4 Eigenvalue repulsions (also known as level repulsions or avoided crossing) The analytical analyses of Section 3 allowed us to find corners of the 2-parameter space of KN where we can obtain good analytical approximations for the QNMs of KN. Importantly, they gave us evidence for the existence of not one but two main families of QNMs: the photon sphere and near-horizon families. These are distinct families because the analytical analyses reveal different origins: the PS family is associated with properties of null orbits in the eikonal limit, while the NH family is related to modes whose wavefunction is very localized about the horizon near extremality. Our numerical search of QNMs, whose findings will be presented in this section and in Sections 5−6, confirm that KN indeed has two families of QNMs and not more, and all the numerical QNM frequencies are well approximated by (3.16) and/or (3.40) in the regimes where the latter are valid. However, the distinction between the two QNM families of KN becomes very fuzzy as we move along the 2-parameter space of KN. In the RN limit (a = 0, Q = 0) this distinction is very sharp: one of the families is the PS family well approximated by (3.16) and the second one is the NH family well described by (3.40); recall Fig. 1. But when we switch on the rotation and allow it to increase we find that the PS and NH families lose their individual identities. Instead branches of these two families combine with each other to produce a combined family that we can appropriately call PS−NH modes (and their radial overtone families). This occurs because the KN spectra has a novel phenomenon that is special to the KN QNM system (i.e. present neither in Kerr nor RN), namely eigenvalue repulsion between QNM families. In subsection 4.2 we will describe in detail this phenomenon in the KN QNM spectra. Although, in the context of black hole QNMs eigenvalue repulsions are particular to KN (see also footnote 3), such a feature is common in some eigenvalue problems, notably: 1) in solid state physics where e.g. it is responsible for energy bands/gaps in the spectra of electrons moving in certain Schrödinger potentials, and in 2) in quantum mechanical eigenvalue systems with the so-called avoided crossing phenomenon. Therefore, before discussing eigenvalue repulsions in KN, in subsection 4.1 we will present a simple textbook example of eigenvalue repulsions that will allow us to understand from first principles what occurs in the KN eigenfrequency spectra.

Complexified eigenvalue repulsion
Eigenvalue repulsion is a phenomenon that occurs in simple quantum mechanical models (albeit it can also occur in classical physics, most notably when two levels of a classical harmonic oscillator are coupled). In quantum mechanics this phenomenon is also known as the Wigner-Teller effect, avoided crossing or level repulsion [64,65]).
To explain the similarities and differences between what is observed in standard quantum mechanics textbooks and the phenomenon that we observe numerically in the QNM spectra of KN, we will start by reviewing the simplest textbook example exhibiting avoided crossing (see for instance §79 of [64] and/or §IV.C of [65])).
For concreteness consider a two-level system with Hamiltonian H 0 , orthonormal eigenstates |ψ i and energy levels E i , so that Let us imagine perturbing H 0 with an interaction W , such that the full Hamiltonian is given by H = H 0 + W . Here W can be thought as coupling the two eigenstates ψ i . In the {|ψ 1 , |ψ 2 } basis, the coupling is given as a 2 × 2 matrix W with entries W ij ≡ ψ i |W |ψ j .
In the {|ψ 1 , |ψ 2 } basis the perturbed Hamiltonian matrix H ij = ψ i |H|ψ j can be written as Self-adjointness of the perturbed Hamiltonian then demands H to be Hermitian, and thus W 21 = W 12 , where the bar denotes complex conjugation.
It is a rather standard exercise to diagonalise H given in (4.2) and find that the eigenvalues of the perturbed Hamiltonian are: Einstein summation convention on the last term). Eigenvalue crossing (i.e. E − = E + ) will only occur if the argument of the square root vanishes. Since the argument of the square root is given by a sum of two positive definite terms, we must demand each to be zero separately: Let us now imagine that W is a function of a number of real parameters, say N . Since we have two conditions to be satisfied in order for crossing to occur, we expect that crossing can only happen over a subspace of the N real variables parametrised by N − 2 real variables. 19 Except at this special subspace, (4.3) predicts that eigenvalues do not cross under the effect of perturbations W (since E − < E + for W 12 = 0). This is known as avoided crossing.
However, the case at hand (QNMs of KN), is more complicated than this standard textbook example because the perturbation operator is not self-adjoint. However, we shall see that progress can nevertheless be made to understand the properties of its intricate QNM spectra in terms of avoided crossing. Let us denote by L 0 the operator whose eigenspectrum yields the QNM spectrum of a RN black hole, which is non-degenerate: see Fig. 1 20 . Let us label such QNMs as {ψ i , ω i } with i = 1, 2 (in the simplest case, we should regard these as the two slowest decaying QNMs for a given value of Q/M as shown in Fig. 1). We would like to investigate what will happen to these two QNMs as we turn on J/M 2 . The operator governing the eigenspectrum will change to L = L 0 + K, so that K = 0 at J = 0.
For quasinormal modes, L is not self-adjoint, but one can nevertheless introduce a nondegenerate bilinear form ·|· with respect to which the ψ i are orthogonal [104]. However, in general, ·|· will be complex. We will choose the normalisation of the ψ i to be such that ψ i |ψ j = δ ij 21 . Note that we cannot choose ψ i |ψ j = δ ij since it could well be that ψ i |ψ i = 0.
As with the Hermitian case, we define L ij = ψ i ||ψ j , which leads to the perturbed matrix with K ij = ψ i |K|ψ j . L can also be straightforwardly diagonalised as whereω i = ω i − K ii (with no Einstein summation convention on the last term). Eigenvalue crossing will only occur if the argument of the square root vanishes. Unlike the Hermitian case, this time this gives only one condition Let us now imagine that K depends on N real parameters. Since the condition (4.7) is in general complex, it provides a restriction on two of the N parameters. This means eigenvalue crossing can only occur on a N − 2 subspace, just as in the Hermitian case. This is the reason why we need also at least two real parameters to see avoided crossing in the non-Hermitian case. In the black hole context, this justifies why we can see this phenomenon in Kerr-Newman [3], RN-dS [63], Myers-Perry-dS [62], but not in RN or Kerr black holes. The analysis above also shows that level crossing (in the complex frequency plane) will only occur at most at a point in the full Kerr-Newman space of parameters (which has N = 2 adimensional parameters, namely Q/M and J/M 2 ). Our numerical analysis of the KN QNM spectra (mainly of of Sections 3.2 and 4.2) provide us with strong evidence to conjecture that this level crossing point lies precisely at extremality when the PS modes reach Im(ω) = 0. This is the point in Fig. 12 (in the case of n = 0 PS and NH modes). This conjecture is backed up not only by our numerical studies, but also by our approximate analytic form of the near-horizon matching asymptotic expansion frequency (3.40), which has the same elements as (4.6), with −λ 2 (m,ã) playing the role of (ω 1 −ω 2 ) 2 4 + K 12 K 21 . The study of level crossing for non-Hermitian systems remains an active topic of research particularly when more than two-levels are considered (see [105] for an excellent topical review on the subject). For instance, in (4.7) we could have demanded that the real (imaginary) part vanishes, but let the imaginary (real) part be arbitrary. This would lead to avoided crossing in the imaginary (real) part, but would allow for crossing in the real (imaginary) part. This example shows that avoided crossing for non-Hermitian matrices can indeed be a richer phenomenon than its Hermitian cousin.

Eigenvalue repulsions in the frequency spectra of KN
Perhaps surprisingly at first sight, but certainly not after the discussions in subsection 4.1, in the numerical search of the KN QNM spectra we find eigenvalue repulsions between the two distinct families of QNMs of the system. These eigenvalue repulsions are unique to the KN QNM system since they are not observed in the spectra of Schwarzschild, RN nor Kerr black holes (for reasons that were understood in subsection 4.1). Our strategy to describe and discuss further these eigenvalue repulsions is the following. In Figs. 13−14 we display a series of panels. Each one of them plots the imaginary part of the dimensionless frequency, Im(ωr + ), as a function of the dimensionless charge,Q = Q/r + , at fixed a/a ext . We choose to use units of r + since some curves change too much in a small range of charge if we use units of M . Different plots of this series are for different values of fixed a/a ext . Namely, moving from top-left into bottom-right panels of Figs. 13−14 we have fixed a/a ext = 0, 0.38, 0.39, 0.5, 0.8, 0.86 and 0.96 (see legend on top of each panel). So we start at the RN family with a = 0 and progressively increase a/a ext till we reach a KN black hole family parametrized by 0 ≤Q ≤ 1 where the whole KN family is at 96% of extremality. We have chosen these particular a/a ext cases because they are representative of what happens to the system in a window of a/a ext centred at the given a/a ext . When we move to the next panel a new major feature appears that justifies introducing a new plot to illustrate it. We only display the imaginary part of the frequency. This is because the plots for the real part of the frequency are not very illuminating since the curves for the different modes quickly become very close to each other as we approach extremality i.e, as a/a ext increases. We will add comments about the real part of the frequency whenever appropriate and at the end of this section.
We can now describe in detail the content of each plot in Figs. 13−14. In the topleft panel of Fig. 13 we start with the RN black hole (a = 0). This describes what happens to the system with a = 0 but it is also representative of small rotation cases with a/a ext below 0.38. We plot the first two overtones (n = 0, 1) of the PS QNM family (that we denote by PS 0 and PS 1 or, more generically, as PS n modes) and the first two overtones (n = 0, 1) of the NH QNM family (denoted as NH 0 and NH 1 or simply as NH n modes). The PS 0 and PS 1 curves are described by orange diamonds and dark-red triangles, respectively. In the Schwarzschild limit (Q = 0), the PS 0 family reduces to the dark-red disk ω r + = 0.74734337 − 0.17792463 i while the PS 1 curve reduces to red square ω r + = 0.69342199 − 0.54782975 i, first computed by [40,41]. On the other hand, the NH 0 and NH 1 families are the green circle and blue square curves, respectively (not shown: forQ < 0.85 these curves plunge quickly to lower Imω). Note that this plot contains the same PS 0 and NH 0 information as the one of Fig. 1 (although here we use units of r + instead of M ). Moreover, w.r.t. Fig. 1, in the top-left panel of Fig. 13 we also   In the RN case, there is an unambiguous QNM family classification: the orange diamond (dark-red triangle) curve is the n = 0 (n = 1) PS family which reduces to the dark-red disk ω r + = 0.74734337 − 0.17792463 i (red square ω r + = 0.69342199 − 0.54782975 i) in the Schwarzschild limit [40,41]. The green circle (blue square) curve is the n = 0 (n = 1) NH family (not shown: forQ < 0.85 these curves extend to lower Imω). In the middle panels one observes eigenvalue repulsions unique to the KN QNM spectra. In the RN case, we also show the frequencỹ ω MAE given by (3.40) for n = 0 (black curve) and for n = 1 (magenta curve).  display the near-horizon matched asymptotic expansion frequencyω MAE as given by (3.40) for n = 0 (solid black curve) and for n = 1 (solid magenta curve); these are better seen in the inset plot where one finds that (3.40) gives the correct slopes near-extremality at Q 1. As emphasized already in subsection 3.2, these analyticalω MAE are in excellent agreement with the numerical NH QNM frequencies, as long as we are near extremality (which for RN occurs atQ = 1). Actually this time we demonstrate that (3.40) is an excellent approximation (near extremality) not only for the first overtone NH 0 but also for NH 1 (and higher overtones n although not shown). A major feature of this a = 0 plot is that the PS n and NH n curves are very well defined and clearly distinct from each other, with the PS n frequencies well approximated byω eikn PS in (3.16), and the NH n frequencies in excellent agreement withω MAE as given by (3.40). It is also important to emphasize that the imaginary part of the PS 0 and NH 0 curves (in particular) cross each other but, as best displayed in the right panel of Fig. 1, this is not the case for the real part of the frequency. This will be a common feature in all the cases displayed in Figs. 13−14: whenever we see crossing between two curves describing the imaginary part of the frequency there is no crossing between the curves that represent the real part of the frequency.
As far as it is possible, we will keep the same colour/shape code for the PS n and NH n QNM families displayed in the top-left panel as we move to the other plots with increasing a/a ext . However, at a certain point we will no longer be able to assign the PS or NH nomenclatures to the QNM curves of the system.
As we switch on a and increase a/a ext , the QNM spectra remains similar to the one on the top-left panel but the PS 1 (dark-red triangles) and NH 0 (green circles) curves start getting deformed in the region where they intersect as a simple crossover in the imaginary part. It is as if each of these curves starts feeling the presence of the other and they start interacting. This is particularly seen in the top-right panel of Fig. 13 for a/a ext = 0.38. Then, increasing a little bit the rotation, at a/a ext = 0.39 (bottom-left panel of Fig. 13) a dramatic new feature occurs. The 'old' PS 1 (by 'old' we mean w.r.t. the previous plot or, ultimately, w.r.t. the a = 0 plot) dark-red triangle curve breaks into two pieces, and the same occurs for the 'old' NH 0 curve. This occurs forQ ∼ 0.875 as best seen in the inset plot. Not less remarkably, the left-branch (Q 0.875) of the 'old' PS 1 curve merges with the right-branch (Q 0.875) of the 'old' NH 0 curve. That is to say, the PS 1 and NH 0 families lose their individual identity and they combine into what we now can call the PS 1 −NH 0 family of QNMs. Similarly, the left-branch of the 'old' NH 0 (Q 0.875) curve joins with the right-branch of the 'old' PS 1 (Q 0.875) curve to form together a new QNM family that we denote as the NH 0 −PS 1 family of QNMs. These breakups and subsequent mergers are even more surprising because they glue two sub-families that were, for lower rotations, assigned different radial overtones n. At this rotation parameter we can say that we have 4 families of QNMs (from top-left to bottom-right): the PS 0 , the PS 1 −NH 0 , the NH 0 −PS 1 and the NH 1 .
Altogether, these features and frequency gaps are characteristic of the phenomenon of eigenvalue repulsion that we discussed in subsection 4.1. In particular, in the breakup region, there is a a 'frequency gap' between the new PS 1 −NH 0 and NH 0 −PS 1 curves. This 'frequency gap' is zero exactly at the breakup rotation (somewhere in the window a/a ext ∈ [0.38, 0.39]), and then it grows as a/a ext increases. This is what is seen e.g. when we move to a/a ext = 0.5 case shown in the bottom-right panel of Fig. 13. In this plot we see that a further eigenvalue repulsion episode happened in the window a/a ext ∈ [0.39, 0.5]. Indeed, the NH 0 −PS 1 curve (green circles plus dark-red triangles) broke up aroundQ ∼ 0.91 and the same happened to the NH 1 curve (blue squares). The left-branch of the 'old' NH 0 −PS 1 curve is now merged with the right-branch of the 'old' NH 1 curve to form what we can call a NH 0 −PS 1 −NH 1 family of QNMs. Simultaneously, the left-branch of the 'old' NH 1 curve (blue squares) is now merged with the right-branch of the 'old' NH 0 −PS 1 curve (or with a portion of the even 'older' PS 1 curve since it only contains dark-red triangles) to form what we can call a NH 1 −PS 1 curve.
So far the original PS 0 family escaped eigenvalue repulsion phenomena, but this changes when we keep increasing a/a ext even further as seen in Fig. 14 (in this figure we drop the subdominant NH 1 −PS 1 curve). Indeed, at a/a ext = 0.8 we already notice that the PS 0 curve (orange diamonds) and NH 0 portion (green circles) of the PS 1 −NH 0 curve are getting deformed by each other in the region where they intersect as a simple crossover. Again, it is as if each of these curves feels the presence of the other and reacts to the interaction (see the inset plot). This is similar to the eigenvalue repulsion observed before between the PS 1 and NH 0 modes and, inevitably, the PS 0 and the PS 1 −NH 0 curves break up in the window a/a ext ∈ [0.8, 0.86]. Indeed, in the top-right panel of Fig. 14, we see that at a/a ext = 0.86 these two curves break atQ ∼ 0.93. The left-branch (orange diamonds) of the 'old' PS 0 curve merges with the right-hand branch (green circles) of the NH 0 portion (green circles) of the PS 1 −NH 0 curve to produce what we denote as the PS 0 −NH 0 family of QNMs. At the sameQ ∼ 0.93, the left-branch of the 'old' PS 1 −NH 0 is now merged with the right branch of the 'old' PS 0 curve to give birth to what we call a PS 1 −NH 0 −PS 0 family. Similar eigenvalue repulsions keep occurring when we increase a/a ext towards extremality. For example, already very close to extremality, namely at a/a ext = 0.96, the two most dominant QNM families are shown in the bottom panel of Fig. 14 (we do not show data for even higher overtones). Here, we identify the PS 0 −NH 0 curve already observed in the previous plot. This is the family that has the lowest |Imω| for allQ. Additionally, we see that the 'old' PS 1 −NH 0 −PS 0 curve of the a/a ext = 0.86 broke again (aroundQ ∼ 0.94) and merged with the right branch of the 'old' NH 1 (blue squares) to form a four colour PS 1 −NH 0 −PS 0 -NH 1 curve (see inset plot).
To conclude by summarizing the key aspects of our findings, the first plot of Fig. 13 together with the last plot of Fig. 14, are those that probably best illustrate the main conclusion of our study. There is no doubt that a = 0, the RN black hole, has two clearly distinct families of QNMs: the PS and NH families, together with overtones for each of them (first plot of Fig. 13). Here, the PS n frequencies are well approximated byω eikn PS in (3.16), and the NH n frequencies are in excellent agreement withω MAE as given by (3.40). However, as the rotation increases, several eigenvalue repulsions progressively appear that increasingly break and combine pieces of the 'old' PS n and NH n curves. Very close to extremality, we end up with a QNM landscape that is definitely very different from the RN one. Indeed, as best illustrated in the last plot of Fig. 14, instead of having the PS n and NH n curves, one now has what we can simply call the 'PS−NH' family and its radial overtones (with higher |Imω|). Interestingly, the near-horizon matched asymptotic expansion frequencyω MAE given by (3.40) describes accurately this PS−NH family (and its overtones) for the whole range ofQ at a fixed a/a ext that is close to extremality. Indeed, in the bottom panel of Fig. 14, the solid black curve describes (3.40) with n = 0 and the solid magenta line represents (3.40) with n = 1. And these match very well the numerical frequencies for the n = 0 and n = 1 PS−NH modes, respectively. This is a conclusion that we had already reached when discussing (3.40) and Figs. 9-10 of subsection 3.2. Notice, that this matching betweenω MAE and the numerical data includes the region of the QNM curve that we can trace back as descending from the RN PS modes (i.e. the orange diamond section of the n = 0 PS−NH curve), in agreement with the discussion of the extension of (3.40) to negative values of λ 2 and associated Fig. 10 that we had in subsection 3.2. In particular, this means that the PS−NH overtone curves (including the two shown in the in bottom panel of Fig. 14) approach Im ω = 0 and Re ω = mΩ ext H as a/a ext → 1 for any value ofQ.
In the analysis of this section we have not discussed much the behaviour of the real part of the frequency. This is because nothing of relevance happens to this quantity as we evolve though the 2-parameter space of KN black holes away from the level crossing point that occurs in the imaginary part. Take for example the n = 0 PS and NH modes. As we move away from the level repulsion point at {â,Q}| ext = {â ,Q } {0.360, 0.932}, during a good neighbourhood the real part of these two modes is very similar (parallel to each other) but they do not cross. Then, sufficiently far away from the level crossing point the two Reω surfaces become clearly distinct. These properties will be observed in the right panel of Fig. 15. It turns out that the eigenvalue repulsions induce strong effects at the level of the imaginary part of the frequencies but leave no (notably visible) imprint on the real part of the frequencies. In more detail, whenever there is crossing between two curves describing the imaginary part of the frequency there is no crossing between the curves that represent the real part of the frequency; that is, the crossing in the imaginary part of the frequency never extends to the full complex frequency plane, with one exception. For each pair of modes, this exception occurs when we approach the particular extremal KN black hole with {â,Q}| ext = {â ,Q }.
To summarise, for definiteness consider again the PS 0 and NH 0 pair of modes. In this case, the star point at extremality has {â ,Q } {0.360, 0.932} and is represented with a in Fig. 12. At this point, the PS and NH modes both have Im ω = 0 and Re ω = mΩ ext H . That is, they have the same complex frequency and, as discussed in subsection 4.1, this is the only level crossing point of the system. As we move away from this point, avoided crossing effects emerge and Figs. 13−14 illustrate that these repulsion effects can be strong and induce intricate features in the behaviour of the Imω curves (but not in the Reω curves) in a neighbourhood of the level crossing point but they become unnoticed far away from this point.

Full frequency spectra of the QNMs with slowest decay rate
We have done a fairly good survey (having in mind the associated computational cost; more in Section 6) of several gravito-electromagnetic QNMs of KN and we conclude that, as expected, the modes that have the slowest decay rate are those that are the {Q, a} = {0, 0} extension of the Schwarzschild mode that Chandraseckar classifies as the Z 2 , = m = 2, n = 0 mode; see Table V, page 262 of [40] and associated discussion. These are also the modes we discussed in Section 4.1, together with the n = 1 overtone of the same family (which was first studied by Leaver [41]).
Therefore, before doing a general survey of other modes of interest in Section 6, in this section we display the QNM spectra of the Z 2 , = m = 2, n = 0 modes and the n = 1 overtones (the latter will allow us to complement or even complete the analysis of Section 4.2). Unlike in Section 4.2 where we made a judicious choice of 2-dimensional plots at fixed a/a ext to exhibit and explain eigenvalue repulsions, in this section we plot the QNM frequencies as a function of the full 2-dimensional parameter space of KN. As discussed previously, we can take these 2 dimensionless parameters to be {ã,Q} ≡ {a/M, Q/M } or {â,Q} ≡ {a/r + , Q/r + }. From an astrophysical perspective, it is appropriate to work in units of M and this is how we present many of our physical results, in particular the frequency ωM . However, in practice we have scanned the 2-dimensional parameter space in units of r + : typically (except when we needed a finer grid to study a particular feature), we divided our numerical grid for {Q,â} ≡ {Q/r + , a/r + } with 100×100 points with 0 ≤Q ≤ 1 and 0 ≤â ≤â ext withâ ext = 1 −Q 2 . This is because some features of the QNM spectra (e.g. the crossovers or eigenvalue repulsions between modes) occur in small windows of (Q/M, a/M ) which translate into wider windows of (Q/r + , a/r + ). For this reason, some of the fine details of the frequency spectra that we discuss in this section are better displayed if we present our results in 3-dimensional plots {Q/r + , a/r + , ωM }. In the figures of this section, the left panel always displays the imaginary part of the frequency, Im(ωM ), while the right panel plots the real part of the frequency, Re(ωM ). In Fig. 15 we present the raw data that we collected for the two QNM families that have the slowest rate, namely the families that we identify with the PS 0 and NH 0 modes in the RN limit and their n = 1 overtone cousins PS 1 and NH 1 (for the Z 2 = m = 2 modes). Namely, the orange and the green surfaces are the n = 0 PS (PS 0 ) and n = 0 NH (NH 0 ) families, respectively. On the other hand, the dark-red surface and the blue surface describe the n = 1 PS (PS 1 ) and n = 1 NH (NH 1 ) families, respectively. Thus, we are using the same colour code that was employed in Figs. 13−14 of section 4.2. The solid brown curves are at extremality. They are parametrized byâ = a ext = 1 −Q 2 and have Imω = 0 and Reω = mΩ ext H (we will use the same colour code for this extremal curve in all other 3-dimensional plots where this curve plays a relevant role). Note that the NH 0,1 surfaces have a very large slope and plunge into very large negative Im(ωM ) as we move away from theâ =â ext (Q) extremal curve or, in the RN case, away from theQ = 1 extremal solution. Therefore we only plot these families at largeQ (say, forQ 0.8) where they can have |Im(ωM )| of the order as or smaller than those for the PS 0,1 surfaces.
From the analysis of Fig. 15, several properties emerge. First, as expected, the n = 0 overtones always have the slowest decay rate of their families. Namely, in the left panel, the orange PS 0 surface is above the dark-red PS 1 surface and the green NH 0 surface is above the blue NH 1 surface. The PS 0 and PS 1 surfaces reduce to the Schwarzschild QNMs (red points) at {Q, a} = {0, 0}, whose frequencies where first computed in [40,41]. The plane with a = 0 in Fig. 15 coincides with the RN plots of Fig. 1  Naively, a "bird's-eye" view of the left panel of Fig. 15 seems to suggest that the four surfaces intersect each other with simple crossovers. For example, the orange PS 0 and the green NH 0 surfaces seem to intersect along a curveQ =Q c (â). In the RN limitâ → 0, this intersection curve gives the RN intersection point, i.e.Q c (â = 0) =Q RN This intersection curveQ =Q c (â) between the PS 0 and NH 0 surfaces indeed is well defined for 0 ≤ a/a ext 0.82 but, a fine-tuning or zoom-in analysis of the left panel proves that this is definitely no longer the case for 0.82 < a/a ext ≤ 1. Indeed, this fine-tuned analysis was already performed in Figs. 13−14: in all plots of Fig. 13 and in the top-left panel of Fig. 14 the PS 0 and NH 0 families indeed intersect with a simple crossover (but only the imaginary part of the frequency cross). However, between the top-left (for a/a ext = 0.8) and top-right (for a/a ext = 0.86) panels of Fig. 14 we concluded that the PS 0 and NH 0 families, instead of intersecting, suffer eigenvalue repulsions that effectively destroy their individual identities and leads to the formation of PS-NH families of modes. Coming back to the left panel of Fig. 15, these eigenvalue repulsions occur roughly for 0.82 < a/a ext ≤ 1 and in the charge window 0.928 Q 0.960. Again, the eigenvalue repulsions in this region are not visible in the "bird's-eye" view of the left panel of Fig. 15; we need to zoomin to make this noticeable very much like we did in Figs. 13−14. However, in Fig. 15 there is a particular point that stands-out. The is the level crossing point located at the extremal brown curve with {â ,Q } {0.360, 0.932} and Im ω = 0 and Re ω = mΩ ext H where the shows that the dark-red PS 1 and green NH 0 surfaces intersect with simple crossovers in the window 0 ≤ a/a ext 0.38 (but only the imaginary part of the frequency cross), but this is replaced by eigenvalue repulsions between the two families roughly in the window 0.38 < a/a ext ≤ 1 and 0.870 Q 0.885. Finally, other eigenvalue repulsions, e.g. between the PS 1 and NH 1 surfaces, also occur in the left panel of Fig. 15 as identified in Figs. 13−14.
The evolution and intersections of the four QNM surfaces is much simpler and much less dramatic at the level of the real part of the frequency, which is plotted in the right panel of Fig. 15. We see that the Re(ωM ) of the orange PS 0 and dark-red PS 1 families is very similar and the same happens for Re(ωM ) of the green NH 0 and blue NH 1 families. So much so that one barely distinguishes the PS 0 and PS 1 surfaces and, even less, the NH 0 and NH 1 surfaces. Moreover, nothing special happens to the real part of the frequency in the regions where the eigenvalue repulsions happen in the imaginary part of the frequency (see further discussions about this in the end of subsection 4.2).
From the analysis of both plots in Fig. 15 we see that the NH 0,1 surfaces always approach the solid brown curve with Im ω = 0 and Reω = mΩ ext H at extremality. On the other hand, the PS 0,1 curves approach this solid brown curve if and only ifâ ext (Q) >â where the point was introduced in the discussion that leads to (3.41), pinpointed in Fig. 12, and identified as the level crossing point of the system in subsection 4.1. For a ext (Q) <â which happens forQ <Q ≤ 1 this is no longer the case, in agreement with the discussions of (3.41) and of Figs. 11−12.
To explicitly demonstrate/justify why we have chosen to display many of our plots in units of r + in the plots of Figs. 13−15, in Fig. 16 we plot the n = 0 PS and NH families, so the same as in Fig. 15, but this time without including the n = 1 families and all quantities in units of M , i.e. the plot {Q/M, a/M, ωM }. In the left panel, the slope of the green NH 0 surface is now even more vertical than in Fig. 15, which indicates that the eigenvalue repulsions occur in windows of Q/M that are much narrower than in Q/r + . Furthermore, in the right panel the NH 0,1 surfaces exist in such a narrow region that they are not visible: they are too close to the extremal solid brown line with a width extremely small and invisible to the naked eye. Finally note that if we are interested on the numerical value of the frequency of the slowest decaying mode of KN, we simply need to take the mode with minimum |Im(ωM )| for a given {Q/M, a/M } in Fig. 16. For completeness, we display the result of this operation in Fig. 17, which was first presented in [3]. The Z 2 = m = 2, n = 0 KN modes with slowest decay rate always terminate at extremality along the extremal solid brown curve, with the frequencies off extremality well approximated by (3.40) as best illustrated in Figs. 9−10 and in the bottom panel of Fig. 14. The red surface family, continuously connected to the Schwarzschild mode (dark-red point [40,41]), is the PS 0 QNM family as we unambiguously identify it in the RN limit. It dominates the spectra for most of the parameter space. However, for largeQ it is instead the green surface NH 0 QNM family (as clearly identified in the RN limit) that has the lowest |Imω|. In between these orange/green regions there is a yellowish zone. This is where either simple crossovers (that trade mode dominance) or eigenvalue repulsions between the PS 0 and NH 0 modes occurs. These were analysed in the discussion of Figs. 13−14 where we also found that as we approach extremality it is appropriate to drop the PS and NH classifcation and adopt the nomenclature PS−NH families and their overtones.
In the three Figs. 15−17, at very large charge, namely forQ > 0.99 there is a gap between the last green NH line (withQ = 0.99) and the extremal solid brown curve. We have not collected data in this region because we already know (see Fig. 9 and the bottom panel of Fig. 14) that in this region so close to extremality, the analytical near-horizon MAE frequencyω MAE − as given by (3.40) − provides an excellent approximation that can be used for any physical application where such high charge values might be needed.

QNM spectra: a survey of key modes
So far we have been assuming that the least damped gravito-electromagnetic QNMs of the KN black holes are the Z 2 = m = 2 modes with n = 0. But we have not yet provided evidence that this is the case. It is certainly the case for the RN black hole subfamily (a = 0) and for the Kerr black hole (Q = 0) since several { , m} modes have already been computed in the literature for these cases. But strictly speaking the Z 2 = m = 2 does not necessarily need to be the mode with slowest decay in the whole parameter space {Q/M, a/M } of the KN black hole away from the RN and Kerr sufamilies. Therefore, in this section we do a survey of what should be the QNMs of KN that could eventually challenge the dominance of the Z 2 = m = 2 mode. Mainly these are all families with = 1, 2 and |m| ≤ that are not pure gauge modes and the Z 2 families with = m = 2, 3, 4, 5. For each mode, we only  [40,41]. In the right panel, the orange and green regions are so close to the extremal brown curve that they are not visible. display the data for the first radial overtone (n = 0) because higher overtones always have larger |Im(ωM )| that the n = 0 one.
To classify and identify more precisely the QNMs families that we will study, note that for Q, a → 0 we must recover the Schwarzschild QNMs. In this limit, it is well known that there are two families of QNMs, namely the Regge-Wheeler (aka odd or axial) modes [30] and the Zerilli (aka even or polar) modes [31,32] . These families are isospectral, i.e. they have exactly the same spectrum [40]. Ultimately, we only need to distinguish the gravitational modes of Schwarzschild (described in Table V [40] by the eigenfunction Z 1 ). In recent decades, these QNMs were computed more accurately as detailed in the review [47]. Each of these Z 2 and Z 1 modes in Schwarzschild can be found by solving a single pair of ODEs that constitute an eigenvalue problem for the angular separation constant and frequency [30][31][32]. The Schwarzschild modes are specified by the harmonic number = 0, 1, 2, 3, · · · that essentially fixes the separation constant of the problem after requiring regularity of its spherical harmonic eigenfunctions (Z 2 perturbations with = 0 and = 1 are modes that change the mass and the angular momentum of the black hole, respectively; thus we do not discuss these further). When the black hole has charge and rotation, we have to scan a two parameter space in {Q/M, a/M }. The above two families become coupled gravito-electromagnetic QNMs and the Schwarzschild eigenvalue does not appear explicitly in the KN PDEs (2.23). However, we can still count the number of nodes along the polar direction of the eigenfunctions of (2.23) and this gives . So, when Q = 0 and a = 0, we can still assign to a given mode the value of that the mode has when we trace it back continuously to the Schwarzschild limit. This is what we will do to catalogue the modes we study. In Table 1 we give the list of all modes we present. The first table is for Z 2 modes while the second is for Z 1 modes. In both tables the first column specifies { , m}, the second column gives the value of the frequency for the Schwarzschild case. It matches the frequencies first computed and listed in Table V of page 262 or in Table IV of page 202 of [40]. Finally, in the third column we identify the figure of our manuscript where the QNMs of KN with the given { , m} of the first column are displayed. In the plots of all these figures, the QNM surfaces reduce to the values of the second column of Table 1 in the Schwarzschild limit (see red points at Q = a = 0 in our figures).    [40].
In the last column of each Table we indicate the figure that extends the Schwarzschild result to thẽ Q = 0,ã = 0 case. Note that for a given , modes with |m| ≤ are degenerate in the Schwarzschild limit but this degeneracy is broken once we switch onQ andã.   [40,41]. The orange surface describes PS 0 modes while the green surface corresponds to the NH 0 modes.
We start by analysing what happens to the QNM Z 2 spectra with = m when = m progressively grows from = m = 2 (Fig. 16), to = m = 3 (Fig. 18), to = m = 4 ( Fig. 19) and, finally, to = m = 5 (Fig. 20). As for the = m = 2 case of Fig. 16, the solid brown curves at extremality are parametrized byâ = a ext = 1 −Q 2 and have Imω = 0 and Reω = mΩ ext H . We see that the main features of Figs. 18−20 for = m = 3, 4, 5 are very similar to those of Fig. 16 for the = m = 2 mode that was already analysed in much detail in Sections 4.2 and 5. 22 In particular, we identify the PS 0 and NH 0 surfaces (as unambiguously identified in the RN limit) and a zoom in of Figs. 18−20 (not shown in 22 As with the = m = 2 case of Fig. 16, note that in the right panels of Figs [40,41]. The red surface describes PS 0 modes while the green surface corresponds to the NH 0 modes. our figures) shows that these surfaces intersect with simple crossovers or with eigenvalue repulsions very much similar to those detailed in Figs. 13−14 for the = m = 2 case. Therefore, the key features of Figs. 18−20 are as discussed before. However, we highlight three features. First, the = m = 3, 4, 5 Im(ωM ) surfaces are always below and thus more damped than the = m = 2 one, and the damping increases as = m increases. Second, the NH 0 surfaces only dominate the spectra for very large Q/M and close to extremality, again very much like in the = m = 2 case. In fact, since black holes with very large charge are not expected to have any astrophysical interest, in the plots for the other modes listed in Table 1 we will no longer display the NH families (when they exist; this is certainly the case for the Z 1 = m = 1, 2 modes). Finally, note that for = m ≥ 3 it is still true that the PS 0 frequencies are well approximated by (3.16) (in fact the approximation gets better as m increases and we approach the eikonal limit; see also the = m = 6 case in Fig. 3) and the NH 0 frequencies are in excellent agreement with (3.40) near extremality. Complementing the analysis reported here, note that in [1] we have reported our findings for Z 2 = 3 with m = −3, −2, −1, 0, 1, 2, 3 and we have concluded that their |Imω| is aways higher than the Z 2 = m = 2 modes.
Next, we consider the several cases of Z 2 modes with |m| ≤ = 2 in Fig In the Schwarzschild limit all these = 2 modes are degenerate withω 0.45759551 − 0.09500443 i, but this degeneracy is broken once we switch onQ andã. The figures speak for themselves and we refrain from describing them further. We note simply that the surfaces for positive m have a qualitative shape that is significantly distinct from the ones for negative m (notably, m ≥ 0 cases have a monotonic behaviour that is not observed in the m < 0 cases) and, as expected, further note that for m = the PS 0 surfaces no longer approach Imω = 0 and Reω = mΩ ext H at extremality (hence we do not display these solid brown curves in the associated plots). This sequence of figures demonstrates, as previously claimed, that Z 2 modes with = m = 2 are the dominant ones among the |m| ≤ = 2

families (and all others).
We can now consider the Z 1 modes which are purely electromagnetic modes in the Schwarzschild (and Kerr) limit. 23 In Figs. 24, 25, 26, 27 and 28, we display the = 2 PS surfaces of this family for m = 2, 1, 0, −1, −2, respectively. Moreover, in Figs. 29, 30 and 27, we display the Z 1 = 1 PS surfaces for m = 1, 0, −1, respectively. Comparing Z 1 modes with the same { , m} as Z 2 modes, we see that the qualitative shape of the surfaces is similar but Z 1 modes are typically more damped than the Z 2 modes. Moreover, Z 1 modes with = m also approach Imω = 0 and Reω = mΩ ext H at extremality if and only ifâ ext (Q) >â (see Fig. 24 for = m = 2 and Fig. 29 for = m = 1) where the point was defined in the discussion leading up to (3.41). Forâ ext (Q) <â which occurs forQ <Q ≤ 1 this is no longer the case, very much like in the Z 2 = m discussions of (3.41) and of Figs. 11−12. For a given = m, the value ofâ for Z 1 modes tends to be higher than the one for Z 2 modes.