Exact solution for wave scattering from black holes: Formulation

We establish an exact formulation for wave scattering of a massless field with spin and charge by a Kerr-Newman-de Sitter black hole. Our formulation is based on the exact solution of the Teukolsky equation in terms of the local Heun function, and does not require any approximation. It serves as simple exact formulae with arbitrary high precision, which realize fast calculation without restrictions on model parameters. We highlight several applications including quasinormal modes, cross section, reflection/absorption rate, and Green function.


I. INTRODUCTION
Observational efforts to prove black holes (BHs) finally began to bear fruit in the last few years: direct detection of gravitational waves (GWs) emanating from a merger of binary BHs [1], and electromagnetic observations with very-long-baseline interferometry (VLBI) [2]. The growing global network of ground-based GW interferometers and VLBI multi-wavelength observations at higher resolution in the near future will be powerful tools to unveil the nature of BHs or to test the Kerr hypothesis. By virtue of the uniqueness of the Kerr solution in General Relativity, they allow us unprecedented tests of gravity in the strong-field regime. Theoretical prediction of the propagation of fields with different spins on BH geometry is thus important.
While in the short-wavelength regime one can rely on the geometrical optics approximation, in the long-wavelength regime the approximation breaks down and one needs to take into account wave optics. Most of the physically interesting cases of the wave equations on BH geometries can be solved by separation of variables. The separability of the Klein-Gordon equation for the Kerr-Newman family with a cosmological constant was clarified by Carter [3][4][5]. This result, with the aid of the Newman-Penrose formalism [6], was generalized to higher-spin wave equations for the Kerr background [7][8][9][10] and for the Kerr-de Sitter background [11,12]. With the separated master equation known as the Teukolsky equation, one can investigate wave propagation, for which scattering analysis is a powerful approach [13]. Observables in the wave optics have been commonly evaluated in the literature with certain approximations such as the WKB approximation.
In [14], Suzuki, Takasugi, and Umetsu (STU) showed that both angular and radial parts of the Teukolsky equations for a massless field with spin and charge on the Kerr-Newman-de Sitter (KNdS) spacetime can be transformed into the Heun equation. *1 This result was further generalized to Petrov type-D vacuum backgrounds with a cosmological constant [33]. The Heun equation is a second-order linear homogeneous ordinary differential equation with four regular singular points [34][35][36][37][38]. There are several types of exact solutions for the Heun equation, depending on the analyticity around singular points. Among them, in a series of works [14,39,40], STU adopted a series of the hypergeometric functions to construct an exact solution for the Teukolsky equation on the KNdS background, along the same lines as the Mano-Suzuki-Takasugi formalism [18][19][20] for the asymptotically flat background. They derived an exact formula for the absorption rate in terms of an infinite series. Their formalism was also applied to the calculation of the quasinormal mode (QNM) frequencies for the Kerr-de Sitter black hole [41], generalizing Leaver's method [42]. On the other hand, in recent work [43] on the Kerr-de Sitter background, Hatsuda employed a simpler exact solution known as the local Heun function or simply the local solution, and obtained a compact formula for the QNM frequencies with arbitrary high precision.
In this paper, we consider a massless test field with spin and charge on the KNdS background, and establish an exact formulation of the scattering problem using the local Heun function. The formulation based on the local Heun function is transparent and provides us with concise formulae for black hole physics such as the greybody factor and Green function. One can evaluate specific values of the local Heun function by using a modern technical computing system, Mathematica, which implemented the various Heun functions as built-in functions in the version 12.1 update in 2020.
The rest of the paper is organized as follows. In §II, we transform the Teukolsky equation for a massless field on the KNdS background to the Heun equation, and provide the exact solution in terms of the local Heun function. In §III, we consider the boundary condition at the horizons, and obtain the connection coefficients, which allow us to solve the scattering problem with the exact solution. In §IV, we highlight several applications of our formulation such as QNMs, S-matrix, cross section, reflection/absorption rate, and Green function. §V is devoted to the conclusion.

II. EXACT SOLUTION
In this section, following [14,43], we present the exact solution of the Teukolsky equation on the KNdS background. After summarizing our notation of the KNdS metric in §II A, we review the transformation of the angular and radial parts of the Teukolsky equation into the Heun equation in §II B. In §II C we provide the exact solution in terms of the local Heun function. We consider the boundary condition of the angular solution in §II D, and deal with the radial solution in §III.
Here, Λ is the cosmological constant, and M, aM , and Q are respectively the mass, angular momentum, and charge of the black hole. The electromagnetic field caused by the charge of the black hole is given by One can consider several limiting cases. For instance, Q = 0 reproduces the Kerr-de Sitter spacetime, whereas Q = 0 and a = 0 reproduce the Schwarzschild-de Sitter (SdS) spacetime. Throughout the paper, we assume Λ > 0, and focus on the case where ∆(r) = 0 has four distinct real roots under the condition [44] where We denote the four roots of ∆(r) = 0 as r ± , r ± . We can then factorize ∆(r) as We set the ordering of the four roots as r − < 0 ≤ r − < r + < r + , where r − , r + , and r + are the inner (Cauchy) horizon, outer (event) horizon, and cosmological horizon, respectively. We are interested in the scattering problem in the range r + ≤ r ≤ r + . Comparing (6) with (2), it holds that r − + r − + r + + r + = 0.
Note that so long as ∆(r) = 0 has four distinct roots, our arguments in §II apply to the asymptotically AdS geometry with Λ < 0, except for the ordering of the four roots. For the KNdS spacetime with ΛM 2 1, we have For the SdS case, if 0 < ΛM 2 < 1/9, there are four real roots, which can be expressed in a simple expression where ξ = 1 ΛM 2 For ΛM 2 1, and hence For instance, for the SdS with ΛM 2 = 10 −3 , we have r − /M = −55.75, r + /M = 2.005, r + /M = 53.74. The tortoise coordinate r * is defined by or where yields the surface gravity at the horizons.

