Teukolsky master equation and Painlev\'e transcendents: numerics and extremal limit

We conduct an analysis of the quasi-normal modes for generic spin perturbations of the Kerr black hole using the isomonodromic method. The strategy consists of solving the Riemann-Hilbert map relating the accessory parameters of the differential equations involved to monodromy properties of the solutions, using the $\tau$-function for the Painlev\'e V transcendent. We show good accordance of the method with the literature for generic rotation parameter $a<M$. In the extremal limit, we determined the dependence of the modes with the black hole temperature and establish that the extremal values of the modes are obtainable from the Painlev\'e V and III transcendents.

M and a = J/M are the mass and angular momentum parameter of the black hole, whereas ω, m and s are the frequency, azimuthal angular momentum parameter and spin of the perturbation. The Teukolsky Master Equation has a long history of study [2]. From its discovery, it has been crucial in the study of linear stability of black hole backgrounds. It has also been a driving force behind early numerical studies of differential equations. Analogues of the equations (1) and (2) for other black hole backgrounds, when they can be derived, are an invaluable tool to classical studies of general relativity, supergravity and string theory [3]. After the discovery of gravitational waves [4], the theory of linear perturbations of Kerr black holes became a fundamental tool to analysis of the ringdown phase post black hole collision. Scattering coefficients, on the other hand, have a fundamental role in the study of "black hole engines" of extreme astrophysical phenomena -see [5] for a review.
Fast algorithms exist for the calculations of solutions of (1) and (2) for real frequencies [6], as well as a very fast method to compute eigenvalues for the angular equation [7]. Implementations of the functions involved -the confluent Heun functions -exist in the major computer algebra systems and tables of quasi-normal modes frequencies have been compiled [8]. On the analytical part of the studies, the study of scattering and normal modes is heavily dependent of the expansion of solutions in terms of hypergeometric and confluent hypergeometric functions [9][10][11], as well as asymptotic series that permit expansions for the angular eigenvalues and quasi-normal modes [12,13].
The asymptotic nature of these expansions, however, brings up technical problems, both in the expansions and in the numerical analysis. Quasi-normal modes present such a challenge: the non-local boundary conditions involved requires consideration of general complex frequencies. As soon as one leaves the real values, the Stokes' phenomenon [14,15] presents challenges for the estimation of frequencies. The situation is aggravated in the near-extremal limit [16,17], where the analytical structure of the solutions is complicated.
In [18,19], an alternative scheme was put forward that allowed for calculation of both the eigenvalue problem as well as scattering coefficients, building in previous work [20,21]. The method related these physical properties of the solutions of Fuchsian equations, and confluent limits, like (1) and (2), to the monodromy properties associated to the equation, relying on the Riemann-Hilbert map between the parameters of the equation and the monodromy data. In [22,23], the map was cast in terms of some of the Painlevé (isomonodromic) τ -functions, whose generic expansion was proposed in [24], using methods from equivariant localization and conformal field theory. Later developments [25,26] provided with a comparatively fast method to compute those τ -functions involved.
The isomonodromic method for computing greybody factors has since been successfully used in a variety of black hole backgrounds [22,27,28]. In [23], the authors studied the eigenvalue problem of (1) and (2), giving explicit representations of the Riemann-Hilbert map in terms of the Painlevé V τ -function. The expansions of the angular eigenvalue at small rotation parameter a were recovered, and a procedure for the calculation of quasi-normal modes' frequencies was derived. As the analytical structure of the τ -function is well understood, imparting a lot of structure of the confluent conformal blocks they are derived from [29], the problems alluded to above are evaded. We also cite an alternative but similar approach based on the semiclassical analysis of the conformal blocks developed in [30,31].
In this paper, we carry on the study of [23], analyzing numerically the solution provided, comparing it with the literature for generic a < M before focussing in the near-extremal limit. In Sec. II, we review the generic Riemann-Hilbert map and the relation between the boundary conditions in the eigenvalue problem and monodromy data. In Sec. II A, we outline the method used for numerical solution of the transcendental equations involved in the eigenvalue problem. In Sec. III, we compute the relevant quantities for the equations of interest (1) and (2). In Sec. IV, we present the numeric solution for generic rotation parameter a < M , and in Sec. V we discuss the extremal a → M limit. We found from the studies that some modes display a finite behavior, analyzed in Sec. V A, while the rest display a double confluent limit, studied in Sec. V B, being given in the extremal case a = M by the Painlevé III τ -function. We close in Sec. VI by discussing the results and prospects. Finally, we include in the Appendix VI the frequencies and eigenvalues found for some modes in the extremal case.

