{'fraction_non_alphanumeric': 0.0501466213792958, 'fraction_numerical': 0.0177211451235206, 'mean_word_length': 4.181679055531264, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Eigenvalues and eigenfunctions of spin-weighted spheroidal harmonics in four and higher dimensions

Spin-weighted spheroidal harmonics are useful in a variety of physical situations, including light scattering, nuclear modeling, signal processing, electromagnetic wave propagation, black hole perturbation theory in four and higher dimensions, quantum field theory in curved space-time and studies of D-branes. We first review analytic and numerical calculations of their eigenvalues and eigenfunctions in four dimensions, filling gaps in the existing literature when necessary. Then we compute the angular dependence of the spin-weighted spheroidal harmonics corresponding to slowly-damped quasinormal mode frequencies of the Kerr black hole, providing numerical tables and approximate formulas for their scalar products. Finally we present an exhaustive analytic and numerical study of scalar spheroidal harmonics in (n+4) dimensions.

have innumerable applications in physics (see eg. [13] for an application to the analysis of anisotropies in the cosmic microwave background). The ordinary spherical harmonics are spin-weighted spherical harmonics with s = 0.
Describing the angular dependence of a spin-s perturbing field in a Kerr black hole background, SWSHs and their eigenvalues find application in many theoretical studies of black hole physics.
An important application concerns quasinormal modes of Kerr black holes [14]. The damped oscillation frequencies of Kerr black holes, ω, are uniquely determined by the black hole's mass M and specific angular momentum a. The calculation of quasinormal frequencies reduces to the solution of a coupled system of differential equations. One equation belongs to the class of generalized spheroidal wave equations, and describes the radial dependence of the perturbations; the other is a SWSH with c = aω, describing the angular dependence [15,16,17,18,19,20]. For this reason, a detailed knowledge of the numerical and analytic properties of SWSHs is necessary to compute black hole quasinormal modes. Recent conjectures suggest a relation between highly damped quasinormal modes and black hole area quantization [21,22]. Verifying these conjectures for Kerr black holes calls for a calculation of the eigenvalues s A lm corresponding to values of c with imaginary part much larger than the real part, |c I | ≫ |c R |. More phenomenologically, an important astrophysical problem is the determination of black hole parameters from gravitational wave observations [23,24]. Gravitational waves emitted by the oscillations of a rotating black hole can be described by a superposition of quasinormal modes. The Laser Interferometer Space Antenna (LISA) has the potential to detect these waves with large signal-to-noise ratio [25]. Given a detection, we would like to extract as much information as possible about the source. Ideally, we would like to determine M , a, the source location and the black hole's spin orientation from the observed waveform. An investigation of this issue requires the calculation of "scalar products" between different quasinormal modes, and in particular between the SWSHs describing their angular dependence [24].
SWSHs also find application in quantum field theory in curved space-time. This theory is plagued with ultraviolet divergences: in particular, the expectation value of the stress-energy tensor, which is the crucial quantity that generates space-time curvature via the semiclassical Einstein's equation, suffers from ultraviolet divergences and must be renormalized. This is usually a very difficult endeavour, particularly in the case of non-spherically symmetric space-times, which have recently been studied within this framework in [26]. An understanding of the high-frequency behaviour (both of the radial and angular parts) of the matter fields is essential in order to renormalize the various quantities in the theory [27].
Many theoretical scenarios rely on our Universe being (4+n)-dimensional, with the n extra dimensions compactified on some small scale [28]. The study of higher-dimensional black hole geometries, such as the Myers-Perry solution [29], requires the introduction of higher-dimensional spheroidal harmonics (HSHs), that we shall define in Sec. III [30,31,32,33]. TeV-scale gravity scenarios allow for the possibility that rotating, higher-dimensional mini-black holes are produced in particle accelerators such as CERN's Large Hadron Collider. These higher-dimensional, rotating black holes should evaporate by emission of Hawking radiation [30,34], and both SWSHs and HSHs are useful for a quantitative study of the evaporation process on the (3 + 1)-brane and in the (4 + n)-bulk, respectively. HSHs also show up in recent studies of smooth geometries in the D1-D5 system [35]. Separating the Klein-Gordon equation leads again to five-dimensional spheroidal harmonics (even though they were not identified as such in [35]).
There are no in-depth investigations of HSHs to date. On the other hand, the properties of SWSHs were (at least partially) investigated as soon as they were introduced. Press and Teukolsky [36] provided a polynomial fit in c of the eigenvalues s A lm , which is valid up to c ∼ 3. A formal perturbation expansion in powers of c was carried out by Fackerell and Crossman [37] (see also [38], where some typos were corrected). An expansion for large real values of c has long been known for scalar spheroidal harmonics [1,6,7]. For general spin-s fields, the first attempt to find a large-c expansion was made by Breuer, Ryan and Waller [39,40]. Their analysis, which was partially flawed and incomplete, has recently been revisited and corrected by Casals and Ottewill [27]. An expansion of the eigenvalue for large, pure-imaginary c is known for scalar spheroidal harmonics [1,6,7], and was partially studied in [19,39] for general spin-s fields. However, as we already mentioned some of the studies for general spin-s were flawed, and we still lack a unified picture of the situation.
The main purpose of this paper is to provide a complete analysis of general (four-dimensional) SWSHs, and to extend this understanding to HSHs. The plan of the paper is as follows.
Sec. II contains analytic and numerical results for four-dimensional SWSHs. In Sec. II A we define SWSHs and present a series solution first obtained by Leaver [15]. In Sec. II B we explain how this series solution provides us with a simple numerical algorithm to compute eigenvalues and eigenfunctions, and discuss alternative methods we used to check our results. In Sec. II C we give the series expansion of the eigenvalues for small c. In Sec. II D and II E we present analytic results for large values of the argument in the oblate and prolate cases, respectively, supplementing these analytic expansions by numerical calculations. Sec. II F contains a numerical calculation of the prolate eigenfunctions. In Sec. II G we compute numerically the SWSHs at the slowly damped quasinormal frequencies of a Kerr black hole, providing numerical tables and a simple analytic approximation of their scalar products. While Secs. II A, II B, II C and II D are mostly review material, the results in Secs. II E, II F and II G are new.
Sec. III is devoted to analytic and numerical calculations of the eigenvalues and eigenfunctions for HSHs. In Sec. III A we review a series representation due to Ida et al. [31], illustrating the corresponding computational algorithm for the eigenvalues and eigenfunctions. The results in the following Sections are entirely new. In Sec. III B we give a series expansion of the eigenvalues for small c. Secs. III C and III D present series expansions for large values of the argument in the oblate and prolate cases, respectively, and back them up by numerical calculations. We conclude with a summary of relevant results and open problems.
Other original results are contained in the Appendices. In Appendix A we give an exact analytic solution of the HSH equation for c = 0, and use it to determine the number of zeros of the eigenfunctions in the region of physical interest. Appendix B provides a simple analytic solution of the HSH equation for special values of the parameters.

A. Series solution
Using the Kinnersley tetrad and Boyer-Lindquist coordinates, the angular equation definining SWSHs results from the separation of the equations describing propagation of a spin-s field in the Kerr background [11]: where x ≡ cos θ and θ is the Boyer-Lindquist polar angle. The angular separation constant s A lm and the SWSHs s S lm are, in general, complex. They take on real values only in the oblate case (c ∈ R) or, alternatively, in the prolate case (c = ic I pure-imaginary) with s = 0. In the limit c → 0 the angular separation constant can be determined analytically: Some useful symmetry properties hold (see eg. [15]): (i) Given eigenvalues for (say) positive m, those for negative m are readily obtained by complex conjugation: (ii) Given eigenvalues for (say) negative s, those for positive s are given by Exploiting these symmetries, in our numerical calculations we only consider s ≤ 0 and m ≥ 0.
(iii) Let us define ρ ≡ ic. If ρ and −s A lm corresponds to a solution for given (s, l, m), then another solution can be obtained by the following replacements: m → −m, ρ → ρ * , −s A lm → −s A * l−m . (iv) In the prolate case, if s S lm is a solution with eigenvalue s A lm and c I > 0, then s S * lm is a solution with eigenvalue s A * lm and c I < 0. Leaver found the following series solution for the angular eigenfunctions [15]: where k ± ≡ |m ± s|/2. The expansion coefficients a p are obtained from the three term recursion relation with α p = −2(p + 1)(p + 2k − + 1) , Given a (generally complex) argument c, the separation constant s A lm can be obtained solving numerically the continued fraction equation 9) or any of its inversions [15].

