Holographic thermalization, quasinormal modes and superradiance in Kerr-AdS

Black holes in anti-de Sitter (AdS) backgrounds play a pivotal role in the gauge/gravity duality where they determine, among other things, the approach to equilibrium of the dual field theory. We undertake a detailed analysis of perturbed Kerr-AdS black holes in four- and five-dimensional spacetimes, including the computation of its quasinormal modes, hydrodynamic modes and superradiantly unstable modes. Our results shed light on the possibility of new black hole phases with a single Killing field, possible new holographic phenomena and phases in the presence of a rotating chemical potential, and close a crucial gap in our understanding of linearized perturbations of black holes in anti-de Sitter scenarios.

horizon and can trigger superradiant instabilities at the linear level, and even induce other nonlinear phenomena. Additionally, BHs in AAdS play a central role in the formulation and applications of the AdS/CFT correspondence. This correspondence [1,2] provides a remarkable framework for studying certain strongly coupled gauge theories in d dimensions by mapping them to weakly coupled quantum gravitational systems in d + 1 dimensions. In a certain limit (namely in the large 't Hooft coupling and planar limit), quantum gravity in the bulk reduces to classical general relativity. Within this holographic framework, a black hole is dual to a thermal state and the question of thermalization in the boundary gauge theory translates, in the gravitational bulk, to understanding the "return to equilibrium" behavior of perturbed black holes [3,4,5,6,7,8,9,10].
Here, we will be interested in the original gauge/gravity duality, namely the AdS d+1 /CFT d correspondence (for the d = 3, 4 cases for which Super-Yang-Mills theory is dual to string theory on AdS 4 × S 7 and AdS 5 × S 5 , respectively). Moreover, we are particularly interested in systems with a rotating chemical potential. This requires looking to CFTs formulated on a sphere (since a rotational shift is a pure gauge transformation on a plane), i.e. for bulk solutions that asymptote to global AdS (which is conformal to the static Einstein Universe R t × S d−1 ). Henceforward, when we refer to AAdS spacetimes it is implicitly assumed that we mean global AdS (although some of the discussions are also valid for planar AdS i.e. the Poincaré patch of AdS that asymptotes to R t × R d−1 ). We will often use the notation D = d + 1 for the bulk spacetime dimension; Greek indices denote the bulk dimensions while Latin indices describe boundary coordinates.
Certainly, important headways on thermalization can be made by studying the behavior of perturbed black hole spacetimes at a linearized level. Naturally, the applicability of such analysis depends on the strength of the perturbation off the stationary black hole and the behavior obtained can hint of possible instabilities [11,12,13,14].
The analog problem in asymptotically flat spacetimes is, to a large degree, understood. The approach to equilibrium depends sensitively on the character of the perturbation: massless fields (scalar, vector or tensor type) die off through their quasinormal modes (QNMs), with a time dependence of the form e −iωt with Im(ω) < 0 (for a review see [15]) 1 . Massive fields on the other hand, have a much richer phenomenology, tied to the fact that they can be trapped inside a cavity with size of the order of the Compton wavelength. This trapping causes the field to decay much slower, or can even become unstable for large black hole rotation (see [16] and references therein). The linear behavior of massive fields around rotating black holes is not fully understood yet (and certainly not the nonlinear regime), but it is akin to that of massless fields in AAdS backgrounds in that both can develop trapping potentials. However, an important difference is that the height of the potential barrier is unbounded in the AdS case while it is finite for a massive scalar in a flat background.
Accurate expressions for the QNMs for generic black holes in the asymptotically flat case are known for both static and stationary black holes (see the review [15]); remarkably this is not the case in the AdS background as they are not known for the Kerr-AdS BH. This status of affairs is, at first sight, surprising given the central role they play in holographic dualities as well as in studies of AAdS black hole stability.
It is thus worth discussing in detail the reason for this gap in our knowledge. Since an AAdS spacetime is a non-globally hyperbolic spacetime (i.e., spatial infinity is a timelike boundary in the Carter-Penrose diagram and thus null rays can reach it in finite time), in order to predict the future time evolution of the system we need to give not only initial data but also to specify the choice of boundary conditions (BCs). At the inner boundary (origin or horizon) regularity fixes the BC. However, at the asymptotic boundary this choice isà priori arbitrary, being fixed by a physical motivation. From a pure gravitational perspective it is often stated in the literature that one is interested in "reflecting BCs" which suggests the idea that we want vanishing flux of energy and angular momentum across the asymptotic boundary. On the other hand, in the context of the AdS/CFT duality we typically want to choose BCs that preserve the asymptotic boundary (conformal) metric. Next, and in Appendix A, we emphasize that these two perspectives require exactly the very same BCs. Formally, the discussion of the asymptotic BCs is more clear if we write the total metric (background plus perturbations, if present) in Fefferman-Graham coordinates (this frame is defined such that g zz = L 2 /z 2 and g zb = 0, where z is the radial distance with boundary at z = 0, and x b are the coordinates on the boundary), and looking at its boundary expansion (see [17] and references therein). For odd boundary dimension d this reads 2 ds 2 = L 2 z 2 dz 2 + g ab (z, x)dx a dx b , ab (x) + · · · + z d g ab (x), (1.1) where g (0) (x) and g (d) (x) are the two integration "constants" of the expansion; the first dots include only even powers of z (smaller than d) and depend only on g (0) (thus being the same for any solution that asymptotes to global AdS d ) while the second dots depend on the two independent terms g (0) , g (d) (we will fix Newton's constant as G N ≡ 1). Within the AdS/CFT duality we are (typically) interested in Dirichlet BCs that do not deform the conformal metric g (0) . Indeed, this defines the gravitational background where the CFT is formulated and we want to keep it fixed; in our case this is the static Einstein Universe. Stated in other words, we allow perturbations in the bulk that only deform the expectation value of holographic stress tensor T ab (x) (that specifies and describes the boundary CFT) 3 but that preserve the asymptotic structure of the original background that we perturb. 4 As discussed in Appendix A these BCs do not allow asymptotic dissipation of energy or angular momentum. In other words, everything that hits the asymptotic boundary is reflected back to the bulk core allowing for a non-trivial interplay between the asymptotic and inner (e.g. horizon) boundaries. We have now growing evidence that these conditions favour the development of instabilities. For instance, it has been shown recently that even arbitrarily small perturbations can trigger black hole formation in global AdS [18,19], indicating that global AdS is nonlinearly unstable to a weak perturbative turbulent mechanism (note however the existence of "islands" of stability [20,21,22]). Additionally, it has recently been shown that turbulent behavior 5 (akin to the one displayed by hydrodynamics) arises in long-wavelength perturbations of 3+1 Kerr-AdS [23,24]. The BCs we need to impose to study AAdS perturbations of global AdS BHs are therefore well known. Yet, we still need to understand why the study of QNMs and superradiant instabilities of global AdS BHs is not a closed chapter. For that, we need to look to the perturbation 2 For even d, the asymptotic expansion (1.1) contains also a logarithmic term z d log z 2g(d) ab and the holographic stress tensor has an extra contribution proportional to the conformal anomalies of the boundary CFT [17]. These details are not essential for the present discussion. 3 Note that g (d) is an integration "constant" but not a free function; it is fixed solving the Einstein equation in the bulk subject to regular BCs at the horizon or radial origin. 4 Other BCs that might be called asymptotically globally AdS (and that promote the boundary graviton to a dynamical field) were proposed in [96]. However, they turn out to lead to ghosts (modes with negative kinetic energy) and thus make the energy unbounded below [97]. 5 As well, in planar AdS backgrounds, turbulent behavior of gravity has also been uncovered for (sufficiently) long-wavelength perturbations of black holes in [52,23].
equations. Studying linearized gravitational perturbations requires solving the linearized Einstein equations for the metric perturbation. For generic perturbations this is a coupled nonlinear system of PDEs. Solving this PDE system directly with the above BCs is a hard problem, even numerically. In certain special cases, however, drastic simplifications occur. Fortunately, and quite remarkably, in four dimensions it has been shown that if we use certain gauge invariant scalar variables we can reduce the problem of looking for the most generic perturbations of the above AAdS BHs to solving a single PDE. Moreover, using the harmonic decomposition of the system, the later reduces to solving two ODEs. This remarkable reformulation of the linearised perturbation problem is known as the Regge-Wheeler−Zerilli or Kodama-Ishisbashi formalism for perturbations of Schwarzchild BHs [25,26,27], and Teukolsky formalism for perturbations of Kerr BHs [28,29]. We ask the reader to see the companion paper [30] for a detailed discussion of these two formalisms and for the map relating them when the background rotation vanishes. Once the solution for the gauge invariant scalars is known a simple differential map generates the corresponding metric perturbation tensors (in a particular gauge).
At this point, to find QNMs or instability timescales of AAdS BHs we just need to take the above BCs, discussed for the metric perturbations, and translate them to get the corresponding BCs that need to be imposed on the gauge invariant scalars. On general grounds we should expect the Dirichlet BCs on the metric to translate to Robin BCs (which relate the field with its derivative) on the gauge invariant scalars. In the static background case, this dictionary was found by [9,10]. There are two families of perturbations: scalar (also called even or polar) and vector (odd or axial) sectors. The associated QNMs of the global AdS Schwarzchild BH were then computed [9,10,30]: vector QNMs agree with those first computed in [31,32,33] (the scalar modes of [31,32,33] do not preserve the asymptotic AdS structure). In the stationary case, the BC map was constructed only recently in the companion paper [30]. With it at hand, we can finally compute the gravitational QNM spectrum and superradiant instability growth rates of the Kerr-AdS BH. This is one of main aims of the work here reported. (Previous work on gravitational QNMs [34] and superradiant instability of Kerr-AdS [35] imposed BCs that do not keep the boundary metric fixed). While many of the methods presented here are readily applicable to arbitrary dimensions we concentrate in dimensions d = 3 and d = 4 because of their interest for the AdS d+1 /CFT d holographic dualities.
The interest on the superradiant instability is not restricted to its growth rate. Indeed, the onset curve of this instability (where the imaginary part of the frequency vanishes) is an exact zero mode that is invariant under the horizon-generating Killing field of Kerr-AdS. Therefore we will argue that, in a phase diagram of stationary solutions, this onset curve signals a bifurcation curve to a new family of BHs that have a single Killing vector field (KVF), i.e. that are periodic but not time dependent neither axisymmetric. A far reaching consequence of this statement is that Kerr-AdS BHs are not the only stationary BHs of Einstein-AdS gravity. These BHs can exist because they evade a main assumption of the rigidity theorems [36,37,38]. We will give the explicit perturbative construction of the leading order thermodynamics and properties of these BHs. These ideas were first proposed in [39] and further developed in [40,19]. Now that we have the precise onset curve of superradiance, we have the opportunity to expand their discussion.
Another aim of the present work is to confirm that long wavelength gravitational QNM frequencies agree with the hydrodynamic relaxation timescales that we obtain when we consider the near-equilibrium and long wavelength effective description of the CF T d . This will provide the first explicit confirmation that the match between the QNM spectrum and the CFT thermalization timescales also holds in the presence of a rotating chemical potential. Incidentally, it provides the first non-trivial confirmation that the Robin boundary conditions for the Teukolsky gauge-invariant variable derived in the companion paper [59] are indeed the ones that we must impose if we want the perturbations to preserve the conformal metric.
This work is divided as follows. Section 2 reviews relevant properties of D = 4 Kerr-AdS spacetime and the equations of motion and the BCs [30] governing the behavior of perturbations at the linear level. Section 3 describes the numerical methods employed to solve them. One of these numerical approaches is novel and we expect it to be of interest for other applications. Section 4 presents our results for the full spectrum of gravitational QNMs and superradiant instability timescales of the Kerr-AdS BH. In Section 5 we construct and discuss the novel single Killing vector field BHs that merge with the Kerr-AdS BH at the onset of superradiance. In section 6 we use the fluid/gravity duality to confirm the match between the gravitational long-wavelength QNM spectrum and the CF T 3 hydrodynamic modes even in the presence of a rotating chemical potential. Section 7 repeats the previous section computations and discussions but this time for the D = 5 rotating system that is of interest for the AdS 5 /CFT 4 duality. It will also contribute to identify universal properties of systems with a rotating chemical potential. We work in a particularly clean environment where we study perturbations around the equal angular momentum Myers-Perry BH. Indeed, this background has enhanced symmetry − it only depends non-trivially on the radial direction − and its perturbations have an exact harmonic decomposition. The present study fills important gaps in our knowledge but confirms and opens some interesting questions. In Section 8 we discuss these open questions in what can be viewed as a roadmap in the subject from our viewpoint.

Gravitational perturbations & boundary conditions of Kerr-AdS black hole
In this section we review the basic properties of Kerr-AdS black holes and their gravitational perturbations.

Kerr-AdS black hole
The Kerr-AdS geometry was originally discovered by Carter in the Boyer-Lindquist coordinate system {T, r, θ, φ} [41]. For our purposes, it is convenient here, to follow Chambers and Moss [42] and introduce the new time and polar coordinates {t, χ}, which are related to the Boyer-Lindquist coordinates {T, θ} by where a is the rotation parameter of the solution and Ξ is to be defined in (2.3). In this coordinate system the Kerr-AdS black hole line element reads [42] 3) The Chambers-Moss (CM) coordinate system {t, r, χ, φ} has the appealing property that the line element treats the radial r and polar χ coordinates on an almost equal footing. This property extends to the radial and angular equations describing gravitational perturbations in the Kerr-AdS background. In this frame, the horizon angular velocity and temperature are given by Ω H = a r 2 + + a 2 , The Kerr-AdS black hole obeys R µν = −3L −2 g µν , and asymptotically approaches global AdS space with radius of curvature L. This asymptotic structure is not manifest in (2.2), one of the reasons being that the coordinate frame {t, r, χ, φ} rotates at infinity with angular velocity Ω ∞ = −a/(L 2 Ξ). However, if we introduce the coordinate change we find that as r → ∞ (i.e. R → ∞), the Kerr-AdS geometry (2.2) approaches which we recognize as the line element of global AdS. In other words, the conformal boundary of the bulk spacetime is the static Einstein universe R t × S 2 : This is the boundary metric where the CFT lives in the context of the AdS 4 /CFT 3 correspondence.
The ADM mass and angular momentum of the black hole are related to the mass M and rotation a parameters through M ADM = M/Ξ 2 and J ADM = M a/Ξ 2 , respectively [43,44]. The horizon angular velocity and temperature relevant for the thermodynamic analysis are the ones measured with respect to the non-rotating frame at infinity [43,44] and are given in terms of (2.4) by T h = Ξ T H and Ω h = Ξ Ω H + a L 2 . The event horizon is located at r = r + (the largest real root of ∆ r ), and it is a Killing horizon generated by the Killing vector K = ∂ T + Ω h ∂ Φ . We discuss our results in terms of the horizon radius and rotation parameter, which uniquely determine a given Kerr-AdS black hole. The mass parameter is given in terms of these by M = r 2 + + a 2 r 2 + + L 2 / 2L 2 r + . All regular black hole solutions must obey T H ≥ 0 and a/L < 1. This translates into the following conditions for r + /L and a/L: The first inequality is saturated for a degenerate extremal regular horizon. On the left panel Fig. 1, we show the allowable domain for a/L and r + /L. Further properties of the Kerr-AdS spacetime are discussed in Appendix A of [45]. We will find it useful to parametrize the black hole in variables that are naturally related to the onset of superradiance, and that are gauge invariant. Here we choose the pair (R + , Ω h ), with R + given by: Extremality is attained at (2.9) Note that R + is just the square root of the area of the spatial section of the event horizon, divided by 4π, often denominated areal horizon radius. The allowed values of R + /L and Ω h L are depicted on the right panel of Fig. 1.

Gravitational master equation and global AdS boundary conditions
In the Newman-Penrose−Teukolsky formalism, all the information about (non-trivial) gravitational perturbations with spin s = −2 is encoded in the single variable δΨ 4 which describes the perturbation of the Weyl scalar Ψ 4 = C abcd n amb n cmd . The equation of motion for this perturbation δΨ 4 is described by the s = −2 Teukolsky master equation [28,29]. Introducing the separation ansatz the spin s = −2 Teukolsky master equation separates into angular and radial equations [42,30]: where we have defined ω m (χ) are the spin-weighted s = −2 AdS-spheroidal harmonics. The positive integer specifies the number of zeros of the eigenfunction along the polar direction which are given by −max{|m|, |s|} (so the smallest is = |s| = 2). The associated eigenvalues λ are functions ofω, , m which can be computed numerically. Regularity imposes the constraints that − ≤ m ≤ must be an integer and ≥ |s|. This equation has been studied in [30] in the limit where the rotation vanishes.
If we solve (simultaneously) the angular and radial equations, which are coupled through the two eigenvaluesω and λ, we get information about the most general perturbation of the Kerr-AdS black hole. In particular, the Teukolsky equation and its solution for the spin s = +2 perturbations, described by the variable δΨ 0 = (r − iχ) −2 e −iωt e imφ R (2) ω m (r) S (2) ω m (χ), follow trivially from the spin s = −2 solution. Namely, R ω m (−χ). The later statement implies that the separation constants are such that λ ω m ≡ λ. The only exceptions to the above are the trivial perturbations, the " = 0" and " = 1" modes, which shift, respectively, the mass and angular momentum of the solution along the original Kerr-AdS family, and to which the Teukolsky formalism is blind [46,27,47,30].
Quasinormal modes and unstable modes of the Kerr-AdS black hole are solutions of (2.11)-(2.12) obeying physically relevant boundary conditions (BCs) [30]. At the horizon, the BCs must be such that only ingoing modes are allowed. A Frobenius analysis at this boundary gives two independent solutions, (2.14) 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 BC, we introduce the ingoing Eddington-Finkelstein coordinates {v, r, χ, φ}, which are appropriate to extend the solution through the horizon. These are defined via The BC is determined by the requirement that the metric perturbation is regular in these ingoing Eddington-Finkelstein coordinates, where the metric tensor is constructed applying a differential operator to δΨ 4 (this is known as the Hertz map; see the companion paper [30]). It follows that the metric perturbation is regular at the horizon if and only if R(r)| H behaves where R IEF (r)| H is a smooth function 6 . Therefore, the appropriate BC at the horizon demands we set A out = 0 in (2.14): is what we might call the superradiant factor. Less formally, but perhaps more intuitively, wheñ ω is real and non-zero we can understand this horizon BC by noting that the wave solution e −iωt (r − r + ) −i = e −i(ωt+ ln(r−r + )) is the one that describes ingoing modes at the horizon 6 This analysis misses the special case in which 2iω −mΩ H 4πT H is a positive integer. For this special value, our boundary conditions still allow for outgoing modes at the horizon. However, by inspecting our numerical data we can a posteriori test if this condition is satisfied. In all our simulations, this never seems to be the case. since r must decrease when t grows to keep the phase constant (classically, we cannot have outgoing modes leaving the horizon).
Consider now the asymptotic boundary. A Frobenius analysis of the radial Teukolsky equation(2.12) at infinity yields the two independent asymptotic decays for arbitrary amplitudes B ± . We want the perturbations to preserve the asymptotic global AdS structure of the background Kerr-AdS black hole, i.e. we want the deformation to preserve the asymptotic line element (2.6). In the companion paper [30] we found that this requirement yields the following Robin BC, B with two possible solutions for β, that we call β s and β v , where we have introduced Perturbations obeying the BCs (2.19)-(2.20) preserve the asymptotically global AdS behavior of the background. These are also natural BCs in the context of the AdS/CFT correspondence: they allow a non-zero expectation value for the CFT stress-energy tensor while keeping fixed the boundary metric. The BC (2.19), (2.20) generates what we might call the "rotating sector of scalar modes", in the sense that when the rotation vanishes, these perturbations reduce continuously to the Kodama-Ishibashi scalar modes [30]. 7 By a similar reasoning the BC (2.19), (2.21) selects the "rotating vector modes" of the spectrum. Having this in mind we will often use the nomenclature "scalar/vector" modes when discussing our results [30].
As discussed previously, the Chambers-Moss coordinate system {t, r, χ, φ} rotates at infinity. However, the coordinate transformation (2.5) introduces the coordinate frame {T, R, Θ, Φ} appropriate to discuss the asymptotic global AdS 4 structure of the geometry and the boundary metric where the dual CFT 3 and its hydrodynamic limit are formulated. Consider a generic linear perturbation in Kerr-AdS written in the Chambers-Moss frame {t, r, χ, φ}. Since ∂ t and ∂ φ are isometries of the background geometry we can Fourier decompose the perturbation in these directions as e −iωt e imφ as we did in (2.10). The frequencyω measured in the frame {t, r, χ, φ} differs from the frequency measured in the frame {T, R, Θ, Φ}. It follows from the coordinate transformation (2.5) that they are related by with ω ≡ω Ξ + m a L 2 .
(2.23) 7 The Kodama-Ishibashi vector master equation is the Regge-Wheeler master equation for odd (also called axial) perturbations [25], and the Kodama-Ishibashi scalar master equation is the Zerilli master equation for even (also called polar) perturbations [26].
The quantity ω can be viewed as the natural or fundamental frequency since it measures the frequency with respect to a frame that does not rotate at infinity. This is also the natural frequency measured by a CFT 3 and associated fluid rest frame observer. Therefore, although we will use the frame {t, r, χ, φ} andω to discuss many of our results, we choose to plot our results in terms of ω. Note that the superradiant factor defined in (2.16) can equally be written as = ω−mΩ h 4πT h where the angular velocity Ω h and temperature T h as measured in the {T, R, Θ, Φ} frame are given below (2.6).

Numerical procedures
In this section we discuss the numerical procedures applied to solve for the characteristic frequency ω and separation constant λ. We present three such methods based on: shooting, series expansion, and Newton-Raphson. The first two methods are typically used in studies of QNMs and the latter we introduce here and have found it to be the most robust when exploring limiting cases. As a powerful check, we find excellent agreement between different methods when more than one is applicable.
Shooting. The first method "shoots" for the correct answer in both the angular and radial component. Regularity of the angular eigenfunctions require that they admit the following expansion at the left-and right-boundaries respectively. The coefficients B L n , B R n can be extracted from the angular equation and are functions of the frequency and the separation constant. We typically keep the first six terms in the expansion, numerically integrate the solutions towards each other where we match the logarithmic derivative at an intermediate point. We proceed identically with the radial equation, by imposing conditions (2.16) and (2.19) at the boundaries. Due to well-known divergences of QNMs at the horizon (stable modes diverge exponentially), we use an analytical, series expansion close to the horizon and a similar expansion close to spatial infinity. An example notebook of how the radial equation is dealt with can be found online [15]. The method gives stable, convergent results for small black holes, but becomes less accurate for large black holes.

Series expansion.
A powerful alternative is based on a series solution of the radial equation which avoids the divergent nature of QNMs at the horizon altogether by factoring the relevant terms [3,15]. For simplicity let us focus on non-rotating BHs in this brief description, the extension to rotating BHs is straightforward. Let us start by re-expressing the boundary condition (2.19) as −r(r/LR ω m , where primes denote derivative with respect to r and all quantities are evaluated at spatial infinity. Redefine the wavefunction to R (−2) ω m = ∆r r 5 e −iωr * Z(r), with dr/dr * = ∆ r /r 2 . Then, make the variable change z = 1/r and re-write the radial equation as and the boundary conditions as where primes now denote derivative with respect to z and z = 1/r + . The idea is now to look for a series solution, Z = a n (z − z ) n , where the coefficients a n are found through the recurrence relation where P n = n(n − 1)s 0 + n t 0 and where s, t, u have been expanded in Taylor series around the horizon. The boundary condition then translates into a n (−z ) where β is given by either Eq. (2.20) or Eq. (2.21). Extension to rotating geometries is obtained simply by replacing ω with the corresponding superradiant factor.
Newton-Raphson. We have also developed a novel numerical procedure based on the Newton-Raphson root-finding algorithm that searches for specific quasinormal modes, once a seed solution is given. In order to proceed we first need to recast Eq. (2.11) and Eq. (2.12) in a different form. Let us introduce the following auxiliary functions: where we have implicitly introduced two new compact coordinates y = 1−r + /r andx = 1+χ/a, which map the problem to the unit square: (x, y) ∈ (0, 1) × (0, 1). The boundary conditions on the q I simply arise from regularity, and translate into four Robin boundary conditions at each integration boundary, i.e. We are now ready to introduce the new numerical procedure that determines {q 1 , q 2 , ω, λ}. For the sake of presentation we will only discuss below the case in which we have a single differential equation to solve. The extension to a coupled system like the one above is straightforward.
Consider the following "nonlinear Stürm-Liouville" problem in {f,λ}: where H(λ) is nonlinear function inλ, and a linear differential operator in f and both {a 0 , b 0 } are constants. In many circumstances H takes the following simple form: where each of the H i is a second order differential operator independent ofλ. The former differential equation is often called a quadratic eigenvalue problem, so long as the constants {a 0 , b 0 } admit a similar expansion. The method we describe below allows for any dependence inλ.
We discretize our Eq. (3.9) by introducing a spatial grid {y i }, with N + 1 grid points. Because we are solving for manifestly analytic functions q I , we can readily use a pseudospectral collocation discretization scheme. We choose the Gauss-Chebyshev-Lobatto grid as our collocation points. The nonlinear Stürm-Liouville problem (3.9) reduces to a nonlinear eigenvalue problem of the form: where D i,j is a Chebyshev differentiating matrix and H i,j is the discretization of the operator H. We now introduce a normalization for the eigenvector {f i }, using an auxiliary constant vector {v i }, such that v i f i = 1. In all cases, we choose {v i } to have only one nonzero component, which without loss of generality we choose to be the horizon and the south pole located at y = 0 and x = 0, respectively. The procedure is now clear: we promoteλ to be a parameter to be determined via the Newton-Raphson method. Recall that we have to solve whereĤ i,j is obtained from H i,j by removing its first and last lines, and substitute them by the last two conditions in Eq. (3.10). The Newton-Raphson method states that the correction to our initial guess for ({f (0) i },λ (0) ) can be determined by inverting the following linear system of equations: We then iterate this procedure until |δf j | and |δλ| are below some tolerance, which in this manuscript we take to be 10 −30 . All computations using this method were performed with octuple precision, which is particularly relevant for small black holes. Our results have been benchmarked using previous results in the literature, specifically for scalar field perturbations of Kerr-AdS BHs [48,35,49]. In particular, we recover to all significant digits the numerical results reported in Ref. [49]. Furthermore, we recover all known results from gravitational perturbations of Schwarzschild-AdS BHs with the same boundary conditions [31,50,10,30].
Finally, we note that an important symmetry of the relevant perturbation equations and boundary conditions for QNMs is that if (ω, λ) is a solution for a given m then (−ω * , λ * ) is a solution for −m. As such, we will only discuss positive real part modes, with the understanding that they come in complex conjugate pairs.

QNMs and superradiance in Kerr-AdS: results
In this section we present the numerical results obtained, make contact with some analytical results, and discuss implications with the phenomena of superradiance.

Comparison between analytical and numerical results
The angular (2.11) and radial (2.12) equations constitute a system of ordinary differential equations coupled through the frequencyω and angular λ eigenvalues that cannot be solved analytically when M, a = 0. For this reason, we solve these equations using the numerical methods outlined in Section 3. There is however a regime where we can use a matched asymptotic expansion procedure to get an approximate analytical solution for the QNM and superradiant instability frequency spectra. This perturbative analytical computation provides useful physical insights about the system and is valuable to check our numerical results. We leave the details of this analytical construction to Appendix C and present here only the main outcome of the computation and its comparison with the numerical results.
As justified in Appendix C, the perturbative analytical results are valid in the regime of parameters where i.e., for Kerr-AdS black holes with small horizon radius in AdS radius units and even smaller rotation parameter, and for perturbations whose wavelength is much bigger than the black hole lengthscales. . This is for = 2 modes with no radial overtone. These are an example of the QNM spectrum in the regime a/L < r + /L 1 where the analytical matching analysis is valid and its approximated results can be used to both test our numerical code (valid in any regime), and estimate more precisely the regime of validity of the analytical approximation. The red dots are the exact results from our numerical code. The green curve is the numerical solution of the matching transcendental equation (4.2), while the dashed black curve is the approximated analytical solution (4.3) or (4.4) of (4.2). In both figures there is a critical rotation where Im(ωL) = 0 and Re(ω) − mΩ H 0 to within 0.01%. For lower rotations the QNMs are damped and with Re(ω) − mΩ H > 0, while for higher rotations we have unstable superradiant modes with Re(ω) − mΩ H < 0.
In Appendix C we find that the matched asymptotic expansion analysis indicates that the frequency spectrum is quantized by the condition (C. 19), for a generic mode with quantum numbers and m. This frequency quantization condition simplifies considerably when we choose a particular harmonic . For instance, for the lowest harmonic = 2, the condition (C.19) reads where the superradiant factor is defined in (2.17), and ε j = 1 describes scalar modes with the BC (C.14) while ε j = −1 represents vector modes with the BC (C.15). We can find the frequency that solves this transcendental equation numerically using a standard rootfinder routine (for instance Mathematica's built-in FindRoot routine). Alternatively we can also provide an approximate analytic solution, still in the limit of a/L r + /L 1, assuming that the frequency has a double expansion in the rotation and in the horizon radius, ω(a, r + )L = n 2) Vector modes: Im(ωL) where γ 0.577216 is the Euler-Mascheroni constant. For both scalar and vector modes the imaginary part of the frequency starts negative for a = 0, consistent with the fact that QNMs of Schwarzschild-AdS are always damped. However, as a/L increases, Im(ωL) increases. A good check of our analytical matching analysis is that we find that at the critical rotation where the crossover occurs, i.e. Im(ωL) = 0, one has Re(ω) − mΩ H 0 to within 0.01%. For smaller rotations one has Re(ω) − mΩ H > 0 and for higher rotations one has Re(ω) − mΩ H < 0 and Im(ωL) > 0. Therefore, the instability which is triggered at large rotation rates has a superradiant origin since the superradiant factor becomes negative, < 0 precisely when the QNMs go from damped to unstable. These analytical matching results provide also a good testbed check to our numerics. Indeed we find that our analytical and numerical results have a very good agreement in the regime of validity of the matching analysis. This is demonstrated in Fig. 2 where we plot our numerical and analytical results for the fundamental = 2 scalar and vector modes. As a rough reference we can take this to be r + /L < 5 × 10 −3 and a/L < 10 −4 . (A similar analysis that lead to the results (4.2)-(4.4) can be repeated for any other harmonic starting from (C.19)).

Properties of superradiant unstable modes and QNMs
We are now ready to present the properties of the superradiant unstable modes and QNMs for generic solutions in the parameter space. We use the numerical methods described in Section 3 to find the solution of the coupled ODE angular (2.11) and radial (2.12) equations that describe the most general linear perturbation of a Kerr-AdS BH. We first present the gravitational scalar perturbations that obey the BCs (2.20), and then the gravitational vector perturbations that obey the BCs (2.21).
Consider a Kerr-AdS BH parametrized by particular values of the gauge invariant parameters {R + /L, Ω h L} described in the end of Section 2.1. A generic perturbation can have a frequency with negative, positive, or vanishing imaginary part. Quasinormal modes are damped, Im(ω) < 0, whereas unstable modes grow exponentially in time, Im(ω) > 0. Thus, a particularly important set of modes, if present, are the marginal modes that define the stability boundary in a phase diagram. The marginal mode (or onset mode) curve is defined to be the locus of points in the parameter space (R + /L, Ω h L) for which a mode with Im(ω) = 0 exists. There will be a marginal mode curve for each distinct pair of wave numbers { , m} resulting in an instability. To understand the nature of this instability it is useful to look into another useful characterization of linear perturbations. It comes from considering the difference between the real part of the frequency and mΩ h , which determines the sign of the the energy and angular momentum fluxes the perturbation carries through the future horizon; see Appendix A 8 . Modes with Re(ω) > mΩ h carry positive flux through the horizon, whereas modes with Re(ω) < mΩ h carry negative flux across the horizon, and are called superradiant. Vanishing flux at the horizon requires Re(ω) = mΩ h . We find that Re(ω) = mΩ h whenever Im(ω) = 0 and that Re(ω) < mΩ h when Im(ω) < 0. Therefore, unstable modes in Kerr-AdS are always associated to the superradiant instability.
As important illustrative examples, in the left panel of Fig. 3 we identify the superradiant onset curves (OC) for = m scalar modes (with vanishing radial overtone) in the phase diagram of Kerr-AdS BHs. The axes are given by the gauge invariant horizon radius R + /L and the horizon angular velocity Ω h L (for the frame that does not rotate at infinity), as previously introduced in Fig. 1. Regular Kerr-AdS BHs exist in the blue shaded area, starting at Ω h L = 0 and all the way up towards the black curve where extremality is attained. We identify the OC for the scalar modes with = m = 2, 3, 4, 5. BHs that are above a particular = m OC are superradiantly unstable to modes with those particular values of = m, while BHs below a particular OC are stable to the associated modes. For completeness, in the right panel of Fig.  3 we plot the angular eigenvalue λ along the superradiant OC. Since Im(ω) = 0 along this OC, it follows from the mathematical structure of the coupled equations that we must also have Im(λ) = 0.
The OCs have some properties that merit a detailed discussion. First, in both plots of Fig. 3 the large black points on the left at R + /L = 0 are computed analytically and serve as additional checks for the numerical code. They describe the scalar normal mode frequencies and the associated angular eigenvalues of global AdS given by [19,21], where p = 0, 1, 2, · · · is the radial overtone (number of radial nodes). In more detail, to get the black points in the left panel of Fig. 3 we use the superradiant onset condition to find Ω h R + =0 = ω AdS s /m and we set p = 0, = m, i.e.
Note that given a { , m} pair there is an OC for each radial overtone p, but p > 0 curves always lie above the p = 0 curve, and therefore p = 0 modes are the first to go unstable as the rotation is increased. For this reason only the p = 0 curves are plotted. The OCs always have Ω h L > 1, monotonically approaching Ω h L → 1 (from above) asymptotically as R + /L → ∞, where all the scalar superradiant OCs pile up. This means that only Ω h L > 1 BHs can be unstable to superradiance, a property that was first proven in [51].
Finally, note that for small BHs (say with R + /L 0.45) as = m increases the corresponding superradiant OC lowers. This means that, e.g. we can have small BHs (those in the triangle-like region between the = m = 2 and = m = 3 curves) that are stable to = m = 2 modes but unstable to all other = m ≥ 3 modes, or e.g. BHs that are stable to = m = 2 and = m = 3 but always unstable to all other = m ≥ 4 modes. However, as the areal radius grows we find that the OCs start crossing each other. For example, the = m = 2 curve crosses the = m = 3 curve at R + /L ∼ 0.45 and for higher radius it crosses the = m = 4 and then the = m = 5 curve. So, e.g. at R + /L = 1 the = m = 2 OC is below the three OCs = m = 3, 4, 5. This means that at this radius we can have Kerr-BHs that are unstable to = m = 2 modes but not to = m = 3, 4, 5 modes. At first sight, this is of course exciting as it seems to indicate that there is a region of parameter space where Kerr-BHs are unstable to = m = 2 modes but stable to any other superradiant modes, with obvious consequences for the endpoint of the superradiant instability. However, this is not the case. Indeed, first notice that as = m → ∞ the corresponding OC still starts precisely at the point defined by (4.6). Thus, as = m grows large, its threshold modes are described by an OC that progressively approaches the line Ω h L = 1, becoming horizontal in the limit = m → +∞. Therefore as the BH rotation is increased, the first modes that become superradiantly unstable are the m → ∞ modes. The conclusion that m → +∞ modes are the "first" to become unstable was first presented in the equal angular momenta Myers-Perry BHs in [39]. Furthermore, as we shall discuss later, all vector modes will be superradiantly unstable.
As stated previously, in the left panel of Fig. 3, BHs that are above a particular = m OC are superradiant unstable to those particular = m modes. That is, their perturbations have frequencies with Im(ω) > 0 and Re(ω) < mΩ. On the other hand, BHs below a particular OC are damped and thus stable (when perturbed these BHs return to equilibrium via the emission of QNMs with Im(ω) < 0 and Re(ω) > mΩ).
Having studied the OCs for scalar modes with = m, we now turn to consider one particular mode throughout a region of the parameter space to gain more insight into the stability properties of these black holes. A natural mode to consider is the = m = 2 one, as this is the mode with the largest value of the growth rate Im(ω) found in our study. The imaginary and real parts of the = m = 2 scalar mode frequencies are plotted in Fig. 4, and the imaginary and real part of the associated angular eigenvalues is shown in Fig. 5. These quantities are plotted as a function of the dimensionless horizon radius r + /L and rotation a/L and they define a 2-dimensional surface. To extract more efficiently the relevant physics, we plot in the right panel of Fig. 4 is the real part of the superradiant factor Re( ) = (Re(ω) − mΩ h ) /(4πT h ), as introduced in (2.17). In all these plots the blue curve is the = m = 2 OC already identified in the phase diagram of Fig. 3. To guide the eye (when appropriate) we draw an auxiliary plane with a grid that intersects the physical 2-dimensional surface along the OC and that has Re( ) = 0, Im(ω) = 0, and Im(λ) = 0. We also plot some black curves at constant radius r + /L. In the left panel of Fig. 4, modes that are above the auxiliary plane grid are superradiant unstable modes. In the right panel of Fig. 4 and in the left panel of Fig. 5 they correspond to the surface region below the auxiliary plane grid. Finally, in the right panel of Fig. 5 these unstable modes are described by the surface region "below" the blue line. In the four plots, the superradiant unstable surface region is a 2-dimensional surface bounded by the superradiant OC (blue line) and by the extremality curve (where the black curves at constant radius end). 9 In all these plots, the surface region that starts at the blue OC that is complementary to the unstable region describes the QNMs of the Kerr-AdS BH. An important feature of the gravitational scalar superradiant instability concerns the order of magnitude of its timescale τ ∼ 1/Im(ω). Inspecting the data we find that the maximum growth rate of the instability is reached in a neighborhood of the point {r + /L, a/L} max {0.445 ± 0.020, 0.589 ± 0.020} where the frequency is given by ωL ∼ 1.397 + 0.032 i. So, the maximum growth rate for the scalar superradiant instability and the gauge invariant properties of the BH where it is attained are This maximum is denoted with a large red dot in the plots of Fig. 4 and Fig. 5. Note that this maximum occurs close to extremality but not at it. In particular, if we plot the instability growth rate as a function of the rotation parameter a/L at fixed radius (e.g. r + /L = 0.445), we find that, typically, starting from the onset the instability timescale first increases, reaches a maximum for a/L close to extremality, and then decreases as we approach the T h = 0 Kerr-AdS BH.
Consider now the gravitational vector modes which obey the BCs (2.21). The left panel of Fig. 6 displays the phase diagram of Kerr-AdS BHs with the OCs for the = m = 2, 3, 4, 5 vector modes displayed (again, only curves with vanishing radial overtone are shown). As in the scalar case, BHs that are above a particular = m vector OC are superradiantly unstable to modes with those particular values of = m, while BHs below a particular OC are stable to the associated modes. In the right panel of Fig. 6 we plot the angular eigenvalue λ along the OC for vector modes.
The large black points at R + /L = 0, in both plots of Fig. 6, describe the vector normal modes of global AdS, namely [19,21], Together with the superradiant onset condition (with p = 0 and = m) these normal modes give the black points of Fig. 6, As in the scalar case, the vector OCs always have Ω h L > 1 but contrary to the scalar case, these curves always end at extremality and the OCs for different = m never cross each other. In particular, this means that a BH that is unstable to = m = 2 modes must also be unstable to all = m ≥ 3 modes. As = m grows, the curves hit extremality at a higher areal radius R + /L and they approach the ΩL = 1 line. Modes with m → +∞ reach extremality only in the limit R + /L → +∞.
To discuss details of the superradiant and quasinormal modes of the vector sector, we focus again our attention in the = m = 2 case. The superradiant and QNM properties can be read from the plots of Fig. 7 (imaginary and real part of the frequencies) and in Fig. 8 (imaginary and real part of the angular eigenvalues). We use a similar color coding and visualization angle as the ones used in the scalar case. Therefore, in all these plots the blue curve is the OC already studied in Fig. 6; again the auxiliary plane with a grid intersects the physical 2-dimensional surface along the OC and helps visualizing the separation between unstable superradiant modes (Im(ω) > 0 and Re(ω) < mΩ) and damped QNMs (Im(ω) < 0 and Re(ω) > mΩ); and we plot some black curves at constant radius r + /L. It follows that in the left panel of Fig. 7 the unstable modes are in the upper region between the blue OC and extremality, while in right panel they are in the lower region (that we do not show it in all its extension). The upper region of the left panel of Fig. 8 shows the imaginary part of the eigenvalues of the QNMs (we do not show the upper surface in its full extension but its completion should be clear from the continuation of the interrupted black curves with constant r + /L = 0.5 and r + /L = 0.6).
In the plots of Fig. 7 and Fig. 8   instability and the gauge invariant properties of the BH where it is achieved are In general, e.g. moving along a constant r + /L, we find that the maximum of the vector superradiant instability is achieved much closer to extremality than in the scalar case. This property is probably related to the fact that the vector OC ends at extremality, as opposed to the scalar OC.
Comparing the properties of the maximum unstable cases (4.7) and (4.10), we see that the instability growth rate of the scalar and vector sectors is of the same order, with the maximum growth rate in the vector sector being approximately twice stronger than in the scalar sector. Moreover, the most unstable case in the vector case occurs for a Kerr-AdS BH that is smaller (i.e. with smaller gauge invariant areal radius R + /L) but rotates faster than the Kerr-AdS BH where the scalar instability is highest.
Finally, note that the strength of the scalar or vector gravitational instabilities can can be orders of magnitude higher than the strength of the same superradiant instability sourced by a scalar field perturbation [48,49].

Large AdS limit and comparison with special QNMs in asymptotically flat cases
As we will discuss in section VI, the slowly decaying QNMs in Kerr-AdS play a key role in the fluid/gravity correspondence. These modes have a particularly appealing interpretation in terms of a relativistic hydrodynamic problem naturally induced at the AdS boundary. This correspondence also indicates that rich and complex hydrodynamic phenomena have counterparts in the gravitational theory, as recently demonstrated in [23,52,24]. Such a remarkable, and previously unexpected, phenomena displayed by gravity in the AdS context raises the question of what analogues to hydrodynamic behavior arise in general scenarios. Studying such question is beyond the scope of this work (for recent works related to the gravity/hydro connection in AF settings see e.g. [53,54]); however, as we here are concerned with QNMs we can explore the connection of hydrodynamical modes in AdS with relevant ones in AF spacetimes. To this end, we examine in particular the purely-imaginary QNM mode (often called "shear mode") in the limit r + /L → 0 for the non-spinning case, see left panel of Fig. 9. In this limit one makes contact with its possible asymptotically flat counterpart describing QNMs of a Schwarzschild black hole. Interestingly, we find the result obtained coincides with the "algebraically special" QNM mode. Furthermore, we can look at the profile of this mode, as we change the cosmological constant. It turns out it is very localized around the horizon (becoming more and more localized as we lower the cosmological constant), perhaps indicating that the dynamics involved here does not feel the boundary in any special way, see right panel of Fig. 9. At this stage we stress this does not necessarily imply complex hydrodynamic phenomena has a gravitational analogue in AF cases as has been shown to be the case in the AdS case. Nevertheless this is certainly a tantalizing observation deserving further exploration.

Superradiance and black holes with a single Killing field
In the previous sections we confirmed that Kerr-AdS BHs with Ω h L > 1 are unstable to superradiance. An interesting observation is that at the onset of the superradiant instability there is an exact zero mode with ω = mΩ h and Im ω = 0. This zero mode is special because it is invariant under the horizon-generating Killing field Consequently it is regular on both the past (H − ) and future (H + ) horizons (generic perturbations can be made regular on the future or past horizons, but not both). In these conditions and for a given m, [39] proposed that, in a phase diagram of stationary solutions, the OC of the instability should signal a bifurcation or merger of the Kerr-AdS BH with a new family of BH solutions that are stable to superradiant modes with the given m and that preserve the same isometry of the superradiant onset mode (see also the nice discussion in [98]). That is, these new BHs have a single Killing vector field (KVF); the helical Killing field K = ∂ T + Ω h ∂ Φ . In the context of superradiance of a scalar field, BHs with a similar helical single KVF that merge with the Kerr-AdS family have scalar hair orbiting around the central core. Examples of such hairy BHs were explicitly constructed perturbatively and non-linearly in [40] 10 . Given this explicit proof of existence in the scalar field case, it is natural to expect that a similar new family of single KVF BHs with "lumpy gravitational hair" merge with the Kerr-AdS BH at the OC of gravitational superradiance. The existence of such purely gravitational single KVF BHs was first proposed in [39] and contact between these BHs and geons was made in [19]. In this section we will give the explicit construction (omitted in [19]) that leads to the leading order thermodynamics and properties of these BHs. Perhaps the most important consequence of this study is that Kerr-AdS BHs are not the only stationary BHs of Einstein-AdS gravity [40,19]. 11 We can discuss some of the main properties of the single KVF BHs [40,19] in terms of general arguments. Recall again the main properties of superradiance in global AdS. A mode e −iωT +imΦ can increase its amplitude by scattering off a rotating BH with angular velocity Ω h satisfying ω < mΩ h . In asymptotically global AdS spacetimes, the outgoing wave is reflected back onto the BH and scatters again further increasing its amplitude. This multiple amplification/reflection leads to an instability. The process decreases Ω h and eventually results in a BH with "lumpy hair" rotating around it. Such a BH is invariant under just a single Killing field which co-rotates with the hair, K = ∂ T + Ω h ∂ Φ . Thus, the BH is stationary (periodic) but not time symmetric nor axisymmetric. However, it does not violate the rigidity theorems [36,37,38]. Indeed, these theorems assume the existence of a Killing vector, typically ∂ T , that is not normal to the horizon, and prove that a second Killing field ∂ Φ must then be present. Such a BH is thus time symmetric and axisymmetric. The single KVF BHs evade the primary assumption of the rigidity theorem because in this case K = ∂ T + Ω h ∂ Φ is normal to the Killing horizon.
As stated previously, single KVF BHs and horizonless boson star solutions of this type with scalar hair have been constructed perturbatively as well as numerically at the full nonlinear level in [40]. Alternatively, the leading order description of these BHs can also be found using a thermodynamic analysis [55,56,57,40] similar to the one done below. The full nonlinear result confirms that this thermodynamic construction gives accurate leading order results 12 . For small charges the single KVF hairy BHs exist in a region of the phase diagram that is bounded by the OC of scalar superradiance and by the boson star curve.
In the purely gravitational sector of Einstein-AdS theory that we discuss here, the gravitational analogue of the horizonless boson stars are the geons constructed in [19]. Using the aforementioned thermodynamic model we will conclude that single KVF BHs exist in a region of the phase diagram that is bounded by the OC of gravitational superradiance and by the geon curve. 13 We are ready to start the leading order thermodynamic construction of the single KVF BHs. We first review the geon and Kerr-AdS solutions, then we construct the single KVF BHs by placing a small Kerr-AdS BH on the top of a geon.
Geons are classical lumps of gravitational energy, with harmonic time dependence e −iωT +imΦ , in which the centrifugal force balances the system against gravitational collapse [19]. They are horizon-free, nonsingular, asymptotically globally AdS, and can be viewed as gravitational analogs of boson stars. Each geon is specified by , which gives the number of zeros of the solution along the polar direction, and azimuthal quantum number m. It is a one-parameter family of solutions parametrized e.g. by its frequency. At linear order, a geon is a small perturbation around the global AdS background and its possible frequencies are given by the AdS normal modes, namely (4.5) in the scalar sector, and (4.8) in the vector sector. The energy and angular momentum of the geon are related by E g = ω m J g + O(J 2 g ); they have zero entropy S g = 0 and undefined temperature, and they obey the first law of thermodynamics, dE g = ω m dJ g . 14 Consider now the Kerr-AdS BH. For small E and J (i.e small r + /L expansion), the leading original definition of stationarity to accommodate these novel periodic BHs as members of the stationary class of solutions. 12 A similar thermodynamic model was introduced and proved to be correct, when compared with the exact non-linear results, also in the charged superradiant systems discussed in [55, 56, 57] 13 The hairy BHs of [40] could be constructed non-linearly because they depend non-trivially only on the radial direction while the gravitational single KVF BHs we discuss here have an additional non-trivial dependence on the polar angle. It is challenging to solve the associated coupled system of PDEs and we leave its construction for future work. 14 Back-reacting to higher order each of the individual normal modes of global AdS we approach the full nonlinear geon, but we do not need this knowledge for our argument [19]. and next-to-leading order thermodynamics of this solution is which obeys the thermodynamic first law, dE K = Ω h dJ K + T h dS, up to next-to-leading order.
We can now construct perturbatively the single KVF BH of the theory by placing a small Kerr-AdS BH at the core of the geon. The associated single KVF of the solution is inherited from the geon component of the system. To argue for the existence of this solution and to find its thermodynamic properties we can use a simple thermodynamic model where the leading order thermodynamics of the single KVF BH is modeled by a non-interacting mixture of a Kerr-AdS BH and a geon Absence of interaction between the two components of the system means that the charges E, J of the final BH are simply the sum of the charges of its individual constituents: In this mixture, the Kerr-AdS component controls the entropy and the temperature of the final BH (since by definition the geon has no entropy and has undefined temperature). The single KVF BH chooses the partition of its charges between the geon and the Kerr-AdS components in such a way that the total entropy S of the system is extremized. Indeed, maximizing S = S K (E − E g , J − J g ) with respect to J g and using the first laws for the geon and for the Kerr-AdS, we find that the partition is such that the angular velocities of the two components are the same, Ω h = ω m , i.e. the two phases are in thermodynamic equilibrium. Actually, there is a much simpler way to derive this result. Since the geon has only one Killing field, K = ∂ T + (ω/m)∂ Φ , and we place a Kerr-AdS BH with a Killing horizon at its centre, the geon's Killing field must coincide with the horizon generator of the single KVF BH.
The non-interacting and equilibrium conditions together with the leading order thermodynamics of the Kerr-AdS BH and of the geon yields that the final distribution of the charges among the system's constituents and the entropy and temperature of the single KVF BH are, respectively, Moving down from this curve, the Kerr-AdS contribution weakens and the leading order thermodynamics of the system is increasingly dominated by the geon component. In the limit where r + → 0, the lower bound curve of single KVF phase is expected to be the geon curve. This discussion is best illustrated in Fig. 10, where we represent the phase diagram associated to the = m = 2 solutions of the scalar sector with frequency ω = ω s given by (4.5).
Note that there is a region in the phase diagram (the blue/gray shaded region in Fig. 10) where the Kerr-AdS and single KVF BH families coexist, i.e. the present system provides the first example of non-uniqueness in Einstein gravity in four dimensions. The two families of BHs can have the same mass and angular momentum but different entropy.
As emphasized previously, in the scalar field superradiant system of [40], and in the charged superradiant system of [55,57], the available full non-linear results confirm that the thermodynamic model we use also here gives the correct leading order thermodynamic properties of the system. We leave for the future the explicit non-linear construction of the single KVF BHs.
We postpone for Sec. 8 the discussion of the stability properties of the single KVF BHs and the role they might have in the time evolution and endpoint of the superradiant instability of the Kerr-AdS BH.

Hydrodynamic thermalization timescales in the AdS 4 /CFT 3 duality
In the context of the gauge/gravity duality, a black hole is dual to a thermal state in the holographic quantum field theory (QFT). Moreover, QNMs are fundamental entities in this correspondence since the QNM frequencies in the bulk black hole describe the thermalization or relaxation timescales in the dual QFT. This map was first proposed in [3,4] and later it was understood and established that the QNM spectrum of a given field perturbation coincides with the poles of retarded correlation functions of the gauge theory operator that is dual to the perturbation at hand [5,6,7]. This was done in the framework of linear response theory appropriate for describing linearized fluctuations of any wavelength about AdS backgrounds as long as the perturbation amplitude is small. A particularly relevant family of perturbations are the lowest QNMs, i.e. those with small frequency whose wavelength is large compared to the thermal scale of the field theory. The relaxation timescales of these modes have a hydrodynamic description and can be computed studying perturbations of the Navier-Stokes equation that describes the hydrodynamic regime of the holographic QFT [7,8,9,10]. These hydrodynamic modes are also captured by the fluid/gravity correspondence which is a formal one-to-one map between Einstein's equations in AdS and non-linear hydrodynamic equations [58,59]. It follows from a perturbation theory analysis where the small expansion parameter is the ratio of the mean free path of the theory (i.e. the thermal scale) over the typical variation wavelength of the fluid variables and gravitational field. With respect to linear response theory it has the advantage that it captures also non-linear physics but it is restricted to long wavelength physics. The two regimes therefore complement each other and intersect in a corner of the phase space corresponding to linearized long wavelength perturbations [59]. These are precisely the hydrodynamic QNMs that we want to study in this section.
A particular example of a gauge/gravity duality is the AdS 4 /CFT 3 correspondence, whereby supergravity on the Kerr-AdS×S 7 background is dual to a thermal conformal field theory (CFT) on the holographic boundary of the global AdS geometry. In this case the Kerr-AdS black hole is dual to a thermal state with a rotational chemical potential in the CFT 3 that is formulated on a sphere.
In this section we aim to compare the long wavelength gravitational QNMs of Kerr-AdS with the hydrodynamic relaxation timescales of the dual CFT 3 . First, in Section 6.1 we compute the hydrodynamic modes both perturbatively and numerically and later, in Section 6.2, we compare them with with the long wavelength gravitational QNMs. The excellent match that we find provides a further confirmation of the holographic interpretation of the QNM spectrum, of the shear viscosity to the entropy density bound, η/s = 1/(4π), and ultimately of the correspondence itself. Not less importantly, it provides the first non-trivial confirmation that the Robin boundary conditions for the Teukolsky gauge-invariant variable derived in the companion paper [30] are indeed the ones that we must impose if we want the perturbations to preserve the asymptotic global AdS structure of the background. Indeed, had we chosen different BCs, e.g Dirichlet or Neumann BCs, and the QNM spectrum would not match the hydrodynamic timescales.

Hydrodynamic thermalization timescales
The conformal boundary of the Kerr-AdS geometry is the static Einstein universe R t × S 2 with line element that this time we write as where X is related to the standard polar angle on the sphere introduced in (2.6) by X = cos Θ. The CFT 3 is described by an holographic stress tensor T bc which can be found using, e.g. the formulation of Haro, Skenderis and Solodukhin [17]. We first introduce the Fefferman-Graham coordinate frame {T, z, X, Φ} whereby the Kerr-AdS geometry can be recast in an asymptotic expansion around the holographic boundary z = 0 (r = ∞) as with ds 2 ∂ defined in (6.1). The coordinate transformation that takes Kerr-AdS in the Chambers-Moss frame into the FG frame is obtained as an expansion in z, with the successive terms of the expansion being fixed by requiring that g zz = L 2 /z 2 and g zb = 0 (b = T, X, Φ) at all orders. Up to the order relevant for our analysis, this FG coordinate transformation is explicitly given by 3) The leading terms in these expansions are fixed by our choice of conformal frame, namely we want the normalization where g T T = −1 and the sphere has radius L 2 in the boundary metric ds 2 ∂ . On the other hand the azimuthal coordinate transformation guarantees that the conformal frame does not rotate.
The holographic stress tensor can be read from the h 3 contribution of the expansion (6.2) via [17] T bc = 3h 3 16πG 4 , (6.4) where b, c run over the boundary metric coordinates {T, X, Φ}. This stress tensor has the form of a perfect fluid with energy density ρ, pressure p, and fluid velocity u given by where are the angular velocity Ω ∞ of the fluid, and the ratio γ −1 = T T between the fluid temperature T and the local temperature T (this gives the redshift factor relating measurements done in the laboratory and comoving frames), and ∂ T , ∂ Φ are the Killing vectors corresponding to the isometries of the boundary background (6.1),. Further, u 2 = −1 and the equation of state ρ = 2p follows from the fact that the holographic QFT and its fluid are conformal which implies that the stress tensor is traceless. The stress tensor is conserved with respect to (6.1), ∇ b T bc = 0, since there are no sources (e.g., scalar or Maxwell) in our system. Our bulk background is stationary and therefore the boundary fluid is also in stationary equilibrium fluid configuration with rigid roto-translational motion. Our choice for the fluid velocity definition is such that it obeys the Landau gauge condition This condition guarantees that the stress tensor components longitudinal to the velocity give the local energy density, in the local rest frame of a fluid element [60]. A generic perturbation of the stationary fluid configuration will drive the system away from equilibrium and dissipation must be included to study the evolution of the system. This dissipative contribution to the total holographic stress tensor is encoded in the term Π bc (this follows from a gradient expansion of Einstein equations around AdS in the regime where the thermodynamic variation lengthscales are much larger than the thermal scale of the stationary background [58,59]), are the shear viscosity tensor σ bc , the fluid expansion ϑ, and the is the projector P bc onto the hypersurface orthogonal to u. The quantity η is the shear viscosity. Since the fluid is conformal, its stress energy tensor must be traceless (i.e. the conformal anomaly is proportional to T c c and vanishes 15 ). Consequently the fluid must have vanishing bulk viscosity. Also, the Landau frame condition (6.6) implies u (0) a Π bc = 0 which discards a possible heat diffusion contribution to the first order dissipative stress tensor (i.e. in this frame all the dissipative contributions are orthogonal to the velocity field) [60].
We recall that a precise statement for the validity of the hydrodynamic regime of dual system can be made as follows. The mean free path of a theory is typically given by the ratio of the shear viscosity to the energy density, mfp ∼ η ρ . We are working with a conformal theory so the associated fluid equation of state is ρ = 2p and the viscosity to entropy bound is saturated, η = s/(4π) [61]. For any fluid we also have the Euler-Gibbs relation ρ + p = T s, where the local temperature is related to the fluid temperature (dual to the black hole temperature T = T h ) by the Lorentz factor. Therefore we can write mfp ∼ η ρ ∼ 3 The hydrodynamic approximation is valid for when the thermodynamic quantities of the fluid and of its perturbations vary on lengthscales that are much larger than mfp , namely r + L 1 , and a L 1 , (6.8) where to get the first relation we used the fact that 0 < γ −1 ≤ 1 and that the temperature scales as T h ∼ r + for large radius black holes − see (2.4) − while the second relation follows from the fact that the background pressure p (0) is not a constant and its gradient scales with the rotation parameter in AdS units. According to the holographic dictionary, the fluid temperature is identified with the Hawking temperature T h of the black hole and it follows from the previous discussion that the angular velocity of fluid Ω ∞ = a/L 2 is precisely the shift in the azimuthal coordinate such that the (non-dynamical) background on which the fluid flows is static. On the other hand, the viscosity is given in terms of the horizon radius of the bulk black hole as This is a universal relation for any fluid that is holographically dual to a black hole of Einstein-AdS 4 theory. It follows from the celebrated viscosity to entropy density ratio of the theory namely η = s/(4π) [61]. This is a constitutive relation that is independent of the rotation of the fluid since it follows from measuring quantities in the rest frame of the fluid. Namely we can write the entropy density as s = S/V = Sρ (0) /E = 2p (0) S/E which yields (6.9) after using the relations for the static black hole entropy and energy, S = πr 2 + and E = r + 1 + r 2 + /L 2 /2, and taking the hydrodynamic limit r + /L → ∞.
The hydrodynamic equations of motion for the perturbed fluid, that will ultimately quantize the relaxation timescales of the system, follow from the conservation of the total stress tensor, These equations can be written as a set of two family of equations, namely the relativistic continuity and Navier-Stokes equations, 16 To study the perturbations of these fluid equations, we use the fact that ∂ T and ∂ Φ are isometries of the background to write the most general perturbations for a conformal fluid as a sum of the following Fourier modes The velocity normalization Plugging these fluctuations in the linearized version of the hydrodynamic equations (6.11) we get the equations of motion (EoM) that the fluid perturbations {δP, δu X , δu Φ } have to obey. We solve these equations exactly using numerical methods like those we use to solve the gravitational equations. In addition, to get extra physical insight and check the numerics, we also find a perturbative explicit analytical expression for the fluid quantities of interest.
To solve the linearized hydrodynamic equations using a perturbative method [55,56,40,57], we assume a double expansion in the shear viscosity and in the rotation, both for the fluid perturbations introduced in (6.12), {Q (f ) (X)} = {Q (1) , Q (2) , Q (3) } ≡ {δP/L, δu X , δu Φ }, and for the perturbation frequency ω: and solve progressively (6.10) or (6.11) in a series expansion in η/L 3 and a/L. For our purpose it will be enough to go up to third order (p = 3) in the rotation expansion.
Inspecting the EoM at leading order O η 0 , a 0 , we immediately conclude that we have to split our analysis into two family of modes, namely the scalar and vector modes. The latter have ω 0,0 = 0 and perturb the fluid velocity but not the pressure, while the former have ω 0,0 = 0 and perturb all fluid variables. At this order rotation is absent and the hydrodynamic modes have an expansion in terms of the scalar S i and vector V i Kodama-Ishibashi harmonics (which are both related to the associated Legendre polynomials [27,10,30]). This is in agreement with the fact that the gravitational QNMs split also into two families as dictated by the two possible global AdS boundary conditions (2.19)-(2.21). As rotation and/or viscosity are turned on these two families naturally continue to follow different paths.
Our main goal is to find the characteristic damped oscillation frequencies of the fluid. We leave the details of our computation to Appendix B and give here only its relevant outcome, namely the hydrodynamic CFT thermalization frequencies that can propagate in the CFT 3 . The frequencies of the hydrodynamic scalar modes are: 17  while the frequencies of the hydrodynamic vector modes are: (6.16) In these expansions (and associated Figs. 11, 12 below) we assume the relation (6.9) for the viscosity. When the rotation vanishes, (6.15) and (6.16) reduce to the hydrodynamic frequencies first computed in [10].   As illustrative examples, Figs. 11 and 12 show the regime of validity of the perturbative expressions (6.15) and (6.16) by comparing them against the exact numerical solutions of the linearized hydrodynamic equations (6.11) for the = m = 2 and = m = 3 harmonics in both the scalar and vector sectors. As is evident from the figures, the match is excellent in the small rotation regime as expected. Fig. 11 describes the hydrodynamic scalar modes. For each harmonic there is a pair of solutions, one with positive and the other with negative real part of the frequency. At zero rotation and only in this case, the background has the t − φ symmetry and thus the two  solutions are physically the same: they form a pair {ω, −ω * } related by complex conjugation. Rotation breaks this degeneracy. Fig. 12 describes the hydrodynamic vector modes. These are characterized by having vanishing frequency real part when the rotation vanishes, so there is only one family of solutions for each harmonic.

Long wavelength QNMs and hydrodynamic modes
In the last subsection we computed analytically and numerically the hydrodynamic relaxation timescales. In this section we compare these timescales with the long wavelength gravitational QNMs.
To perform the comparison we recall that the hydrodynamic and gravitational modes are expected to match in the regime of parameters (6.8), namely r + /L 1 and a/L 1. We thus consider a Kerr-AdS black hole with radius parameter r + /L = 100 to do the comparison. A measure of the deviation between the numerical hydrodynamic frequencies (call them ω hydro ) and the numerical gravitational QNM frequencies (call them ω) is given by |1 − ω hydro /ω|. In Fig. 13 we plot this deviation measure as a function of the rotation parameter a/L (a/L < 1 for regular black holes) for a Kerr-AdS BH with r + /L = 100 1. The brown curve (disks) is for scalar modes, while the green curve (squares) is for vector modes. We see that the match between the hydrodynamic and long wavelength QNM frequencies is very good even when the rotation grows large and thus moves away from the hydrodynamic validity regime a/L 1: for scalar (vector) modes the maximum deviation is below 10 −4 (2 × 10 −3 ).
This perfect match when the rotating chemical potential is present is a further confirmation of the holographic interpretation of the QNM spectrum, of the shear viscosity to the entropy density bound, η/s = 1/(4π), and ultimately of the AdS/CFT correspondence itself. 1 Ω hydro Ω Figure 13: Comparison between the long wavelength gravitational QNMs and the hydrodynamic modes (for the later we use the exact numerical results) for r + /L = 100. The brown (disks) curve describes the scalar modes while the green (squares) curve is for the vector modes.

QNMs and superradiance in 5 dimensions
In this section we extend the study of thermalization, quasinormal modes, and superradiance to five dimensions. There are many motivations for doing so. Firstly, there is currently a general interest in studying gravity in higher dimensions, for a modern review see [62]. It is interesting to ask how the solutions to the Einstein equations and their properties vary with D. As the dimension increases, the types of black hole solutions increase dramatically. Some examples of the non-standard black holes possible in higher dimensions are black rings, black Saturns, and black branes. Many of these solutions challenge the intuition gained from studying four dimensions by exhibiting non-uniqueness and a variety of interesting instabilities, some of which likely lead to topology-changing transitions once quantum effects are included. In addition to this very general motivation, we shall see that superradiance in five dimensions is qualitatively very similar to four dimensional Kerr-AdS case. Thus we expect that the intuition we gain from studying superradiance in four and five dimensions will be useful for thinking about other dimensions. Additionally, because of the properties of the particular class of black holes we chose to study, certain aspects of the problem will turn out to be more tractable than the Kerr-AdS case. A second motivation for studying five dimensional asymptotically AdS black holes is that they play an important role in understanding strongly coupled field theories in four dimensions via gauge/gravity duality. In particular, the most well-developed example of this duality is Type IIB string theory on AdS 5 × S 5 spacetimes, which is dual to N = 4 Supersymmetric Yang-Mills. The physics of five dimensional AdS black holes can thus lead to an improved understanding of this specific duality, and more generally about the physics of four dimensional field theories at finite temperature. As in the Kerr-AdS case, large black holes will be particularly interesting in this regard as they will be dual to field theories admitting a hydrodynamic description.

Myers-Perry−AdS black holes with equal angular momenta
The generalization of the Kerr metric to higher dimensions was found by Myers and Perry [63]. It was then further generalized to include negative cosmological constant for the case of five dimensions in [64], and then to arbitrary dimensions by [65,66]. Although our numerical results are for five dimensions only, in the presentation that follows we will keep the dimension general whenever possible. We therefore refer to these black holes as Myers-Perry−AdS (MP-AdS) black holes.
In higher dimensions there are more planes for an object to rotate in than in four dimensions, and therefore these black holes are described by n = (D − 1)/2 angular momenta parameters. For generic choices of the angular momenta, the symmetry of the MP-AdS black hole is R × U (1) n . The first factor is due to time translations, and the n U (1)'s are due to the n independent planes of rotation. For certain choices of the angular momenta, this symmetry can be increased. For odd dimensions, a particularly dramatic enhancement occurs when all angular momenta are equal. In this case, the isometry becomes R × U (1) × SU (N + 1), where D = 2N + 3. For this case, the line element becomes cohomogeneity-1, which is to say that it depends non-trivially only on the radial coordinate. This feature makes the study of linear stability particularly tractable for these black holes, as the linearized perturbation equations can be easily separated and reduced to ODE's. Therefore, in what follows, we shall restrict ourselves to odd dimensions and the equal angular momenta sector of the full parameter space.
Here we introduce the MP-AdS black holes for odd dimension D = 2N + 3 and with all angular momenta equal. The line element is with the metric functions defined as follows: Hereĝ ab is the Fubini-Study metric on CP N . We will adopt the convention that lowercase latin indices run over CP N coordinates, and that hatted tensors are associated with this space. The Fubini-Study metric is Einstein, with the following proportionality constant R ab = 2(N + 1)ĝ ab . (7.4) When the angular momenta are set to zero, this metric reduces to the usual Schwarzschild-AdS metric with the unit sphere written in terms of the Hopf fibration: where here A is related to the Kähler form by 2J = dA. Constant t, r slices have the geometry of homogeneously squashed spheres, with the amount of squashing determined by h(r). The energy, angular momenta, and angular velocity of the horizon are [44]: The horizon-generating Killing field is K = ∂ t + Ω H ∂ ψ . The physics of perturbations of this black hole will depend crucially on Ω H . For Ω H L < 1, K is timelike everywhere outside the horizon.
If Ω H L > 1, then it becomes spacelike sufficiently far away from the horizon, and in particular is spacelike at the conformal boundary. For Ω H L = 1, K is exactly null at the conformal boundary. The metric is described by three dimensional parameters, (L, r M , a). We will find it useful to use instead the parameters (L, r + , Ω H ), where r + is the horizon radius, to describe the solution. Note that these parameters are defined independently of the parameters that appear in the Kerr-AdS case, and also note that the boundary metric is not rotating in these coordinates.
For these equal angular momenta black holes the angular velocity cannot be arbitrarily large and must obey the extremal bound (7.8) In Fig. 14 we plot the domain of the parameter space for the case D = 5. We also plot the lines Ω H L = 1 and r + /L = 1. These divide the sub-extremal parameter space into four distinct regions. This division comes from the analysis of Hawking and Reall [51]. The importance of the line Ω H L = 1 is that black holes with Ω H L < 1 are expected to be stable, whereas those with Ω H L > 1 are expected to be susceptible to instabilities. Additionally, the partition function of these black holes in a grand canonical ensemble becomes ill-defined for Ω H L ≥ 1, as the dual CFT will be rotating faster than the speed of light. The importance of r + /L = 1 is that black holes larger than this are thermodynamically preferred over thermal, rotating AdS in the grand canonical ensemble, whereas smaller black holes are not.

Scalar-gravitational perturbations
We now review the problem of linear perturbations of the above metric. The decomposition we employ was first utilized in [67], where scalar-gravitational perturbations of asymptotically flat MP black holes with equal angular momenta were studied. Metric perturbations may be decomposed according to how they transform under the isometries of the CP N base space. There are three sectors of perturbations to consider: scalar, vector, and tensor. Tensor and scalar field perturbations of the MP-AdS black holes were studied in [39]. A major simplifying feature of five dimensions is that vector and tensor perturbations do not exist, because the associated vector and tensor harmonics do not exist on CP 1 [39,68]. Thus, we need only consider the scalar sector of perturbations in five dimensions. We now briefly review charged scalar harmonics on CP N following [69,67]. First, introduce a charged covariant derivative: That this is the natural derivative operator to consider can be seen from the dimensional reduction of the fibre coordinate in the Hopf fibration. Charged scalar harmonics (with charge m) are then those functions of the CP N coordinates that satisfy Here the eigenvalue is a function of two quantized parameters, (κ, m): where κ = 0, 1, 2..., and m ∈ Z. Charged scalar-derived vectors can be obtained by differentiating, These can be further decomposed into Hermitian and anti-Hermitian parts Lastly, the scalar-derived tensors are given by (7.14) In order to implement the harmonic decomposition of the perturbation, it will be useful to introduce the 1-forms, e 0 = f (r)dt, e 1 = g(r)dr, e 2 = h(r)(dψ + A a dx a − Ω(r)dt). (7.15) The CP N scalar sector of metric perturbations can then be written as With the above parametrization, the CP N dependence in the Einstein equations will separate, resulting in a system of coupled ODE's. These equations are rather lengthy, and depend nontrivially on the particular harmonic (κ, m) under consideration, and we therefore omit their presentation here. The non-trivial way in which the equations depend on the CP N harmonic in question is as follows: For certain values of (κ, m), some of the harmonic tensors do not exist and therefore their coefficient functions are zero. As an example, for N = 1 (D = 5), and m > 0, Y + a = Y ++ ab = Y +− ab = 0, and so therefore the functions f + A , H ++ , H +− do not enter into the perturbation.

Boundary conditions
We now turn to a discussion of the boundary conditions. Boundary conditions must be supplied at the horizon and at the conformal boundary. The appropriate boundary conditions for quasinormal modes are those that correspond to an ingoing perturbation on the future horizon H + , and a normalizable perturbation at the boundary. For the case of asymptotically flat equal angular momenta MP black holes, the requirement of being ingoing at the horizon was translated into boundary conditions in Ref. [67]. The same method applies here, and so we omit a detailed discussion of the horizon boundary conditions.
The boundary conditions at infinity are less straightforward. Here we describe a general method for finding the boundary conditions of a normalizable gravitational perturbation of asymptotically locally AdS spacetimes. Recall that in the Kerr-AdS case the perturbation equations were reduced to a single gauge invariant equation, the Teukolsky equation. There were two steps in the process of finding the boundary conditions at infinity. First, a Frobenius expansion analysis yielded the allowed fall-off's, and then the requiring that the perturbation be normalizable fixed a certain linear combination of the two solution branches.
We wish to generalize this method to a system of n 2nd order ODE's for n functions f i . In will be convenient to convert to a new radial coordinate, z = 1/r. First we demonstrate that z = 0 is a regular singular point of the equations. To do so, change to a new basis of functions q i = z α i f i , and write the equations as If some choice of the α i and overall multiplicative factors can be made such that the coefficient matrices approach finite and non-zero constant matrices in the limit z → 0, then we will have demonstrated that z = 0 is a regular singular point of this system of equations. We do not know of an algorithmic way to determine the α's, but for the problem at hand we have demonstrated that the equations can be put into this form. Then, after changing coordinates to ∂ t = z∂ z , the equations become: This is a 2nd order system of ODE's with constant coefficients, which can be solved via the standard method of writing it as a system of coupled 1st order ODE's: (7.23) and a = 1, ..., 2n. The generic solution (excluding possible logarithmic terms in z) is then a are the eigenvalues and vectors of M . There are 2n coefficients c b , which is expected: for n 2nd order ODE's there should be 2n constants of integration. This is the generalization of the first step in the Frobenius method, in which the two independent branches of a solution to an ODE around a regular singular point are determined. The next and last step is to use physical considerations to determine the c b that correspond to the type of perturbation being studied. For normalizable metric perturbations, the boundary metric is held fixed and the boundary stress tensor is varied. For the purpose of translating this condition into choices of the constants of integration c b , it is useful to put the full metric (background plus perturbation) into Fefferman-Graham gauge. In this gauge the metric takes the form Here D = d + 1 and the lower case Latin indices run over all but the radial coordinate. In order to make the perturbation normalizable we require that it only affects the terms g d and higher in the above expansion. This corresponds to holding fixed the boundary metric and only allowing the metric perturbation to affect the expectation value of the stress tensor of the dual field theory. This requirement will fix n of the constants c b . This method generalizes the usual Frobenius method for finding normalizable fall-off's for fields in AdS.
As an explicit example, we display the fall-off's for D = 5 and for the (0, m) harmonic (for m > 0. For this mode, the non-zero perturbation functions are Here we have introduced the quantities , ∆(r) = g(r) −2 . (7.33) Once these q i functions are known, one can easily find the boundary conditions by expanding the equations near the endpoints.

Numerical results
Here we present our numerical results. We calculated scalar-gravitational QNM frequencies for the five dimensional, equal angular momenta MP-AdS black hole. The possible perturbations are parametrized by two integers (κ, m). We are particularly interested in the onset of superradiant instabilities, which was discussed earlier in Sec. 4.2. We remind the reader that superradiant modes are characterized by the condition Re(ω) < mΩ H , and linear instabilities by the condition that Im(ω) > 0. For superradiant instabilities in AdS, these two conditions are satisfied simultaneously as the mode becomes unstable. In considering the onset of these instabilities, we will find it useful to consider the function Ω H,onset (r + /L), which for a given (κ, m), is the value of the rotation such that ω = mΩ H,onset . Before presenting our results on scalar perturbations, we briefly review the results for tensor perturbations [39], which exist for D odd and D ≥ 7. There an infinite number of superradiant instabilities were found, all for Ω H L > 1. Lower m modes become unstable for larger values of the rotation, and in the limit m → ∞, the critical rotation approaches Ω H,onset → 1. Also, Ω H,onset depends only very weakly on the size of the black hole, and in particular the ordering of Ω H,onset 's for different m's is independent of the size of the black hole. For a given m the instability terminates once the black hole is taken to be sufficiently large. We will find that some scalar perturbations behave qualitatively differently than these tensor perturbations.
In Fig. 15 we plot the onset of the superradiant instability for (0, m) modes, for a range of m. We calculated these threshold unstable modes using the Newton-Raphson method presented in Sec. 3. As a check on our results, it is useful to first consider small black holes. For small black holes, the onset of the instability can be very easily predicted by first setting ω = mΩ H for the onset of superradiance, and then also setting ω = ω AdS , where ω AdS is the normal mode frequency of AdS. For (0, m) modes, this results in Ω H,onset L = 1 + 2 m , (7.34) Notice that for m → ∞, Ω H,onset → 1. As the size of the black hole is increased, the onset curves for different m begin to cross. This behaviour was also observed for scalar perturbations of Kerr-AdS in Sec. 4.2, and is qualitatively different from the behavior of the tensor instabilities discussed above. Although the onset curves cross as the size is increased, the fact that the m → ∞ instability hugs the line Ω H L = 1 indicates that the "m = ∞ mode will never be crossed", i.e. arbitrarily large m-modes will be the first ones to go unstable as Ω H L is increased, even for arbitrarily large black holes. This same phenomenon occurs for scalar perturbations of Kerr-AdS, and is discussed in more detail in Sec. 4.2. Another difference between these scalar instabilities and the tensor ones is that, for a given m, the scalar onset curves extend for arbitrarily large black holes, whereas the tensor curves terminate. For very large black holes, we can examine the approach of the onset curve to it's limiting value of Ω H L → 1. Interestingly, the approach is power-law, with the exponent independent of m: In Fig. 16 we plot the onset of the superradiant instability for (1, m) modes. These instabilities are evidently qualitatively different from the (0, m) instabilities in that a) the onset curves do not cross, and b) for a given m, the instabilities do not persist for arbitrarily large black holes. In fact, these modes appear to be qualitatively very similar to the tensor modes on CP N for N ≥ 2.
In Fig. 17 we plot contour plots of the real and imaginary parts of ω for the (0, 2) mode. The black dashed line corresponds to the onset mode, ω = mΩ H , and the black dot corresponds to the data point with the largest positive value of Im(ω). This point is near the extremality bound, and lies at the end of the grid used to scan the parameter space, and so it is likely that true maximum either lies along or very near the extremality curve. This point is given by (r + /L, Ω H L) = (0.579, 1.548), ωL = 2.481 + 0.039i. (7.36) Lastly, we remark on the importance of the onset modes, ω = mΩ H , on the phase diagram of black hole solutions in five dimensions. In Sec. 5 it was discussed how threshold unstable modes with ω = mΩ H signified new branches of single KVF black hole solutions. We expect much of this analysis to carry over to five dimensions. It would be interesting to study these putative solutions further. Also, as in the Kerr-AdS case, the endpoint of the superradiant instability remains a very interesting open question, especially given the crossing of the onset curves observed in the (0, m) scalar sector of perturbations.

Hydrodynamic thermalization timescales in the AdS 5 /CFT 4 duality
We now turn to study the hydrodynamic QNM's of the five dimensional MP-AdS black hole. As discussed in Sec. 6, this serves as a powerful check on our numerics, the hydrodynamic approximation, and more generally, gauge/gravity duality itself. Compared to the Kerr-AdS case, the hydrodynamic approximation for the cohomogeneity-1 MP-AdS black holes is conceptually simpler and has a wider range of validity. Recall that in the Kerr-AdS case, the pressure was a function of the angular coordinate θ, and this introduced another length scale into the approximation which limited it's domain of applicability (although the agreement turned out to be excellent even for large rotations). In contrast, the pressure for these five dimensional black holes is constant, and the approximation is valid for all rotations Ω H L < 1. As is often the case, we will find that the hydrodynamic approximation agrees excellently with our numerical for reasonably large black holes. The general theory of hydrodynamic modes was reviewed in Sec. 6, and so here we only remark on the differences that occur for these five dimensional equal angular momenta MP-AdS black holes. The boundary metric is now h µν dx µ dx ν = −dt 2 + L 2 (dψ + A a dx a ) 2 +ĝ ab dx a dx b . (7.37) The energy density, pressure, and fluid velocity are related to the black hole parameters via We will also need to use the viscosity to entropy relation, η = s/4π, which in our conventions yields η = r 3 + /(2L). Hydrodynamic modes are obtained through perturbations of the stress tensor that are traceless and divergenceless. The CP N dependence of these equations can again be separated using the charged harmonics introduced above. Since our numerical data is for scalar perturbations, we will restrict our attention to hydrodynamic scalar modes. The decomposition of the fluid variables is: The perturbed quantities {δp, δu t , δu ψ , δu + δu − } and ω can then be solved for by using the conservation of the stress tensor. As the expressions are rather lengthy, and depend nontrivially on the harmonic under consideration, we omit their full presentation here. Given the difficulty of obtaining analytic predictions for black hole QNM's, we will however include the expressions in an expansion about Ω H = 0. For D = 5 and the (0, m), harmonics (with m ≥ 2), and to first order in both L/r + , Ω L = Ω H L, the hydrodynamic modes are Next, we compare our numerical data for large black holes with the hydrodynamic approximation. In Fig.'s 18 and 19, we fix r + /L = 100 and plot the analytic prediction of the hydro modes against our numerical data for two choices of CP 1 harmonics: (0, 2) and (1, 1). We work to leading order in L/r + and find excellent agreement for all Ω H L < 1. The typical error of the hydro approximation is 10 −4 , which is exactly what we'd expect as the next term in the expansion comes in at O(L 2 /r 2 + ) ∼ 10 −4 . Although we find excellent agreement between the hydrodynamic approximation and our numerical data, the approximation fails to capture the superradiant instabilities. As mentioned earlier, all superradiant instabilities necessarily have Ω H L > 1. For these black holes, the boundary is rotating faster than the speed of light, and the hydrodynamic approximation shouldn't be expected to be valid. So, while the hydrodynamic approximation has again proven to be an excellent approximation scheme, in this case it fails to capture the instabilities. Im Ω L Figure 19: Comparision of the leading order hydrodynamic approximation for (κ, m) = (1, 1) modes. Here r + /L = 100. Once again, we find excellent agreement between our data and the hydrodynamic prediction.

Discussion and open problems
We studied the most general linear perturbations of the Kerr-AdS BH and of the equal angular momentum Myers-Perry-AdS BH in D = 5. We imposed asymptotic BCs that preserve the conformal metric [17,30]. These BCs also guarantee that the energy and angular momentum fluxes across the asymptotic boundary vanish. Using a novel numerical approach, which we believe might also be useful for other applications, we computed the QNM spectrum of these BHs and the growth rate of their instabilities. The only linear instability that we find in the D = 4 and D = 5 stationary BHs have a superradiant nature and they appear only in BHs with Ω h L > 1. We focused on these spacetimes because of their interest for AdS 4 /CF T 3 and AdS 5 /CF T 4 dualities formulated on the static Einstein Universe, i.e. on the sphere. Higher dimensional stationary BHs with D ≥ 6 were not considered here but they should have novel features that might be worth investigating. Indeed, it is established that D ≥ 6 stationary BHs are also unstable to the ultraspinning instability, whose onset was identified in [70] (this instability was first studied in asymptotically flat stationary BHs in D ≥ 6 [71,72,73,67,74]). However, it is still an open question whether another instability that is present in D ≥ 6 vacuum stationary BHs, namely the bar-mode instability [71,75,76,77], is present in AAdS rotating BHs. The onset of the superradiant instability is an exact zero mode that is invariant under the horizon-generating Killing field of the Kerr-AdS (MP-AdS) BH. On the shoulders of an idea originally proposed in [39], we argued that, in a phase diagram of stationary solutions, the superradiant onset curve is a bifurcation line to a new family of BH solutions with a single Killing field that span a region that is further limited by the geon family constructed in [19]. We have constructed perturbatively the leading order thermodynamics of these novel BHs using a simple thermodynamic model [55,56,57,40]. In the future, it is important to construct explicitly (numerically) these single KVF BHs and geons at full nonlinear level to confirm the ideas here discussed (in the context of scalar superradince, similar single KVF BHs and boson stars have already been constructed nonlinearly in [40]. Their properties are in agreement with the thermodynamic model we use here, in the regime of small charges). It is worth emphasizing that these single KVF BHs are periodic but not time symmetric neither axisymmetric and their existence shows that the Kerr-AdS and MP-AdS BHs are not the only stationary BHs of Einstein-AdS theory (as discussed in detail before, their existence is not in conflict with the rigidity theorems).
An interesting open question concerns the time evolution and endpoint of the superradiant instability. Before addressing this issue for rotating systems, it is useful to discuss first the situation for global AdS Reissner-Nordström BH (RN-AdS BH) with chemical potential µ that are unstable to charged superradiance. This is the case if the RN-AdS BH is scattered by a charged scalar wave with frequency ω and charge e that obeys ω ≤ eµ. Here, the marginal mode with ω = eµ signals a bifurcation curve, in a phase diagram of static solutions, to a new family of charged BHs with scalar hair that have been explicitly constructed (perturbatively and nonlinearly) in [55,56,78,57]. When they coexist, the entropy of the hairy BHs is always higher than the entropy of the RN-AdS BH with same mass and charge, for fixed e. Moreover, the hairy BHs are not unstable to superradiance since for a given mass and charge (and fixed e) the chemical potential of the hairy BH is smaller than the chemical potential of the RN-AdS BH; therefore the would-be superradiant modes of hairy BHs no longer fit inside the global AdS box. So far we have just discussed the phase diagram of solutions but said nothing about the time evolution of the original RN-AdS BH. The expectation, to be confirmed by a full time evolution, is that the endpoint of the charged superradiant instability in the RN-AdS BH is one of the hairy BHs constructed in [55,56,57]. This follows from the fact that for a given mass and charge (and fixed scalar charge e) the hairy BH has higher entropy and lower chemical potential than the RN-AdS BH. Therefore a time evolution towards the hairy BH is compatible with the second law of thermodynamics and the endpoint would be stable (to superradiant modes with the given fixed e). Given the properties of this charged system, and the obvious similarities with the rotating superradiant system, it is often assumed that we can use the charged system to extrapolate on evolution properties of the rotating system. However, we next argue that such an extrapolation for time evolution properties is not appropriate. To begin, notice a fundamental difference between the charged and rotating systems. The charge of the scalar field e, that enters the superradiant condition ω ≤ eµ, is fixed. However, the azimuthal quantum number m, that enters the superradiant condition ω ≤ mΩ h , is not fixed since the nonlinearities of Einstein equation will excite other m modes during a time evolution. This means that a given single KVF BH, constructed in association with a given m mode, can be at most just a metastable state but never the endpoint of the superradiant instability. This is because the single KVF BH is stable to the particular m-mode but not to other m superradiant modes that are inevitably excited in a time evolution. Therefore the endpoint of the superradiant instability in rotating BHs is not known at all and finding it is one of the most interesting open questions in BH perturbation physics. Not much can be said about it without performing the full time evolution but it is interesting to observe that, typically, stable BHs to a given m-mode are nevertheless unstable to higher m-modes. So one possibility is that the system will evolve to configurations with higher and higher m structure. Another important observation is that only BHs with angular velocity Ω h L < 1 are stable to superradiance, as first proved in [51]. So a natural expectation for the endpoint of the superradiant instability would be a (single KVF) BH with Ω h L < 1. Finding whether such a BH exists requires constructing the single KVF BHs at full nonlinear level. However, in the similar scalar superradiant system of [40], where much of the present discussion about the time evolution also applies, the single KVF BHs of the theory have been explicitly constructed nonlinearly but none of them has Ω h L < 1.
Within the gauge/gravity correspondence, black hole QNMs are dual to thermalization timescales in the dual CFT. An explicit check of this statement is possible in the regime where the CFT admits a near equilibrium, long wavelength effective hydrodynamic description. We have explicitly checked that for BHs with radius much larger than the AdS length (and small rotations in the four dimensional case), the long wavelength gravitational QNMs of stationary BHs match the hydrodynamic relaxation timescales of the dual CFT. This confirms that the holographic interpretation of the QNM spectrum extends to systems with a rotating chemical potential. It is also a further check of the validity of the shear viscosity to the entropy density bound, η/s = 1/(4π), and a non-trivial confirmation that the global AdS BCs derived in the companion paper [30] preserve the conformal boundary.
The damped QNM modes of a BH have a well defined dual CFT interpretation, but not much is known about the holographic interpretation of the superradiant instability (for discussions in this direction see [51]). From the gravity side it is clear that the superradiance instability has to do with a quenched cooling of the system, since increasing the angular momentum very rapidly, cools the system down. It would be interesting to connect this interpretation with a simple CFT model for such a phenomena, where one could perhaps understand the final state of the system. In addition, it would also be important to understand the novel holographic phases or states that are dual to the single KVF BHs and geons that appear in the superradiant context. From a different perspective, in the bulk we have discussed BHs only at the classical level. However, when quantum effects are included, Hawking radiation is also present and it is entangled with spontaneous superradiant emission. These phenomena should have a microscopic or statistical description. Within string theory, certain BHs can be described by a configuration of D-branes. In this context, Hawking radiation can be microscopically understood as the emission of a closed string off the D-branes as a result of the collision of two open strings that are attached to the D-branes (see [79] and references therein). Similarly, superradiant emission has a microscopic description in terms of collisions of fermionic (spinning) left and right moving string excitations, and the superradiant condition ω ≤ mΩ h follows from the Fermi-Dirac statistics for the fermionic open strings [79]. It would be interesting to extend this microscopic description to BHs (like Kerr-AdS) that do not have a D-brane description.
The properties of spheroidal wavefunctions and eigenvalues are known analytically for some time in AF spacetimes [80]. Our (numerical) analysis leaves the corresponding analytical analysis of spheroidal harmonics in AAdS unexplored, but clearly a compelling topic. Specially interesting is the extremal regime a = L, which might be amenable to a full analytic treatment, both in the angular eigenvalue and in the eigenfrequency.
In a similar vein, a detailed analysis of superradiance in extremal, AF and AAdS geometries is seemingly lacking. Quasinormal mode results for the extremal, AF Reissner-Nordstrom geometry uncovered an interesting symmetry between different perturbations [81] which might propagate to superradiant amplification factors and to other geometries. Note that our analysis does not apply to the extremal Kerr-AdS BH because it has a double horizon and thus our BCs are not appropriate. However, as one approaches the extremal BH, superradiance emission persists as first observed for the Kerr BH in [82]. An interesting observation is that in the AF case, when the Kerr BH is extremal and the perturbations have a frequency that saturates the superradiant bound, i.e. ω = mΩ ext h , the radial Teukolsky equation has an exact solution in terms of hypergeometric functions [83]. However, for the Kerr-AdS BH we can no longer solve the radial equation analytically, even in the above particular conditions [84]. Extremal Kerr(-AdS) BHs are also interesting because they have a near-horizon limit where a Kerr/CFT correspondence can be formulated (see e.g. [85,84] and references therein). The study of gravitational perturbations in the Kerr(-AdS) near-horizon geometries was done in [85,84]. It is interesting to note that all frequencies in the near-horizon geometry correspond to the single frequency ω = mΩ ext h in the original full geometry. Probably this is the reason why no signature of superradiance is found in the near-horizon geometry. It might be useful to explore further this system since its radial solutions are analytical.
We conclude this discussion section with some important general remarks concerning perturbations of AAdS spacetimes. One might wonder whether Robin boundary conditions, such as the ones used throughout this paper, lead to a well defined initial value problem for fields propagating in arbitrary asymptotically AdS backgrounds. This has been shown to be the case for the propagation of a real scalar field in [86,13], where no assumption about the stationarity of the background or separability of the wave equation was made. The proof given there can be readily extended to complex of multi-component fields, including the gravitational perturbations discussed here. In the absence of a linear instability, one might think that linear perturbations about Kerr-AdS will decay exponentially with time, in a manner dictated by the QNM spectrum. However, it turns out that this is not the generic case, and indeed, depending on the smoothness of the initial data, the decay might be a lot slower than that. A simple argument suggest logarithmic decay: modes with very large angular momentum have a very large timescale, their growth rate can be shown (using for instance the WKB approximation) to decay exponentially with increasing angular quantum number , i.e. τ ≡ Im(ω) −1 ∼ exp(α ), where α is independent of . This suggests that very long timescales can be achieved if the initial data contains support in very large angular momentum quantum numbers, that is to say ∼ log τ . Since the initial data has to live in a Sobolev space of sufficiently high order, we conclude that ||ψ|| ∼ (log τ ) −p , where p is related to the Sobolev norm we are considering. Note that even real analytic data might not decay exponentially, this would correspond to taking the limit p → +∞, which would lead to ||ψ|| ∼ τ −β , for some constant β. The logarithmic behavior has been rigorously shown to be sharp in [14,12]. Perhaps more worrying, the long time behavior of generic perturbations about black hole in global AdS might not even be related with quasinormal modes at all! In [87], it was shown, by counterexample, that quasinormal modes do not form a complete basis and that some perturbations can never be described by their dynamics. However, we should stress that in the presence of a linearly unstable mode, such as a superradiant instability, the linear spectrum of perturbations does provide an accurate description of the dynamics at early times.

A Fluxes across the horizon and asymptotic boundary
In this Appendix we explicitly show that the energy and angular momentum fluxes across the asymptotic boundary vanish if we impose boundary conditions (BCs) that preserve the conformal metric. We also review why the flux across the horizon is proportional to the superradiant factor.
The energy and angular momentum fluxes of gravitational perturbations are calculated using the Landau-Lifshitz "pseudotensor" whose definition we review next (see e.g. [85]). Consider metric perturbations h µν around a backgroundḡ µν up to second order in the amplitude, The linearized Einstein equation reads At second order, the Einstein equation relates terms linear in h (2) to terms quadratic in h (1) : where the RHS is quadratic in h (1) . Written out explicitly, for generic perturbations it reads (here, we use the notation h µν ≡ h (1) µν ; if we choose the traceless-transverse gauge this is known as the Landau-Lifshitz "pseudotensor"): We can now define the fluxes associated with the first order perturbation. Let ξ be one of the Killing vector fields ξ = ∂ t or ξ = −∂ φ of (Kerr-)AdS, that are conjugate to the energy (E) and angular momentum (J) of the solution, respectively. Conservation of the "pseudotensor" T µν , ∇ µ T µν = 0, and the Killing equation, ∇ (µ ξ ν) = 0, imply that the 1-form J µ = −T µν ξ ν is conserved, d J = 0, where is the Hodge dual. We can then define the energy or angular momentum flux across a hypersurface Σ (like the horizon or the asymptotic boundary) as where n ν is the normal vector to Σ and dV Σ is the induced volume on Σ. Consider first the asymptotic boundary Σ = Σ ∞ which is the timelike hypersurface defined by z = 0 (where z is the FG radial coordinate). This has unit normal n = z/L dz. As discussed in association with the FG expansion (1.1), AAdS backgrounds start differing from each other only at order O g (d) z d−2 . This in particular also implies that the most general perturbation of a global AdS background that preserves the asymptotic structure of the background has an asymptotic expansion around z = 0 that starts at order O z d−2 , i.e. h zµ = 0 and h ab = e −iωt e imφ f (X)z d + · · · (the Fourier decomposition in t, φ follows from the fact that these directions are isometries of the background). Inserting this general perturbation in the "pseudotensor" (A.4) and computing the fluxes (A.5) we find that they vanish because the integrand of the fluxes has a polynomial expansion that starts at O z d (for both Killing fields): That is, perturbations that preserve the conformal metric (the static Einstein Universe) have vanishing energy and angular momentum fluxes at the asymptotic boundary. Take now the Killing horizon (null) hypersurface, Σ = Σ H , defined by r = r + . To find the flux across the horizon we work with the ingoing Eddington-Finkelstein coordinates {v, r, χ, φ}, introduced in (2.15), that extend the solution through the horizon. The horizon generator is by definition normal to the horizon, i.e. n ≡ K = ∂ v +Ω h ∂ φ. The metric perturbation h µν ≡ h (1) µν is constructed applying a differential operator to the Teukolsky variable δΨ 4 (this is known as the Hertz map; see the companion paper [30]). This yields long expressions for the components of h µν that are not at all illuminating. The keypoint is that inserting them in the "pseudotensor" (A.4) we find that the fluxes across the horizon are proportional to the superradiant factor (C ξ are positive constants if {ω > 0, m > 0}), 18 That is, these fluxes are negative (inwards the BH) if ω > mΩ h for which we perturbations are damped (QNMs); positive (outwards the BH) if ω < mΩ h in which case this energy and angular momentum fluxes feed the superradiant instability growth; and finally they vanish when ω = mΩ h , i.e. at the onset of superradiance.

B Details of the hydrodynamic QNM computation (D = 4)
In this appendix we give details of the hydrodynamic computation that leads to the frequency quantization (6.15) and (6.16). Our starting point is the double expansion (6.14) in the shear viscosity and in the rotation, both for the fluid perturbations introduced in (6.12), {Q (f ) (X)} = {Q (1) , Q (2) , Q (3) } ≡ {δP/L, δu X , δu Φ }, and for the perturbation frequency ω. These expansions are inserted in the hydrodynamic equations of motion (6.10) or (6.11) that are then solved progressively in a series expansion in η/L 3 and a/L. For our purpose it will be enough to go up to first order in the viscosity (n = 1) and up to second order in the rotation (p = 2) expansions. There are two families of modes, namely the scalar and the vector modes.

B.1 Scalar modes
Consider first the scalar modes. At leading order in the aforementioned expansion, the viscosity and rotational effects are absent, and we are interested in finding the quantities S (f ) 0,0 and ω 0,0 (for scalar modes we use the notation S j,i ). In these conditions, the pressure perturbation is proportional to the Kodama-Ishibashi scalar harmonic S(X, Φ) ∼ e imΦ P m (X), where P m (x) is the associated Legendre polynomial, while the velocity perturbation is proportional to the vector derived scalar harmonics obtained by taking angular derivatives of the scalar harmonic S i ∝ D i S (where D j is the covariant derivative associated to the unit radius metric on S 2 ). We thus have S (1) 0,0 e imΦ = A 1 e imΦ P m (X), S 0,1 e imΦ = A 2 e imΦ P m (X) and S (3) 0,1 e imΦ = i m A 2 e imΦ P m (X), for arbitrary amplitudes A k . Inserting these expressions in the equations of motion (EoM) we fix the ratio A 1 /A 2 and quantize the frequency ω 0,0 . This yields (we introduce the notation and ω 0,0 that can be read from (6.15). This conclusion agrees with the static results first derived in [10]. Still at leading order in the viscosity, we now consider the first order correction introduced by the rotation. It follows from two of the EoM at this order that the perturbations Q 0,1 and/or its derivative. Plugging these relations in the third EoM we fix the frequency correction ω 0,1 as written in (6.15) (this is done doing the procedure exemplified below for the ω 0,2 conribution) and the differential equation for S where B 0 is a new arbitrary amplitude that is introduced at this order. We can improve our approximation by finding the correction up to second order in the rotation (at this point still at vanishing viscosity). This requires looking to the EoM at order O η 0 , a 2 that involve the unknown quantities S (f ) 0,2 and ω 0,2 . We use this case to exemplify in detail how we typically solve equations of our problem to get the perturbative frequency corrections. Two of the EoM at order O η 0 , a 2 yield two algebraically equations for S Note that the ODE for S Namely, we can subtract the vanishing total divergence contribution dX∂ X P m f 1 S (1) 0,2 +ŝ 2 P m + s 1 P m = 0, where we have redefined the coefficients f 2 → f 2 and s 2 →ŝ 2 to absorb the new contributions arising from the integration by parts. We use again a similar approach, namely we subtract the total divergence term dX∂ X P m f 2 S (1) 0,2 = 0 and use integration by parts to get dXP m (s 2 P m + s 1 P m ) = 0 where we made the redefinition s →s 2 and a would be S (1) 0,2 contribution is absent since f 3 −f 2 = 0. Subtracting the total divergence dX∂ X s 2 (P m ) 2 and a third final integration by parts finally yields dXP m ŝ 1 P m = 0 withŝ 1 =s 2 + s 1 . Explicitly, this final condition is − ( + 1) ( + 1) 4 √ 2 ( + 1)ω 0,2 + 2 + − 3 − 16 +m 2 ( + 1) 2 + + 4 + 12 To proceed we use the integrals P m (X)P m (X)dX = 2 (2 + 1) to rewrite (B.4) as This condition finally quantizes the frequency contribution ω 0,2 as To include the effects of dissipation we now consider the linear order contribution in the viscosity, while still doing also an expansion in the (adimensional) rotation parameter, i.e. we solve the perturbative EoM at order O η, a 0 , O η, a 1 , O η, a 2 . The technical analysis proceeds in a way that is very similar to the procedure already outlined for the zero-order contribution in the viscosity so we now omit further details and just give the final results for the frequencies ω 1,i and for the perturbation eigenfunctions S 1,0 = K 2 P m (X) , S 1,0 = i m K 2 P m (X). (B.7) where K 2 is a new arbitrary amplitude, and the frequency ω 1,0 is written in (6.15).
At order O η, a 1 the eigenfunctions are 1,1 = P m (X) where C 0 is a new arbitrary amplitude, and the frequency contribution ω 1,1 can be found in (6.15). Finally, to get the frequency correction at order O η, a 2 we use again the integration by parts procedure that we already described to get the O η 0 , a 2 contribution. Going through this procedure we find the frequency ω 1,2 that can be read from (6.15) and we omit here the associated long expressions for S (f ) 1,2 . The frequency contribution ω 1,3 written in (6.15) is computed in a similar way.

B.2 Vector modes
Consider now the vector modes. These distinguish from the scalar modes because at leading order in the viscosity and rotation they have vanishing pressure perturbation and vanishing frequency: V (0) 0,0 = 0 and ω 0,0 = 0 (for scalar modes we use the notation V 0,2 and its first and second derivatives, in addition to two source contributions proportional to the associated Legendre polynomial and its derivative. It is used to find the frequency contribution ω 0,3 as given in (6.16), after using several integrations by parts as exemplified in the previous scalar mode treatment. The long relations associated to this discussion are omitted here.
We now consider the viscosity contributions. Using the EoM at order O η, a 0 and O η, a 1 we find the eigenfunctions where K 2 is an arbitrary amplitude, and fix the frequency contribution ω 1,0 as written in (6.16) and ω 1,1 = 0. The EoM at order O η, a 1 also determine V 1,1 and V 1,1 that we write below, while EoM at order O η, a 2 1 find algebraic relations for V 1,2 (that we do not write here) and a 19 To be more clarifying, the first order EoM determine ω0,1 = 0 but leave V   1,1 = eigenvalues and quantum numbers are quantized as λ = ( − 1)( + 2) − 2m 2 + + 4 + 1 aω + O a 2ω2 , a 2 L 2 , with = 2, 3, 4, · · · , |m| ≤ , (C.1) where the azimuthal quantum number m is an integer, and we have introduced the quantum number with properties discussed after (2.13). This fixes the angular eigenvalue spectrum and we just need to solve now the radial equation in a regime of parameters that is consistent with the approximation where (C.1) is valid.
We follow a standard matching asymptotic expansion analysis whereby we divide the exterior spacetime of the Kerr-AdS black hole into two regions; a near-region where r − r + 1 ω and a far-region where r−r + r + . In each of these regions, some of the terms in radial equation make a sub-dominant contribution and can be consistently discarded. We will find that if we further require r + L 1, this procedure yields an equation with an (approximate) analytical solution in both spacetime regions. The next important step is to restrict our attention to the regime r +ω 1. In this regime, the far and near regions have an overlaping zone, r + r − r + 1 ω , where the far and near region solutions are simultaneously valid. In this matching region, we can then match/relate the set of independent parameters that are generated in each of the two regions. We will also find that if we further restrict our analysis to the regime a r + 1, it is sufficient to work only with the the leading order contribution for the angular eigenvalues in (C.1), λ ∼ ( − 1)( + 2).
The regime of validity of the matching analysis can be expressed in a much simplified form.
Indeed, the rotation parameter is constrained by the extremity condition a ≤ r + 3r 2 + +L 2 L 2 −r 2 + for r + < √ 3L and by a < L for r + > √ 3L (see e.g. [45]). For the regime we are interested, r + L 1, we thus have a ≤ r + 3r 2 + +L 2 L 2 −r 2 + = r + + O(r 3 + ). Thus r + L 1 automatically implies a L 1. Moreover, for r + L 1 the (real part) of the QNM frequencies of the BH do not differ much from the normal mode frequencies of global AdS that are orderωL ∼ O(1). Therefore r + L 1 also implies r +ω 1 and aω 1. To sum, our analysis will be valid in the regime of parameters (4.1). We discuss the different regions separately and discuss how to match the solutions obtained next.

C.1 Near-region equation and regular solution at the horizon
The near-region is defined by r − r + 1 ω . Introducing the wave function Φ = ∆ r R (−2) ωm , the radial equation (2.12) reads If we further restrict to the regime r + L 1, the cosmological constant contribution can be neglected. Specifically, in the radial equation (C.2) the following approximations are valid ∆ r r∼r + r 2 + a 2 − 2M r + · · · (r − r + )(r − r − ) , with r − a 2 r + , 6r 2 L 2 + 4iK r r∼r + 6r 2 + L 2 − 8ir +ω 1 − a 2 L 2 ∼ −8ir +ω + · · · , (C.3) K 2 r − 2i K r ∆ r ∆ r r∼r + Ξ 2 r 2 + + a 2 2 (4πT H ) 2 ( + 2i) ∆ r + 8ir +ω + · · · (r + − r − ) 2 ( + 2i) (r − r + )(r − r − ) + 8ir +ω + · · · , where Ω H , T H are the angular velocity and temperature defined in (2.4), and motivated by the BC (2.16) we have introduced the superradiant factor, With these near-region approximations the radial equation (C.2) is then where, in the approximation regime (4.1), we replaced the eigenvalue λ by its leading contribution in (C.1) (the requirement a/r + 1 is fundamental here since we neglect a contribution proportional to msa/r + when compared with λ ∼ ( − 1)( + 2)). Introducing a new radial coordinate z and wavefunction F defined as This is a standard hypergeometric equation [95], z(1−z)∂ 2 z F + [c − (a + b + 1)z]∂ z F − abF = 0, whose most general solution in the neighborhood of z = 0 is A in z 1−c F (a − c + 1, b − c + 1, 2 − c, z) + A out F (a, b, c, z). Using (C.6), one finds that the most general solution of the near-region equation is therefore F (a, b, c, z) , The first term represents an ingoing wave at the horizon z = 0, while the second term in (C.8) represents an outgoing wave which we set to zero, A out = 0, to guarantee that no perturbations come off the horizon. For the matching we need the large r (i.e z → 1) behavior of the ingoing near-region solution. To get this, we use the z → 1 − z transformation law for the hypergeometric function [95],

C.2 Far-region wave equation and global AdS solution
The far-region is defined by r−r + r + . In this region the effects induced by the black hole mass and angular momentum can be neglected to leading approximation. The far-region background where the gravitational perturbation propagate is then simply global AdS spacetime. Our approximations then yield ∆ r r 2 1 + r 2 L 2 and, in the regime where the eigenvalue λ is given by the leading contribution in (C.1), the radial equation (2.12) boils down to ∂ r ∆ r ∂ r R B 0 F − 1, + 1 + Lω; 2( + 1); 2r r + iL where B 0 , B 1 are at this point arbitrary amplitudes whose ratio will be constrained by the asymptotic global AdS BC. 21 Asymptotically the solution decays as In the regime we are working one has a 0 and the BC expressions (2.20)-(2.22) for β simplify considerably reducing to 1) β = β s = −Lω 1 + λ λ − 2 (L 2ω2 − 1) , (C.14) 2) β = β v = λ 2Lω − Lω , (C. 15) for scalar and vector modes, respectively. Here, λ = ( −1)( +2). Going through this asymptotic matching we find how the amplitudes B 0 B

C.3 Matching. QNM and superradiant frequencies
In the regime r +ω 1, the near and far regions have an overlaping zone, r + r − r + 1 ω , where both are simultaneously valid. The requirement that the solutions can be matched across the overlapping zones related the amplitudes A in , B (−2) + and quantizes the frequencyω. In particular, the frequencies that are allowed to propagate in the Kerr-AdS black hole are found matching the large r behavior (C.10) of the near-region solution with the small r behaviour (C.16) of the far-region solution. This yields two conditions, one following from the matching of the r +2 coefficients and the other from the matching of the r 1− coefficients. One of these constraints is used to find the ratio between the near and far region amplitudes where the superradiant factor was introduced in (2.17), and the asymptotic BC parameter β is given by (C.14) for scalar, and by (C.15) for vector perturbations. Recall that this expression is valid in the approximation regime (4.1). This frequency quantization condition simplifies considerably when we choose a particular harmonic . In particular, for the lowest harmonic, = 2, it reduces to (4.2). We leave the detailed discussion of the solution of this frequency quantization condition and the comparison with the associated exact numerical results to subsection 4.1.