II. MONODROMY PROPERTIES AND BOUNDARY CONDITIONS
The angular (1) and radial (2) equations are stances of the confluent Heun differential equation which is characterized by its 3 singular points, two of them regular at z = 0 and z = t 0 and an irregular point of Poincaré rank 1 at z = ∞. Following the treatment of linear differential systems in isomonodromic deformations, we call θ = {θ 0 , θ t , θ } the single trace monodromy parameters, c t0 the accessory parameter and t 0 the conformal modulus of the equation (4). In previous work, the authors have described a map from c t0 and t 0 to the monodromy data {σ, η} associated with the solutions of (4). The Riemann-Hilbert map from c t0 and t 0 to {σ, η} is defined implicitly from the equations where θ − = {θ 0 , θ t −1, θ +1}. The τ -function τ V is defined in the theory of the isomonodromic deformations, through the embedding of the equation (4) into the matricial system where t is not necessarily equal to t 0 . The function τ V is defined, up to a multiplicative constant, by the equation The function τ V is associated to the fifth Painlevé transcendent [32].
In order to make practical use of the map (5) we need a procedure to evaluate the τ -function τ V . The small t expansion for the generic fifth Painlevé transcendent was proposed in [24], using c = 1 conformal blocks. One can find the parallel discussion for the Nekrasov-Shatashvilii (semiclassical) limit of the conformal blocks in [31]. In [26], the authors made use of Riemann-Hilbert problems methods in the theory of integrable systems and gave an expression for τ V based on Fredholm determinants (see also [33]). The Fredholm determinant expression allows for efficient computation of τ V , and it is given by where 1 F 1 the confluent (Kummer's) hypergeometric series. Finally, From the definition (8), one can expand the determinant and work out the small-t expansion of the τ -function for Painlevé V, recovering the results found in the literature [23,26,34,35] whereτ V comprises the expansion of the Fredholm determinant in (8). When interpreted in terms of two variables µ = κ V t σ and t, the series obtained from the small t expansion is analytic in t and meromorphic in µ: The structure of the expansion (17) imports a great deal of structure from the conformal block expansion it is derived from [24]. Although it is not clear from (8), the τ V is almost periodic in σ, where the extra factor f is a function of the monodromy variables but not of t. This feature is more clearly stated in the Nekrasov expansion of the τ -functions associated to the Painlevé VI, V, and III transcendents [24]. Assuming σ in the fundamental domain σ ∈ (0, 2), the terms displayed in (17) are actually the most relevant ones. Each coefficient in the expansion ofτ V is a rational function on the monodromy parameters where θ enters in the numerator only and only σ enters in the denominator. The term of order µ m t n has single or double poles at integer values of −n < σ < n.
The σ parameter in τ V has an interpretation in terms of solutions of the associated differential equation (4). As seen in [36], σ parametrizes the Floquet solutions of (4), which converges in an annulus t 0 < |z| < R < ∞. Substituting this solution into the equation, the 3-term recurrence relation is where The recurrence equation can be solved using continued fractions. The result can be written as which, when truncated by the convergents of order N , and using Euler's formula for continued fractions, gives an expansion for the accessory parameter which agrees to order t N 0 with the logarithm derivative of τ V expansion of c t0 one derives from (5). The 3-term recurrence equation (20) and its solution in terms of continued fractions (25) are the basis for the so-called continued fraction, or Leaver's method [7] to compute angular eigenvalues and quasi-normal modes for the Teukolsky Master Equation. This method allows for fast convergence, but suffers near the extremal regime a → M , and the analytic structure of the accessory parameter c t0 defined implicitly by (25) is unclear. Finally, the symmetry σ → 2 − σ, which is clear from the term-by-term expansion of c t0 , and a direct consequence of the reflection property of semiclassical conformal blocks, is completely missing from (25).
The relevance of the Riemann-Hilbert map {t 0 , c t0 } → {σ, η} given by (5) relies on the fact that connection and eigenvalues problems for the equation (4) can be stated in terms of the monodromy data. For instance, one can show that the existence of Frobenius solutions with prescribed asymptotic behavior near z = 0 and z = t 0 : which do not exist for generic {t 0 , c t0 }, can be cast as simpler condition on the monodromy parameter σ: which will be relevant to the angular eigenvalue problem. By the same structure, the conditions involved in the radial problem can also be phrased in terms of monodromy data. Consider the set of local solutions of (4) We can see that y t,± consist of a basis of solutions near z = t 0 and y ∞,± are a basis near z = ∞. By algebraic manipulation of the monodromy matrices, one can show that the connection matrix C t between these local solution has the form where and ρ t ,ρ t , ρ ∞ ,ρ ∞ are normalization constants. These entries of C t can be used to relate the scattering coefficients to {σ, η}, provided there is a symmetry connecting the solutions, relating ρ i andρ i -usually time reversal [19]. For quasi-normal modes, the condition is of no energy flux out of the black hole outer horizon -corresponding to z = t 0 , and no energy flux out of infinity requires C t to be lower triangular. In turn, this requires η = η 0 with A fairly comprehensive formulation of the boundary problems involving the Teukolsky master equation in terms of monodromy data can be found in [11]. Given these monodromy conditions, the system (5) will become overdetermined, and will only allow solutions for particular discrete values of {t 0 , c t0 }. In our application, t 0 and c t0 are known once the quantum numbers s, , m the black hole parameters r + and r − and the frequency ω is set. The transcendental equations (5) will then be used to determine the σ variable in terms of these parameters through (26) and the condition that τ V = 0 will then allow for a solution only for discrete values of ω, corresponding to the quasi-normal modes' frequencies. The solution thus obtained is not unique: given the quasi-periodicity of τ V (18), to any such value of σ there is a family σ + 2n, n ∈ Z. Note that this periodicity is manifest in the continuous fraction method, where the shift in σ by an even integer in (25) is compensated by an integer shift in n.
Given the meromorphic expansion (17), we can invert the expansion for the zero of the τ -function in (5) and derive an equation for µ = κ V t σ , or e iπη as a series in t. Writing σ = 2n +σ, and supposing 0 < σ < 1, we have that τ ( θ; σ, η; t 0 ) = 0 implies where and χ V ( θ;σ; t 0 ) is analytic near t 0 = 0: For −1 < σ < 0, one can obtain a similar expression Either of the expressions (34) or (37) will hold in the fundamental branch ofσ ∈ [0, 2]. Lastly, by making use of the identity Γ(z)Γ(1 − z) = π/ sin πz, we can show that, with the quantization condition for η (33), we have Taking these relations into account, the equation for the zero of τ V turns into As a final remark, the successive terms in χ V (36) can be computed from the Fredholm expansion, but it is computationally more efficient to solve for τ V = 0 with (33) directly.
A. Evaluating the τ -function The Fredholm determinant formulation of the τ -function (8) allows us with a method to compute it numerically up to order t N , in polynomial time O(N α ). There is a variety of methods to numerically compute Fredholm determinants [37]. One such method is to truncate the kernels of the operators A and D c defined in (9) through simple Riemann quadratures: with the l'Hôpital's rule for the diagonal terms. This method allows for fast evaluation in cases where |t| R, but relies on implementations of the hypergeometric and confluent hypergeometric functions which may not be compatible with arbitrary-precision arithmetic.
On the other hand, in order to recover the Nekrasov expansion in [24], we must use a different basis for truncation. First we expand the parametrices and compute the matrix elements associated to A and D c in the Fourier basis g(z ) = n g n (z ) n . The matrix-valued coefficients G n (σ, θ t , θ 0 ) and G c,n (−σ, θ ) can be computed from the expansion of the Gauss hypergeometric series: and where (z) n = Γ(z + n)/Γ(z) is the Pochhammer symbol. The kernels A(z, z ) and D c (z, z ) can be suitably expanded: The resulting matrices are semi-infinite, and truncation to order N gives an approximation to the Painlevé V τfunction of order O(t N , |t| (1± σ)N ). We refer to [25] for the corresponding calculation for the Painlevé VI τ -function 1 . From the computational point of view, this expansion is costlier, but has the bonus of not requiring pre-existing implementations of special functions, except the gamma function.

III. ANGULAR AND RADIAL SYSTEMS
We are now ready to present the solution of the problem as a series. The angular equation (1) can be brought to the standard confluent Heun form (4)θ in (1) being the independent angular variable -by the change of variables where one can read the monodromy parameters and the modulus and accessory parameter can be obtained directly after some manipulations t Ang = −4aω, t Ang c Ang,t = λ + 2aωs + a 2 ω 2 .
As derived previously [23], the τ -function expansion of the accessory parameter (5) can be used, along with the quantization condition (28) can be used to derive an expression of the angular eigenvalue s λ ,m s λ ,m (aω) = ( − s)( + s + 1) − 2ms 2 ( + 1) aω which agrees with the value found in the literature [12]. Along the same lines, the radial equation can be brought to the canonical form by changing variables where We can now define the modulus and accessory parameter for the radial equation For the subsequent analysis, let us define in terms of which the parameters (51) and (52) are given by and We remark that θ + + θ − = 2s + 4iM ω has no explicit dependence on ν. This ν parameter defined by (53) correlates with the black hole temperature as ν → 0: Finally, for real ω, the imaginary part of the monodromy parameter associated to the outer horizon θ Rad,t where δS has the interpretation, given by the first law of black hole thermodynamics, as the entropy increase of the black hole when it absorbs a quantum of field with energy ω and angular momentum m. In [23] the authors expanded in this relation and the implications to any putative underlying quantum theory. We can now proceed to the study of the equations (58).

IV. RESULTS FOR GENERIC a
The solution of the eigenvalue problem for the radial equation is written as (5) where the value η 0 is given in terms of θ Rad and σ by the quantization condition (33). By substituting the expansion of the angular eigenvalue (48), the system (58) can be seen as transcendental equations determining M ω and σ. By the procedure outlined in Sec. II A, these equations can be solved numerically. We are going to restrict ourselves to the fundamental mode throughout. The implementation of the τ -function was done in Julia language using ArbNumerics, a port of the Arb C language library for arbitrary-precision ball arithmetic, with 192-digit accuracy. The determinant in (8) was truncated at N f = 36, and the angular eigenvalue was computed the continuous fraction method (25), capped at N c = 64 convergents. The roots of the transcendental equations (58) were found using a simple Newton method. In Fig. 1, we compare the results with those using the continued fraction method [8,38] for s = −2 and = 2. In all light modes studied, there was excellent agreement with the modes studied and the literature.
The theorem on the Painlevé property of the isomonodromic τ -functions, of which τ V is an example, assures that the roots sought for in (58) are isolated [39] away from the essential singularity at t 0 = 0. Near that point, there is an accumulation of roots, which complicates the numerical analysis in the extremal limit. In addition, the equations (58) are satisfied by any mode, be it fundamental, excited or even non-normalizable. In order to make sure we are following the right root, an initial guess based on the results in the literature is made at reasonable values of a/M and then the root is followed by small variations of the a/M factor. The difference between the roots found in the literature are shown in Fig. 2, with excellent agreement throughout, and the discrepancies increasing as one approaches the extremal limit, where one expects the method based solely on the continued fraction expansion of the accessory parameter to fail. We will see below in more detail the dependence of frequencies as r − → r + , where the continuous fraction method of the literature is not valid.

V. THE a → M LIMIT
The extremal, or r − → r + , limit can be studied using the τ -function by taking the appropriate limit of the conditions (5). The limit evades somewhat the numerical hurdles alluded to in the last Section, as well as provides analytic tools to study the extremal case.
Let us first note that the angular equation has a smooth limit of the parameters (46), (47), and then the eigenvalue expansion will essentially be the same as above (48). We will assume that one can in principle compute s λ ,m as a function of the frequency near the extremal value. We will argue that this is in principle an amenable task.
The radial equation, on the other hand, will have a more complex limit depending on the particular mode. Given the ν parameter defined above (53), we now define in terms of ν the confluence parameter and the new isomonodromic variable We observed two distinct behaviors for M ω as we approach the extremal limit ν → 0: A. M ω converges to m/2 with ν or higher order as ν → 0. In this case the confluence parameter has a finite limit, and the system is actually well described by (58); B. M ω does not go to m/2 as ν → 0. In this case, the confluence parameter Λ will diverge and one has to consider the confluent limit of the equations for the radial system (58).
We now proceed to analyze each case separately.

A. The Finite Λ Limit
We observed numerically that for the modes = m, with m = 0, the eigenfrequencies tend to m/(2M ) in the extremal limit with ν or higher. To describe the behavior of the solutions of (58) in this limit we propose the expansion M ω = m 2 + β 1 ν + β 2 ν 2 + . . . , σ = 1 + α, α = α 0 + α 1 ν + α 2 ν 2 + . . . , where the coefficients β i and α i can be computed recursively from the equations (26) and (39). The consideration of the series is simplified from the fact that the expansion parameter z 0 of the expressions for the accessory parameter (26) and (36) is small and Λ defined through (59) is finite. We then have that each term of the expansions (26) and (36) is finite and the term of order t n 0 is of order ν n . In terms of ν, the accessory parameter for the radial equation is Finally, we will also expand the angular eigenvalue (48) where we note that, despite being an expansion in aω = M ω cos ν, it can have odd terms in ν through its dependence with M ω. We also remark that λ 0 = s λ ,m (m/2) is the extremal value of the angular variable. Substituting the series for M ω and σ into the expansion for the accessory parameter (26), we can relate the coefficients α i with β i and λ i , and compute them recursively. The first terms are: The sign of α 0 chosen will actually depend on the mode. We can now substitute the expansions into (39) and find a transcendental equation for the β i . Supposingσ = 1 +α, with −1 < α 0 < 1, the first non-trivial term from (34) is The expansion of χ V in (34) is analytic in ν, whereas the expansion of the tσ −1 and Θ V will include non-analytic terms like ν log ν. As one takes ν to zero, the term να above will go to zero if α > 0, and the only way to satisfy the equation will be if one of the gamma functions' arguments in the numerator becomes very close to zero. Then So, if α 0 is real, which should be expected if m is small, then one picks the sign so that α 0 > 0. This behavior seems to be restricted to the mode s = 0, = m = 1 and it is shown in the top part of Fig. 3. Note that in this case the correction to the real part of the frequency is of higher order in ν. When m is large enough, one expects α 0 to be purely imaginary, and then the ν-dependent term in (67) will oscillate logarithmically. This correction term will still be small if α 0 < 0, which selects the negative root in (65). The behavior is much more common and it is represented in the bottom part of Fig. 3, where both real and imaginary parts of the frequency display the linear behavior with ν. Furthermore, in the case whereα 0 is purely imaginary, the imaginary slope is approximately −ν/4. These modes are perturbatively stable, but their decay time diverges as one approaches the extremal limit.
Finally, we note that the equation (66) admits infinite roots of the sort β 1 → β 1 − in/2, corresponding to the poles of the gamma function at negative integral values of the argument. This facts leads to the accumulation of zeros of the τ -function alluded to in Sec. II. Similar remarks in the same context were made from the continuous fraction expansion in the excellent papers [16] and [17], albeit the last one with a slightly different value for α 0 . Some of the results in this section were anticipated by [40].

B. The Confluent Limit and the Third Painlevé Transcendent
All modes with m = , including those with negative m, will not tend to M ω ext = m/2 in the extremal limit. In this case, the parameter Λ (59) goes off to infinity, and the equations (5) will undergo a confluent limit, where Λ → ∞ with u = Λt finite. In order to write the extremal version of (58), we first have to take the confluent limit of the τ -function (8).
As the calculation is relatively short and to our knowledge not present in the literature we include it here. We start by considering (8) and deform the circle C multiplying its radius by t. This has the effect of shifting the t dependence in the argument from the kernel of D to A, so that With this provision, we can now implement the confluent limit on the parametrix Ψ where θ • = θ t +θ 0 is fixed and Ψ c (σ, θ • ; uz) is the same confluent parametrix as above (14) and the first Λ −1 correction is also given in terms of confluent hypergeometric functions (71) Given the expansion of the parametrix, the expansion of the kernelÃ(z, z ) then follows whereÃ We now turn into the confluent limit of the monodromy parameters. We will assume that η has a well-defined limit, which is certainly the case in our application. The parameter κ has a well-defined function in terms of Λ, provided Λ is not close to the negative real axis. Expanding the gamma functions in (15), we find For convenience, we will refer to κ III = Π III e iπη and µ = κ III u σ in the following expressions. The first term in the A kernel expansion (72) gives the Painlevé III τ -function This definition can be compared to the expansion given in [24] by comparing the first terms, see (80) below. It will be interesting to consider the first order term in Λ −1 of the expansion of (8). Using well-known properties of the determinant, we have where, again, µ = e iπη Π III u σ . We note that the correction is well-defined even when the determinant vanishes. Generically, for finite-dimensional matrices, is the adjugate to M, which is the transpose of the cofactor matrix. The calculation of the τ III from (76) follows the same strategy of Sec. II A, expanding the parametrices where the coefficients G c,n are the same as (43). The expansion of the confluent kernels A c (z, z ) and D c (z, z ) is analogous to (44), and the expansion of the Painlevé III τ -function (76) gives For small u, the first correction in order Λ −1 is surprisingly simple so, to first order in Λ −1 , the zero of the τ -function does not change from the extremal value. The confluent limit of (5) can now be written explicitly where u 0 k 0 is the confluent limit of t 0 c t0 . We solve these conditions by using the same procedure as the one used with τ V . Inverting the series to find the value of η corresponding to the zero of τ III , we find an expression analogous to (34), where, as before, and As the quasi-periodicity in σ is inherited from τ V , we have the same remarks about the fundamental domain of σ as before. In the expansion of (83), we assumed 0 < σ < 1. Analogous expansions for σ < 0 can be derived following the same procedure that led to (37). The expansion of the accessory parameter k 0 also follow the same lines, and has a structure parallel to (26): Finally, the quantization condition (33) also has a well-defined confluent limit. If Λ goes to ∞ in a ray with argument sufficiently close to π/2, we have e iπη0 = e −2πiσ sin π 2 (θ + σ) sin π 2 (θ − σ)  4. The near-extremal behavior for the fundamental quasi-normal frequency for s = −1, = 2 and m = 1. We note the roughly quadratic dependence for small ν, which is consequence of the first near-confluent correction to the τIII -function involved.
For our application, this limit should hold if M ω > m/2 in the extremal limit. If M ω < m/2, then the argument of the exponential factor should be replaced by +2πiσ.
With the quantization condition, we can use the gamma function reflection formula to simplify (85) into We are now ready to consider the extremal ν = 0 case. Let us define the extremal variables and the results for low-lying modes are shown in Tables I, II and III. As can be checked, none have M ω close enough to m/2 to warrant a meaningful truncation of the expansions (88) and(86) at small order in u 0 . Based on the analysis above, however, we can remark that, given the Λ −1 corrections to (86) and (88) only appear at order ν 2 , we have that the near-extremal corrections to the frequencies below are all of order O(T 2 + ), as illustrated in Fig. 4.

VI. DISCUSSION
In this paper, we applied the isomonodromic method outlined in [19,23], and revised in Sec. II to study quasi-normal modes of the generic rotating Kerr black holes. The procedure outlined in Sec. II A provides us with a numerical procedure which is an alternative to the continuous fraction method which, despite having a slower convergence, has a more controlled behavior at the extremal limit r − → r + . In the method, the eigenvalue problem is reduced to solving two transcendental equations (5), which is numerically amenable given the analytical properties of the functions involved. In Sec. III we translated the conditions satisfied by the angular and radial equations, and in Sec. IV we showed the excellent accordance between the isomonodromic and the continued fraction method for generic rotating parameter.
In Sec. V we turned to the extremal limit a → M limit. Based on the results of the previous section, we found two distinct behaviors as r − → r + . In the = m modes, we observed in Sec. V A that the eigenfrequencies approached ω ext = m/(2M ), in a manner generically linear with the temperature. In the rest of the modes, the limit was seen in Sec. V B to be generically described by a confluent version of (5), involving the third Painlevé transcendent, whose Fredholm determinant formulation for the τ -function derived. By keeping the first near-extremal correction, we could assert that the extremal value for the frequencies approached the extremal value in a manner quadratic with the temperature.
As was anticipated in [23], the method does not converge faster than the standard continuous fraction method, but on the other hand has its analytical properties more transparent. Specifically, one can derive asymptotic formulas and estimate the error terms in expressions like (67). Although we restricted our study to the fundamental mode, the method can be used to study higher modes, as they also satisfy (5). We hope these initial results contained in this text show the promise of the application of isomonodromy to study problems related to accessory parameters of Fuchsian differential equations and their confluent limits. This is a very generic problem in the field of mathematics, which include scattering and eigenvalue problems as particular subcases.
We have found that in general the numerical results fall in excellent accordance with the results in the literature, with 7-8 digit accuracy. While one has, for the numerical methods of choice, that the Stokes phenomenon displayed in the solutions of (4) complicates the enforcement of boundary conditions as soon as one considers non-real frequencies, it may be the case that, in principle, non-analytical (in t) terms contribute to the expansion (8). Either option is fully deserving of separate study. In a less technical matter, it is still not clear which mechanism selects a priori the different behaviors seen at the extremal limit.
We hope that the method proves to be useful in other black hole backgrounds whose perturbations are governed by solutions of the confluent Heun equation (4), specifically the Reissner-Nordström and the Kerr-Newman backgrounds. The more technical relations to conformal field theories and integrable systems outlined in [23] are also a very promising prospect, whose understanding may shed light on a quantum description of the degrees of freedom involved in the perturbation, and the finite, non-linear perturbations. We hope to be able to return to these points in future work. modes are special in that they do not undergo the confluent limit and are best described by the method in Sec. V A.