B. Transformation of the Teukolsky equation into the Heun equation
We consider the propagation of a massless field with spin s and charge e on the KNdS background. In terms of the Newman-Penrose formalism, the master variables ψ s are given by where each case corresponds to the gravitational, electromagnetic, Dirac, and scalar field, respectively. Note that s = 0 corresponds to a conformally coupled massless scalar field, whose equation of motion is given by ( The Teukolsky equations for spin 0, 1 2 , 1, 3 2 , 2 fields on the Kerr-de Sitter background and those for spin 0, 1 2 fields on the KNdS are separable and take the unified form [40]. With and the separation constant λ, the angular and radial parts of the Teukolsky equation are given by *2 where x = cos θ, c = aω, ∆ = d∆/dr, and It was clarified in [14] that the angular and radial Teukolsky equations on the KNdS background can be transformed into Heun equations, and the exact solution was constructed in terms of a series of hypergeometric functions. Regarding the transformation from the Teukolsky equation to the Heun equation, there are 4! = 24 independent transformations depending on how to map the four regular singular points. For the angular part we follow the transformation adopted in [14], whereas for the radial part we follow the transformation adopted in [43] for the Kerr-de Sitter background, so that our parameter regions of interest, −1 ≤ x ≤ 1 or r + ≤ r ≤ r + are mapped to 0 ≤ z ≤ 1, where z is the independent variable after the transformation. Further, it was shown in [33] that, for massless field on Petrov type-D vacuum backgrounds with a cosmological constant, the separated Teukolsky equations can be transformed into Heun equations. While we focus on (18) and (19) on the KNdS background, our analysis can be straightforwardly generalized to such a case.

Angular part
Let us begin with the angular part (18) of the Teukolsky equation. Since (18) does not depend on the charge Q, the argument on the angular part remains the same regardless of the charge. Note also that the cosmological constant Λ enters the equation only via α = Λa 2 /3. Therefore, for the nonrotating limit a → 0, for which α (and c) vanishes, the equation does not depend on Λ.
For simpler geometries, the angular Teukolsky equation (18) allows a simple exact solution. For nonrotating black hole, i.e., Schwarzschild(-de Sitter) or Reissner-Nordström(-de Sitter) black hole, the exact solution is known as the spin-weighted spherical harmonics, S s (θ)e imϕ = s Y m (θ, ϕ), with the eigenvalue λ = ( + 1) − s(s − 1). The explicit form is given by (22) For the Kerr or Kerr-Newman geometry, the exact solution is denoted as the spin-weighted spheroidal function. No analytic expression for the eigenvalue λ is known in this case.
For the more general case of a rotating black hole in the presence of the cosmological constant, the spin-weighted spheroidal function is not the analytic solution. However, we can still derive the exact solution since the angular equation (18) can be transformed into the Heun equation. With a nonzero cosmological constant, the angular Teukolsky equation (18) has four regular singular points at x = ±1, ±i/ √ α after removing a removable singularity at x = ∞.
We transform the independent and dependent variables as *2 There is a typo in the angular equation (3.1) in [40]: The second last term −2m(1 + α)ξ in the second line should be +2m(1 + α)ξ.
to map the four regular singular points (−1, 1, −i/ √ α, i/ √ α) to (0, 1, z a , ∞). Here, the superscript (a) denotes the angular part. Note that the boundaries x = −1, 1 are now mapped to z = 0, 1, respectively. Here, we denote z ∞ = z| x→∞ and z a = z| x→−i/ √ α , namely, and define which satisfy an identity The transformations (23) and (24) allow us to rewrite the angular equation (18) as where Equation (28) is nothing but the Heun equation, at which we shall take a closer look in §II C.

Radial part
Next, we proceed to the radial part of the Teukolsky equation (19). The equation has four regular singular points at r = r ± , r ± after removing a removable singularity at r = ∞. We transform the independent and dependent variables as to map the four regular singular points (r + , r + , r − , r − ) to (0, 1, z r , ∞). Here, the superscript (r) denotes the radial part. To avoid notational complexity, here we use z to denote the independent variable as in the angular part, but no confusion should occur as the arguments on the angular and radial parts are independent of each other. Note that the black hole horizon r = r + and the cosmological horizon r = r + are now mapped to z = 0, 1, respectively. Therefore, again, the parameter range that we are interested in is 0 ≤ z ≤ 1. Here we denote z ∞ = z| r→∞ and z r = z| r→r − , namely, both of which are larger than unity. Also, we define a purely imaginary function and denote which satisfy an identity With the transformations (30) and (31) and the identities (7) and (35), the radial Teukolsky equation (19) can be rewritten as where These expressions are much simpler than those in [14] and a natural generalization of those in [43] for Q = 0.

C. Local Heun function
In §II B 1 and §II B 2, we see that we can transform the angular and radial Teukolsky equations into (28) and (36) respectively, which are the same type of differential equation, as pointed out first in [14]. This type of differential equation, i.e., the second-order Fuchsian equation with four regular singular points on the Riemann sphere, is known as the Heun equation [35][36][37], which is given by with the condition γ + δ + = α + β + 1, a = 0, 1.
The Heun equation has six independent parameters. a is called a singularity parameter, α, β, γ, δ (and ) are called exponent parameters, and q is called an accessory parameter. In §II C only, we use α, β, γ, a to denote the parameters of the Heun equation, rather than the parameters for the KNdS geometry. The angular and radial Teukolsky equations in the forms (28) and (36) are nothing but the Heun equation (40) and respectively. Note that the conditions (41) are satisfied by virtue of the identities (27) and (35). The Heun equation has four regular singular points at z = 0, 1, a, ∞. At the vicinity of each regular singular point, we can construct two linearly independent local solutions, or Frobenius solutions. Following the standard notation, we denote the local Heun function Hl(a, q; α, β, γ, δ; z) as the canonical local solution of the Heun equation at z = 0, namely, where the coefficients c k are defined by the three-term recurrence relation The local Heun function (44) converges for |z| < min (1, |a|). Therefore the maximum of the radius of convergence is unity for |a| > 1. However, the local Heun function Hl can be analytic at z = 0, 1 for some discrete values q = q m (m = 0, 1, 2, · · · ). In this case the function is called the Heun function and is denoted by Hf . Further, it can be analytic at z = 0, 1, a with α = −n (n = 0, 1, 2, · · · ) and q = q m (m = 0, 1, 2, · · · , n). In this case the function becomes polynomial and is called the Heun polynomial Hp. In this paper, we only use the local Heun function (44). The local Heun functions at z = 0, 1 are of special interest to us in discussing scattering from black holes. Two local Heun functions at z = 0 are given by y 01 (z) = Hl(a, q; α, β, γ, δ; z), and two local Heun functions at z = 1 are given by The asymptotic behavior of the exact solutions (46)-(49) is determined by the characteristic exponents The local Heun functions at z = 0 are related to the local Heun functions at z = 1 via linear combinations y 01 (z) = C 11 y 11 (z) + C 12 y 12 (z), (52) y 02 (z) = C 21 y 11 (z) + C 22 y 12 (z).
The connection coefficients are formally given by the ratio of the Wronskians as where W z [u, v] = u dv dz − du dz v. Note that from (36) it holds that, for linearly independent solutions y a , y b , Therefore, while the Wronskian itself is not constant, the ratio between two Wronskians is constant. Conversely, the local Heun functions at z = 1 can be expressed as y 11 (z) = D 11 y 01 (z) + D 12 y 02 (z), (56) y 12 (z) = D 21 y 01 (z) + D 22 y 02 (z), where namely, While the connection coefficients can be formally written down analytically [45], this approach requires the evaluation of the local Heun function on the maximum convergence radius, and in general it is not clear whether it is convergent [46]. Even if it is convergent, it typically requires the analytic continuation of the local Heun function, which has a high computational cost. The expressions (54) or (59) are more practical. To obtain the connection coefficients C ij or D ij , one can evaluate the right-hand sides of (54) or (59) at any z within the overlapping region of the two disks of convergence. The advantage of this formulation is that the scattering problem is defined between z = 0 and 1 and the calculation remains within the circle of convergence of local Heun functions at z = 0 and 1. This situation should be compared with the case where one needs a calculation outside the circle of convergence, for which one needs analytic continuation or other types of exact solutions of the Heun equation valid for a wider range, such as hypergeometric function series. In our case, we can calculate the connection coefficients at some point between z = 0 and 1 without analytic continuation. We shall see in §III that the connection coefficients play a central role for the scattering problem.
For the specific calculations in the present paper, we use the built-in function HeunG implemented in Mathematica 12.1 or later, which yields the local Heun function Hl (44) inside the circle of convergence, whereas it gives an analytic continuation of Hl outside the circle of convergence. The analytic continuation typically takes more computational time, and sometimes causes a multi-value issue. For the radial Teukolsky equation, since a = z r > 1 holds, the radius of convergence for the local Heun functions (46) and (47) at z = 0 is unity. Therefore, there always exists an overlapping region of the two disks of convergence for the local Heun functions at z = 0 and z = 1, where we can use both local Heun functions without analytic continuation. The general solution of the radial Teukolsky equation (36) can thus be written as a linear combination of y Ii,s denotes the radial exact solution, i.e., the exact solution y Ii with the parameter set (43) for I = 0, 1 and i = 1, 2. We define the angular exact solution y (a) Ii,s in the same manner with the parameter set (42). For the scattering problem, we shall focus on two specific radial solutions imposing a certain set of boundary conditions, which we shall discuss in §III. We shall also see that both local Heun functions are useful to see the asymptotic behavior close to the black hole horizon or cosmological horizon.

D. Angular solution
Before proceeding to the scattering problem with the radial solution in §III, let us check the requirement on the regularity of the angular solution in terms of the exact solutions. Since the angular Teukolsky equation (18) does not depend on the charge Q, we can directly apply the argument of the angular part in [43] for the Kerr-de Sitter case. From (50) and (51), we see that the angular function S Ii, The general solution S s (x) is given by a linear combination of S Ii,s (x). To make the angular solution regular at x = ±1, we should respectively choose S 01,s (x) or S 02,s (x) for s − m 0, and S 11,s (x) or S 12,s (x) for m + s 0. For S s (x) to satisfy both regularities at x = ±1, we require linear dependence of the exact solutions, namely, For a nonrotating black hole with a/M = 0, this equation is satisfied by the eigenvalue λ = ( + 1) − s(s − 1). For a rotating black hole, this equation depends on λ and ω implicitly. For a fixed frequency ω, this condition determines λ, which we can obtain by using a root-finding algorithm. On the other hand, to obtain the QNM frequencies, we should solve (62) and a boundary condition on the radial solution to obtain λ and ω simultaneously, as we shall see in §III and §IV. In either case, we need an initial input value sufficiently close to the roots. In Fig. 1, we present the eigenvalue λ for scalar waves on the Kerr-de Sitter background obtained by the above method. We compare our exact results with the analytic expansion formula given by Eq. (4.18) in [14] for small aω and Λa 2 /3. We denote these two results as λ Heun and λ STU , respectively. So long as one considers low-frequency waves scattered by a slowly rotating black hole with a small cosmological constant, the analytic expansion formula works well and the difference between λ Heun and λ STU is negligible. To see its validity and limitation, we consider a rapidly rotating black hole a/M = 0.9 with a small cosmological constant ΛM 2 = 10 −3 . In the left panel of Fig. 1, λ Heun and λ STU are shown by solid and dashed curves, respectively, for m = and = 2, 4, 6. For the calculation of λ Heun , we pick up sampling points with the interval ∆(M ω) = 0.05 for the range 0 ≤ M ω ≤ 3. We take λ = λ STU as the initial input value for the root-finding algorithm FindRoot in Mathematica, and set PrecisionGoal → 15. For the algorithm to work well with this initial input, we need to set PrecisionGoal larger than 12. To get the plots in Fig. 1, we use ParallelTable with 8 cores and get the list of data. The computation time for each curve is about 2.5 sec. In the right panel of Fig. 1, we present the relative errors between λ Heun and λ STU . As expected, the relative error increases as the frequency increases. In this setup, we see that for M ω ≤ 1 and ≥ 2, the error remains O(10 −1 )%, so it is reasonable for this parameter range to use the analytic expansion formula. On the other hand, for low-multipole and high-frequency waves, the error of the analytic expansion formula becomes large, and hence one should use the exact formula.
L X c q g 8 e p / P u / q i r N P n a / V H 9 6 9 r G N u d C r I O 9 O y A S 3 M J v 6 + u H p a 3 4 + N 9 G Y Z F f s h f x f s i d 2 T z e w 6 m / m 9 S r P X S B O H 6 D 9 f O 5 W U J x O a z N p b T W j Z h e j r + j G K M Y x R e 8 9 i y y W s Y I C n b u H E 5 z h P P a i J J W U M t J M V W K R Z h j f Q l E / A d 0 t i 8 I = < / l a t e x i t >`= 4 < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 f 0 q g i m n k T V f I X 4 m x u t N y c e W U i U = " As a test of the application range of the present method, we check the case where λ = 0 is adopted as the initial input value. For this initial value, the root-finding algorithm requires a longer computational time and larger value of PrecisionGoal since the initial values for larger are far from the true value. For example, if we set PrecisionGoal smaller than 15 for = 2, the method does not work well with the initial input λ = 0. In practice, one can also adopt the eigenvalue λ = ( + 1) − s(s − 1) for the nonrotating case (a = 0) as a simpler initial input value than λ STU , while in that case the computation time becomes about 1.5 times as long as the case with λ STU . However, the precision reaches, e.g., 20 digits so long as one requires the option PrecisionGoal → 20.

III. SCATTERING PROBLEM
In this section we focus on the radial solution and provide the exact solution for the scattering problem. In §III A we consider the asymptotic solution at the black hole and cosmological horizons, respectively. We shall see that the asymptotic solutions correspond to in/outgoing waves and are consistent with the asymptotic form of the exact solution in terms of the local Heun function obtained in §II. In §III B we exploit the asymptotic solution as the boundary condition, and write down the coefficients for in/outgoing waves in terms of the connection coefficients for the local Heun function.
For the following we omit the superscript (r) from the radial solution y (r) Ii,s for simplicity. Since we do not discuss the angular solution y (a) Ii,s below, no confusion should occur.

A. Asymptotic behavior
We can obtain the boundary condition by considering the asymptotic behavior of the radial equation at the black hole and cosmological horizons, for which the Schrödinger form is useful. We employ the tortoise coordinate r * defined in (13) as an independent variable, and transform the dependent variable as We can then rewrite the radial Teukolsky equation (19) in the Schrödinger form with the potential The potential depends on the spin only via s 2 and is, except s(1 − α) − λ, apparently. Actually, the combination s(1 − α) − λ is invariant under s → −s [40]. Hence, the potential (65) has a symmetry V * −s (r) = V s (r), where z * is a complex conjugate of z. This implies that Y s (r * ) and Y * −s (r * ) are two linearly independent solutions of the same differential equation (64). Therefore, if R s = ∆ −s/2 (r 2 + a 2 ) −1/2 Y s is a solution of the radial Teukolsky equation, ∆ −s R * −s = ∆ −s/2 (r 2 + a 2 ) −1/2 Y * −s is the solution linearly independent to R s . The potential asymptotically approaches a constant value where we denote r 1 = r + , r 2 = r + , and f h = f (r h ) for h = 1, 2. Consequently, the asymptotic behavior of two independent solutions is given by From (14), at the vicinity of the horizon r + or r + , the tortoise coordinate behaves as Using (68), we obtain Plugging (69) into (67) and multiplying ∆ −s/2 (r 2 +a 2 ) −1/2 , we obtain the asymptotic solutions of the radial Teukolsky equation where we have omitted proportional constants. One can check that for the SdS case the asymptotic solutions (70) are e iωr * and ∆ −s e −iωr * , respectively.

B. Scattered waves
In general, the asymptotic behavior of a general solution R s (r) is given by a linear combination of the two asymptotic solutions (70). For the scattering problem, we focus on two independent solutions R in (r) and R up (r) that satisfy the following asymptotic behaviors [47] The physical meaning is transparent once combined with the time-dependent part e −iωt . The "in" solution is defined by the boundary condition that there is no wave coming out from the black hole horizon. On the other hand, the "up" solution is defined by the boundary condition that there is no incoming wave from the cosmological horizon. Both boundary conditions are appropriate for the classical picture of the horizons. Combined with two other solutions defined by any two solutions among the four solutions (71)-(73) are linearly independent solutions for the same radial Teukolsky equation. For the scattering problem, we mainly use R in,s and R up,s . In the definition of R in,s in (71) and R up,s in (72) there are six coefficients. Not all the coefficients are independent. Clearly, one can omit the overall factors as the degrees of freedom for the normalization, but here we keep them for later convenience. On the other hand, we can derive relations between coefficients for R in,s and R up,s as follows. From (19), for a set of two linearly independent solutions R 1 , R 2 , it holds that Plugging in (R 1 , R 2 ) = (R in,s , R up,s ) and (R out,s , R up,s ), we obtain where Note that F * −s = F s holds. The ratios between the coefficients C s , D s determine the scattering problem and yield the S-matrix, reflection/transmission rate, and so on. Our aim in this section is thus to write down the coefficients C s , D s using the exact solution in terms of the local Heun function given in §II.
As we shall see below, the asymptotic behavior suggests that R in,s (r), R up,s (r) respectively corresponds to y 02,s (z), y 11,s (z), namely, R up,s (r) = D 11,s R 01,s (r) + D 12,s R 02,s (r), (r → r + ), where each R Ii,s is defined by (31) with the corresponding solution y Ii,s , with I = 0, 1 and i = 1, 2. Note that here we are not using any approximation but using the exact relations (53) and (56). R in,s , R up,s are given exactly by the local Heun functions at z = 0 and z = 1, and each two expressions coincide with each other for the region where two disks of convergence overlap. Using r − r h ∆(r)/∆ h for r → r h , we obtain z A∆(r), (z → 0; r → r + ), where With these relations and the asymptotic expansions (50) and (51), we see that the solutions (78), (79) indeed satisfy the boundary conditions given in (71), (72), respectively. Hence, we can express the coefficients C s , D s in (71) and (72) as For the scattering problem, the "squares" of the coefficients are important; these take the following form without A and A : In addition, from (75) and (76) we obtain the following relations: Here we implicitly assume that C 22,s and C * 21,−s are nonvanishing. If we consider C 22,s = 0 for instance, then we should go back to (58) and see D 11,s = 0.
To summarize, we have solved the scattering problem exactly in the sense that we have expressed the coefficients for the "in" and "up" solutions in terms of the connection coefficients between the local Heun function, which is the exact solution of the Teukolsky equation. Our calculation does not rely on any approximations such as the high-/low-frequency limit or slow-rotation limit. A specific example is the WKB or eikonal approximation for the high-frequency regime, which is commonly used in the literature. Such approaches with approximations are helpful to extract a simple intuitive picture and formulae for a limited setup. On the other hand, our exact formulation actually provides a simple expression without restriction of the parameter set. This allows us to use a simple and fast computation to understand black hole physics, which we shall explore in §IV.

IV. APPLICATIONS
In this section, we highlight several applications of the exact formulation of the scattering problem in §III. We discuss the quasinormal modes in §IV A, the S-matrix and cross section in §IV B, the reflection and absorption rates, greybody factor, and superradiant scattering in §IV C, and the Green function in §IV D. Our exact formulation serves as a simple and fast computational method with arbitrary high precision compared to the direct numerical integration of the Teukolsky equation.

A. Quasinormal modes
We can obtain QNM frequencies by requiring the regularity condition (62) on the angular part, as well as the boundary condition on the radial part with a purely ingoing wave at the black hole horizon r → r + and a purely outgoing wave at the cosmological horizon r → r + . Specifically, the condition on the radial part is given by C 22,s = 0, i.e., In parallel to the angular condition (62), the radial condition (95) also depends on ω and λ implicitly. For a nonrotating black hole, we can plug in the eigenvalue λ = ( + 1) − s(s − 1) and solve (95) only to obtain the QNM frequencies with a root-finding algorithm. For a rotating black hole, we obtain ω and λ by solving (62) and (95) simultaneously. The Wronskian is given by the exact solution in terms of the local Heun function, which in practice we can calculate by using the built-in function HeunG implemented in Mathematica 12.1 or later. This method gives us an arbitrary-precision arithmetic for the QNM frequencies. As already shown in [43] for the Kerr-de Sitter black hole, this method is quite fast, typically within O(1) second, and yields QNM frequencies that are consistent with the results in the literature. Therefore, we do not repeat the calculation of the QNM frequencies here. A caveat is that, to numerically find out the correct root, one needs to set an initial value sufficiently close to the root.
Let us note some technical details. To optimize the calculation, one can choose the evaluation point of the Wronskians either within the overlapping region of both disks of convergence for the local Heun functions at z = 0 and 1, or outside but still near the overlapping region. While the radius of convergence for y 02,s is always 1 and does not depend on the parameters of the KNdS geometry or scattered waves, this is not the case for y 11,s . Specifically, as one takes smaller ΛM 2 , the radius of convergence for y 11,s becomes smaller. In such a case we find that z = 0.9 is a convenient choice that yields a short computation time.
Let us note some differences from Leaver's method [42]. Leaver's method is one of the most successful algorithms to calculate the QNM frequencies, and it is implemented in the Kerr-de Sitter black hole in [41]. Both Leaver's method and the above method yield QNM frequencies with arbitrary high precision without any approximations, and there are the following qualitative differences. In Leaver's method, one solves three-term recurrence relations associated with the angular and radial equations in terms of infinite continued fractions. The boundary conditions determine the eigenvalues and the QNM frequencies implicitly as roots of two equations containing infinite continued fractions. One should then truncate the continued fractions appropriately, and use a root-finding algorithm. In principle this procedure allows one to obtain an analytic expansion formula. However, for high-precision computation one should take care with the convergence of the truncation. On the other hand, for the root finding procedure for (62) and (95), the ambiguities of the truncation do not appear. One can control the precision of the QNM frequencies solely by the precision of the calculation of the Wronskians. However, the Wronskians are calculated numerically and their analytic expressions are unclear. Therefore, the two arbitrary-precision arithmetics have complementary advantages.

B. S-matrix and cross section
By definition (71), the "in" solution R in,s is the solution satisfying the boundary condition with the purely ingoing boundary condition at the black hole horizon. We can then define the S-matrix S ,s (ω) as a ratio between the coefficients for the ingoing and outgoing waves at the cosmological horizon: Again, this expression reduces the number of Wronskians and minimizes the computational cost. This formula does not require numerical integration. We can obtain the S-matrix by calculating the ratio of the Wronskians of the local Heun function. Given the S-matrix, we can write down the differential cross sections and the scattering amplitudes, which are defined through an infinite series of the partial wave expansion. In practice, one needs to truncate the infinite series at some finite max . However, it is known that a naive truncation of the partial wave expansion introduces a numerical error, and hence special care is required [48]. The situation is analogous to the computation of the Coulomb scattering series [49]. Since our main goal in the present paper is to establish the analytic formulation of the wave scattering from black holes, here we do not address this issue further.

C. Reflection and absorption rates
We can express the conserved current of the scattered wave [50] in terms of the exact solution. We shall see below that the exact formulation provides a simple formula for the reflection rate and absorption (transmission) rate or the greybody factor. As we explained above, R s and ∆ −s R * −s are linearly independent solutions of the same differential equation (19). Plugging (R 1 , R 2 ) = (R in,s , ∆ −s R * in,−s ) and (R up,s , ∆ −s R * up,−s ) into (74) and evaluating it at r → r + and r → r + , we obtain These relations imply energy conservation [50]. While we do not specify the relative normalization between the s and −s solutions, the normalization degrees of freedom do not enter if we write down the energy conservation in the form The physical meaning is transparent. The first term on the right-hand side of (100) indicates the probability of the incoming wave being reflected by the black hole, whereas the second term means the probability of the incoming wave transmitting the effective potential and falling into the black hole. Therefore, the first and second terms yield the reflection rate R s and transmission rate T s , respectively. Note that, in the context of black hole scattering, T s is also called the absorption rate since transmission through the effective potential means absorption by the black hole. Similar logic also holds for each term on the right-hand side of (101). In particular, the second term on the righthand side of (101) is the greybody factor Γ s , which is the probability of the outgoing wave reaching the cosmological horizon. With (75) and (76), we can see that (100) and (101) are equivalent. Namely, we can rewrite them as where The relation (104) guarantees that the absorption rate T s coincides with the greybody factor Γ s . Using the exact solution, we can derive the following simple expression: where we have used (96) and (97). The absorption rate is then given by T s = 1 − R s . Alternatively, from (87) and (89), we obtain By virtue of (104), these formulae also allow us to calculate the greybody factor Γ s . In particular, for the scalar wave s = 0, the absorption rate can be written as Therefore, the absorption rate can be negative if F 0 = J(r + )/J(r + ) < 0 is satisfied. This is nothing but superradiant scattering. The condition J(r + )/J(r + ) < 0 can be rewritten as am + eQr + 1+α Here we stress that we have obtained these formulae exactly without any approximations, e.g., the high-/lowfrequency limit. Recalling that z ∞ and F s defined in (77) can be algebraically obtained, the only necessary calculation that one needs to perform to obtain (105) or (106) is the evaluation of the Wronskians between two local Heun functions at z = 0 and 1, which is achievable within the overlapping region of the two disks of convergence. Furthermore, compared to the calculation for the QNM frequencies in §IV A, the calculation for the reflection/absorption rate does not require an initial value close to the solution.
In Fig. 2, we present the reflection rate R s of the scalar wave with = 2, 4, 6 by the SdS black hole with ΛM 2 = 10 −3 as a function of M ω. We calculated R s for 0 ≤ M ω ≤ 1.7 with the sampling mesh size ∆(M ω) = 0.02. The solid curves are the results obtained by the exact formula (105). As a consistency check, we also calculated the reflection rate numerically using Mathematica. First, using NDSolve with the method "StiffnessSwitching", we numerically integrate the radial Teukolsky equation in the Schrödinger form (64) with the ingoing boundary condition near the BH horizon for Y s , which is the negative sign of (67). We then fit the behavior of the obtained wave function for the large-r * region, where the effective potential converges as (66), by the asymptotic solutions and read off the coefficients for the in/outgoing waves. Here we use the asymptotic solutions (67) in terms of the tortoise coordinate r * rather than r since the frequency of oscillation diverges in r space. Note that the numerical calculation is done with the default machine precision. For the sake of clarity, let us note the specifications of our computer for the numerical calculation. We use a Mac Pro with a 3 GHz, 8-core processor and the command ParallelTable is used to get the list of data. For all the numerical calculations for the reflection rate in the present paper (Figs. 2 and 3), we adopt the above method.
The results of the numerical integration are shown by dashed curves in Fig. 2 and are in good agreement with the results of the exact formula shown by solid curves. While analytic calculations known in the literature are valid under certain approximations such as high/low multipoles, our exact formula is based on the exact solution without approximation and hence can be used for a wider range of multipoles . The result in Fig. 2 is also consistent with physical intuition since partial waves with c are absorbed by black holes, where c = 3 √ 3M ω is the critical angular momentum. This is because partial waves with the impact parameter ∼ b c ≡ c /ω = 3 √ 3M are marginally scattered at the vicinity of the peak of the effective potential.
There are several differences between the exact formula (105) and the numerical calculation. First of all, the exact formula allows us to obtain the reflection rate with arbitrary high precision. One can easily improve the precision by requiring higher precision for the root-finding algorithm. On the other hand, for the numerical integration, it is more difficult to control the precision since one needs to take account of several processes to solve the differential equation. For the above calculations, the computation time for the exact formula is comparable to the numerical calculation. Specifically, to obtain each curve in Fig. 2, it takes about 10 sec for the numerical calculation without setting PrecisionGoal in Mathematica (the machine precision) and about 20 sec for the exact formula with PrecisionGoal → 15, respectively. For scalar waves (s = 0) scattered by the Kerr-de Sitter black hole (Q = 0, a/M = 0.9, ΛM 2 = 10 −3 ), we plot the reflection rate R s of m = ± modes with = 2 in Fig. 3. For the eigenvalue λ, instead of solving the angular part, we use the analytic expansion formula [14] since the error remains O(10 −1 )% for this setup, as shown in Fig. 1. We obtain the solid curves by using the exact formula (105), and the dashed curves by numerical integration. The sampling mesh size is ∆(M ω) = 0.01; i.e., we take 120 points in the range 0 ≤ M ω ≤ 1.2. We used PrecisionGoal → 15 for the exact formula, and the default machine precision for the numerical calculation. Using Table to list the data, the computation time to get each curve in Fig. 3 is 100 sec for numerical calculation and 240 sec for the exact formula. We see that the two results are in good agreement. The blue (red) curve depicts m = (m = − ), corresponding to the case where the angular momentum of the black hole and incident wave are (oppositely) aligned. The difference in the alignment causes the difference in the critical impact parameter b c ≡ c /ω at which a transition occurs from absorption to reflection.
In the right panel of Fig. 3, we show a closer look of the left panel to confirm the superradiant scattering for m = + . Indeed, we can see that the reflection rate exceeds unity, shown by the horizontal dashed line. For this parameter set, the condition (108) on the superradiant frequency reads 6.23 × 10 −4 < M ω < 6.25 × 10 −1 . The vertical dashed line corresponds to the upper bound of the superradiant freqeuncy, which is consistent with our calculations.

D. Green function
In this section, we construct the Green function for the wave scattering problem by a KNdS black hole in terms of the local Heun function. We choose the KNdS black hole as the origin of the spherical coordinate system (r, θ, ϕ) with the rotation axis at θ = 0. We assume a stationary point source, whose spatial location is denoted by x s = (r s , ϑ s , ϕ s ), and the observing point at x = (r, ϑ, ϕ), where ϑ is related to the polar angular variable θ of the spherical coordinates as ϑ = π/2 − θ. Therefore, ϑ = 0 is the equatorial plane of the KNdS black hole. The relationship between these points and the black hole is shown in Fig. 4. For the case where a spin-s wave is emitted by a stationary point source, the spatial part of the Green function G(x, x s ) can be expanded with the partial waves as where s S m (θ, ϕ) is the modified spin-weighted spheroidal harmonics due to the presence of the cosmological constant, which can be expressed by the local Heun function as given in Eq. (24). The differential equation that the radial partG(r, r s ) obeys can be derived from the master equation for the spin-s wave with Dirac's δ function as the source term as As we discussed in §III, the homogeneous equation (19) can be exactly solved in terms of the local Heun function, and satisfies the relation (74). Following the standard prescription, we can construct the Green function by using the two linearly independent solutions and the constant (74), which is given bỹ where Θ(r) is the unit step function. Here, we have chosen R in and R up given in (78) and (79), respectively, as a suitable pair of independent solutions to the radial Teukolsky equation by considering the boundary condition of the wave scattering problem by black holes. Note that we omit the subscript s of R in,s and R up,s for the spin-s wave for simplicity. The denominator ∆ s+1 W r [R in , R up ] is the constant given in (74), which should be evaluated at some r between r + < r < r + , or 0 < z < 1. Plugging (111) into (109), we obtain the spatial part of the Green function as (112) In particular, as mentioned in §II B 1, for the scattering of scalar waves by a nonrotating black hole, the angular solution is given by the spherical harmonics. In this case we can use the addition theorem for the spherical harmonics: where the variable γ represents the angle between source and observer, which is defined by cos γ = cos θ cos θ s + sin θ sin θ s cos (ϕ − ϕ s ). We then arrive at Let us check whether this formula reproduces the formula derived with asymptotic solutions of scalar fields (s = 0) in Schwarzschild spacetime [51]. We use the rescaled radial function Y s (63) with a = 0 and s = 0, and the tortoise coordinate r * . The relationship between the Wronskian for R s with respect to r and that for Y s with respect to r * is then given by Plugging this into the Green function (114) yields which amounts to the Green function derived in [51]. Let us highlight several differences between the exact Green function (112) and the analysis performed in [51]. First, in [51] the Green function (116) was evaluated by substituting the asymptotic forms of the radial function corresponding to (71) and (72). However, in that case the sum over the partial waves does not converge due to 1/r behavior of the gravitational potential. This issue originates from the use of the asymptotic solutions. Indeed, in [51], the convergence issue was circumvented by adding a finite-distance correction to the asymptotic solutions, which play the role of regulator for the partial wave sum. In contrast, for the exact Green function (112), there is no convergence issue intrinsically. This is because the exact Green function does not rely on the asymptotic solution or approximations, and inherits finite-distance effects at nonlinear order.
Second, in [51], it is assumed that both the source and observer are located at a sufficiently distant r/M, r s /M 1, but are finite points so that one can substitute the asymptotic solution with the correction term into the Green function. Furthermore, the wavelength of the scalar wave is restricted to the short-wavelength case M ω 1 to evaluate the phase shift within the WKB approximation. A small deflection angle (ϑ ∼ 0, ϕ ∼ 0) was additionally assumed to obtain a simple formula. In contrast, in our formulation there is no approximation and no restriction on the scattered wave and the configuration of the source and the observer since the radial functions R in/up represented in terms of the local Heun function are the exact solution to the radial Teukolsky equation. Moreover, our Green function (109) applies to a more general case, i.e., spin 0, 1 2 , 1, 3 2 , 2 massless fields on the Kerr-de Sitter background and those for spin 0, 1 2 massless fields on the KNdS background. Therefore, the Green function (112) is the most general exact formula for wave scattering by a KNdS black hole.
In Fig. 5, we present the power spectrum, i.e., the absolute square of the Green function (114) measured at r = 20M for forward scattering by the SdS black hole with ΛM 2 = 10 −3 of scalar waves emitted from the source located at (r s , ϑ s , ϕ s ) = (6M, 0, π). To obtain the power spectrum, we evaluate R in/up in the following two ways: First, we employ the exact solution (78), (79) in terms of the local Heun function with PrecisionGoal → 25, which is shown as solid red curves in Fig. 5. On the other hand, the power spectrum obtained by the numerical integration is shown by dashed blue curves. For the numerical integration, here we improve a similar calculation performed in [52]. In [52], the WKB approximation was partially employed, but here we do not use the approximation. Here, to obtain R in/up , we numerically integrate the differential equations (64) and (13) with the boundary conditions (71) and (72). We choose the location to impose the boundary conditions sufficiently close to each horizon, and confirm that the results are almost unaffected by some change of the location. Specifically, we start the numerical integration with the purely ingoing boundary condition as the initial condition at a nearby point of the black hole horizon r i /M = r + /M + 10 −6 ∼ 2.0027 and solve the radial equation towards the source point r s /M = 6, whereas, for R up , the radial equation is solved with a purely outgoing boundary condition from a point near the de Sitter horizon r f /M = r + /M − 0.74 ∼ 53 to r s /M = 6. Then, substituting these solutions into (114), the Green function is obtained. To get both results, we truncate the partial wave sum at max = 8M ω + 6, which we find yields a good convergence. As shown in Fig. 5, the two curves are in good agreement. The computational time to obtain the curves in Fig. 5 is about 30 min for the exact formula and 1 min for the numerical integration. While the exact formula takes longer, note that the numerical integration here does not have high precision. Actually, the numerical result matches the exact result up to 3 digits only. To improve the numerical result, one needs to choose the location for the boundary condition closer to each horizon, and to require higher precision for the root-finding algorithm and the differential equation solver. For instance, if we take r f /M = r + /M − 10 −6 , the numerical result matches the exact result up to 5 digits, and the computational time is 3 min in this case. It would thus be fair to say that the exact formula serves a simple calculation method with high precision. It allows an arbitrary high-precision calculation and it is easier to control the precision without numerically solving the differential equation.
For exact forward scattering with (r, ϑ, ϕ) = (20M, 0, 0) in the left panel of Fig. 5, the behavior that |G| 2 increases linearly stems from the property of the caustics at the forward position of the present scattering problem, which will diverge for M ω → ∞. The period of oscillation on the linear growth reflects the scale of the peak of the effective potential. This corresponds to the position of the unstable circular photon orbit in the geometrical optics limit, and is evaluated as M ∆ω ∼ 1/(3 √ 3) ∼ 0.2 for ΛM 2 1. On the other hand, for the case of the slightly off forward scattering with (r, ϑ, ϕ) = (20M, 0, π/10) in the right panel of Fig. 5, there is one more oscillating scale with a longer period. This originates from the breaking of the symmetry of the relative relation of the source-black hole-observer positions, which causes interference due to the difference of light ray paths in the limit of the geometrical optics. As another demonstration, we present the angular dependence of the absolute square of the scattered scalar wave for fixed frequency M ω = 7, 4, 1 in Fig. 6. As expected, it shows a peak at ϕ = 0 and decays with oscillations depending on the fixed frequency. We see that our exact formula is valid for a wide range of the azimuthal angle. We have provided several examples of scalar wave scattering by the SdS black hole with a small cosmological constant and compared the results with previous works. Since the main goal of the present paper is to establish the formulation, we have avoided to present too many specific calculations. However, our formula (112) is quite general and applies to the wave scattering of the spin-s field from the KNdS black hole. We will investigate the details of several observables in wave optical gravitational lensing for a more general case in a future work.

V. CONCLUSION
In this paper we have established the exact formulation for the wave scattering problem by the KNdS black hole. We consider the propagation of a massless field with spin and charge on the KNdS background. The Teukolsky equations for spin 0, 1 2 , 1, 3 2 , 2 fields on the Kerr-de Sitter background and those for spin 0, 1 2 fields on the KNdS are separable and take the unified form. Here, the spin 0 field corresponds to a scalar field conformally coupled to gravity. Transforming the angular and radial Teukolsky equations into Heun equations, we can write down the exact solution in terms of the local Heun functions at regular singular points. For the angular solution, we can impose the regularity condition by requiring the linear dependence of the local Heun functions. For the radial equation, with the appropriate transformation, we can respectively map the black hole horizon r = r + and the cosmological horizon r = r + into z = 0 and 1, and discuss the scattering problem within the range 0 ≤ z ≤ 1. For this setup, there exists an overlapping region of the two disks of convergence of the local Heun functions at z = 0 and 1, and we can discuss the scattering problem. We can write down the "in" and "up" solutions, which satisfy certain boundary conditions, in a fully analytic way without any approximations. We have expressed the coefficients for the asymptotic in/outgoing waves exactly in terms of the connection coefficients for the local Heun functions at z = 0 and 1, which are given as the ratio of the Wronskians of the local Heun function. Once the coefficients are obtained exactly, we can write down various important quantities for black hole scattering exactly.
We have highlighted several applications of our exact formulation. It has already been shown in [43] that the local Heun function is a powerful tool to calculate the QNM frequencies for the Kerr-de Sitter black hole with arbitrary high precision. Given a sufficiently close initial value as an input, one can obtain the QNM frequencies very quickly. We have generalized this result to the KNdS geometry and provided the arbitrary-precision arithmetic for the KNdS QNM frequencies for spin 0, 1 2 massless fields for the first time. We can also write down the S-matrix, with which the differential cross sections and the scattering amplitudes can be written down. Further, we have explored the conserved current for the scattering problem in terms of the exact solution, and derived simple formulae for the reflection/absorption rate and the greybody factor (see also [53] for a recent study on the greybody factor and Hawking radiation for the Kerr-de Sitter black hole in this context). We have checked the consistency between the results obtained by our exact formula and numerical integration, and clarified the efficiency of our formula for the reflection rate in comparison with the numerical integration. Finally, we have constructed the Green function for the wave scattering from the KNdS black hole. We have calculated the power spectrum as the absolute square of the Green function to see the frequency dependence of the forward and slightly off forward scattering, as well as the angular dependence for the fixed frequency waves. They are consistent with the numerical results as well as the previous results in the literature, where some approximations were employed.
Our exact formulation of the wave scattering from the KNdS black hole provides simple and practical formulae, which are arbitrary-precision arithmetics. Unlike known (semi-)analytic calculations in the literature, we do not use any approximations. There is no restriction on parameters such as the frequency of the scattered waves, or the relative locations of the source of the waves, black hole, and observer. The exact formulae predict the scattering problem with arbitrary high accuracy. While we have highlighted several specific applications, it would be intriguing to apply our formulation to more general cases or other observables. We leave these topics for future work.