Kerr-de Sitter Quasinormal Modes via Accessory Parameter Expansion

Quasinormal modes are characteristic oscillatory modes that control the relaxation of a perturbed physical system back to its equilibrium state. In this work, we calculate QNM frequencies and angular eigenvalues of Kerr--de Sitter black holes using a novel method based on conformal field theory. The spin-field perturbation equations of this background spacetime essentially reduce to two Heun's equations, one for the radial part and one for the angular part. We use the accessory parameter expansion of Heun's equation, obtained via the isomonodromic $\tau$-function, in order to find analytic expansions for the QNM frequencies and angular eigenvalues. The expansion for the frequencies is given as a double series in the rotation parameter $a$ and the extremality parameter $\epsilon=(r_{C}-r_{+})/L$, where $L$ is the de Sitter radius and $r_{C}$ and $r_{+}$ are the radii of, respectively, the cosmological and event horizons. Specifically, we give the frequency expansion up to order $\epsilon^2$ for general $a$, and up to order $\epsilon^{3}$ with the coefficients expanded up to $(a/L)^{3}$. Similarly, the expansion for the angular eigenvalues is given as a series up to $(a\omega)^{3}$ with coefficients expanded for small $a/L$. We verify the new expansion for the frequencies via a numerical analysis and that the expansion for the angular eigenvalues agrees with results in the literature.


Introduction
Quasinormal Modes (QNMs) are free oscillatory modes of an open physical system. Contrary to normal modes, QNMs oscillate but also decay exponentially in time. The exponential ringdown of these modes is thus characterized by complex frequencies ω QN M = ω R + i ω I , where ω R ∈ R corresponds to the frequency of the oscillation and ω I < 0 sets the timescale of the decay. The study of QNMs and their exponentially-growing counterparts (ω I > 0) is important for, e.g., the linear stability analysis of background spacetimes (e.g., [1]), to characterize the ringdown stage of a waveform from a black hole binary inspiral [2] and for setting the thermalization scale in holographic gauge theories [3]. For reviews on QNMs, see, e.g., [4,5].
Black holes are compact objects characterized by the existence of an event horizon. Among the class of stationary and axisymmetric solutions of General Relativity, we can find exact black hole spacetimes. They are classified by their mass M , angular momentum J, charge Q and by the cosmological constant Λ. In the following, we use the more convenient angular momentum parameter per unit mass, a ≡ J/M , and the de Sitter radius L ≡ 3/Λ (Λ > 0).
In this paper, we obtain QNMs and angular eigenvalues for the spin-s perturbation equations of four-dimensional Kerr-de Sitter black holes (Λ > 0). These black holes contain four horizons and we focus on the QNMs obtained from solutions which are valid between the outer event horizon (at radius r + ) and the cosmological horizon (at radius r C ≥ r + ).
A standard approach for studying massless spin-field perturbations of black hole spacetimes is to write the spin-s field equations using the Newman-Penrose formalism. This formalism enables the decoupling and simplification of the linear field perturbation equations into one single ("master") partial differential equation for any Petrov type D spacetime. In the absence of sources, this master equation is separable into two ordinary differential equations (ODEs) for modes with fixed frequency ω and azimuthal angular number m. This approach was first pursued by Teukolsky for a Kerr black hole [6]. More generally, one can find a master equation for spin s = 0, ±1/2, ±1, ±3/2, ±2 massless field perturbations of any Petrov type D spacetime [7][8][9][10], including Kerr-de Sitter black holes. The associated ODEs are all reducible to Heun equations [10,11]. In this work, we exploit an integrable structure behind Heun equations, deeply related to 2D conformal field theory (CFT), to find QNMs and angular eigenvalues of Kerr-de Sitter black holes.
Heun equations are ODEs with four regular singular points [12]. We can fix these points to be at the values z = 0, 1, x, ∞ of the independent variable, without loss of generality. These equations can then be characterized by the moduli parameter x, its local monodromy data and the so-called accessory parameter. The monodromy data governs the behaviour of the analytically continued solutions around the singular points, while the accessory parameter contains global information of the solutions. Heun equations can be formally extended to have an extra apparent singularity, with trivial monodromy. The position of the apparent singularity can then be continuously changed in a precise way to leave unaffected the monodromy data. This is called the isomonodromic deformation of the Heun equation [13,14]. Perhaps surprisingly, this deformation is Hamiltonian and integrable. One can use this structure to construct an expansion of the accessory parameter in x with coefficients given in terms of the local and global monodromy data [15,16]. We call this the Accessory Parameter Expansion (APE), revised in section 4.1. This approach was successfully applied in the context of 2D conformal field theory in relation to the classical limit of conformal blocks [16] and for 2D conformal mappings from a simply connected planar region to the interior of a circular arc quadrilateral [17]. For other approaches to obtaining the APE, we refer the reader to [15,18,19].
In this work, we use the APE to calculate Kerr-de Sitter QNMs and angular eigenvalues for the first time. In particular, we apply the APE to the case of Kerr-de Sitter black holes near the extremal limit r C → r + , also called the rotating Nariai limit, about both the upper and lower superradiant bounds. 1 More specifically, we obtain an analytical expansion for the QNM frequencies of the form: s ω m = mΩ + + ω 0 (a) +ω 1 (a) 1 +ω 2 (a) 2 + O( 3 ) , (1.1) in the extremality parameter ≡ (r C −r + )/L, whereω 0,1,2 are coefficients we determine.
Here and in what follows Ω k , k = +, C denotes the angular velocity of the horizon r k , k = +, C. We provide this expansion: (i) up to third order in with the coefficients expanded up to third order in a/L; (ii) up to second order in for general a. The expansion in Eq. (1.1) is about mΩ + (unexpanded in ). The expansion about mΩ C can be obtained by letting Ω + → Ω C , → − and taking the complex conjugate of the right hand side of Eq. (1.1). We thus show, and check numerically, that the lower superradiant bound ω = mΩ C is an accumulation point of QNMs, similarly to the upper bound ω = mΩ + .
To the best of our knowledge, only the termω 0 in Eq. (1.1) for the specific case of s = 0 and coupling constant with the Ricci scalar ξ = 0 had so far been obtained in the literature. It was obtained in [20] by calculating QNMs directly on the rotating Nariai geometry, which is found by a near horizon limit of the Kerr-de Sitter spacetime. Apart from the APE derivation of Eq. (1.1), we also calculate QNMs directly in the rotating Nariai geometry in the case of s = 0 and general ξ, thus obtaining an alternative derivation ofω 0 for this case. Overall, we believe that the expressions we derive forω 0 for non-zero spin (as well as for s = 0 and ξ = 0) and forω 1,2 for any spin are new. We have checked our high-order general-spin Eq. (1.1) against a numerical analysis that we have carried out using Leaver's method to high accuracy. This method was first used in [21] to obtain Kerr-de Sitter QNMs and angular eigenvalues and our results are consistent with that work.
Kerr-de Sitter angular eigenvalues were first obtained numerically in [9]. Our expansion for the eigenvalues for small α and aω agrees with that already obtained in [11] using a completely different method (and with that in [22] in the Kerr limit L → ∞). Our eigenvalue derivation here is a proof of principle of the APE, which is extendable to higher orders than what we present here.
Summing up, with our method, we obtain: (i) an angular eigenvalue expansion for small aω and a/L consistent with [11,22], and (ii) a new expansion for the QNM frequencies near the rotating Nariai limit. For the reader who is only interested in the final result, the coefficients in the QNM expansion (1.1) are given in Eqs. (4.15) and (4.18) and in appendix D, and the coefficients in the eigenvalue expansion (1.3) are given in appendix E.
There exist various other methods for calculating QNM frequencies, from asymptotic analyses to numerical techniques. Ref. [21] provides a numerical calculation of QNMs of Kerr-de Sitter black holes, for modes between r C and r + , by using a method originally developed by Leaver [23,24]. Their analysis suggests that the QNM frequencies accumulate towards ω = mΩ + in each of the extremal limits r + → r − and r + → r C , where Ω + is the angular velocity of the event horizon r + . The frequency ω = mΩ + is the upper bound of the superradiant regime. For non-extremal (i.e., r + = r C ) Kerr-de Sitter black holes, the superradiant regime is |m|Ω C < |ω|< |m|Ω + , where Ω C is the angular velocity of the cosmological horizon r C . Other works we found in the literature on QNMs in Kerr-de Sitter are the following: [25] provides an analytic expansion of the QNM with the lowest overtone (N = 0) in the double approximation of small a (to linear order) and large-( + 1/2) (up to order ( + 1/2) −6 ), where is the multipolar number; similarly, [26] provides a Bohr-Sommerfeld type of condition which QNMs for small a must satisfy and which could be solved numerically, being better suited for larger ; within the context of the so-called strong cosmic censorship conjecture, [27] obtained an analytic expression for the QNMs for = |m| 1 and complement it with a numerical study.
We finish the introduction by outlining the contents of this paper. We describe the geometry and configuration space of the Kerr-de Sitter black hole in section 2. In section 3, we review the spin-field master equation in Kerr-de Sitter and the scattering problem. We express the product of scattering coefficients, the so-called greybody factor, in terms of the monodromies of the radial ODE. We then describe the monodromy conditions for angular eigenvalues and QNMs. We finish section 3 by obtaining QNMs directly in the rotating Nariai geometry [20] and obtaining the same result from the full Heun equation in Kerr-de Sitter under the appropriate limits. In section 4, we review how to obtain the APE via the isomonodromic τ -function, first proposed in [28]. In Subsecs. 4.2 and 4.3, we apply the APE around = 0 to obtain expansions for the angular eigenvalues and QNMs. In section 5, we compare our analytic results with a numerical calculation using Leaver's method. In section 6 we present our conclusions and potential future developments.
We also include complementary sections with more details on our results. In appendix A, we introduce the gauge transformations of Fuchsian equations used to reduce the master equations to Heun equations. In appendix B, we present an alternative derivation of the leading order QNM equation. In appendix C, we review how to obtain the APE from the isomonodromic τ -function as in [16], with the crucial modification of expanding the monodromies in . Finally, in appendices D and E we display the coefficients of the QNM expansion and angular eigenvalues respectively. In this paper, we choose Planck unitsh = c = G = 1.