B. Numerical calculation of the eigenvalues and eigenfunctions
Flammer's classical reference [1] is largely dedicated to tabulating eigenvalues and eigenfunctions of scalar spheroidal harmonics. The main limitation of Flammer's impressive work is that his tables only deal with the prolate (pure imaginary frequency) and oblate (pure real frequency) cases. The general case of complex frequencies (which is of interest, for example, in the calculation of the angular dependence of quasinormal modes) is not covered. A pioneering numerical work in this sense was carried out by Oguchi [7], who computed angular eigenvalues for complex values of c and s = 0. Quite recently this topic received more attention. Li et al. published a useful review of numerical methods to compute eigenvalues and eigenfunctions for s = 0 [8]. Falloon et al. developed a Mathematica notebook to compute s = 0 harmonics for general complex values of the frequency [9]. Finally, Barrowes et al. provided a rapid and accurate method to calculate the prolate and oblate scalar spheroidal wave functions and their eigenvalues for complex frequencies in the limit |c| → ∞ [10].
In comparison, numerical calculations of the eigenfunctions and eigenvalues for general spin s are scarce. To make things worse, until recently the few analytic predictions were contradictory (the situation for large real frequencies has finally been clarified in [27], while large pure imaginary frequencies were considered in [19]).
Leaver's solution gives a simple and practical algorithm for the numerical calculation of eigenvalues and eigenfunctions. Start from the known analytic eigenvalue for c = 0, Eq. (2.2). Use this as an initial guess, increase the value of c and solve numerically Eq. (2.9) to get the eigenvalue for c = 0. Once the eigenvalue is known, compute any number of series coefficients a p using the recursion relation, and plug them into the series solution (2.5) to get the corresponding eigenfunction to any required precision. In our numerical calculations we truncate the series at p = p max = 1000 (but usually an excellent approximation to the eigenfunction can be obtained keeping only ∼ 10 terms).
This algorithm only determines the eigenfunction up to a normalization constant, which can easily be fixed by imposing the normalization condition (2.10) We independently checked our numerical results for the eigenvalues and eigenfunctions using different methods. These methods are variants of those described in detail in [27], where they were used to deal with oblate SWSHs, so we only describe them briefly. To check the eigenvalues from Leaver's continued fractions we adapt the procedure derived by Sasaki and Nakamura [41] for s = −2 and c ∈ R. The idea is to approximate the angular differential equation by a difference equation, which can be written in matrix form. The eigenvalue is then obtained by imposing the determinant of this (tridiagonal) matrix to be zero. The Sasaki-Nakamura method faces a numerical difficulty as the discretization grid is refined: the value of the determinant becomes very large, making it increasingly difficult to obtain the eigenfunctions. Therefore, in order to find the eigenfunctions we complement the Sasaki-Nakamura analysis by an adaptation of the shooting method introduced in [42] for the special case of (scalar) spheroidal harmonics with real frequency. We start the integration of the angular equation with arbitrary values for (i) the eigenvalue, (ii) the angular function at x = −1. The eigenvalue results from imposing the derivative of s S lm near x = +1, as found by numerical integration, to agree with a certain analytic approximation of the derivative, coming from a power series expansion around x = +1. In [27], for large real frequency it was sometimes necessary to look for an extreme (rather than a zero) of either the determinant or the difference of derivatives near x = +1. That is because in the oblate case the eigenvalues for some modes "pair up" for large frequency. This degeneracy does not occur for large and pure-imaginary frequency (prolate case). On the downside, when using the "shooting" method we must now find the zero of a complex function. We locate this zero using Newton's method [42], which is globally convergent.
As a first example of the application of Leaver's method, in Fig. 1 we show some prolate and oblate eigenvalues for s = 0. Progressively increasing the inversion index of the continued fraction [15] we can track the eigenvalues (which in this special case are real) up to quite large values of c. Our numerical code passed a number of tests. We first verified that our results are in agreement with Tables 10-12 of [1] in the prolate case and Tables 130-132 of [1] in the oblate case. For complex c we were able to reproduce Table 2 in [7]. As a final sanity check, we verified in a few representative cases that the Mathematica notebooks presented in [8,9] are also in agreement with our code. Notice that prolate and oblate eigenvalues coincide when c → 0, but their behavior for large frequencies is markedly different. Prolate eigenvalues grow linearly with |c I | for large |c I |, while oblate eigenvalues are proportional to −c 2 for large c. In the following Section we derive analytically asymptotic expansions in the limit of small c, large real c and large pure-imaginary c.