Kerr-de Sitter Black Holes
The Kerr-de Sitter metric [29] can be written in Chambers-Moss coordinates [9] as where and α ≡ a/L. The coordinate ranges are t, r ∈ R, θ ∈ [0, π], φ ∈ [0, 2π). The solutions are labeled by the parameters (a, M, L), where M is the mass, J = M a is the angular momentum and L ≡ 3/Λ is the de Sitter radius, where Λ is the cosmological constant. In particular, the Ricci scalar is equal to R = 12/L 2 . When L → ∞, we recover the Kerr metric in Boyer-Lindquist coordinates and when a = 0 we recover the Schwarzschild-de Sitter metric. The polynomial ∆ r has four roots: r i , i ∈ I, where I = {C, −, +, −−}, with r −− = −(r − +r + +r C ). It is easy to check that, in order to have a black hole in a de Sitter Universe, i.e, for the two largest roots to be positive, 0 < r + ≤ r C , then the other two roots must be real and must satisfy r −− < 0 ≤ r − ≤ r + < r C , and thus r − is the inner Cauchy horizon. This implies that r C , r + and r − are the radii of, respectively, the cosmological horizon, the (outer) event horizon and the (inner) Cauchy horizon. The negative horizon r −− is reachable if one goes inside the inner horizon and, avoiding the ring-like singularity (located at r = 0, θ = π/2), goes to negative values of r, thus reaching another asymptotically de Sitter region. The global extension of the Kerr-de Sitter black hole is described in [9,[29][30][31] and is shown in figure 1.
The limit r + → r − corresponds to an extremal black hole (r + = r − is an event horizon), whereas r + → r C corresponds to an extreme naked singularity (r + = r C is a cosmological horizon) [31] and it is called the rotating Nariai limit [20,32]. 2 The extremal black hole is a maximally-rotating black hole, and it is interesting to note that it can have a > M (a = M being the maximal rotating limit in Kerr). . Kerr-de Sitter configuration space in terms of ( , a) with L = 1. The green region corresponds to the same one in figure 2 of Kerr-de Sitter black holes. The = 0 line corresponds to the rotating Nariai limit r + = r C and the point = 1, a = 0 corresponds to empty de Sitter spacetime. The blue line corresponds to the extremal limit r + = r − , while the orange line corresponds to the equal temperature condition T C = T + for r + = r C (lukewarm solutions).
Two or more roots of ∆ r coincide if and only if the discriminant D = L −12 i =j∈I (r i −r j ) is zero. It is possible to write the discriminant in terms of only M/L and a/L. One can show that D ≥ 0 and M, a ≥ 0 if and only if all the roots of ∆ r are real and r C ≥ r + ≥ r − ≥ 0 > r −− . This condition thus defines the region in blue in figure 2, corresponding to the parameter space where black hole solutions exist [9].
The QNM expansion that we calculate later in this paper is most easily organized as a series in a and in the extremality parameter ≡ (r C − r + )/L, instead of a and M . The horizons are given by Each horizon r k , k ∈ I, has an associated temperature T k and an angular velocity Ω k given by , Ω k = a r 2 k + a 2 , where the prime denotes derivative with respect to the argument of ∆ r (r). The configuration space in terms of ( , a) is given in figure 3. The extremal limit corresponds to the blue line, while the rotating Nariai limit is clearly the = 0 line. Given the temperatures of the event and cosmological horizons, respectively, T + and T C , the orange line in figure 3 corresponds to the condition T + = T C for r + = r C , called lukewarm solutions in [20]. As noted in [20], these solutions are not in thermal equilibrium, as the two horizons can still exchange angular momentum while keeping the temperature fixed. Finally, it is worth mentioning that, recently, it was shown that the region between the event and cosmological horizons of Kerr-de Sitter is non-linearly stable for small angular momentum [33].
3 Scattering on Kerr-de Sitter Black Holes

Linear Field Perturbations
The analysis of linear field perturbations is important in order to address, for example, the linear stability of the spacetime and wave scattering around black holes. The standard formalism describing spin-field perturbations of the Kerr black hole was set up by Teukolsky [6]. Teukolsky used the so-called Newman-Penrose formalism, whereby the various field components are projected onto a null tetrad. Teukolsky decoupled the equations for the spin-field perturbations in Kerr, obtaining one single master equation which is obeyed by all possible spins of the field. Furthermore, in Boyer-Lindquist coordinates, the Teukolsky master equation essentially separates into two ODEs, one for the radial coordinate and the other one for the polar angular coordinate. This approach was later generalized to a Kerr-de Sitter background spacetime in [8,9,11]. The two ODEs are actually of Heun type [12] for any Petrov type D spacetime with a cosmological constant [10]. We note that there are different tetrad normalization choices in the literature, but the resulting equations are easily obtained by simple transformations that do not change the linear analysis. Here we follow the conventions of [11,21], which are different from those in [7,9]. We next review the radial and angular equations for linear spin-field perturbations of Kerr-de Sitter.
Let ψ s,ω,m (t, φ, r, θ) = e −iωt e imφ R s,ω,m (r)S s,ω,m (θ) be a mode solution of the master equation for the Kerr-de Sitter metric in Chambers-Moss coordinates (2.1). 3 Here, ω ∈ C is the mode frequency and m ∈ Z its azimuthal number. The spin-s master equation separates into the following angular and radial ODEs [8,9,11]: and with s λ m = s λ m (aω, α) being the angular eigenvalue. 4 The prime in ∆ r denotes derivative with respect to r and in ∆ u with respect to u. Notice that the s = 0 radial equation above corresponds to the Klein-Gordon equation with conformal coupling [21], i.e., with Ricci scalar coupling parameter ξ = 1/6. The solutions of the angular equation (3.1a), assuming regularity at u = −1 and +1, are usually called (Anti-)de Sitter spheroidal harmonics [7] and its eigenvalues s λ m = s λ m (aω, α) are only known numerically [21] or as an expansion to order O(α 2 , (aω) 2 ), as derived in [11].

Reduction of Radial and Angular Equations to Heun Equations
In the following, we want to analyze the connection problem for the angular and radial equations (3.1). By the connection problem we mean rewriting the solution near a certain singular point of the differential equation in terms of a linear combination of solutions near another singular point. The tasks of finding the eigenvalues of the angular equation and the QNMs of the radial equation can both be described in terms of the connection problem [28,34,35], as we show in the next subsection. It is thus helpful to reduce the radial and angular equations to some known special form, so that we can understand the structure of the solutions. In this section, we are going to apply Möbius transformations to the radial and angular variables r and u. This will map the roots of the polynomials ∆ r (r) and ∆ u (u) to (z k ) k∈J ≡ (0, 1, x, ∞), with J ≡ {0, 1, x, ∞}.
Equations (3.1a) and (3.1b) naively possess five singular points each. In the radial case, these points are the four roots of ∆ r (r), i.e., r k , k ∈ I = {C, −, +, −−}, together with the point r = ∞. In the angular case, the singular points are the four roots of ∆ u (u), i.e. u = ±1, ± i α , together with the point u = ∞. For Λ = 0, all singular points of the radial and angular equations are regular and thus both equations are Fuchsian. Surprisingly, as firstly shown in [11], both equations can be reduced to Heun equations [12], i.e., Fuchsian equations with four regular singular points. After a series of transformations, the points r = ∞ and u = ∞ cease to be singular points of the corresponding ODEs, thus being characterized as removable singularities [36]. We rederive this result in a systematic way in appendix A. On the other hand, for Λ = 0 there is a confluence of singular points, producing irregular singular points at r = ∞ and at u = ∞. Both radial and angular equations then become confluent Heun equations [34,37]. In the following, we derive the Heun equations from the angular (3.1a) and radial (3.1b) equations.

Generic Heun Equation.
The basis of our method is to classify the Heun equation in terms of its monodromy data. For a function y = y(z), the Heun equation can be written in the canonical form (see appendix A or [12,36]) where the singular points are (z k ) k∈J ≡ (0, 1, x, ∞), with x ∈ C. The coefficients θ ≡ (θ 0 , θ 1 , θ x , θ ∞ ), also known as the local monodromy coefficients (to which we shall simply refer as monodromies), correspond to half of the difference of Frobenius exponents of Frobenius series solutions y (j) k , j = 1, 2 close to the singular points. That is, In their turn, b 1 and b 2 are the Frobenius exponents at z = ∞: The regularity of z = ∞ in (3.5) implies the Fuchs condition The two relations in Eqs. (3.8) and (3.9) can be used to write the Frobenius exponents b 1,2 in terms of the θ's as The accessory parameter K x is the only parameter which, for generic monodromy data, cannot be written only in terms of θ's, as it also depends on global information of the solutions. Therefore, the Heun equation is completely determined by five parameters for fixed moduli parameter x: M Heun ≡ (θ 0 , θ 1 , θ x , θ ∞ ; K x ).

Angular Equation.
The Frobenius exponents of local solutions of (3.1a) around (u k ) k∈J ≡ −1, i α , 1, − i α are given by Explicitly, the difference of Frobenius exponents is given by Starting with Eq. (3.1a), we make the coordinate transformation and We also make the s-homotopic transformation [36] S s,ω, which, as proven in appendix A.2, reduces (3.1a) to the Heun equation The accessory parameter is , (3.18) where we used the fact that (3.19) as the sum of residues k H k = k res u=u k (H(u)/∆ u (u)) is zero from the Cauchy theorem. It should be clear that the extra property (3.19) is not generally true for an arbitrary Heun equation; it is a special condition for the black hole case and reduces the number of parameters from five to four, as one of the monodromies is determined in terms of the others. The x → 0 (α = a/L → 0) limit of (3.17) is regular for L fixed, as neither the θ k 's nor "x(x − 1)K x " diverge, whereas the coefficients diverge when taking the same limit but with a fixed, instead of L. The first limit is the Schwarzschild-de Sitter limit, whereas the second one corresponds to the Kerr limit. In the limit a → 0, we have s λ m → ( + 1) − s 2 .

Radial Equation.
The Frobenius exponents of the radial solutions are θ k = iW k + s/2, where W k ≡ res r=r k W (r) ∆r(r) , k ∈ J , are residues (see appendix A.1). Explicitly, The plus or minus sign in (3.20) arises from the sign of ∆ r (r k ), chosen in such a way that the temperatures T k are positive [9,38]. From the residue theorem, we have that k∈J θ k = 2s. (3.21) which implies, from Eq. (3.10), that b 2 = 1 − 2s and b 1 = 1 − 2s + θ ∞ . In order to obtain the Heun equation from Eq. (3.1b), we apply the coordinate transformation corresponding to the mapping and the s-homotopic transformation The transformations (3.22) and (3.24) result in the Heun equation According to (A.40), the accessory parameter in (3.25) is given by Notice that we use some of the same symbols in the radial and in the angular equations, but it should be clear from the context whether such symbol corresponds to the radial or to the angular case.
The limits x → 0 and x → 1 correspond to the two possible extremal limits: r + → r C and r + → r − , respectively. If we do not make any assumptions on ω, two of the θ k 's diverge in each of these limits and so the resulting confluence yields an irregular singular point (confluent Heun equation). However, if we expand ω = mΩ + + xω + O(x 2 ) (or starting the series with mΩ C ), with some coefficientω, and take the limit x → 0 in Eq. (3.25), then the θ k 's are all finite and we obtain a regular singularity from the coalescence of singularities. We can also expand ω = mΩ + + (1 − x)ω + O((x − 1) 2 ) (or starting with mΩ − ), take the limit x → 1 to reach the same conclusion. 5 As we show below, this procedure for x → 0 gives the near-horizon, near-extremal limit associated to the rotating Nariai limit.

Kerr-de Sitter Connection Coefficients and Greybody Factor
A standard approach for analytically calculating (approximations to) scattering coefficients is to match series expansions near each singular point of the radial ODE in an intermediate region. In many cases, the radial ODE can be solved exactly in the asymptotic regions in terms of known special functions (such as hypergeometric functions) -see, e.g., [39][40][41] for QNMs in near-extremal Kerr. However, when trying that approach for the Kerr-(Anti-)de Sitter family of solutions, one obtains a Heun or confluent Heun equation in the simplest of the cases. Therefore, in order to reduce the ODEs in these cases to the hypergeometric equation, one needs to make some extra assumption, such as low frequency [42] or large (Anti)-de Sitter radius [43]. Another approach for calculating scattering coefficients in the case of Heun equations is the so-called Mano-Suzuki-Takasugi (MST) method [12,[44][45][46]. This method consists of expressing the solutions as infinite series of hypergeometric functions, and it was used in [11,47,48] in Kerr-de Sitter. Given a certain augmented convergence condition, these series solutions converge in an ellipsis with two singular points as foci.
Here we take a different path and use the monodromy approach to scattering amplitudes, first discussed in [34] and later generalized in [28,35,37]. An alternative monodromy approach for calculating QNMs with large imaginary part has also been developed in, e.g., [49,50].
Ingoing and Outgoing Solutions. Let us first focus on the radial equation (3.1b) supposing that we know the angular eigenvalue s λ m . For definiteness, we consider the problem of connecting a local solution near the black hole horizon r + to a local solution near the cosmological horizon r C . Let us define two linearly independent solutions R (j) ± (r) by their local behaviour near r = r j , for j = +, C, as and κ j is the surface gravity (see Eq. (2.4) for the angular velocity and temperature definitions and Eq. (3.20) for the θ's).
The notion of ingoing and outgoing solutions is essential to define a scattering problem. Here these notions should match the ingoing and outgoing null coordinates in Kerr-de Sitter, (v, u) respectively: which in turn motivate the definition of the radial tortoise coordinate r * as This implies that r * → +∞ as r → r C and r * → −∞ as r → r + . Partly in terms of r * , the above local solutions may be written as Therefore, travelling waves have the form exp(−i(ωt ±ω j r * )) and the phase velocity depends on the relative sign of the real part of the frequencies, i.e., − is outgoing at r = r j . (3.32) In our case, as r C ≥ r + > 0, we have that Ω + ≥ Ω C . This allows for the possibility |m|Ω + > | (ω)|> |m|Ω C , which for ω ∈ R corresponds to the superradiant regime, valid only for bosons [48,51]. The physical definition of superradiance is that the net flux of radiation with respect to each horizon should have alternate signs. Given that the notion of ingoing and outgoing is equivalent to the sign of the radiation flux, for ω ∈ R, we havẽ For further details, we refer the reader to [48]. Notice that definitions (3.32) and (3.33) are invariant under (ω, m) → (−ω, −m).

Connection Coefficients for Radial Solutions.
We now wish to obtain formulas relating the local solutions R (j) ± , j = +, C. This means writing a set of local solutions at a singular point as a linear combination of a set of local solutions at another singular point in terms of certain connection coefficients: transmission coefficients T s and reflection coefficients R s . These coefficients become scattering coefficients if the frequency is real, in which case we can consistently define the radiation flux at each one of the horizons. We classify the scattering problem between r + and r C into two possibilities with respect to the boundary conditions: purely ingoing or purely outgoing solutions at one of the horizons.
We define the IN mode by imposing the boundary conditions [40], where T s = T s (ω, m) and R s = R s (ω, m) are complex connection coefficients. If (ω) (ω j ) > 0, this corresponds to a purely ingoing solution at r + . We can obtain another linearly independent solution from (3.34) via the discrete transformation and multiplying the solution by ∆ −s r , i.e., If ω ∈ R, then (ω, m) → (−ω, −m) is equivalent to taking the complex conjugate solution, and so R This implies that R s = R * −s and T s = T * −s when ω ∈ R. The "Wronskian" between two solutions is It is a simple fact that this current is constant in r, as long as u 1 and u 2 obey the same radial equation (3.1b). Evaluating the quantities at r = r + , we obtain the following expression: 39) whose real part represents the flux of radiation transmitted into the black hole when ω ∈ R.
Evaluating the quantities at r = r C , yields for the total flux, Because the flux is conserved, W hor and thus .
The greybody factor G s (ω, m) is defined to be the ratio between the horizon flux and the incoming flux For ω ∈ R and outside the superradiant regime, the greybody factor must be real and positive. This has been proven in [48] using the Teukolsky-Starobinsky identities, which relate T −s with T s . This procedure also makes the superradiance regime explicit. We refer the reader to [48] for details. We can do the same exercise for the UP mode, defined via the following boundary condition at the cosmological horizon: 44) and the reflected mode, which is defined via the boundary condition: We then obtain similar results for the flux as above for the purely in/outgoing solutions on the event horizon: Now, writing τ s = |τ s |e iφ , φ = arg(τ s ), we define the vectors from the above linearly independent solutions. Then it is straightforward to show that [34,40] with the connection matrix given by The condition of unit determinant for M +C is equivalent to the current conservation (3.42) R sR s +T sT s = 1. (3.53)

Monodromy Approach to Angular Eigenvalues and Quasinormal Modes
Let us start by reviewing how we can write greybody factors in terms of monodromies [28,34,35,37]. They informed reader eager to skip this technical discussion can find the main results for the Kerr-de Sitter greybody factor in (3.61), the quasinormal mode condition (3.63) and the angular eigenvalue condition (3.71).
A general solution of a second order linear ODE can be written as a linear combination of local solutions, such as those in (3.27) for the radial equation. Thus, we can define a vector of two general solutions by with g k ∈ SL(2, C) constant matrices and Φ k denotes a vector of local Frobenius solutions, with k = 1, . . . , n referring to any of the n singular points r = r k of the ODE. If we go around a positively-oriented loop γ k in the complex-r plane enclosing only one of the singularities (r k ) k=1,...,n (see figure 4), then Φ transforms to  where the initial value Φ of the solution at the loop is multiplied by a constant monodromy matrix M k ∈ SL(2, C).
In principle, there should be an overall phase e iπs in (3.56), obtained from the common factor (r − r k ) −s/2 in the definition (3.27), but we absorb it, as a common overall factor, in the solution Φ to obtain SL(2, C) representations of the monodromy group. Notice that Given two singular points r k and r j with j, k = 1, . . . , n and j = k, we wish to study the connection matrix where the last equality follows from Eq. (3.54). Thus, the connection matrix is constant and determined by the g k , k = 1, . . . , n. Choosing a proper normalization for the local solutions, in terms of the g k , such that det M jk = 1, we can write this connection matrix in the form (3.51).
Let us now choose a path γ jk which encloses two singular points r j and r k , as in figure 4, in such a way that the general solution changes to Φ γ jk (z) = Φ(z)M j M k . We define the composite monodromy coefficient σ jk by The reason for the extra minus sign here, in comparison with the definition of the single traces, is explained in appendix B. Applying Eqs. (3.56), (3.58) and (3.59) to the case of (3.51), we find .
This is the greybody factor for a connection matrix normalized as per Eq. (3.51).
Quasinormal Modes. For the radial equation (3.1b), according to (3.27) and (3.49), we can rewrite (3.60) as , (3.61) where σ ≡ σ +C . We wish to calculate QNM frequencies in Kerr-de Sitter. These frequencies are defined by the condition that the modes are purely ingoing at the black hole horizon and purely outgoing at the cosmological horizon. 6 According to Eq. (3.34), the QNM condition corresponds to (T s ) −1 = 0, i.e., poles of T s , which appears in Eq. (3.61). The poles ω = ω p on the right hand side of Eq. (3.61) are However, these poles may correspond to poles of τ s or T s = T −s (−ω, −m) as well. Not being able to distinguish between poles of T s and those of T s or τ s is a shortcoming of the monodromy formula above, as it only gives information about the product of these coefficients. However, below we recover the exact greybody factor in the rotating Nariai limit r C → r + by using Eq. (3.62). This will also allow us to identify the poles of T s (ω, m) with N ∈ Z + , and the poles of T −s (−ω, −m) with N ∈ Z − . While the former corresponds to the actual QNMs, the later corresponds to the reflected QNMs, and thus, in principle, they are not QNMs. In order to assess the possible poles coming from τ s , one needs to generalize the Teukolsky-Starobinsky identities obtained in Ref. [47] for complex ω. We leave this for future work.
In section 4, we are going to apply the APE derived in Ref. [16] to find QNMs. Because the APE has the symmetry σ → −σ, which is also compatible with Eq. (3.60), changing the sign in front of (θ C − θ + ) in Eq. (3.62) is equivalent to letting N → −N − 1. Therefore, from the APE point of view, there is no distinction between the positive and negative N possible poles in (3.62) and, from here on, we choose the QNM condition to be As mentioned above, the condition N ≥ 0 will be justified in the rotating Nariai limit in the next subsection. Finally, it is worth mentioning that, in this paper, we focus on the radial equation for which 2θ k / ∈ Z, and so we do not cater for any potential QNM frequencies that could correspond to 2θ k ∈ Z. This can introduce logarithm singularities in the solutions and thus change the QNM condition above.
Accumulation points of Quasinormal Modes. As a function of complex frequency, the transmission coefficient may have not only poles (QNMs) but also branch cuts. Apart from spurious (unphysical) branch cuts arising from the angular eigenvalue [52,53], the transmission coefficient may also have cuts arising from the behaviour of the radial potential which may have physical consequences. Let us first write the radial equation (3.1b) as ∂ 2 r * + V(r) R(r) = 0, where R(r) ≡ ∆ s/2 r (r 2 + a 2 ) 1/2 R s,ω,m (r) and the radial potential V(r) is straightforwardly read-off from Eq. (3.1b). It is easy to check that, if r C = r + , then V(r) → V +,C + O(e α +,C r * ), as r → r +,C , for some constants V +,C and α + > 0 and α C < 0. Based on the arguments of [54] (see also [55]), this asymptotic behaviour does not lead to branch points in the complex ω-plane for the radial solutions which obey purely-outgoing (or purely-ingoing) boundary conditions at either r + or r C . On the other hand, in the case r + = r C , this singular point becomes an irregular singular point and it is easy to show that, then, V(r) →V C + O(1/r * ) as r → r + = r C , for some constantV C . This is essentially due to the fact that r * = O(1/(r − r + )) as r → r + when r + = r C , as opposed to r * = O (ln(r − r + )) as r → r +,C when r + = r C . Furthermore, V C possesses a zero at the superradiant bound frequency (namely, at the frequency such thatω C = 0; recall that, in this extremal case, it isω + =ω C ). According to Ref. [54], this asymptotic behaviour means that radial solutions obeying purely-outgoing (or purely-ingoing) boundary conditions at r → r + (see Eq.(3.31)) possess a branch point atω C = (ω + =)0 when r + = r C . 7 Such branch point is carried over to the corresponding transmission coefficient (T s , T s ,T s ,T s ).
The above relationship between asymptotic behaviour of the radial potential and branch points in the complex ω-plane is already well-known in Kerr. Indeed, the transmission coefficient in Kerr possesses a branch point at ω = 0, associated with the irregular singular point at r = ∞. This branch point gives rise to a power-law tail decay at late-times [52,56,57]. More relevantly to the study here, in extremal (r + = r − ) Kerr, the transmission coefficient has an extra branch point atω + = 0, which gives rise [58,59] to the so-called Aretakis phenomenon at the horizon [60]. It has been shown [39,41,58,61] that, in near-extremal Kerr, QNMs accumulate nearω + = 0, effectively leading to a branch cut originating atω + = 0 in (exactly) extremal Kerr spacetime. We expect that a similar accumulation of QNMs occurs in Kerr-de Sitter both nearω + = 0 and nearω C = 0 as the extremal limit r + → r C is approached. In section 4.3, we show that this is indeed the case.
Renormalized Angular Momentum Parameter. In Ref. [11], the authors use a series of hypergeometric functions to express the solution of Heun equations. The hypergeometric functions depend on the so-called renormalized angular momentum parameter ν. As discussed in [35], ν parametrizes the monodromy at infinity of the hypergeometric functions. Because the hypergeometric differential equation possesses only three regular singular points, we see that its monodromy at infinity is equivalent to its composite monodromy. In its turn, this composite monodromy of the hypergeometric differential equation is equal to the composite monodromy σ 01 of the full Heun equation modulo 2n, where n ∈ Z. That results in the following relationship: σ 01 = 2(ν − θ 0 + θ 1 ) + 1 (mod 2n), where n ∈ Z.
Angular Eigenvalues. In order to find the angular eigenvalues s λ m , we need to impose a boundary condition on the solutions of the spin-s angular equation (3.1a). We are interested in solutions which are regular in the domain u ∈ [−1, 1], corresponding to z ∈ [0, x] in the Heun equation (3.17). It is well-known in the literature [11] that (s, , m) are either all integers or half-integers and that ≥ max(|s|, |m|). (3.64) Therefore, 2θ 0 = m−s and 2θ x = −m−s are integers and so the difference between Frobenius exponents for solutions at either z = 0 or z = x is an integer. This opens up the possibility of a logarithmic singularity in one of the local solutions at each singular point. According to standard Frobenius analysis [14,36], the local solutions in this case are for some coefficients A k . If we define the vectors irr ), we find the following form for their monodromy matrices (recall Eq. (3.55)) at, respectively, z = 0 and x: We thus define a general solution of the master angular equation (3.1a) as irr , u → 1, and for some coefficients A reg = A reg (ω, m) and A irr = A irr (ω, m), and similarly for the tilde connection coefficients. Therefore, in this angular case, we have a connection problem similar to the QNM one above. Given that the monodromy matrices of the general solution (S A irr = ± 2(1 + (−1) |m−s|+|m+s| cos(2πσ)). (3.69) The regularity condition corresponds to choosing A irr = 0 in Eq. (3.67), which is equivalent to In the next section, we shall use the constraints (3.63) and (3.71) on the composite monodromy parameter σ 0x in order to find, respectively QNM frequencies and angular eigenvalues.

Quasinormal Modes in the Rotating Nariai Limit
In this subsection, we focus on the region between r + and r C (thus discarding the other horizons) and study the extremal limit r + → r C , which is also called the rotating Nariai limit [20,32,62]. The rotating Nariai limit is applied directly on the metric and it yields a new geometry which is topologically equivalent to dS 2 S 2 . 8 Here we review this limit by finding the QNM frequencies of a scalar field which possesses general coupling with the curvature of this new extremal-limit geometry. This slightly extends the results in Ref. [20] for minimal coupling to general coupling. We also show how to obtain the scalar radial equation of the rotating Nariai metric directly from the spin-s master radial equation (3.1b) of Kerr-de Sitter. The result for QNMs that we obtain in this limit here will be used as a check of the leading-order term in the higher-order expansion in the extremality parameter = (r C − r + )/L that we obtain in the next section.
(3.76) 8 The semi-direct product here means that the isometries of dS 2 must be accompanied by S 2 isometries so that they are isometries of the full spacetime. 9 The definitions used in [20] can be recovered by letting → r C L andκ → L Let us now consider a scalar field, with general coupling parameter ξ, propagating on the rotating Nariai geometry of Eq. (3.74). The corresponding Klein-Gordon equation separates by variables and, for a field mode e −iωt+imφ S m (θ)R ωm (ȳ), withω ∈ C and m ∈ Z, the resulting angular and radial equations are, respectively, and λ m is the eigenvalue of the angular equation (3.77). Notice that (3.77) does not depend on the frequencyω. Therefore, the angular eigenvalue λ m does not depend onω either. The radial equation (3.78) has three regular singular points atȳ = 0, 1 and ∞ and can be transformed to the hypergeometric equation [20] and The Frobenius exponents of (3.78) θ 0,1,∞ at, respectively,ȳ = 0, 1 and ∞ can be found from these parameters to be given via Let R in be a solution of (3.78) defined via the boundary condition for some coefficients A m . Assuming that ω ≥ 0, Eq. (3.84) represents scattering in the regionȳ ∈ (0, 1) with a purely ingoing condition atȳ = 0, while the notion of ingoing and outgoing atȳ = 1 depends on the sign of ω + mΩ . According to the asymptotic behaviour of hypergeometric functions [63], we have . (3.86) We find the QNM frequencies by requiring that A (in) m = 0 in the case ω + mΩ > 0. 10 This condition straightforwardly yields that the QNM frequenciesω =ω ± N m are given by The greybody factor for real frequencies follows readily from Eq. (3.85): Eqs. (3.87) and (3.88) agree with [20] for ξ = 0 (Eq. (3.87) also agrees with the QNM frequencies of the geometry resulting from the extremal Nariai limit r C → r + of Schwarzschild-de Sitter spacetime [64] by setting a = 0, ξ = 0, r − = 0 and λ m = ( + 1)). 11 Eqs. (3.87) and (3.88) are new results for the nonminimal coupling case (ξ = 0). Alternatively, the result (3.88) can be derived easily if we use the monodromy technique [34,35] reviewed in the previous section. The greybody factor is given by Eq. (3.60) with σ replaced by θ ∞ (as explained below): G(ω, m) = 2 sin(2πθ 0 ) sin(2πθ 1 ) cos (2π(θ 0 − θ 1 )) + cos(2πθ ∞ ) = 2 sinh(2πω) sinh 2π(ω + mΩ) cosh 2π(2ω + mΩ) + cosh(πβ) , (3.89) which is exactly the same result as Eq. (3.88). Therefore, we can avoid the calculation of the exact connection coefficients using the general formula (3.60). The poles of the greybody factor Eq. (3.89) are given by (3.87) but with N ∈ Z instead of N = 0, 1, 2, . . . . This is an intrinsic limitation of finding QNMs from the greybody factor, as the two possibilities of N ∈ Z + and N ∈ Z − appear here because we take the ratio of connection coefficients in (3.88). From the exact connection coefficients (3.85), we see that positive N corresponds toω, m > 0 and negative N corresponds toω, m < 0. The transformation {ω, m} → {−ω, −m} is a symmetry of the radial equation (3.78) but changes the 10 For ω + mΩ < 0, we should instead look for zeros of A (out) m , but it does not possess any non-trivial zeros.
11 When comparing, one should taking into account the different coordinate times.
notion of ingoing and outgoing solutions. We can decide which branch to take by going back to the solutions (3.84).

Spin-s Quasinormal Modes from the Kerr-de Sitter Master Equation.
The formula (3.89) deserves further explanation. Firstly, notice that (3.60) is valid for a generic connection problem between two singular points of a Fuchsian equation which possesses any number of singular points, as long as the monodromy matrices are diagonalizable. Secondly, as the hypergeometric equation is a Fuchsian equation with three regular singular points, the composite monodromy aroundȳ = 0 andȳ = 1 must be equal to the monodromy θ ∞ at y = ∞. This can be made even more precise by noticing that we can recover the rotating Nariai radial equation This means that, with respect to the monodromies in Eq. (3.20), we make the replacement Both ζ ∞ and x diverge in the extremal limit r + → r C as −1 , and this implies that θ 0 and θ 1 diverge as x → ∞. We can regularize the monodromies in Eq. (3.20) by defininḡ which gives, for small and fixedω, where the leading order terms should be calculated at = 0. The accessory parameter K x in (A.40), with the appropriate identifications of Eqs. (3.90) and (3.91), goes to zero as x → ∞, but x(x − 1)K x /(z − x) is finite. Applying this limit (i.e., x → ∞ -or, equivalently, → 0 -with z finite) to the Kerr-de Sitter Eq. (3.25) yields and sλ m ≡ where we remind the reader that s λ m is the Kerr-de Sitter angular eigenvalue. One can check that by taking the rotating Nariai limit on the master angular equation ( That is what we achieve in the next section.

Accessory Parameter Expansion
In the previous section, we considered an arbitrary second order Fuchsian ODE and we expressed the greybody factor (written as a product of connection coefficients) in terms of the monodromies of local solutions -see Eq. (3.60). By looking for poles of this formula we found condition (3.62) for σ +C , which sets the boundary condition for QNMs (3.63)of the Kerr-dS master radial equation for N ∈ Z + . We also similarly discussed the monodromy condition (3.71) for angular eigenvalues. We thus mapped special boundary conditions for a solution of the radial or angular equation to constraints on its monodromy exponents. This approach was proposed for the confluent Heun equation in [34,37] and for the Heun equation in [28,35]. We have in fact shown above that this method works in the particular case of the rotating Nariai limit. This limit yields a simpler case since the composite monodromy then reduces to a local monodromy of the hypergeometric equation.
In this section, we apply the composite monodromy constraints to solutions of the Heun equation (3.5) obtained from the Kerr-de Sitter master equation. Below we present our results using the APE and the numerical checks done using Leaver's method.

How to obtain the APE
Here we review the mathematical background behind the APE and its derivation done in [16]. The monodromy group of functions on the complex plane with four regular singular points are labeled by six parameters [16]. Apart from the four monodromies θ = (θ 0 , θ 1 , θ x , θ ∞ ) at the singular points, there are two extra composite monodromies, (σ 0x , σ 1x ). We thus define the monodromy data of the four-punctured sphere as the six parameters M SL(2,C) ≡ (θ 0 , θ 1 , θ x , θ ∞ ; σ 0x , σ 1x ). This means that the Heun equation is not the most general Fuchsian equation with four regular singular points, as it has only 5 parameters M Heun = (θ 0 , θ 1 , θ x , θ ∞ ; K x ) (see below Eq. (3.10)). Instead, we can say that M Heun corresponds to a reduced case of M SL(2,C) with only one composite monodromy, say, M red ≡ (θ 0 , θ 1 , θ x , θ ∞ ; σ 0x , σ 1x = σ 1x (σ 0x )). Then, M red can be mapped to M Heun via the relation K x = K x (θ, σ 0x , x). For fixed x, the APE provides the Riemann-Hilbert map between M red and M Heun . The small-x expansion of this relation is called the accessory parameter expansion (APE). In Ref. [16], the APE was obtained using isomonodromic deformations, as we review below. For other approaches to obtaining the APE, we refer the reader to, for example, [15,18,19]. For recent numerical results on finding 2D conformal maps from the APE using the approach of [16,28], see [17].
The APE obtained in [16] depends on monodromy data M red which is independent of x. This is contrary to what happens in the case of the ODEs satisfied by perturbations of Kerr-de Sitter black holes. In the case of the radial ODE in Kerr-de Sitter, x ∝ and the monodromies diverge for small . However, as discussed in section 3, if we fix the frequencȳ ω defined in Eq. (3.92), the monodromies converge as an expansion in and the accessory parameter is analytic for small . Therefore, by setting σ ≡ σ 0x and expanding all parameters (x( ), θ( ), σ( )) in terms of as we obtain the APE for some coefficients x j and θ k,j and functions K x,n = K x,n (x j , θ k,j , σ j ), j ≥ 0, n ≥ −1, k = 0, 1, x, ∞.
In the case of the angular ODE, we have x ∝ α and we set = α. The monodromies in Eq. (3.12) converge for small α and fixed L and we obtain the APE for the angular eigenvalue from Eq. (4.2). In the case of QNMs, we instead set = (r C − r + )/L and use the frequency definition in Eq. (3.92), which gives a regular expansion for the monodromies as → 0 for fixedω. Given the APE expansions in both cases, we can equate it with the exact accessory parameters coming from the angular and radial ODEs (see, respectively, Eqs. (3.18) and (3.26)), using the exact monodromies θ and the QNM condition Eq. (3.63) and the angular eigenvalue condition Eq. (3.71). This procedure, which we shall present in more detail below, will give the expansions in Eqs. (1.1) and (1.3). First, however, we shall show how to obtain the APE from a technique called isomonodromic deformation [65][66][67].
The most general monodromy data M SL(2,C) for an ODE with four singular points includes two composite monodromies, (σ 0x , σ 1x ). However, as we commented above, the Heun equation (3.5) has only one extra parameter K x . To account for these two composite mon-odromies, we write the most general Fuchsian equation as and µ(t) are some t-dependent functions with t ∈ C and 12 The form of (4.4) follows by the requirement that λ(t) is an apparent singularity, i.e., the solution has no branch cut at this point. The ODE parameters can be mapped to (σ 0x , σ 1x ) via initial conditions for the isomonodromic flow (λ(t), µ(t)), whose Hamiltonian is given by Eq. (4.4). However, there is a special set of initial conditions at t = x [28,35] that reduce the general Eq. (4.3) to the canonical Heun equation (3.5) with the same monodromy data M Heun . These initial conditions are given by In this case, we reduce the number of parameters of the ODE from two (λ, µ) to only one (K x ), and similarly for the composite monodromies, since we have σ 1x = σ 1x (σ 0x , x). This is explained in full detail in Ref. [16]. Here we just report on the solution for the accessory parameter.
The main difference between the current work and [16] is that here we want K x as an expansion in using the expansion of the monodromy data in (4.1), as described above. This changes the expansion obtained in [16], which becomes more complicated. It is more convenient to carry out the expansion for the accessory parameter of the ODE in the normal form (see appendix A for definitions) This is the form most straightforwardly connected to 2D CFT. Using the procedure outlined in appendix C, we find where we made use of the expansions in Eqs. (4.1). We have calculated this series up to order 3 in order to compare with the numerical results, but the higher-order terms are particularly long so we do not display them here.
The key point is that now we can equate two small-expansions for H x . One is the APE (4.7). The other expansion is obtained from Eq. (4.6) with the accessory parameter given by Eq. (3.26) or (3.18) in, respectively, the radial or angular cases, and expanded for small . We call the equality between these two small-expansions for H x the "APE equation". We shall next proceed to do this, first for the angular case and then for the radial case.

Angular Eigenvalues
Here we wish to apply the APE in Eq. (4.7) to find the angular eigenvalues s λ m in Eq. (3.1a). The expansion parameter in APE here is = α = a/L and we expand in α only after replacing a by αL in the monodromies in Eq. (3.12). Since a in Eq. (3.12) only appears in the combination aω/α, the mentioned expansion procedure is equivalent to expanding for both small α and small aω (without the replacement of a by αL), both quantities being of the same order in smallness. From the APE point of view, we expand it up to order α 3 , which, from the previous observation, gives all powers α n (aω) m for n + m = 0, . . . , 3 of s λ m (aω, α), as in (1.2). The monodromies (3.12) are polynomial in α (after replacing a by αL) and we expand x to one order higher in α (4.8) Using Eqs. (3.18), (4.6) and (4.8) we obtain the following expansion for the accessory parameter: Equating this with the APE (4.7), using the monodromies (3.12) and the angular eigenvalue condition (3.71) yields the APE equation in the angular case. We then use this APE equation to isolate for s λ m , thus yielding an expression which we expand for small α. We thus obtain Eq. (1.3), where the explicit form of the coefficients λ ω,k expanded for small α is given in appendix E. We have checked that this result matches (in the orders contained in both expansions) the angular eigenvalue expansion up to order (α 2 , ω 2 ) calculated in [11]. We have also checked that, in the Kerr limit L → ∞, all the orders that we give agree with the small-aω expansion given in [22]. Thus, the eigenvalue expansion that we have obtained already exists in the literature (split between [11] and [22]) -the derivation here serves to validate the APE and provides a way of extending the expansion to higher orders.

Quasinormal Mode Frequencies
Let us now calculate the QNM frequencies in Kerr-de Sitter around the rotating Nariai limit using the APE. This calculation is more complicated than in the angular eigenvalue case, as both sides of the APE equation depend explicitly on powers of the frequency ω. Therefore, we cannot just isolate ω as we did for the angular eigenvalue in the angular case. However, as mentioned above, the QNM frequency can still be calculated order-by-order in as displayed in Eq. (1.1). We show here the explicit calculation near the rotating Nariai limit and just display the result up to order . We, in fact, obtained the expansion up to order ( 3 , a 3 ) and provide it in appendix D. We numerically-checked the higher order expansion, which is discussed in the next section.
In this section, we will use the change of coordinates given by Eq. (3.22), where, in particular, r = r C is mapped to z = 0 and r = r + is mapped to z = x. The expansion parameter here is = (r C − r + )/L. Therefore, L as r C → r + . (4.10) So the extremal limit r C → r + corresponds to x → 0. We expand the monodromies (3.20) by expressing ω = mΩ + + κω and expanding for small withω fixed: (4.11) whereΩ is defined in Eq. (3.76) andκ in Eq. (3.73). We now set out to obtain the APE equation. First, we expand for some coefficientsω 0,1,2 which we set out to determine. In order to obtain one of the sides of the APE equation, we use Eq. (4.6) for H x with the accessory parameter given by Eq. (3.26) consistently expanded for small . For this expansion, we use Eq. (4.11) for the monodromies and Eq. (4.12) for the frequency, which in its turn is also used for expanding the angular eigenvalue in .
The other side of the APE equation is the APE (4.7). This expression depends on the coefficients of the expansion of the monodromies as well as on σ, for which we use the QNM condition given by Eq. (3.63). In order to obtain the final expansion for H x from the APE, because of the replacement in Eq. (3.63), the monodromies need to be expanded as per Eq. (4.11) and the frequency needs to be expanded as per Eq. (4.12).
We have now explained how to obtain the two sides of the radial APE equation. This APE equation gives polynomial equations for the coefficients ofω at each order in . The first coefficientω 0 obeys a quadratic equation and the other ω 1,2,... obey linear equations, depending on the previously obtained ones in a recursive way. We next proceed to give the explicit expression forω to leading order in , which coincides with the rotating Nariai limit result.
The first side of the APE equation, which is obtained from Eqs. (3.26), (4.6), (4.11) and (4.12), reads, to leading order in , where sλ m is given in Eq. (3.96). Equating the expansion for H x coming from (4.13) with the one coming from the APE (4.7) as explained above, yields, to leading order, the following quadratic equation forω 0 : sλ m − 2imΩ(N + 1/2) + N (N + 1) − s(s + 1) = 0, (4.14) with N ∈ Z (which comes from the QNM condition in Eq. (3.63)). The two solutions of this equation are This expression agrees with Eq. (3.87) for s = 0 and ξ = 1/6. As dictated by Eq. (3.87), from now on N is restricted to be a non-negative integer -we note that the APE method does not determine the sign of N , this is just imposed via an independent calculation (namely, the one that led to Eq. (3.87)). 13 Putting together Eqs. (3.92), (4.12) and (4.15), we have that the near-extremality expansion for the QNM frequencies in Kerr-de Sitter is: where N = 0, 1, 2, . . . . We remind the reader that Ω + ,Ω,κ and sλ m are given in, respectively, Eqs. (2.4), (3.76), (3.73) and (3.96). To the best of our knowledge, this result is new in the literature. We emphasize that, even though we write Eq. (4.16) as a small-expansion, Ω + in it is exact and so not expanded for small (the reason being that Ω + became buried withinω in Eq. (3.92), which we kept fixed when carrying out the expansions).
We remind the reader that, up to this point, all of our expressions in this subsection were considering the limit r → r + (together with the small-expansion). By doing this, we obtained the modes that approach mΩ + in Eq. (4.16). In order to obtain the modes that, instead, approach mΩ C , all we have to do is, in the quantities in Eq. (4.16), swap r + ↔ r C (which, in particular, implies → − ) and take the complex conjugate. Under this transformation, we obtain: Eqs. (4.16) and (4.17) show that, as advanced in section 3.4, both the lower and the upper superradiant bound limits (i.e., mΩ C and mΩ C , respectively) are accumulation points of QNMs as the extremal limit r C → r + is approached. In order to obtain the next order correctionω 1 , we need to expand all quantities in Eqs. (4.10), (4.11) and (4.13) to one order higher in . We then obtain: The -derivative of s λ m at = 0 can be written as (4.19) To the best of our knowledge, no orders higher than leading-order in Eq.(4.12) have been previously obtained in the literature for any spin. Above we have given explicit expressions for general a for the two leading-order terms in Eq.(4.12). Although one could, in principle, obtain higher order terms, we found it difficult to obtain them if keeping α = a/L general. On the other hand, by expanding the various quantities for small α, we managed to reach up to one higher order in . Specifically, what we did is the following. We obtained both Eqs. (4.7) and (4.13) up to including O ( 3 ). We then expanded the coefficient of each -order (including theω 0 ,ω 1 , . . . which appear within these coefficients, as well as the eigenvalue as given by (1.3) with Eqs. (E.1)-(E.4)) for small α up to including O (α 3 ). We then equated the left-and right-hand-sides at each order in this double series in α and , thus obtaining an expansion for the QNM frequencies s ω m valid up to order ( 3 , α 3 ). This expansion is the one in Eq. (1.1), where the explicit form of the coefficients is given in appendix D.

Comparison with Numerical Calculation
We have verified our expansions of the QNM frequencies and the angular eigenvalues with a numerical calculation. In this section, we first briefly sketch the approach used -we note that the method will be described in detail in [68]. We then present the numerical results as a check of the previous analyical expansions: Eq. (1.1) (with the coefficients given in appendix D) for the QNM frequencies and Eq. (1.3) for the angular eigenvalues.
The method we used for the numerical calculation is a direct implementation of the one described in [21], which, in its turn, is a generalization to Kerr-de Sitter of the technique introduced by Leaver in Kerr spacetime [23,24]. Essentially, the method consists of expressing a solution of the radial ODE (3.1b) as a powers series in "(r C − r −− )(r − r + )/((r − r −− )(r C − r + ))". By construction, this series satisfies the QNM boundary condition at the event horizon (which is the same as the ingoing boundary condition (3.34) there). However, this infinite series only satisfies the QNM boundary condition at the cosmological horizon (which is the same as the upgoing boundary condition (3.44) there) if it converges at that horizon. This convergence condition leads to an equation in terms of a 'radial' continued fraction which a frequency must satisfy for it to be a QNM frequency.
The method for calculating the angular eigenvalue is similar. One essentially expands a solution of the angular ODE (3.1a) about θ = π, thus obeying the regularity condition there by construction. The regularity condition at the other endpoint, θ = 0, is only obeyed if s λ m satisfies a certain equation in terms of an 'angular' continued fraction.
We numerically solve simultaneously for both the 'radial' continued fraction equation and the 'angular' continued fraction equation. The solution of this system of equations is then a QNM frequency ω = s ω m and an angular eigenvalue s λ m .
In Figs. 5 and 6 we plot token results of our numerical check. The plots show that each order in our small-expansion of the QNM frequencies consistently improves the expansion relative to the numerical results. The plots are for s = −2, but we also carried out numerical checks for s = 0, −1/2 and −1 and found similarly good agreement with our analytical expansions.

Conclusions
Quasinormal Modes are an interesting topic of research from theoretical as well as observational points of view. In this work, we provided a new method for calculating QNMs frequencies and angular eigenvalues for Kerr-de Sitter black holes. This method consists of using the accessory parameter expansion, given in terms of monodromies of the associated Fuchsian ODE. In particular, we have extended the APE described in Ref. [16] to the case where the monodromies are expanded in a small parameter . This allowed us to obtain analytical expansions for both the angular eigenvalues and the QNM frequencies. Our expansion for the eigenvalues agrees with results in the literature, whereas, to the best of our knowledge, our high-order expansion for the frequencies is new. Finally, we provided numer- Relative error of (ω) Relative error of (ω) ical evidence that these analytical formulas are a very good approximation for reasonably small .
There are a couple of directions for improving or extending our results. Firstly, as the τ -function used here has a complete expansion, there is a possibility that a deeper study of its properties and symmetries can provide even more compact formulas or with faster convergence than the ones used here. Our implementation of the algorithm provided in [16] is limited as the expressions of order higher than 3 are too cumbersome to be simplified with our available computing power. We are confident that simplification of these expressions is possible, but we leave this optimization for future work. Another possibility is to implement the APE numerically, as it was done in [17]. There are also other approaches to obtain the APE, like the ones in [18,19,34,69] or the direct application of CFT operator product expansions, as a variation of what was done in [15,16]. Secondly, although our method works very well for a large region of configuration space ( , a), we were limited to study the APE expansion close to the extremal limit → 0. We need the other extremal limit → 1 (r + → r − ) in order to try to cover the full configuration space for Kerr-de Sitter black holes. This could in principle be implemented with an expansion of the τ -function as t → 1, but, in this case, the expansion would be in terms of the other composite monodromy σ 1x . One would therefore need to study the relationship between σ 1x and σ 0x [15,70], which sets the relevant boundary conditions on the modes.
There are also other possible generalizations and applications of the APE method itself. For example, in [71], the scalar-field perturbation equations of a 5-dimensional Kerr-anti-de Sitter black hole were written as Heun equations and the eigenvalue and QNM conditions were written in terms of the monodromies (similarly to section 3.4 here in the case of Kerr-de Sitter). The APE method applied to this black hole setting, verified numerically, is available in [72]. Moreover, the connection between the AdS/CFT approach to QNMs and the CFT structure behind isomonodromic method still remains to be completely understood.
One could also try to extend the APE to the case of an ODE with either less or more singular points. In particular, the APE could be used to study confluent cases of the Heun equation, corresponding to two regular and one irregular singular points. The linear perturbations of Kerr and Schwarzschild black holes are described by confluent Heun equations and the monodromy setup for these cases has been studied in [37]. However, this example still lacks a systematic expansion using the Painlevé V τ -function [73,74], like in [16] for the Painlevé VI case. The same principle of finding the APE could also be used for an ODE possessing more than four singular points, using the respective τ -function [75,76]. It is known that linear perturbation equations for Kerr-NUT-(Anti-)de Sitter black holes are also reducible to separable ODEs for dimensions higher than four [77]. In particular, the number of singular points of these equations increases as the dimension increases.
Another potentially interesting, but speculative, line of investigation is to try to relate the isomonodromic approach with the Kerr/CFT correspondence [78], which so far only deals with the (near-)extremal (r + = r − ) case [79]. The near-extremality expansion used here might be helpful to understand whether the CFT picture in the extremal case could in some sense be deformed to cater for the non-extremal case. In this sense, it is curious that the τ -function can be understood as a c = 1 CFT chiral correlator (even in the non-extremal case). However, the connection with the Kerr/CFT description, if any, still remains to be understood. Discussions along these lines can be found in [28,71,80]. In particular, there is also room to explore the Kerr-dS/CFT correspondence [62], which has been less explored than Kerr/CFT. The approach of [81] to obtain the full solutions of Heun's equation from CFT is also an interesting related direction.
Finally, the APE can also be used outside the realm of black hole physics. Particular examples are the Rabi model in quantum optics and its extensions [82,83] and finding conformal mappings in 2-dimensions [17], as well as potential applications in condensed matter systems.
The APE method has opened new directions for studying spectral problems of Fuchsian ODEs -here, in particular, within the realm of black hole physics and Heun equations. We hope that the study presented here may be helpful for a better understanding of the integrable structure behind black holes, their QNMs and their connections to CFT.
We consider the case that U and V are rational functions of r and (A.1) has n regular singular points at r = r k , k = 1, . . . , n, with r n = ∞. This equation can be rewritten as As in electromagnetism, we shall call A(r) the gauge connection of these Fuchsian equations. What happens to this system after performing a homographic transformation and a s-homotopic transformation? First, let us apply the following homographic transformation where d 1,2,3,4 are complex constants. This maps the points r = r k to new points z = z k , k = 1, . . . , n, and leads to a new gauge connection with Ψ = (ψ(z),Ũ (z)∂ z ψ(z)) T . Now, let us apply the following s-homotopic transformation to the solution Thus, and the new gauge potential isĀ =Û −1ÃÛ −Û −1 ∂ zÛ . Calculating all terms, we have that This can be simplified further by writing the linear system in terms of (ψ, ∂ zψ ) T , removing U by makingΨ Removing the tildes and primes in the new gauge connection and linear system solution (i.e., relabelĀ by A,Ψ by Ψ andψ by ψ), we are left with a system of the form where Ψ = (ψ, ∂ z ψ) T . This results in the following ODE for the transformed ψ

A.1 Radial Master Equation
The master radial equation (3.1b) is of the form We now apply the Möbius transformation with i, j = 1, 2, 3, 4, and the inverse transformation We also apply the homotopic transformation for some constant β, which generates a new ODE of the form (A.7). Notice that we omitted the indices s, ω, m in the solution, as they do not play any role in this appendix. To make (A.7) more explicit, first, we havẽ Then, we can find that and It is helpful to expand Eq. (A.16) into partial fractions. Essentially, we have the following two types of terms: for some constants A, B, C, D. The main difference between these two terms is that the first one has no pole at z = ζ ∞ , whereas the second one has a pole there due to the F (z) factor. Thus, we have Then, going back to Eq. (A.16), we get Let us first check that the singularity z = ζ ∞ is removable. We start by demanding that res z=z k P (z) = 2 (β − (2s + 1)) = 0, (A.27) Both equations have a solution for β = 2s + 1 and we will assume this value for β from now on in this subsection. Now, quite surprisingly, substituting β in Eq. (A.18), we obtain (1 + s)(2s + 1) z i∞ (A.29) = (2s + 1)(2s + 2) where in the second line we used z i∞ = −r 14 ζ ∞ /r i4 , coming from Eq. (A.12), and in the last line we used the fact that there is no third order term in the polynomial ∆ r (r). Therefore, the singularity z = ζ ∞ is completely removable from the ODE. This fact has been proved for any type D vacuum solution with cosmological constant in [10]. Our derivation is an improvement to that work, as it gives the explicit dependence on the ρ i and a more compact and suggestive final result in terms of the monodromies.
At this stage, we have We now have two possibilities to simplify the Heun equation (A.7): the canonical form and the normal form.
Canonical Form. In order to find the canonical form of the Heun equation, we choose ρ k so as to cancel the (z − z k ) −2 term in Q(z) in (A.18). In order to find the Frobenius exponents ρ k , we take Its roots are given by Then, if we define θ k = iW k + s/2 and choose ρ k = ρ (+) k = s + 2θ k we obtain We can rewrite the last term in two convenient ways Comparing with the radial Heun equation (3.17), we can extract Heun's accessory parameter in (A.35) by K x = − res z=x Q(z) = −Q 3 and this gives 39) or which is the best form to multiply it by x(x − 1).
Normal Form. In the APE, we actually use the normal form of the accessory parameter, since our formulas for the τ -function are adapted to this form. The normal form is obtained by making ρ k = 1 + s, k = 1, 2, 3 in Eq. (A.32), which yields Then, H x = − res z=x Q(z) = −Q 3 , which is related to (A.40) by (4.6).

A.2 Angular Master Equation
Most of the expressions in this angular subsection readily follow from those in the previous radial case presented in Subsection A.1. Since they are separate subsections, we shall use some of the same symbols to denote quantities here defined as the quantities with the same symbol in section A.1 merely under the change r → u, and so we will not give the definitions of those symbols explicitly here. The angular master equation (3.1a) is of the form (A. 45) we have that h ν (z) is a solution of the hypergeometric equation Now, we use the two solutions of the quadratic equation (B.4) for ν, which we denote by ν ± , to obtain the composite monodromy as x → 0 where σ 0 ≡ lim x→0 σ 0x (see Eq. (4.1a)). Finally, by equating the QNM condition in Eq. (3.62) and the appropriate monodromies (4.11) onto Eq. (B.6), we obtain the correct Eq. (4.14) and the corresponding QNMs. Finally, we see that, since ν ± = θ 0,0 + θ x,0 − 1 2 ± σ 0 , the monodromy matrix of the solution z ν h ν (z) is given by The usefulness of the approach via Eq. (C.1) relies on the complete expansion of the τ -function in terms of monodromy data, first obtained in [84]. Such expansion is where the structure constants C(θ, σ) are given in terms of Barnes functions G as 17 . ( where (i, j) denotes the respective box in a Young diagram n, n i the number of boxes in row i, n j the number of boxes in column j and h n (i, j) = n i + n j − i − j + 1 its hook length.
As mentioned, the two integration constants of the τ -function are (σ, s) and, as shown in [16], the condition λ(x) = x introduces a constraint s = s(θ, σ, x) = s 0 + s 1 x + O(x 2 ), with s i = s i (θ k , σ). In order to find the APE, we just need to plug this constraint into Eq. (C.1) and consistently expand the result for small x.
The main difference between the APE used in the present paper and the one in [16] is that here the monodromies depend on the moduli x = x( ) via the parameter . As explained in the main text, we assume the monodromies and moduli have regular expansions in given by (4.1). If we solve the condition λ(x) = x for s as an expansion in using the expanded monodromies and moduli, we get s =s 0 +s 1 + O( 2 ), withs i =s i (θ k,j , σ j ), and this result can be used in (C.1) to obtain our main result (4.7).

D Coefficients of the Quasinormal Mode Expansion
In this appendix, we provide expressions below are the coefficients in Eq. (1.1) expanded up to O(a 3 ). For simplicity, in this appendix we take L = 1 and we define δ σ ≡ ( + 1) and ζ ≡ √ 12δ σ − 8s 2 + 5.