C. Small-c expansion
The analytic properties of the eigenvalues and eigenfunctions of the SWSH equation (2.1) have been studied by many authors [1,6,19,27,36,37,38,39,40]. Here and in the following Sections we present for the first time a compact and complete summary of the analytic results for the eigenvalues in various asymptotic regimes.
The small-c expansion was worked out many years ago for the scalar case s = 0 (see eg. [1]). The general s case was considered by Press and Teukolsky [36], Fackerell and Crossman [37] and Seidel [38]. Here we present an alternative and simpler derivation of their results, which can easily be extended to the higher-dimensional case (see Sec. III B). Our method is directly based on the continued fraction representation. For c = 0, the eigenvalue s A lm is explicitly determined from the requirement that the series expansion has a finite number of terms, since otherwise it is divergent [15]. Taking c = 0 in (2.9) and imposing β r = 0 for some integer r ≥ 0 we obtain s A lm = (r + |m|)(r + |m| + 1) − s(s + 1) . (2.11) Appendix A gives an alternative derivation of this result in a more general framework. Setting r = l − |m| we recover the result (2.2); furthermore, the condition r ≥ 0 implies the well-known constraint on the angular quantum numbers, l ≥ |m|.

D. Large and real c (oblate case)
For s = 0, analytic power-series expansions for large (pure real and pure imaginary) values of c have long been known [1,6]. In this Subsection and in the next we review and extend known analytic results for general values of the spin-weight s. Large-c expansions of the eigenvalues are usually based on the following idea. First, an asymptotic approximation of the solution is found. An asymptotic expansion of the eigenvalue is then obtained by matching the number of zeros of the asymptotic solution with the (known) number of zeros of the exact solution (as determined, eg., in [40]). This matching relies on the solution being analytic in c.
Breuer, Ryan and Waller [40] perform an asymptotic analysis for large, real c (and fixed m) of the SWSH equation. Their analysis is based on similar analyses in the s = 0 limit [1,6] and on a partially-flawed previous study for general spin [39]. It yields the following asymptotic expansion for the separation constant: Unfortunately, Ref. [40] left the parameter s q lm undetermined. The recent work by Casals and Ottewill [27] determines s q lm through a complete asymptotic analysis for large, real c (and fixed m). First they find an asymptotic expansion for the angular solution valid everywhere in x ∈ [−1, +1]. Then they determine s q lm by imposing the asymptotic angular solution to have the right number of zeros. This could not be done successfully by Breuer, Ryan and Waller as their asymptotic angular solution was not valid near the origin x = 0, where the solution may have one zero. The value of the parameter s q lm obtained in [27] is: where s l m ≡ max(|m|, |s|) + (|m + s| − |m − s|)/2 + s, and the variable gives the number of zeros of the asymptotic angular solution at x = 0. In Fig. 2 we use Leaver's method to give numerical support to the analytic results, computing selected oblate eigenvalues for s = 0. Fig. 2 only shows the lowest radiatable multipoles (that is, the lowest values of l such that l ≥ |s|), which are physically the most relevant. Numerical results fully confirm the analytic expectations. The imaginary part of the eigenvalues is always zero within the numerical accuracy. The small-c and large-c behaviors are both in agreement with the analytic predictions. As an additional code check we also verified that the numerically computed eigenvalues satisfy the symmetry relation (2.4).

E. Large and pure-imaginary c (prolate case)
In the scalar case s = 0, the asymptotic analysis of the eigenvalues for large and pure-imaginary c presents no particular difficulties, and can be found in [1,6]. A very accurate approximation of the numerical eigenvalues in the left panel of Fig. 1, which is valid within a few percent down to c ≃ 2, is: where L = l − |m| is the number of zeros of the scalar harmonics in the interval −1 < x < 1. An expression accurate to O(1/|c I | 6 ) can be found in Eq. (8.1.11) of [1].
Ref. [39] provides, to our knowledge, the first attempt to generalize Eq. (2.21) to general spin s. Unfortunately the number of zeros of the SWSH given on page 115 of [39] is wrong, and was later corrected by Breuer et al. (Theorem 4.1 in [40]). Therefore, although the derivation in [39] is essentially right, the expression for the leading term is not.
Breuer's analysis for general spin was corrected in [19], but two open issues remain. In the scalar case, the asymptotic solution found in [1,6] is a valid approximation to the spheroidal harmonics only in a region far from the end-points x = ±1. We expect the function, which is real, to have most zeros far from the end-points; still, in principle we cannot discard the possibility that the function does have a zero near the end-points. As mentioned in the beginning of Sec. II D, the omission of a zero would lead to the wrong asymptotic expansion for the eigenvalue (this is indeed what happened in [40] in the case of large, real c and general spin s).
For general spin not only we run into the same problem, but we have the additional complication that the solution is complex. However, in [19,39] the number of zeros of the exact, real solution for real c was matched with the number of zeros of the real part of the asymptotic (complex) solution for pure-imaginary c. As indicated in [1,6] for the case s = 0, the number of zeros of the SWSH is the same for c 2 ∈ (−∞, +∞), but this is only true for real functions. Numerical calculations (see Sec. II F below, in particular Figs. 5 and 6) show that the number of zeros of the real part and of the modulus of the SWSHs are both functions of c (when c is pure-imaginary). Therefore it is not completely clear why the asymptotic analysis in [19,39] for the eigenvalue should be correct.
Below we provide analytic and numerical arguments in support of the result obtained in [19,39]. We derive an asymptotic solution for general spin-weight s which is a valid approximation to the SWSH near the end-points. This "outer" solution does not contain any zeros in the case s = 0, and this explains why the asymptotic behaviour of the eigenvalue is correct, at least for scalar perturbations. We could not generalize this argument to electromagnetic and gravitational perturbations, but we give numerical evidence that, in all cases we considered, the asymptotic behaviour of the eigenvalue given in [19] for |s| = 1 and |s| = 2 is correct.
In order to find an asymptotic solution we define a new angular wavefunction s y lm (x) via and change the independent variable setting u ≡ 2|c I |x. Substituting into Eq. (2.1) we get When |c I | → ∞ and in the region 1/|c I | ≪ |x| ≪ 1, the equation satisfied by the real part of s y lm becomes a parabolic cylinder equation. Therefore the real part of the "inner solution" s y inn lm , which is a valid approximation in this inner region, is a parabolic cylinder function D L with L zeros: ℜ s y inn This is essentially how the inner solution was found in [1,6,19,39]. Note, however, that the slightly different change of variable (2.22) means that the inner solution s S inn lm (x) we have obtained is always regular at x = ±1, whereas the one given in [19,39] is not.
It is known [43] that the zeros of the parabolic cylinder function D L 2|c I |x lie within |x| < (2L + 1)/|c I |, so they are all within the inner region. We can determine the value of L by equating the number of zeros of the parabolic cylinder function to the number of zeros of the real part of the exact solution for pure-imaginary c = ic I (i.e., the real part of the prolate SWSH) in the region 1/|c I | ≪ |x| ≪ 1. Unfortunately, the latter number is not known. In [19,39] the number of zeros of the parabolic cylinder function is instead equated to the number of zeros of the exact solution for real c (i.e., the oblate SWSH) in the whole region x ∈ (−1, +1), as given in [40]. The value of L so determined turns out to be Higher order corrections in the asymptotic expansion can be obtained as indicated in [1]. We checked Eqs. (2.24) and (2.25) numerically, solving Leaver's angular continued fraction for s A lm as a function of the complex parameter c. Now we take one more step towards understanding the asymptotic behaviour given by Eqs. (2.24) and (2.25) providing an "outer" asymptotic solution, which is a valid approximation to the SWSH in a region far from the origin. A WKB-type expansion similar to the one in [27] gives the result: where s a lm and s b lm are constants of integration. This solution is valid when |x| >> 1/ |c I | and |x| ∼ 1 − O(|c I | ǫ ), with −1 < ǫ ≤ 0. Matching the outer and inner solutions so that they coincide in their overlap region, we get the following final expression for the outer solution (valid near x = ±1): (2.27) This solution formally diverges at x = ±1, and actually it is not valid at the end-points, but only in a region arbitrarily close to the end points. In such a region the exponential decay with |c I | takes care of the diverging factor 1 − x 2 −1/4 . It can easily be checked that: (1) the outer solution has no zeros when s = 0; (2) its real part has no zeros within x ∈ (−1, +1) when |s| = 1; (3) its real part has two zeros located at x = ±1/ √ 2 when |s| = 2. We finally have a complete account of scalar perturbations (s = 0). The outer solution has no zeros, so all zeros of the asymptotic solution are provided by the inner solution. Equating the number of zeros of the parabolic cylinder function to that of the oblate SWSH is therefore valid in this case, which is why Eqs. (2.24) and (2.25) are correct for s = 0. For s = 0 we are one step closer to understanding the asymptotic behaviour. Now we know all the zeros of the real part of the asymptotic solution, and we "only" need to find out the number of zeros of the real part of the SWSH for pure-imaginary frequency. We have numerically checked in many cases that Eqs. (2.24) and (2.25) give the right behaviour, and we believe them to be correct even for s = 0, but only an analytic proof can settle the issue once and for all.   3 shows prolate eigenvalues s A lm for s = 0 and a few selected values of (l, m), providing numerical support to our analytic conclusions. The angular eigenvalue s A lm is now complex (unless m = 0). We can still limit our calculations to positive m's because of the symmetry property (2.3). From Fig. 3 we see that in the prolate case the real part dominates as c I → ∞. The numerical results for any s are consistent with the linear growth predicted by Eqs. (2.24) and (2.25). A peculiar feature can be seen in the real part of −2 A 22 , which seems to "bend down" at c I ≃ 3. We verified that the presence of this "bending" is not a numerical artifact using both Leaver's method and the shooting method. There is no discontinuity at the bending location, just a smooth change of slope. Ref. [18] found that, in the intermediate damping regime, the real part of the Kerr quasinormal frequencies with l = m = 2 tends to ω R = 2Ω, Ω being the angular velocity of the horizon, unlike the frequency of modes with any other value of l and m. We suspect this anomalous behavior of quasinormal modes with l = m = 2 and s = −2 in the intermediate damping regime could be related with the observed bending of the angular eigenvalue. Finally, Fig. 4 is a computational exercise to illustrate in a specific case the transition from the prolate to the oblate regime. We set c = |c|e iθ and compute the real part of −2 A 22 along trajectories in the complex plane with fixed values of θ ∈ [0, π/2]. A phase angle θ = 0 corresponds to the oblate case (Fig. 2) and a phase angle θ = π/2 to the prolate case (Fig. 3). For intermediate values of θ the function ℜ −2 A 22 (|c|) smoothly deforms from the prolate to the oblate character. Notice also that the bending of the l = m = 2 prolate eigenvalue we see at |c| ≃ 3 (see also Fig. 3) disappears immediately as θ becomes nonzero. In our opinion, this provides some circumstantial evidence that the bending is not due to a branch cut. We refer to [7,10] for an extensive discussion of branch cuts and their effect on the prolate/oblate nature of the eigenvalues when s = 0.

F. Numerical calculation of the prolate eigenfunctions
Once we have obtained the eigenvalues it is a simple matter to compute the eigenfunctions for given values of s, l, m and c: we just need to compute the coefficients a p from the recursion relation (2.6)-(2.8) and plug them into Leaver's series solution, Eq. (2.5). The eigenfunctions for s = 0 are well known [1], and so are the eigenfunctions for oblate harmonics (real c) with general spin s [27]. Therefore here we concentrate on prolate eigenfunctions with different values of c I . In black hole perturbation theory, according to Leaver's convention on the Fourier transform [15] stable perturbations correspond to quasinormal frequencies with ω I < 0. For this reason, below we pick c I < 0. Using the symmetry property (iv) of Sec. II A, eigenvalues and eigenfunctions for c I > 0 are trivially obtained by complex conjugation. The left panels refer to the simple case c I = 0. In this limit the eigenfunctions reduce to spin-weighted (s = −2) spherical harmonics, which can be thought of as being prolate and oblate at the same time [12]. The "spikes" in this logarithmic plot are zeros of the eigenfunctions, and their number is in agreement with Eq. (2.25) for c I = 0. However, as anticipated, Fig. 5 shows that (for s = 0) as soon as c I = 0 the number of zeros in the region −1 < x < 1 changes. Fig. 5 only gives the modulus of the eigenfunctions, to illustrate the point that prolate harmonics with spin-weight s = −2 have no zero for −1 < x < 1 when c I = 0. In Fig. 6   G. Scalar products of the eigenfunctions at the Kerr quasinormal frequencies As a useful application to general relativity of the numerical algorithm described in Sec. II B, here we study the angular eigenfunctions corresponding to different quasinormal mode frequencies and different values of a/M , for l = 2. To our knowledge, these eigenfunctions have never been shown before in the literature. The procedure we use is by now standard and can be found in many papers [15,17,18,19,20], so here we only sketch the basic idea. In the Kerr space-time, linear gravitational perturbations are described by a pair of coupled differential equations: one for the angular part of the perturbations, and the other for the radial part. The radial equation resulting from the separation of the field equations in the Kerr background is given (eg.) in [11,15]. The angular equation is a SWSH equation with c = aω. Boundary conditions for the two equations can be cast as a couple of three-term continued fraction relations of the form (2.9). Computing quasinormal frequencies ω and the related angular eigenvalues is simple. For assigned values of s, l, m, a and ω, first find the angular separation constant s A lm looking for zeros of the angular continued fraction. Then replace the corresponding eigenvalue into the radial continued fraction, and look for its zeros as a function of ω.
Once we know the quasinormal mode frequencies and the angular eigenvalues, we can find the angular eigenfunction following the procedure of Sec. II B. The results are plotted in Fig. 7 for different gravitational (s = −2) quasinormal frequencies with l = 2, |m| ≤ l and different values of a/M . The modulus of the eigenfunction is scarcely affected by rotation and almost unaffected by the overtone index, at least for slowly damped modes. The phase (hence the relative magnitude of the imaginary part of the eigenfunctions) is more sensitive to angular momentum and mode damping. It is generally larger for counterrotating modes, and it grows with the overtone index N . In general, the magnitude of the phase seems to be related with the magnitude of the imaginary part of the quasinormal frequencies.
For corotating modes in the extremal limit the imaginary part of quasinormal frequencies goes to zero, as first shown by Detweiler [44]. This explains what, at first sight, might look like a puzzling feature of the plots: for near-extremal Kerr black holes (a/M = 0.98, green lines) the phase is actually smaller than for black holes with a/M = 0.80 if the modes are corotating. The opposite is true for counterrotating modes, whose imaginary part does not tend to zero in the extremal limit.
In analysing data from gravitational waveforms emitted by an oscillating black hole it is necessary to compute the scalar products [24] Here and in the following we append to the quasinormal frequencies the angular indices (l, m) and a third index (the "overtone" index N = 0, 1, 2 . . . ) that sorts frequencies by the magnitude of their imaginary part: in the limit a/M = 0, modes with larger N have larger imaginary part and damp faster. The angular momentum-dependent quantity α mll ′ N N ′ (a) will be determined below, and the Kronecker symbol δ m,m ′ comes from the e imφ -dependence of the harmonics. For a/M = 0 the spin-weighted spheroidal harmonics reduce to spin-weighted spherical harmonics s Y lm [12], for which We evaluated SWSHs and their scalar products by two independent, approximate calculations.
In the first calculation we use the sixth-order polynomial approximation for the separation constant s A lm computed by Press and Teukolsky (Table I in Ref. [36]). We plug this expansion in the spheroidal wave equation and integrate the equation numerically, imposing the appropriate boundary conditions at one extremum of the integration interval. In principle we should expect the agreement with Leaver's method to be good only for small values of aω. Surprisingly, results for the eigenfunctions and the scalar products turn out to be very close also for large aω, where the polynomial approximation of the "true" separation constant is not so accurate. Even for a/M = 0.98 the approximate scalar products agree with Leaver's method within ∼ 5%.
As a second independent, analytic check we solved the spheroidal wave equation expanding in powers of aω and using standard perturbation theory. For aω = 0 the solutions are the spin-weighted spherical harmonics [12]. The next order correction was found by Press and Teukolsky (Eq. (3.7) of Ref. [36]). Using their approximate solution we can show that, to leading order, where α mll ′ N N ′ (a) has been defined in Eq. (2.28), and j 1 j 2 m 1 m 2 |JM is a Clebsch-Gordan coefficient. For l = l ′ the same first-order perturbative analysis gives α mll ′ nn ′ (a) ≃ 1.
In Table I we compare values of the scalar products obtained using our two approximate methods with the results of Leaver's method. There we focus on the worst possible scenario of a nearly maximally rotating Kerr black hole (a/M = 0.98). Both approximations work really well, providing a powerful consistency check of the numerical calculation.  For reference, the constants α mll ′ nn ′ (a) computed using Leaver's method are given in Tables II and III for fixed values of m = m ′ (for m = m ′ the scalar product is zero). In the context of gravitational wave data analysis, it is most useful to consider scalar products between the quasinormal modes which are more likely to be excited during collapse or merger: say, the first three Kerr overtones (N, N ′ ) = 0, 1, 2 with (l, l ′ ) = 2, 3, 4 and |m| ≤ 2. Tables II  and III show that, to a reasonable level of approximation, we can assume α mll ′ nn ′ (a) ≃ δ ll ′ , that is, The errors on this approximate relation increase with the black hole's angular momentum a/M . They are typically larger for small m, but in most cases the approximate formula is valid to within ∼ 10%, at least for cases of physical interest (say, for l < 5).       To save space, we omit the leading zeros and the i's in the imaginary parts (so, for example, .9997 ± .0213 actually means 0.9997 ± 0.0213i). Entries replaced by an asterisk can be obtained from the symmetric entries by complex conjugation: i.e., α mll ′ NN ′ (a) = α * ml ′ lN ′ N (a). m a/M = .98 In four dimensions there is only one possible rotation axis for a cylindrically symmetric space-time, so there is only one angular momentum parameter. In higher dimensions there are several possible choices of rotation axes. Therefore higher-dimensional black holes will have, in general, several angular momentum parameters, each referring to a particular rotation axis [29]. Here we focus on the simple case where there is only one rotation axis (and correspondingly, only one angular momentum direction). The angular separation constant A kjm now depends on three indices. The integer m = 0, ± 1, ± 2 . . . comes from separation of the angle describing the azimuthal dependence of the perturbations around the symmetry axis. Since (3.1) is invariant under the replacement m → −m we will restrict our calculations to m ≥ 0. The number of extra dimensions n and the parameter j = 0, 1, 2, . . . are related with the eigenvalues of hyperspherical harmonics on the n-sphere, which are given by −j(j + n − 1). The parameter k(= 0, 1, 2 . . . ) plays the same role as r/2 in (2.11) in the four-dimensional case, i.e. it labels the eigenvalues for fixed values of j and m. Eq. (3.1) describes the angular dependence of a scalar field perturbing a (n + 4)-dimensional Myers-Perry black hole with a single rotation parameter, and is a simple generalization of the four-dimensional case [30,31,32,33]. When n = 1 it also describes scalar perturbations of smooth geometries in the D1-D5 stringy system [35].
For finite but small c we proceed as in Sec. II C. When c = 0 the recursion relation (3.7) has a finite number of terms k. Therefore, expanding around c = 0 we find it convenient (although strictly not necessary) to use the k-th inversion of Eq. (3.7):β Now we expand the separation constant as a Taylor series: wheref 0 = l(l + n + 1) , (3.12a) The coefficientsf 1 , . . . ,f 5 were obtained substituting (3.11) into Eq. (3.10) and expanding the resulting expression in powers of c. Higher order coefficients can be obtained easily, but the expressions are too lengthy to reproduce them here. The four-dimensional results can be obtained by setting j = n = 0. They are in agreement with [37,38] and with our Eqs. (2.16a)-(2.16g).
C. Large and real c (oblate case) We now proceed to calculate the asymptotic behaviour of the oblate HSHs (i.e., Eq. (3.1) with c ∈ R) for large c, following the same method as in Sec. II D. In order to find asymptotic solutions valid near the end-points (i.e. "inner solutions") we make the change of independent variable u ≡ 2c(1 − x). The ordinary differential equation (3.3) then becomes u d 2 y kjm du 2 + (|m| + 1) This equation is the same as the corresponding equation in four dimensions where δ kjm plays the role of s A lm , except for the presence of the last term in Eq. (3.13). This last term can be neglected to leading order in c. The asymptotic solution of Eq. (3.13) valid near the end-points is therefore the same as the inner solution in four dimensions (given by Eq. (3.10) of [27] with h = 0), so that δ kjm = −c 2 + 2qc Our task is now to determine the parameter q. We can find an asymptotic solution which is valid far from the end-points (the "outer" solution) as in [27]. The result is: where a kjm and b kjm are constants of integration. This solution is similar to the outer solution in four dimensions, i.e., Eq. (3.26) of [27] with h = 0. Unlike the four-dimensional case, (3.14) is not valid very near the origin x = 0. Despite this inconvenience, the matching of the inner and outer solutions in the overlap region proceeds exactly as in four dimensions: the number of extra dimensions n plays no part in the matching. The next step is to find the number of zeros of the HSHs. Equating this number with the number of zeros of the asymptotic solution determines the value of q. Refs. [1,6] show that scalar spheroidal harmonics (a special limit of HSHs with s = n = j = 0) are entire functions of c, so their number of zeros does not vary with c. In particular, their number of zeros is the same as that of the spherical harmonics (s = n = j = c = 0), which is known. Similarly, Refs. [40,45,46] show that SWSHs (n = j = 0) are entire functions of c. Therefore their number of zeros is easily obtained, since it does not vary with c. The argument used in both cases relies crucially on two properties: (P1) the fact that the indicial equation at any possible regular singular point of the differential equation does not involve c, (P2) the continuity in x of the coefficients in the differential equation when the coefficient of the second order derivative is equal to one.
We must be careful if we wish to use a similar argument to show that the number of zeros of the HSHs does not vary with c. First of all, the characteristic exponentsα 0 andα 1 for the HSH equation do not depend on c, and the differential equation (3.3) does not possess any singular points. Therefore the differential equations for S kjm and y kjm both satisfy the property (P1) above. However, the coefficients of either differential equation are not continuous at x = 0 when the coefficient of the second order derivative is equal to one. Therefore, property (P2) is not satisfied by either differential equation. Nevertheless, we can still apply the argument separately to the regions x ∈ (−1, 0) and x ∈ (0, +1), where both properties (and other necessary ones) are satisfied by both differential equations. Note that no zero of y kjm can "move" from either region x ∈ (−1, 0) or x ∈ (0, +1) to the other as c varies, because this function cannot possess a zero at x = 0. We then conclude that the number of zeros of y kjm in either interval (−1, 0) or (0, +1) [and therefore in the whole region (−1, +1)] is independent of c. In particular, it is equal to the number of zeros of y kjm in the corresponding region for c = 0. So the number of zeros of the HSHs for x ∈ (−1, +1) for any real c can be obtained from the special case c = 0.
In Appendix A we find an exact solution of the HSH equation for c = 0 and show that the number of zeros of the HSHs in x ∈ (−1, +1), excluding x = 0, is equal to 2k. The HSHs possess a zero at x = 0 if j > 0 and do not if j = 0; however, we will not need this property. As a check of the analytic results we performed numerical calculations of the eigenfunctions using Eq. (3.4); some of these numerically computed, normalized eigenfunctions are shown in Fig. 8.
The matching of the number of zeros of the HSH with that of the asymptotic solution proceeds now as in four dimensions. Because of the forms of the inner and outer asymptotic solutions, the number of zeros of the asymptotic solution (excluding x = 0) has the dependence on q given in Eq. (4.4) of [27] with h = 0 (h in that paper is what we call s here). We can equate this number to the number of zeros of the HSH (excluding possible zeros at x = 0), which is 2k = l − |m| − j.
The resulting equation is the same as Eq. (4.4) in [27] with h = 0, and with z 0 replaced by j. The value of the parameter q is therefore q = l − j + 1 = 2k + |m| + 1 . The asymptotic behaviour of the separation variable A kjm is then given by There is an imprint of the extra dimensions (through the parameter j) already at order O(c). Notice that the above asymptotic expansion (3.16) does not agree in the limit n = j = 0 with the four-dimensional result, Eqs.  To conclude this Section, we back up our analytic results by numerical calculations of the oblate eigenvalues obtained solving Eq. (3.7). In Fig. 9 we plot the eigenvalues for n = k = 1, m = 1, 2 and j = 0, 1, 2; these values of (n, k, j, m) have been chosen to match those in Figs. 2-7 of [31]. We checked that the asymptotic behavior of the oblate eigenvalues, both in the limit c → 0 and in the limit c → ∞, is in agreement with the analytic predictions. In this Section we complete the picture computing a series expansions of the eigenvalue A kjm for large and pureimaginary values of c. The present analysis is similar to the one in Sec. II E.
We first find an approximation valid near the origin when |c I | → ∞ by introducing a new independent variable u ≡ 2|c I |x and using again the transformation (3.2) (withk 0 = |m| andk 1 = j). The resulting ordinary differential equation for the function y kjm is the same as Eq. (3.3), but in terms of the new variable u: Now we apply the same approximation as in Sec. II E. In the limit |c I | → ∞ and in the inner region 1/|c I | ≪ |x| ≪ 1, Eq. (3.17) reduces to where we have set δ kjm ∼ 2q|c I | + O |c I | 0 , and q is an unknown parameter that we shall now determine, not to be confused with the parameter introduced in (3.15). With respect to the four-dimensional case, Eq. (3.18) contains an extra term with a first derivative. The presence of this term means that the inner solution y inn kjm is not a parabolic cylinder function, as in the four-dimensional case. By applying the change of variables Eq. (3.18) becomes The solution of this equation can be expressed in terms of a Laguerre polynomial L (α) n [47]. It can be checked that the asymptotic approximation of the prolate HSHs valid in the inner region 1/|c I | ≪ |x| ≪ 1 as |c I | → ∞, which is regular at x = 0, is y inn kjm (u) = e −u 2 /4 L n (x) is equal toñ if (as in our case)α > −1 [47]. It can easily be checked that all these zeros lie within the inner region. The number of zeros of y inn kjm in the inner region [or equivalently in the region x ∈ (−1, +1)] is therefore equal to 2ñ.
For |c I | → ∞, an outer solution which is a valid approximation of the prolate HSH near the end-points can be obtained using a WKB-type expansion, as in [27] and in Sec. II E. The result is where a kjm and b kjm are constants of integration. For n = j = 0, (3.23) reduces to the corresponding outer solution (2.26) in the prolate SWSH case if s = 0 and we replace L + 1/2 there by q. Furthermore This outer solution has no zeros. Like in the prolate SWSH case with s = 0, all zeros of the asymptotic solution for |c I | → ∞ are provided by the inner solution. We want to equate their number to the number of zeros of the exact solution y kjm in the region x ∈ (−1, +1), which is given in Appendix A for c = 0, is independent of c 2 ∈ (−∞, +∞) and is equal to 2k. Therefore 2q = 4k + n + 2j + 1, and the separation constant for large |c I | is A kjm = δ kjm + O |c I | 0 = (4k + n + 2j + 1)|c I | + O |c I | 0 = [2(l − |m|) + n + 1] |c I | + O |c I | 0 . The expansion (3.25) is in agreement with the eigenvalues obtained from a numerical solution of Eq. (3.7) in the prolate case, some of which are plotted in Fig. 10.

IV. CONCLUSIONS
In this paper we provided a complete numerical and analytic analysis of four-dimensional spin-weighted spheroidal harmonics (SWSHs) and of higher-dimensional spheroidal harmonics (HSHs). For reference we list below what we regard as the most useful results, be they known in the existing literature or originally derived in this paper.
For SWSHs, the main results are: In our opinion, the main unsolved theoretical problem is to clarify why Eq. (2.24) is in good agreement with the numerics for |s| = 0, even though the number of zeros of the eigenfunction seems to depend on c when c is not real. Given the excellent accuracy of Eq. (2.30) in the four-dimensional case, it would be desirable to extend the perturbative formalism developed by Press and Teukolsky [36] to higher dimension. It could also be useful to push the expansions (2.24), (3.16) and (3.25) to higher order, and to extend the present understanding of the location of the spheroidal branch points [10] to the case s = 0. Finally, it would be an exciting and interesting challenge to extend our results for higher-dimensional scalar harmonics to higher-dimensional harmonics with s = 0. We leave these issues for future work.