Analytic Black Hole Perturbation Approach to Gravitational Radiation

Abstract

We review the analytic methods used to perform the post-Newtonian expansion of gravitational waves induced by a particle orbiting a massive, compact body, based on black hole perturbation theory. There exist two different methods of performing the post-Newtonian expansion. Both are based on the Teukolsky equation. In one method, the Teukolsky equation is transformed into a Regge-Wheeler type equation that reduces to the standard Klein Gordon equation in the flat-space limit, while in the other method (which was introduced by Mano, Suzuki, and Takasugi relatively recently, the Teukolsky equation is used directly in its original form. The former’s advantage is that it is intuitively easy to understand how various curved space effects come into play. However, it becomes increasingly complicated when one goes to higher and higher post-Newtonian orders. In contrast, the latter’s advantage is that a systematic calculation to higher post-Newtonian orders can be implemented relatively easily, but otherwise, it is so mathematical that it is hard to understand the interplay of higher order terms. In this paper, we review both methods so that their pros and cons may be seen clearly. We also review some results of calculations of gravitational radiation emitted by a particle orbiting a black hole.

1 Introduction

1.1 General

In the past several years, there has been substantial progress in the projects of ground-based laser interferometric gravitational wave detectors, which include LIGO [65], VIRGO [108], GE0600 [50], and TAMA300 [102, 3, 103]. TAMA300 was in operation from 1999 until 2004. LIGO and GE0600 began operating in 2002. LIGO has performed it’s fifth science run from November 2005 to October 2007, and over one year of science-quality data was taken with all of it’s three LIGO interferometers in simultaneous operation. The Virgo detector performed its first science run in 2007 and 4.5 months of joint data-taking with LIGO was done. There are several future projects as well. Most importantly, the Laser Interferometer Space Antenna (LISA) project is in progress [67, 66]. The DECIGO [2] and BBO [36] are more ambitious space interferometer proposals which aim to cover the frequency gap between the ground based interferometers and LISA. For a review of ground and space laser interferometers, see, e.g., the respective Living Reviews article [88].

The detection of gravitational waves will be done by extracting gravitational wave signals from a noisy data stream. In developing the data analysis strategy, detailed knowledge of the gravitational waveforms will help us greatly to detect a signal, and to extract the physical information about its source. Thus, it has become a very important problem for theorists to predict with sufficiently good accuracy the waveforms from possible gravitational wave sources.

Gravitational waves are generated by dynamical astrophysical events, and they are expected to be strong enough to be detected when compact stars such as neutron stars (NS) or black holes (BH) are involved in such events. In particular, coalescing compact binaries are considered to be the most promising sources of gravitational radiation that can be detected by the ground-based laser interferometers. The last inspiral phase of a coalescing compact binary, in which the binary stars orbit each other for ∼ 10⁴ cycles, will be in the bandwidth of the interferometers, and this phase may not only be detectable: it could provide us with important astrophysical information about the system, if the theoretical templates are sufficiently accurate.

Unfortunately, it seems difficult to attain the sensitivity to detect NS-NS binary inspirals with the first generation of interferometric detectors. However, the coalescence of BH-NS/BH-BH binaries with a black hole mass of ∼ 10–20 M_⊙ may be detected out to the distance of the VIRGO cluster if we are lucky enough. In any case, it will be necessary to wait for the next generation of interferometric detectors to see these coalescing events more frequently [73, 82].

To predict the waveforms, a conventional approach is to formulate the Einstein equations with respect to the flat Minkowski background and apply the post-Newtonian expansion to the resulting equations (see the Section 1.2).

In this paper, however, we review a different approach, namely the black hole perturbation approach. In this approach, binaries are assumed to consist of a massive black hole and a small compact star which is taken to be a point particle. Hence, its applicability is constrained to the case of binaries with large mass ratio. Nevertheless, there are several advantages here that cannot be overlooked.

Most importantly, the black hole perturbation equations take full account of general relativistic effects of the background spacetime and they are applicable to arbitrary orbits of a small mass star. In particular, if a numerical approach is taken, gravitational waves from highly relativistic orbits can be calculated. Then, if we can develop a method to calculate gravitational waves to a sufficiently high PN order analytically, it can give insight not only into how and when general relativistic effects become important, by comparing with numerical results, but it will also give us a knowledge, complementary to the conventional post-Newtonian approach, about as yet unknown higher-order PN terms or general relativistic spin effects.

Moreover, one of the main targets of LISA is to observe phenomena associated with the formation and evolution of supermassive black holes in galactic centers. In particular, a gravitational wave event of a compact star spiraling into such a supermassive black hole is indeed a case of application for the black hole perturbation theory.

1.2 Post-Newtonian expansion of gravitational waves

The post-Newtonian expansion of general relativity assumes that the internal gravity of a source is small so that the deviation from the Minkowski metric is small, and that velocities associated with the source are small compared to the speed of light, c. When we consider the orbital motion of a compact binary system, these two conditions become essentially equivalent to each other. Although both conditions may be violated inside each of the compact objects, this is not regarded as a serious problem of the post-Newtonian expansion, as long as we are concerned with gravitational waves generated from the orbital motion, and, indeed, the two bodies are usually assumed to be point-like objects in the calculation.

In fact, Itoh, Futamase, and Asada [57, 58] developed a new post-Newtonian method that can deal with a binary system in which the constituent bodies may have strong internal gravity, based on earlier work by Futamase and Schutz [44, 45, 42]. They derived the equations of motion to 2.5PN order and obtained a complete agreement with the Damour-Deruelle equations of motion [26, 25], which assumes the validity of the point-particle approximation. In the Futamase-Schutz method, each star in a binary is first expressed as an extended object and then the limit is taken to set the radius to zero in a specific manner first proposed by Futamase [42]. At the same time, the surface integral approach (a la Einstein-Infeld-Hoffmann [32]) is taken to derive the equations of motion. More recently Itoh and Futamase [56, 54] derived the 3PN equations of motion based on the Futamase-Schutz method, and they are again in agreement with those derived by Damour, Jaranowski and Schafer [27] and by Blanchet et al. [10] in which the point-particle approximation is used.

There are two existing approaches of the post-Newtonian expansion to calculate gravitational waves: one developed by Blanchet, Damour, and Iyer (BDI) [12, 7] and another by Will and Wiseman (WW) [111] based on previous work by Epstein, Wagoner, and Will [33, 109]. In both approaches, the gravitational waveforms and luminosity are expanded in time derivatives of radiative multipoles, which are then related to some source multipoles (the relation between them contains the “tails”). The source multipoles are expressed as integrals over the matter source and the gravitational field. The source multipoles are combined with the equations of motion to obtain explicit expressions in terms of the source masses, positions, and velocities.

One issue of the post-Newtonian calculation arises from the fact that the post-Newtonian expansion can be applied only to the near-zone field of the source. In the conventional post-Newtonian formalism, the harmonic coordinates are used to write down the Einstein equations. If we define the deviation from the Minkowski metric as

$${h^{\mu v}} \equiv \sqrt { - g} {g^{\mu v}} - {\eta ^{\mu v}},$$

((1))

the Einstein equations are schematically written in the form

$$\square {h^{\mu v}} = 16\pi |g|{T^{\mu v}} + {\Lambda ^{\mu v}}(h),$$

((2))

together with the harmonic gauge condition, ∂_νh^μν = 0, where □ = η^μν∂_μ∂_ν is the D’Alambertian operator in flat-space time, η^μν = diag (−1, 1, 1, 1), and Λ^μν (h) represents the non-linear terms in the Einstein equations. The Einstein equations (2) are integrated using the flat-space retarded integrals. In order to perform the post-Newtonian expansion, if we naively expand the retarded integrals in powers of l1c, there appear divergent integrals. This is a technical problem that arises due to the near-zone nature of the post-Newtonian approximation. In the BDI approach, in order to integrate the retarded integrals, and to evaluate the radiative multipole moments at infinity, two kinds of approximation methods are introduced. One is the multipolar post-Minkowski expansion, which can be applied to a region outside the source including infinity, and the other is the near-zone, post-Newtonian expansion. These two expansions are matched in the intermediate region where both expansions are valid, and the radiative multipole moments are evaluated at infinity. In the WW approach, the retarded integrals are evaluated directly, without expanding in terms of 1/c, in the region outside the source in a novel way.

The lowest order of the gravitational waves is given by the Newtonian quadrupole formula. It is standard to refer to the post-Newtonian formulae (for the waveforms and luminosity) that contain terms up to $\mathcal{O}\left( {{{(v/c)}^n}} \right)$ beyond the Newtonian quadrupole formula as the (n/2)PN formulae. Evaluation of gravitational waves emitted to infinity from a compact binary system has been successfully carried out to the 3.5 post-Newtonian (PN) order beyond the lowest Newtonian quadrupole formula in the BDI approach [12, 7, 18, 19, 15, 16, 11]. Up to now, the WW approach gives the same answer for the gravitational waveforms and luminosity to 2PN order.

The computation of the 3.5PN flux requires the 3PN equations of motion. As mentioned in the above, the 3PN equations of motion have been derived by three different methods. The first is the direct post-Newtonian iteration in the harmonic coordinates [14, 28, 10]. The second employs the Arnowitt-Deser-Misner (ADM) coordinates within the Hamiltonian formalism of general relativity [59, 60, 27]. The third is based on the Futamase-Schutz method [56, 54].

Since the first two methods use the point particle approximation while the third one is not, let us first focus on the first two. In both methods, since the stars are represented by the Dirac delta functions, the divergent self-fields must be regularized. In earlier papers, they used the Hadamard regularization method [59, 60, 14, 28]. However, it turned out that there remains an unknown coefficient which cannot be determined within the regularization method. This problem was solved by Damour, Jaranowski and Schäfer [27] who successfully derived the 3PN equations of motion without undetermined numerical coefficients by using the dimensional regularization within an ADM Hamiltonian approach. Then the 3PN equations of motion in the harmonic coordinates were also derived without undetermined coefficients by using a combination of the Hadamard regularization and the dimensional regularization in [10]. The 3.5PN radiation reaction terms in the equations of motion are also derived in both approaches [76, 62]. See reviews by Blanchet [8, 9] for details and summaries on post-Newtonian approaches.

In the case of Futamase-Schutz method, as mentioned in the beginning of this subsection, the 3PN equations of motion is derived by Itoh and Futamase [56, 54], and the 3.5PN terms are derived by Itoh [55]. See a review article by Futamase and Itoh [43] for details on this method.

There are other methods in which stars are treated as fluid balls [51, 63, 80, 81]. Pati and Will [80, 81] use an method which is an extension of the WW approach in which the retarded integral is evaluated directly. With these method, the 2PN equations of motion as well as 2.5PN and 3.5PN radiation reaction effects are derived.

1.3 Linear perturbation theory of black holes

In the black hole perturbation approach, we deal with gravitational waves from a particle of mass μ orbiting a black hole of mass M, assuming μ ≪ M. The perturbation of a black hole spacetime is evaluated to linear order in μ/M. The equations are essentially in the form of Equation (2) with η_μν, replaced by the background black hole metric g ^BH_μν and the higher order terms Λ(h)_μν, neglected. Thus, apart from the assumption μ ≪ M, the black hole perturbation approach is not restricted to slow-motion sources, nor to small deviations from the Minkowski spacetime, and the Green function used to integrate the Einstein equations contains the whole curved spacetime effect of the background geometry.

The black hole perturbation theory was originally developed as a metric perturbation theory. For non-rotating (Schwarzschild) black holes, a single master equation for the metric perturbation was derived by Regge and Wheeler [87] for the so-called odd parity part, and later by Zerilli [112] for the even parity part. These equations have the nice property that they reduce to the standard Klein-Gordon wave equation in the flat-space limit. However, no such equation has been found in the case of a Kerr black hole so far.

Then, based on the Newman-Penrose null-tetrad formalism, in which the tetrad components of the curvature tensor are the fundamental variables, a master equation for the curvature perturbation was first developed by Bardeen and Press [6] for a Schwarzschild black hole without source (T^μν = 0), and by Teukolsky [106] for a Kerr black hole with source (T^μν ≠ 0). The master equation is called the Teukolsky equation, and it is a wave equation for a null-tetrad component of the Weyl tensor ψ₀ or ψ₄. In the source-free case, Chrzanowski [23] and Wald [110] developed a method to construct the metric perturbation from the curvature perturbation.

The Teukolsky equation has, however, a rather complicated structure as a wave equation. Even in the flat-space limit, it does not reduce to the standard Klein-Gordon form. Later, Chandrasekhar showed that the Teukolsky equation can be transformed to the form of the Regge-Wheeler or Zerilli equation for the source-free Schwarzschild case [21]. A generalization of this to the Kerr case with source was done by Sasaki and Nakamura [92, 93]. They gave a transformation that brings the Teukolsky equation to a Regge-Wheeler type equation that reduces to the Regge-Wheeler equation in the Schwarzschild limit. It may be noted that the Sasaki-Nakamura equation contains an imaginary part, suggesting that either it is unrelated to a (yet-to-be-found) master equation for the metric perturbation for the Kerr geometry or implying the non-existence of such a master equation.

As mentioned above, an important difference between the black-hole perturbation approach and the conventional post-Newtonian approach appears in the structure of the Green function used to integrate the wave equations. In the black-hole perturbation approach, the Green function takes account of the curved spacetime effect on the wave propagation, which implies complexity of its structure in contrast to the flat-space Green function. Thus, since the system is linear in the black-hole perturbation approach, the most non-trivial task is the construction of the Green function.

There are many papers that deal with a numerical evaluation of the Green function and calculations of gravitational waves induced by a particle. See Breuer [20], Chandrasekhar [22], and Nakamura, Oohara, and Kojima [72] for reviews and for references on earlier papers.

Here, we are interested in an analytical evaluation of the Green function. One way is to adopt the post-Minkowski expansion assuming GM/c² ≪ r. Note that, for bound orbits, the condition GM/c² ≪ r is equivalent to the condition for the post-Newtonian expansion, v²/c² ≪ 1. If we can calculate the Green function to a sufficiently high order in this expansion, we may be able to obtain a rather accurate approximation of it that can be applicable to a relativistic orbit fairly close to the horizon, possibly to a radius as small as the inner-most stable circular orbit (ISCO), which is given by r_ISCO = 6GM/c² in the case of a Schwarzschild black hole.

It turns out that this is indeed possible. Though there arise some complications as one goes to higher PN orders, they are relatively easy to handle as compared to situations one encounters in the conventional post-Newtonian approaches. Thus, very interesting relativistic effects such as tails of gravitational waves can be investigated easily. Further, we can also easily investigate convergence properties of the post-Newtonian expansion by comparing a numerically calculated exact result with the corresponding analytic but approximate result. In this sense, the analytic black-hole perturbation approach can provide an important test of the post-Newtonian expansion.

1.4 Brief historical notes

Let us briefly review some of the past work on post-Newtonian calculations in black-hole perturbation theory. Although the literature on numerical calculations of gravitational waves emitted by a particle orbiting a black hole is abundant, there are not so many papers that deal with the post-Newtonian expansion of gravitational waves, mainly because such an analysis was not necessary until recently, when the construction of accurate theoretical templates for the interferometric gravitational wave detectors became an urgent issue.

In the case of orbits in the Schwarzschild background, one of the earliest papers was by Gal’tsov, Matiukhin and Petukhov [47], who considered the case when a particle is in a slightly eccentric orbit around a Schwarzschild black hole, and calculated the gravitational waves up to 1PN order. Poisson [83] considered a circular orbit around a Schwarzschild black hole and calculated the waveforms and luminosity to 1.5PN order at which the tail effect appears. Cutler, Finn, Poisson, and Sussman [24] worked on the same problem numerically by applying the least-square fitting technique to the numerically evaluated data for the luminosity, and obtained a post-Newtonian formula for the luminosity to 2.5PN order. Subsequently, a highly accurate numerical calculation was carried out by Tagoshi and Nakamura [99]. They obtained the formulae for the luminosity to 4PN order numerically by using the least-square fitting method. They found the log v terms in the luminosity formula at 3PN and 4PN orders. They concluded that, although the convergence of the post-Newtonian expansion is slow, the luminosity formula accurate to 3.5PN order will be good enough to represent the orbital phase evolution of coalescing compact binaries in theoretical templates for ground-based interferometers. After that, Sasaki [91] found an analytic method and obtained formulae that were needed to calculate the gravitational waves to 4PN order. Then, Tagoshi and Sasaki [100] obtained the gravitational waveforms and luminosity to 4PN order analytically, and confirmed the results of Tagoshi and Nakamura. These calculations were extended to 5.5PN order by Tanaka, Tagoshi, and Sasaki [105]. Fujita and Iyer [39] extended this work and derived 5.5PN waveforms.

In the case of orbits around a Kerr black hole, Poisson calculated the 1.5PN order corrections to the waveforms and luminosity due to the rotation of the black hole, and showed that the result agrees with the standard post-Newtonian effect due to spin-orbit coupling [84]. Then, Shibata, Sasaki, Tagoshi, and Tanaka [94] calculated the luminosity to 2.5PN order. They calculated the luminosity from a particle in circular orbit with small inclination from the equatorial plane. They used the Sasaki-Nakamura equation as well as the Teukolsky equation. This analysis was extended to 4PN order by Tagoshi, Shibata, Tanaka, and Sasaki [101], in which the orbits of the test particles were restricted to circular ones on the equatorial plane. The analysis in the case of slightly eccentric orbit on the equatorial plane was also done by Tagoshi [95, 96] to 2.5PN order.

Tanaka, Mino, Sasaki, and Shibata [104] considered the case when a spinning particle is in a circular orbit near the equatorial plane of a Kerr black hole, based on the Papapetrou equations of motion for a spinning particle [79] and the energy momentum tensor of a spinning particle by Dixon [29]. They derived the luminosity formula to 2.5PN order which includes the linear order effect of the particle’s spin.

The absorption of gravitational waves into the black hole horizon, appearing at 4PN order in the Schwarzschild case, was calculated by Poisson and Sasaki for a particle in a circular orbit [85]. The black hole absorption in the case of a rotating black hole appears at 2.5PN order [46]. Using a new analytic method to solve the homogeneous Teukolsky equation found by Mano, Suzuki, and Takasugi [68], the black hole absorption in the Kerr case was calculated by Tagoshi, Mano, and Takasugi [98] to 6.5PN order beyond the quadrupole formula.

If gravity is not described by the Einstein theory but by the Brans-Dicke theory, there will appear scalar-type gravitational waves as well as transverse-traceless gravitational waves. Such scalar-type gravitational waves were calculated to 2.5PN order by Ohashi, Tagoshi, and Sasaki [77] in the case when a compact star is in a circular orbit on the equatorial plane around a Kerr black hole.

In the above works the energy and angular momentum flux at infinity or the absorption rate at the horizon were evaluated. In the Kerr case, in order to specify the evolution of particle’s trajectory under the influence of radiation reaction, we need to determine the rate of change of the Carter constant which is not directly related to the asymptotic gravitational waves. Mino [70] proved that the average rate of change of the Carter constant can be evaluated by using the radiative field (i.e., retarded minus advanced field) in the adiabatic approximation. An explicit calculation of the rate of change of the Carter constant was done in the case of a scalar charged particle in [30]. Sago et al. [90] extended Mino’s work and found a simpler formula for the average rate of change of the Carter constant. They derived analytically the rate of change of the Carter constant as well as the energy and the angular momentum of a particle for orbits with small eccentricities and inclinations up to O(v⁵) [89]. In Ref. [48], the method was extended to the case of the orbits with small eccentricity but arbitrary inclination angle, and the rate of change of the energy, angular momentum and the Carter constant up to O(v⁵) were derived.

In the rest of the paper, we use the units c = G = 1.

2 Basic Formulae for the Black Hole Perturbation

2.1 Teukolsky formalism

In terms of the conventional Boyer-Lindquist coordinates, the metric of a Kerr black hole is expressed as

$$\begin{gathered} d{s^2} = - \frac{\Delta }{\Sigma }{(dt - a\;{\sin ^2}\theta d\varphi )^2} + \frac{{{{\sin }^2}\theta }}{\Sigma }{[({r^2} + {a^2})d\varphi - a\;dt]^2} \hfill \\ \;\;\;\;\;\;\;\;\; + \frac{\Sigma }{\Delta }d{r^2} + \Sigma d{\theta ^2}, \hfill \\ \end{gathered} $$

((3))

where Σ, = r² + a² cos² θ and Δ = r² - 2Mr + a². In the Teukolsky formalism [106], the gravitational perturbations of a Kerr black hole are described by a Newman-Penrose quantity ${\psi _4} = - {C_{\alpha \beta \gamma \delta }}{n^\alpha }{\bar m^\beta }{n^\gamma }{\bar m^\delta }$ [74, 75], where C_αβγδ is the Weyl tensor and

$${n^\alpha } = \frac{1}{{2\Sigma }}(({r^2} + {a^2}), - \Delta ,0,a),$$

((4))

$${m^\alpha } = \frac{1}{{\sqrt 2 (r + ia\;\cos \theta )}}(ia\;\sin \theta ,0,1,i/\sin \theta ).$$

((5))

The perturbation equation for φ ≡ ρ⁻⁴ψ₄, ρ = (r - is cos θ)⁻¹, is given by

$${}_s\mathcal{O}\phi = 4\pi \Sigma \hat T$$

((6))

Here, the operator $_s\mathcal{O}$ is given by

$$\begin{gathered} {}_s\mathcal{O} = - \left( {\frac{{{{({r^2} + {a^2})}^2}}}{\Delta } - {a^2}{{\sin }^2}\theta } \right)\partial _t^2 - \frac{{4Mar}}{\Delta }{\partial _t}{\partial _\phi } - \left( {\frac{{{a^2}}}{\Delta } - \frac{1}{{{{\sin }^2}\theta }}} \right)\partial _\phi ^2 \hfill \\ \;\;\;\;\;\;\;\; + {\Delta ^{ - s}}{\partial _r}({\Delta ^{s + 1}}{\partial _r}) + \frac{1}{{\sin \theta }}{\partial _\theta }(\sin \theta {\partial _\theta }) + 2s\left( {\frac{{a(r - M)}}{\Delta } + \frac{{i\;\cos \theta }}{{{{\sin }^2}\theta }}} \right){\partial _\phi } \hfill \\ \;\;\;\;\;\;\;\; + 2s\left( {\frac{{M({r^2} - {a^2})}}{\Delta } - r - ia\;\cos \theta } \right){\partial _t} - s(s\;{\cot ^2}\theta - 1), \hfill \\ \end{gathered} $$

((7))

with s = −2. The source term $\hat T$ is given by

$$\hat T = 2({B'_2} + B{_2^{*^\prime }}),$$

((8))

$$\begin{gathered} {{B'}_2} = - \frac{1}{2}{\rho ^8}\bar \rho {{\hat L}_{ - 1}}[{\rho ^{ - 4}}{{\hat L}_0}(\rho {}^{ - 2}{{\bar \rho }^{ - 1}}{T_{nn}})] \hfill \\ \;\;\;\;\;\;\; - \frac{1}{{2\sqrt 2 }}{\rho ^8}\bar \rho {\Delta ^2}{{\hat L}_{ - 1}}[{\rho ^{ - 4}}{{\bar \rho }^2}{{\hat J}_ + }({\rho ^{ - 2}}{{\bar \rho }^{ - 2}}{\Delta ^{ - 1}}{T_{\bar mn}})], \hfill \\ \end{gathered} $$

((9))

$$\begin{gathered} {{B'}_2}* = - \frac{1}{4}{\rho ^8}\bar \rho {\Delta ^2}{{\hat J}_ + }[{\rho ^{ - 4}}{{\hat J}_ + }({\rho ^{ - 2}}\bar \rho {T_{\overline {mn} }})] \hfill \\ \;\;\;\;\;\;\; - \frac{1}{{2\sqrt 2 }}{\rho ^8}\bar \rho {\Delta ^2}{{\hat J}_ + }[{\rho ^{ - 4}}{{\bar \rho }^2}{\Delta ^{ - 1}}{{\hat L}_{ - 1}}({\rho ^{ - 2}}{{\bar \rho }^{ - 2}}{T_{\bar mn}})], \hfill \\ \end{gathered} $$

((10))

where

$${\hat L_s} = {\partial _\theta } - \frac{i}{{\sin \theta }}{\partial _\varphi } - ia\;\sin \theta {\partial _t} + s\;\cot \theta ,$$

((11))

$${\hat J_ + } = {\partial _r} - \frac{1}{\Delta }(({r^2} + {a^2}){\partial _t} + a{\partial _\varphi }),$$

((12))

and ${T_{nn}},{T_{\bar mn}},\;\text{and}\;{T_{\overline {mm} }}$, are the tetrad components of the energy momentum tensor (T_nn = T_μνn^μn^ν etc.). The bar denotes the complex conjugation.

If we set s = 2 in Equation (6), with appropriate change of the source term, it becomes the perturbation equation for ψ₀. Moreover, it describes the perturbation for a scalar field (s = 0), a neutrino field (|s| = 1/2), and an electromagnetic field (|s| = 1) as well.

We decompose ψ₄ into the Fourier-harmonic components according to

$${\rho ^{ - 4}}{\psi _4} = \mathop \Sigma \limits_{\ell m} \frac{1}{{\sqrt {2\pi } }}\int {d\omega {e^{ - i\omega t + im\varphi }} - 2{S_{\ell m}}(\theta ){R_{\ell m\omega }}(r).} $$

((13))

The radial function R_ℓmω, and the angular function _sS_ℓm(θ) satisfy the Teukolsky equations with s = -2 as

$${\Delta ^2}\frac{d}{{dr}}\left( {\frac{1}{\Delta }\frac{{d{R_{\ell m\omega }}}}{{dr}}} \right) - V(r){R_{\ell m\omega }} = {T_{\ell m\omega }},$$

((14))

$$\left[ {\frac{1}{{\sin \theta }}\frac{d}{{d\theta }}\left( {\sin \theta \frac{d}{{d\theta }}} \right) - {a^2}{\omega ^2}{{\sin }^2}\theta - \frac{{{{(m - 2\cos \theta )}^2}}}{{{{\sin }^2}\theta }} + 4a\omega \cos \theta - 2 + 2ma\omega + \lambda } \right] - 2{S_{\ell m}} = 0.$$

((15))

The potential V(r) is given by

$$V(r) = - \frac{{{K^2} + 4i(r - M)K}}{\Delta } + 8i\omega r + \lambda ,$$

((16))

where λ is the eigenvalue of ₋₂S ^aω_ℓm and K = (r² + a²)ω - ma. The angular function _sS_ℓm(θ) is called the spin-weighted spheroidal harmonic, which is usually normalized as

$${\int_0^\pi {|{}_{ - 2}S} _{\ell m}}{|^2}\sin \theta d\theta = 1.$$

((17))

In the Schwarzschild limit, it reduces to the spin-weighted spherical harmonic with \gl \ar ℓ(ℓ + 1). In the Kerr case, however, no analytic formula for λ is known. The source term T_ℓmω is given by

$${T_{\ell m\omega }} = 4\int {d\Omega dt{\rho ^{ - 5}}{{\bar \rho }^{ - 1}}({{B'}_2} + {{B'}_2}*){e^{ - im\varphi + i\omega t}}\frac{{{}_{ - 2}S_{\ell m}^{a\omega }}}{{\sqrt {2\pi } }},} $$

((18))

We mention that for orbits of our interest, which are bounded, T_ℓmω has support only in a compact range of r.

We solve the radial Teukolsky equation by using the Green function method. For this purpose, we define two kinds of homogeneous solutions of the radial Teukolsky equation:

$$R_{\ell m\omega }^{in} \to \left\{ \begin{gathered} B_{\ell m\omega }^{trans}{\Delta ^2}{e^{ - ikr*}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;for\;r \to {r_ + } \hfill \\ {r_3}B_{\ell m\omega }^{ref}{e^{ - i\omega r*}} + {r^{ - 1}}B_{\ell m\omega }^{inc}{e^{ - i\omega r*\;\;\;\;}}for\;r \to + \infty ,\;\; \hfill \\ \end{gathered} \right.$$

((19))

$$R_{\ell m\omega }^{up} \to \left\{ \begin{gathered} C_{\ell m\omega }^{up}{e^{ikr*}} + {\Delta ^2}C_{\ell m\omega }^{ref}{e^{ - ikr*}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;for\;r \to {r_ + } \hfill \\ C_{\ell m\omega }^{trans}{r^3}{e^{i\omega r*}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;for\;r \to + \infty ,\;\; \hfill \\ \end{gathered} \right.$$

((20))

where k = ω - ma/2Mr₊, and r^* is the tortoise coordinate defined by

$$\begin{gathered} r* = \int {\frac{{dr*}}{{dr}}} \hfill \\ \;\;\;\; = r + \frac{{2M{r_ + }}}{{{r_ + } - {r_ - }}}\ln \frac{{r - {r_ + }}}{{2M}} - \frac{{2M{r_ - }}}{{{r_ + } - {r_ - }}}\ln \frac{{r - {r_ - }}}{{2M}}, \hfill \\ \end{gathered} $$

((21))

where ${r_ \pm } = M \pm \sqrt {{M^2} - {a^2}}$, and where we have fixed the integration constant.

Combining with the Fourier mode e^−iwt, we see that R ⁱⁿ_ℓmω has no outcoming wave from past horizon, while R^up has no incoming wave at past infinity. Since these are the properties of waves causally generated by a source, a solution of the Teukolsky equation which has purely outgoing property at infinity and has purely ingoing property at the horizon is given by

$${R_{\ell m\omega }} = \frac{1}{{{W_{\ell m\omega }}}}\left( {R_{_{\ell m\omega }}^{up}\int_{{r_ + }}^r {dr'} R_{_{\ell m\omega }}^{in}{T_{\ell m\omega }}{\Delta ^{ - 2}} + R_{_{\ell m\omega }}^{in}\int_r^\infty {dr'} R_{_{\ell m\omega }}^{up}{T_{\ell m\omega }}{\Delta ^{ - 2}}} \right),$$

((22))

where the Wronskian W_ℓmω is given by

$${W_{\ell m\omega }} = 2i\omega C_{\ell m\omega }^{trans}B_{\ell m\omega }^{inc}.$$

((23))

Then, the asymptotic behavior at the horizon is

$${R_{\ell m\omega }}(r \to {r_ + }) \to \frac{{B_{\ell m\omega }^{trans}{\Delta ^2}{e^{ - ikr*}}}}{{2i\omega C_{\ell m\omega }^{trans}B_{\ell m\omega }^{inc}}}\int_{{r_ + }}^\infty {dr'} R_{\ell m\omega }^{up}{T_{\ell m\omega }}{\Delta ^{ - 2}} \equiv \tilde Z_{\ell m\omega }^H{\Delta ^2}{e^{ - ikr*}},$$

((24))

while the asymptotic behavior at infinity is

$${R_{\ell m\omega }}(r \to \infty ) \to \frac{{{r^3}{e^{i\omega r*}}}}{{2i\omega B_{\ell m\omega }^{inc}}}\int_{{r_ + }}^\infty {dr'} \frac{{{T_{\ell m\omega }}(r')R_{\ell m\omega }^{in}(r')}}{{{\Delta ^2}(r')}} \equiv \tilde Z_{\ell m\omega }^\infty {r^3}{e^{i\omega r*}}.$$

((25))

We note that the homogeneous Teukolsky equation is invariant under the complex conjugation followed by the transformation m → -m and ω → -ω. Thus, we can set $\bar R_{\ell m\omega }^{\text{in},\;\text{up}}$, where the bar denotes the complex conjugation.

We consider T_μν of a monopole particle of mass μ. The energy momentum tensor takes the form

$${T^{\mu v}} = \frac{\mu }{{\Sigma \sin \theta dt/d\tau }}\frac{{d{z^\mu }}}{{d\tau }}\frac{{d{z^v}}}{{d\tau }}\delta (r - r(t))\delta (\theta - \theta (t))\delta (\varphi - \varphi (t)),$$

((26))

where z^μ = (t, r(t),θ(t),ϕ(t)) is a geodesic trajectory, and τ = τ(t) is the proper time along the geodesic. The geodesic equations in the Kerr geometry are given by

$$\begin{gathered} \Sigma \frac{{d\theta }}{{d\tau }} = \pm {\left[ {\hat C - {{\cos }^2}\theta \left( {{a^2}(1 - {{\hat {\mathcal{E}}}^2}) + \frac{{\hat l_z^2}}{{{{\sin }^2}\theta }}} \right)} \right]^{1/2}} \equiv \Theta (\theta )|, \hfill \\ \Sigma \frac{{d\varphi }}{{d\tau }} = \pm \left( {a\hat {\mathcal{E}} - \frac{{{{\hat l}_z}}}{{{{\sin }^2}\theta }}} \right) + \frac{a}{\Delta }\left( {\hat {\mathcal{E}}({r^2} + {a^2}) - a{{\hat l}_z}} \right) \equiv \Phi , \hfill \\ \Sigma \frac{{dt}}{{d\tau }} = \pm \left( {a\hat {\mathcal{E}} - \frac{{{{\hat l}_z}}}{{{{\sin }^2}\theta }}} \right) + a{\sin ^2}\theta + \frac{{{r^2} + {a^2}}}{\Delta }\left( {\hat {\mathcal{E}}({r^2} + {a^2}) - a{{\hat l}_z}} \right) \equiv T, \hfill \\ \Sigma \frac{{dr}}{{d\tau }} = \pm \pm \sqrt R , \hfill \\ \end{gathered} $$

((27))

where

$$R = {[\hat {\mathcal{E}}({r^2} + {a^2}) - a{\hat l_z}]^2} - \Delta {[\hat {\mathcal{E}}a - {\hat l_z})^2} + {r^2} + \hat C].$$

((28))

and $\hat {\mathcal {E}},\;\hat {\mathcal {l}}_z, \text{and}\;\hat {\mathcal {C}}$ are the energy, the z-component of the angular momentum, and the Carter constant of a test particle, respectively. These constants of motion are those measured in units of μ. That is, if expressed in the standard units, they become $\mathcal{E} = \mu \hat {\mathcal{E}},{l_z} = \mu {\hat l_z}$, and $C = {\mu ^2}\hat C$.

Using Equation (27), the tetrad components of the energy momentum tensor are expressed as where

$$\begin{gathered} {T_{nn}} = \mu \frac{{{C_{nn}}}}{{\sin \theta }}\delta (r - r(t))\delta (\theta - \theta (t))\delta (\varphi - \varphi (t)), \hfill \\ {T_{\bar mn}} = \mu \frac{{{C_{\bar mn}}}}{{\sin \theta }}\delta (r - r(t))\delta (\theta - \theta (t))\delta (\varphi - \varphi (t)), \hfill \\ {T_{\overline {mm} }} = \mu \frac{{{C_{\overline {mm} }}}}{{\sin \theta }}\delta (r - r(t))\delta (\theta - \theta (t))\delta (\varphi - \varphi (t)), \hfill \\ \end{gathered} $$

((29))

where

$${C_{nn}} = \frac{1}{{4{\Sigma ^3}\dot t}}{\left[ {\hat {\mathcal{E}}({r^2} + {a^2}) - a{{\hat l}_z} + \Sigma \frac{{dr}}{{d\tau }}} \right]^2},$$

((30))

$${C_{\bar mn}} = \frac{\rho }{{2\sqrt 2 {\Sigma ^2}\dot t}}\left[ {\hat {\mathcal{E}}({r^2} + {a^2}) - a{{\hat l}_z} + \Sigma \frac{{dr}}{{d\tau }}} \right]\left[ {i\sin \theta \left( {a\hat {\mathcal{E}} - \frac{{{{\hat l}_z}}}{{{{\sin }^2}\theta }}} \right) + \Sigma \frac{{d\theta }}{{d\tau }}} \right],$$

((31))

$${C_{\overline {mm} }} = \frac{{\rho 2}}{{2\Sigma \dot t}}{\left[ {i\sin \theta \left( {a\hat {\mathcal{E}} - \frac{{{{\hat l}_z}}}{{{{\sin }^2}\theta }}} \right) + \Sigma \frac{{d\theta }}{{d\tau }}} \right]^2},$$

((32))

and $\dot t = dt/d\tau$. Substituting Equation (10) into Equation (18) and performing integration by part, we obtain

$$\begin{gathered} {T_{\ell m\omega }} = \frac{{4\mu }}{{\sqrt {2\pi } }}\int_{ - \infty }^\infty {dt\int {d\theta {e^{i\omega t - im\varphi (t)}}} } \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \times \left\{ { - \frac{1}{2}L_1^\dag ({\rho ^{ - 4}}L_2^\dag ({\rho ^3}S)){C_{nn}}{\rho ^{ - 2}}{{\bar \rho }^{ - 1}}\delta (r - r(t))\delta (\theta - \theta (t))} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \frac{{{\Delta ^2}{{\bar \rho }^2}}}{{\sqrt {2\rho } }}\left( {L_2^\dag S + ia(\bar \rho - \rho )\sin \theta S} \right){J_ + }\left[ {\frac{{{C_{\bar mn}}}}{{{\rho ^2}{{\bar \rho }^2}\Delta }}\delta (r - r(t))\delta (\theta - \theta (t))} \right]\;\;\; \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \frac{1}{{2\sqrt 2 }}\left( {L_2^\dag ({\rho ^3}S({{\bar \rho }^2}{\rho ^{ - 4}}){,_r}} \right){C_{\bar mn}}\Delta {\rho ^{ - 2}}{{\bar \rho }^2}\delta (r - r(t))\delta (\theta - \theta (t)) \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \frac{1}{4}\;{\rho ^3}{\Delta ^2}S{J_ + }({{\bar \rho }^4}{J_ + }(\bar \rho {\rho ^{ - 2}}{C_{\overline {mm} }}\delta (r - r(t))\delta (\theta - \theta (t))]\} , \hfill \\ \end{gathered} $$

((33))

where

$$L_s^\dag = {\partial _\theta } - \frac{m}{{\sin \theta }} + a\omega \sin \theta + s\cot \theta ,$$

((34))

$${J_ + } = {\partial _r} + iK/\Delta ,$$

((35))

and S denotes $_{ - 2}S_{\ell m}^{a\omega }(\theta )$ for simplicity.

For a source bounded in a finite range of r, it is convenient to rewrite Equation (33) further as

$$\begin{gathered} {T_{\ell m\omega }} = \mu \int_{ - \infty }^\infty {dt{e^{i\omega t - im\varphi (t)}}{\Delta ^2}\{ ({A_{nn0}}} + {A_{\bar mn0}} + {A_{\overline {mm} 0}})\delta (r - r(t)) \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + {[({A_{\bar mn1}} + {A_{\overline {mm} 1}})\delta (r - r(t))]_{,r}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + [{({A_{\overline {mm} 2}}\delta (r - r(t))]_{,rr}}\} , \hfill \\ \end{gathered} $$

((36))

where

$${A_{nn0}} = \frac{{ - 2}}{{\sqrt {2\pi {\Delta ^2}} }}{\rho ^{ - 2}}{\bar \rho ^{ - 1}}{C_{nn}}L_1^\dag [{\rho ^{ - 4}}L_2^\dag ({\rho ^3}S)],$$

((37))

$${A_{\bar mn0}} = \frac{2}{{\sqrt {\pi \Delta } }}{\rho ^{ - 3}}{C_{\bar mn}}[(L_2^\dag S)\left( {\frac{{iK}}{\Delta } + \rho + \bar \rho } \right) - a\sin \theta S\frac{K}{\Delta }(\bar \rho - \rho )],$$

((38))

$${A_{\overline {mm} 0}} = - \frac{1}{{\sqrt {2\pi } }}{\rho ^{ - 3}}\bar \rho {C_{\overline {mm} }}S\left[ { - i{{\left( {\frac{K}{\Delta }} \right)}_{,r}} - \frac{{{K^2}}}{{{\Delta ^2}}} + 2i\rho \frac{K}{\Delta }} \right],$$

((39))

$${A_{\bar mn1}} = \frac{2}{{\sqrt {\pi \Delta } }}{\rho ^{ - 3}}{C_{\bar mn}}\left[ {L_2^\dag S + ia\sin \theta (\bar \rho - \rho )S} \right],$$

((40))

$${A_{\overline {mm} 1}} = - \frac{2}{{\sqrt {2\pi } }}{\rho ^{ - 3}}\bar \rho {C_{\overline {mm} }}S\left( {i\frac{K}{\Delta } + \rho } \right),$$

((41))

$${A_{\overline {mm} 2}} = - \frac{2}{{\sqrt {2\pi } }}{\rho ^{ - 3}}\bar \rho {C_{\overline {mm} }}S.$$

((42))

Inserting Equation (36) into Equation (25), we obtain ${\tilde Z_{\ell m\omega }}$ as

$${\tilde Z_{\ell m\omega }} = \frac{\mu }{{2i\omega B_{\ell m\omega }^{inc}}}{\int_{ - \infty }^\infty {dt{e^{i\omega t - im\varphi (t)}}W} _{\ell m\omega ,}}$$

((43))

where

$${W_{\ell m\omega }} = {\left\{ {R_{\ell m\omega }^{in}[{A_{nn0}} + {A_{\bar mn0}} + {A_{\overline {mm} 0}}] - \frac{{dR_{\ell m\omega }^{in}}}{{dr}}[{A_{\bar mn1}} + {A_{\overline {mm} 1}}] + \frac{{{d^2}R_{\ell m\omega }^{in}}}{{d{r^2}}}{A_{\overline {mm} 2}}} \right\}_{r = r(t)}}.$$

((44))

In this paper, we focus on orbits which are either circular (with or without inclination) or eccentric but confined on the equatorial plane. In either case, the frequency spectrum of ${T_{\ell m\omega }}$ becomes discrete. Accordingly, ${\tilde Z_{\ell m\omega }}$ in Equation (24) or (25) takes the form,

$${\tilde Z_{\ell m\omega }} = \mathop \Sigma \limits_n \delta (\omega - {\omega _n}){Z_{\ell m\omega }}.$$

((45))

Then, in particular, ψ₄ at r → ∞ is obtained from Equation (13) as

$${\psi _4} = \frac{1}{2}\mathop \Sigma \limits_{\ell mn} {Z_{\ell m{\omega _n}}}\frac{{ - 2S_{\ell m}^{a{\omega _n}}}}{{\sqrt {2\pi } }}{e^{i{\omega _n}(r* - t) + im\varphi }}.$$

((46))

At infinity, ψ₄ is related to the two independent modes of gravitational waves h₊ and h_× as

$${\psi _4} = \frac{1}{2}({\ddot h_ + } - i{\ddot h_ \times }).$$

((47))

From Equations (46) and (47), the luminosity averaged over t ≫ Δt, where Δt is the characteristic time scale of the orbital motion (e.g., a period between the two consecutive apastrons), is given by

$$\left\langle {\frac{{dE}}{{dt}}} \right\rangle = \sum\limits_{\ell ,m,n} {\frac{{|{Z_{\ell m{\omega _n}}}{|^2}}}{{4\pi \omega _n^2}} \equiv } {\sum\limits_{\ell ,m,n} {\left( {\frac{{dE}}{{dt}}} \right)} _{\ell mn}}.$$

((48))

In the same way, the time-averaged angular momentum flux is given by

$$\left\langle {\frac{{d{J_z}}}{{dt}}} \right\rangle = \sum\limits_{\ell ,m,n} {\frac{{m|{Z_{\ell m{\omega _n}}}{|^2}}}{{4\pi \omega _n^3}} \equiv } {\sum\limits_{\ell ,m,n} {\left( {\frac{{d{J_z}}}{{dt}}} \right)} _{\ell mn}} = \sum\limits_{\ell ,m,n} {\frac{m}{{{\omega _n}}}{{\left( {\frac{{dE}}{{dt}}} \right)}_{\ell mn}}} .$$

((49))

2.2 Chandrasekhar-Sasaki-Nakamura transformation

As seen from the asymptotic behaviors of the radial functions given in Equations (24) and (25), the Teukolsky equation is not in the form of a canonical wave equation near the horizon and infinity. Therefore, it is desirable to find a transformation that brings the radial Teukolsky equation into the form of a standard wave equation.

In the Schwarzschild case, Chandrasekhar found that the Teukolsky equation can be transformed to the Regge-Wheeler equation, which has the standard form of a wave equation with solutions having regular asymptotic behaviors at horizon and infinity [21]. The Regge-Wheeler equation was originally derived as an equation governing the odd parity metric perturbation [87]. The existence of this transformation implies that the Regge-Wheeler equation can describe the even parity metric perturbation simultaneously, though the explicit relation of the Regge-Wheeler function obtained by the Chandrasekhar transformation with the actual metric perturbation variables has not been given in the literature yet.

Later, Sasaki and Nakamura succeeded in generalizing the Chandrasekhar transformation to the Kerr case [92, 93]. The Chandrasekhar-Sasaki-Nakamura transformation was originally introduced to make the potential in the radial equation short-ranged, and to make the source term well-behaved at the horizon and at infinity. Since we are interested only in bound orbits, it is not necessary to perform this transformation. Nevertheless, because its flat-space limit reduces to the standard radial wave equation in the Minkowski spacetime, it is convenient to apply the transformation when dealing with the post-Minkowski or post-Newtonian expansion, at least at low orders of expansion.

We transform the homogeneous Teukolsky equation to the Sasaki-Nakamura equation [92, 93] which is given by

$$\left( {\frac{{{d^2}}}{{dr{*^2}}} - F(r)\frac{d}{{dr*}} - U(r)} \right){X_{\ell m\omega }} = 0.$$

((50))

The function F(r) is given by

$$F(r) = \frac{{\eta ,r}}{\eta }\frac{\Delta }{{{r^2} + {a^2}}},$$

((51))

where

$$\eta = {c_0} + \frac{{{c_1}}}{r} + \frac{{{c_2}}}{{{r^2}}} + \frac{{{c_3}}}{{{r^3}}} + \frac{{{c_4}}}{{{r^4}}},$$

((52))

With

$$\begin{gathered} {c_0} = 12i\omega M + \lambda (\lambda + 2) - 12a\omega (a\omega - m), \hfill \\ {c_1} = 8ia[3a\omega - \lambda ](a\omega - m), \hfill \\ {c_2} = 24iaM(a\omega - m) + 12a{}^2[1 - 2{(a\omega - m)^2}], \hfill \\ {c_3} = 24i{a^3}(a\omega - m) - 24M{a^2}, \hfill \\ {c_4} = 12{a^4}. \hfill \\ \end{gathered} $$

((53))

The Function U(r) is given by

$$U(r) = \frac{{\Delta {U_1}}}{{{{({r^2} + {a^2})}^2}}} + {G^2} + \frac{{\Delta {G_{,r}}}}{{{r^2} + {a^2}}} - FG,$$

((54))

where

$$G = \frac{{2(r - M)}}{{({r^2} + {a^2})}} + \frac{{r\Delta }}{{{{({r^2} + {a^2})}^2}}},$$

((55))

$${U_1} = V + \frac{{{\Delta ^2}}}{\beta }\left[ {{{\left( {2\alpha \frac{{{\beta _{,r}}}}{\Delta }} \right)}_{,r}} - \frac{{{\eta _{,r}}}}{\eta }\left( {\alpha + \frac{{{\beta _{,r}}}}{\Delta }} \right)} \right],$$

((56))

$$\alpha = - i\frac{{K\beta }}{{{\Delta ^2}}} + 3i{K_{,r}} + \lambda + \frac{{6\Delta }}{{{r^2}}},$$

((57))

$$\beta = 2\Delta \left( { - iK + r - M - \frac{{2\Delta }}{r}} \right).$$

((58))

The relation between ${R_{\ell m\omega }}$, and ${X_{\ell m\omega }}$, is

$${R_{\ell m\omega }} = \frac{1}{\eta }\left[ {\left( {\alpha + \frac{{{B_{,r}}}}{\Delta }} \right){\chi _{\ell m\omega }} - \frac{\beta }{\Delta }{\chi _{\ell m\omega ,r}}} \right],$$

((59))

where ${\chi _{\ell m\omega }}\Delta /{({r^2} + {a^2})^{1/2}}$. Conversely, we can express ${X_{\ell m\omega }}$, in terms of ${R_{\ell m\omega }}$, as

$${X_{\ell m\omega }} = {({r^2} + {a^2})^{1/2}}{r^2}{J_ - }{J_ - }\left( {\frac{1}{{{r^2}}}{R_{\ell m\omega }}} \right),$$

((60))

where J_ (d/dr) - i(K/Δ).

If we set a = 0, this transformation reduces to the Chandrasekhar transformation for the Schwarzschild black hole [21]. The explicit form of the transformation is

$${R_{\ell m\omega }} = \frac{\Delta }{{{c_0}}}\left( {\frac{d}{{dr*}} + i\omega } \right)\frac{{{r^2}}}{\Delta }\left( {\frac{d}{{dr*}} + i\omega } \right)r{X_{\ell m\omega }},$$

((61))

$${X_{\ell m\omega }} = \frac{{{r^5}}}{\Delta }\left( {\frac{d}{{dr*}} - i\omega } \right)\frac{{{r^2}}}{\Delta }\left( {\frac{d}{{dr*}} + i\omega } \right)\frac{{{R_{\ell m\omega }}}}{{{r^2}}},$$

((62))

where c₀, defined in Equation (53), reduces to ${c_0} = (\ell - 1)\ell (\ell + 1)(\ell + 2) - 12iM\omega$. In this case, the Sasaki-Nakamura equation (50) reduces to the Regge-Wheeler equation [87], which is given by

$$\left( {\frac{{{d^2}}}{{dr{*^2}}} + {\omega ^2} - V(r)} \right){X_{\ell \omega }}(r) = 0,$$

((63))

where

$$V(r) - \left( {1 - \frac{{2M}}{r}} \right)\left( {\frac{{\ell (\ell + 1)}}{{{r^2}}} - \frac{{6M}}{{{r^3}}}} \right).$$

((64))

As is clear from the above form of the equation, the lowest order solutions are given by the spherical Bessel functions. Hence it is intuitively straightforward to apply the post-Newtonian expansion to it. Some useful techniques for the post-Newtonian expansion were developed for the Schwarzschild case by Poisson [83] and Sasaki [91].

The asymptotic behavior of the ingoing wave solution Xⁱⁿ which corresponds to Equation (19) is

$$X_{\ell m\omega }^{in} \to \left\{ \begin{gathered} A_{\ell m\omega }^{ref}{e^{i\omega r*}} + A_{\ell m\omega }^{inc}{e^{ - i\omega r*}}\;\;\;\;\;\;\;\;\;\;for\;r* \to \infty , \hfill \\ A_{\ell m\omega }^{trans}{e^{ - ikr*}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;for\;r* \to - \infty . \hfill \\ \end{gathered} \right.$$

((65))

The coefficients A^inc, A^ref, and A^trans are related to B^inc, B^ref, and B^trans, defined in Equation (19), by

$$B_{\ell m\omega }^{inc} = - \frac{1}{{4{\omega ^2}}}A_{\ell m\omega }^{inc},$$

((66))

$$B_{\ell m\omega }^{ref} = - \frac{{4{\omega ^2}}}{{{c_0}}}A_{\ell m\omega }^{ref},$$

((67))

$$B_{\ell m\omega }^{trans} = \frac{1}{{{d_{\ell m\omega }}}}A_{\ell m\omega }^{trans},$$

((68))

where

$$\begin{gathered} {d_{\ell m\omega }} = \sqrt {2M{r_ + }} [(8 - 24iM\omega - 16{M^2}{\omega ^2})r_ + ^2 \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + (12iam - 16M + 16amM\omega + 24i{M^2}\omega ){r_ + } \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - 4{a^2}{m^2} - 12iamM + 8{M^2}]. \hfill \\ \end{gathered} $$

((69))

In the following sections, we present a method of post-Newtonian expansion based on the above formalism in the case of the Schwarzschild background. In the Kerr case, although a post-Newtonian expansion method developed in previous work [94, 101] was based on the Sasaki-Nakamura equation, we will not present it in this paper. Instead, we present a different formalism, namely the one developed by Mano, Suzuki, and Takasugi which allows us to solve the Teukolsky equation in a more systematic manner, albeit very mathematical [68]. The reason is that the equations in the Kerr case are already complicated enough even if one uses the Sasaki-Nakamura equation, so that there is not much advantage in using it. In contrast, in the Schwarzschild case, it is much easier to obtain physical insight into the role of relativistic corrections if we deal with the Regge-Wheeler equation.

3 Post-Newtonian Expansion of the Regge-Wheeler Equation

In this section, we review a post-Newtonian expansion method for the Schwarzschild background, based on the Regge-Wheeler equation. We focus on the gravitational waves emitted to infinity, but not on those absorbed by the black hole. The black hole absorption is deferred to Section 4, in which we review the Mano-Suzuki-Takasugi method for solving the Teukolsky equation.

Since we are interested in the waves emitted to infinity, as seen from Equation (25), what we need is a method to evaluate the ingoing wave Teukolsky function $R_{\ell m\omega }^{\text{in}}$ or its counterpart in the Regge-Wheeler equation, $X_{\ell m\omega }^{\text{in}}$, which are related by Equation (59). In addition, we assume ω > 0 whenever it is necessary throughout this section. Formulae and equations for ω < 0 are obtained from the symmetry $\bar X_{\ell m\omega }^{{\text{in}}} = X_{\ell - m - \omega }^{\text{in}}$.

3.1 Basic assumptions

We consider the case of a test particle with mass μ in a nearly circular orbit around a black hole with mass M ≫ μ. For a nearly circular orbit, say at r ∼ r₀, what we need to know is the behavior of $R_{\ell m\omega }^{\text{in}}\;\text{at}\;r \sim {r_0}$. In addition, the contribution of ω to $R_{\ell m\omega }^{\text{in}}$ comes mainly from ω ∼ mΩ_ϕ, where Ω_ϕ ∼ (M/r ³₀ )^1/2 is the orbital angular frequency.

Thus, if we express the Regge-Wheeler equation (63) in terms of a non-dimensional variable z ≡ ωr, with a non-dimensional parameter ε ≡ 2Mω, we are interested in the behavior of $X_{\ell m\omega }^{\text{in}} (z)$ at z ∼ ωr₀ ∼ m(M/r₀)^1/2 ∼ v with ε ∼ 2m(M/r₀)^3/2 ∼ v³, where v ≡ (M/r₀)^1/2 is the characteristic orbital velocity. The post-Newtonian expansion assumes that v is much smaller than the velocity of light: v ≪ 1. Consequently, we have ε ≪ v ≪ 1 in the post-Newtonian expansion.

To obtain $X_{\ell m\omega }^{\text{in}}$ (which we denote below by $X_{\ell}$ for simplicity) under these assumptions, we find it convenient to rewrite the Regge-Wheeler equation in an alternative form. It is

$$\left[ {\frac{{{d^2}}}{{d{z^2}}} + \frac{{2d}}{{zdz}} + \left( {1 - \frac{{\ell (\ell + 1)}}{{{z^2}}}} \right)} \right]{\xi _\ell }(z) = \in {e^{ - iz}}\frac{d}{{dz}}\left[ {\frac{1}{{{z^3}}}\frac{d}{{dz}}({e^{iz}}{z^2}{\xi _\ell }(z))} \right],$$

((70))

where ${\xi _\ell }$ is a function related to ${X _\ell }$ as

$${X_\ell } = z{e^{ - i \in \ln (z - \in )}}{\xi _\ell }.$$

((71))

The ingoing wave boundary condition of ${\xi _\ell }$ is derived from Equations (65) and (71) as

$${X_\ell } \to \left\{ \begin{gathered} A_\ell ^{inc}{e^{i \in \ln \in }}{z^{ - 1}}{e^{ - iz}} + A_\ell ^{ref}{e^{ - i \in \ln \in }}{z^{ - 1}}{e^{i(z + 2 \in \ln z)}}\;\;\;\;\;\;\;\;\;\;for\;r* \to \infty , \hfill \\ A_\ell ^{trans}{ \in ^{ - 1}}{e^{i \in (\ln \in - 1)}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;for\;r* \to \infty . \hfill \\ \end{gathered} \right.$$

((72))

The above form of the Regge-Wheeler equation is used in Sections 3.2, 3.3, 3.4, and 3.5.

It should be noted that if we reinstate the gravitational constant G, we have ε = 2GMω. Thus, the expansion in terms of e corresponds to the post-Minkowski expansion, and expanding the Regge-Wheeler equation with the assumption ε ≪ 1 gives a set of iterative wave equations on the flat spacetime background. One of the most significant differences between the black hole perturbation theory and any theory based on the flat spacetime background is the presence of the black hole horizon in the former case. Thus, if we naively expand the Regge-Wheeler equation with respect to e, the horizon boundary condition becomes unclear, since there is no horizon on the flat spacetime. To establish the boundary condition at the horizon, we need to treat the Regge-Wheeler equation near the horizon separately. We thus have to find a solution near the horizon, and the solution obtained by the post-Minkowski expansion must be matched with it in the region where both solutions are valid.

It may be of interest to note the difference between the matching used in the BDI approach for the post-Newtonian expansion [7, 12] and the matching used here. In the BDI approach, the matching is done between the post-Minkowskian metric and the near-zone post-Newtonian metric. In our case, the matching is done between the post-Minkowskian gravitational field and the gravitational field near the black hole horizon.

3.2 Horizon solution; z ≪ 1

In this section, we first consider the solution near the horizon, which we call the horizon solution, based on [85]. To do so, we assume z ≪ 1 and treat ε as a small number, but leave the ratio z/ε arbitrary. We change the independent variable to x = 1 - z/ε and the wave function to

$$Z = {\left( {\frac{ \in }{z}} \right)^3}\frac{{{X_\ell }}}{{A_\ell ^{trans}}} = {\left( {\frac{ \in }{z}} \right)^2}\frac{{ \in {\xi _\ell }}}{{A_\ell ^{trans}{e^{i \in (\ln \in - 1)}}}}.$$

((73))

Note that the horizon corresponds to x = 0. We then have

$$x(x - 1)Z'' + [2(3 - i \in )x - (1 - 2i \in )]Z' + [6 - \ell (\ell + 1) - 5i \in + { \in ^2}(2 - 3x + {x^2})]Z = 0,$$

((74))

where a prime denotes differentiation with respect to x. We look for a solution which is regular at x = 0.

First, we consider the lowest order solution by setting ε = 0 in Equation (74). The boundary condition (72) requires that Z = 1 at x = 0. The solution that satisfies the boundary condition is

$$Z = \sum\limits_{n = 0}^{\ell - 2} {\frac{{{{(2 - \ell )}_n}{{(\ell + 3)}_n}}}{{n!}}{x^n},\;\;\;\;\;\;{{(a)}_n} = \frac{{\Gamma (a + n)}}{{\Gamma (a)}}.} $$

((75))

Thus, the lowest order solution is a polynomial of order ℓ - 2 in x = 1 - z/ε.

Next, we consider the solution accurate to $\mathcal{O}(\epsilon)$. We neglect the terms of $\mathcal{O}(\epsilon^2)$ in Equation (74). Then, the wave equation takes the form of a hypergeometric equation,

$$x(x - 1)Z'' + [(a + b + 1)x - c]Z' + abZ = 0,$$

((76))

with parameters

$$\begin{gathered} a = - (\ell - 2) - i \in = \mathcal{O}({ \in ^2}), \hfill \\ b = \ell + 3 - i \in + \mathcal{O}({ \in ^2}), \hfill \\ c = 1 - 2i \in . \hfill \\ \end{gathered} $$

((77))

The two linearly independent solutions are F(a, b; c; x) and x^1-cF(a+1-c, b+1-c; 2-c; x), where F is the hypergeometric function. However, only the first solution is regular at x = 0. Therefore, we obtain

$${\xi _\ell }(z) = A_\ell ^{trans}{ \in ^{ - 1}}{e^{i \in (\ln \in - 1)}}{\left( {\frac{z}{ \in }} \right)^2}F\left( {a,b;c;1 - \frac{z}{ \in }} \right).$$

((78))

The above solution must be matched with the solution obtained from the post-Minkowski expansion of Equation (70), which we call the outer solution, in a region where both solutions are valid. It is the region where the post-Newtonian expansion is applied, i.e., the region ε ≪ z ≪ 1. For this purpose, we rewrite Equation (78) as (see, e.g., Equation (15.3.8) of [1])

$$\begin{gathered} {\xi _\ell } = A_\ell ^{trans}{ \in ^{ - 1}}{e^{i \in (\ln \in - 1)}}\left[ {{{\left( {\frac{z}{ \in }} \right)}^{\ell + i \in }}\frac{{\Gamma (c)\Gamma (b - a)}}{{\Gamma (b)\Gamma (c - a)}}F\left( {a,c - b;a - b + 1;\frac{ \in }{z}} \right)} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. { + {{\left( {\frac{z}{ \in }} \right)}^{ - \ell - 1 + i \in }}\frac{{\Gamma (c)\Gamma (b - a)}}{{\Gamma (a)\Gamma (c - )}}F\left( {b,c - a;b - a + 1;\frac{ \in }{z}} \right)} \right]. \hfill \\ \end{gathered} $$

((79))

This naturally allows the expansion in ε/z. It should be noted that the second term in the square brackets of the above expression is meaningless as it is, since the factor Γ(a - b) diverges for integer ℓ. So, when evaluating the second term, we first have to extend ℓ to a non-integer number. Then, only after expanding it in terms of ε, we should take the limit of an integer ℓ. One then finds that this procedure gives rise to an additional factor of $\mathcal{O}(\epsilon)$. For ε/z ≪ 1, it therefore becomes $\mathcal{O}(\epsilon^{2\ell+2})$ higher in ε than the first term. Then, we obtain

$${\xi _\ell }( \in \ll z \ll 1) = \frac{{(2\ell )!}}{{(\ell - 2)!(\ell + 2)!}}\frac{{{z^\ell }}}{{{ \in ^{\ell + 1}}}}\left[ {1 + i \in ({a_\ell } + \ln z) - \frac{{(\ell - 2)(\ell + 2)}}{{2\ell }}\frac{ \in }{z} + \mathcal{O}({ \in ^2})} \right],$$

((80))

where

$${a_\ell } = 2\gamma + \psi (\ell - 1) + \psi (\ell + 3) - 1,$$

((81))

and ψ(n) is the digamma function,

$$\psi (n) = - \gamma + \sum\limits_{k = 1}^{n - 1} {{k^{ - 1}},} $$

((82))

and γ ⋍ 0.57721 is the Enter constant.

As we will see below, the above solution is accurate enough to determine the boundary condition of the outer solution up to the 6PN order of expansion.

3.3 Outer solution; ε ≪ 1

We now solve Equation (70) in the limit ε ≪ 1, i.e., by applying the post-Minkowski expansion to it. In this section, we consider the solution to $\mathcal{O}(\epsilon)$. Then we match the solution to the horizon solution given by Equation (80) at ε ≪ z ≪ 1.

By setting

$${\xi _\ell }(z) = \sum\limits_{n = 0}^\infty {{ \in ^n}} \xi _\ell ^{(n)}(z),$$

((83))

each ξ ⁽ⁿ⁾_ℓ (z) is found to satisfy

$$\left[ {\frac{{{d^2}}}{{d{z^2}}} + \frac{2}{z}\frac{d}{{dz}} + \left( {1 - \frac{{\ell (\ell + 1)}}{{{z^2}}}} \right)} \right]\xi _\ell ^{(n)} = {e^{ - iz}}\frac{d}{{dz}}\left[ {\frac{1}{{{z^3}}}\frac{d}{{dz}}\left( {{e^{iz}}{z^2}\xi _\ell ^{(n - 1)}(z)} \right)} \right].$$

((84))

Equation (84) is an inhomogeneous spherical Bessel equation. It is the simplicity of this equation that motivated the introduction of the auxiliary function ξ_ℓ [91].

The wroth-order solution ξ ⁽⁰⁾_ℓ satisfies the homogeneous spherical Bessel equation, and must be a linear combination of the spherical Bessel functions of the first and second kinds, j_ℓ(z) ∼ z^ℓ and n_ℓ(z). Here, we demand the compatibility with the horizon solution (80). Since j_ℓ(z) ∼ z^ℓ and n_ℓ(z) ∼ z^−ℓ−1, n_ℓ(z) does not match with the horizon solution at the leading order of ε. Therefore, we have

$$\xi _\ell ^{(0)}(z) = \alpha _\ell ^{(0)}{j_\ell }(z).$$

((85))

The constant α ⁽⁰⁾_ℓ represents the overall normalization of the solution. Since it can be chosen arbitrarily, we set α ⁽⁰⁾_ℓ = 1 below.

The procedure to obtain ξ ⁽¹⁾_ℓ (z) was described in detail in [91]. Using the Green function G(z, z′) = j_ℓ(z_<)n_ℓ(z_>), Equation (84) may be put into an indefinite integral form,

$$\xi _\ell ^{(n)} = {n_\ell }{\smallint ^z}dz{z^2}{e^{ - iz}}{j_\ell }{\left[ {\frac{1}{{{z^3}}}({e^{iz}}{z^2}\xi _\ell ^{(n - 1)}(z))'} \right]^\prime } - {j_\ell }{\smallint ^z}dz{z^2}{e^{ - iz}}{n_\ell }{\left[ {\frac{1}{{{z^3}}}({e^{iz}}{z^2}\xi _\ell ^{(n - 1)}(z))'} \right]^\prime }.$$

((86))

The calculation is tedious but straightforward. All the necessary formulae to obtain ξ ⁽ⁿ⁾_ℓ for n ≤ 2 are given in the Appendix of [91] or Appendix D of [71]. Using those formulae, for n = 1 we have

$$\begin{gathered} \xi _\ell ^{(1)} = \alpha _\ell ^{(1)}{j_\ell } + \beta _\ell ^{(1)}{n_\ell } \hfill \\ \;\;\;\;\;\;\;\; + \frac{{(\ell - 1)(\ell + 3)}}{{2(\ell + 1)(2\ell + 1)}}{j_{\ell + 1}} - \leqslant \left( {\frac{{{\ell ^2} - 4}}{{2\ell (2\ell + 1)}} + \frac{{2\ell - 1}}{{\ell (\ell - 1)}}} \right){j_{\ell - 1}} \hfill \\ \;\;\;\;\;\;\;\; + {R_{\ell ,0}}{j_0} + \sum\limits_{m = 1}^{\ell - 2} {\left( {\frac{1}{m} + \frac{1}{{m + 1}}} \right)} {R_{\ell ,m}}{j_m} - 2D_\ell ^{nj} + i{j_\ell }\ln z. \hfill \\ \end{gathered} $$

((87))

Here, D ^nj_ℓ and R_ℓ,m are functions defined as follows. The function D ^nj_ℓ is given by

$$D_\ell ^{nj} = \frac{1}{2}[{j_\ell }Si(2z) - {n_\ell }(Ci(2z) - \lambda - \ln 2z)],$$

((88))

where $Ci(x) = - \int_x^\infty {dt\cos t/t}$ and $Si(x) = \int_0^x {dt\sin t/t}$. The function R_m,k is defined by R_m,k = z²(n_mj_k - j_mn_k). It is a polynomial in inverse powers of z given by

$${R_{m,k}}\left\{ \begin{gathered} - \sum\limits_{r = 0}^{\frac{1}{2}(m - k - 1)} {{{( - 1)}^r}} \frac{{\Gamma (m - k - r)\Gamma \left( {m + \frac{1}{2} - r} \right)}}{{r!\Gamma (m - k - 2r)\Gamma \left( {k + \frac{3}{2} + r} \right)}}{\left( {\frac{2}{z}} \right)^{m - k - 1 - 2r}}\;\;\;\;for\;m > k, \hfill \\ - {R_{m,k}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;for\;m < k. \hfill \\ \end{gathered} \right.$$

((89))

Here, we again perform the matching with the horizon solution (80). It should be noted that ξ ⁽¹⁾_ℓ , given by Equation (87), is regular in the limit z → 0 except for the term β ⁽¹⁾_ℓ n_ℓ. By examining the asymptotic behavior of Equation (87) at z ≪ 1, we find β ⁽¹⁾_ℓ = 0, i.e., the solution is regular at z = 0. As for α ⁽¹⁾_ℓ , it only contributes to the renormalization of α ⁽⁰⁾_ℓ . Hence, we set α ⁽¹⁾_ℓ = 0 and the transmission amplitude A ^trans_ℓ is determined to $\mathcal{O}(\epsilon)$ as

$$A_\ell ^{trans} = \frac{{(\ell - 2)!(\ell + 2)!}}{{(2\ell )!(2\ell + 1)!}}{ \in ^{\ell + 1}}[1 - i \in {a_\ell } + \mathcal{O}({ \in ^2})].$$

((90))

It may be noted that this explicit expression for A ^trans_ℓ is unnecessary for the evaluation of gravitational waves at infinity. It is relevant only for the evaluation of the black hole absorption.

3.4 More on the inner boundary condition of the outer solution

In this section, we discuss the inner boundary condition of the outer solution in more detail. As we have seen in Section 3.3, the boundary condition on ξ_ℓ is that it is regular at z → 0, at least to $\mathcal{O}(\epsilon)$, while in the full non-linear level, the horizon boundary is at z = ε. We therefore investigate to what order in ε the condition of regularity at z = 0 can be applied.

Let us consider the general form of the horizon solution. With x = 1 - z/ε, it is expanded in the form

$${\xi _\ell } = \xi _\ell ^{\{ 0\} }(x) + \in \xi _\ell ^{\{ 1\} }(x) + { \in ^2}\xi _\ell ^{\{ 2\} }(x) + ...$$

((91))

The lowest order solution ξ ^{0}_ℓ (x) is given by the polynomial (75). Apart from the common overall factor, it is schematically expressed as

$$\xi _\ell ^{\{ 0\} } = {\left( {\frac{z}{ \in }} \right)^\ell }\left[ {1 + {c_1}\frac{ \in }{z} + ... + {c_{\ell - 2}}{{\left( {\frac{ \in }{z}} \right)}^{\ell - 2}}} \right].$$

((92))

Thus, ξ ^{0}_ℓ does not have a term matched with n_ℓ, but it matches with j_ℓ. We have ξ ^{0}_ℓ = z^ℓε^−ℓ ∼ ε^−ℓj_ℓ. A term that matches with n_ℓ first appears in ξ ^{1}_ℓ . This can be seen from the horizon solution valid to $\mathcal{O}(\epsilon)$, Equation (79). The second term in the square brackets of it produces a term ε(z/ε)^−ℓ−1 = ε^ℓ+2z^−ℓ−1 ∼ ε^ℓ+2n_ℓ. This term therefore becomes $\mathcal{O}(\epsilon^{2\ell+2}/z^{2\ell+1})$ higher than the lowest order term ε^−ℓj_ℓ. Since ℓ ≥ 2, this effect first appears at $\mathcal{O}(\epsilon^6)$ in the post-Minkowski expansion, while it first appears at $\mathcal{O}(v^{13})$ in the post-Newtonian expansion if we note that $\epsilon=\mathcal{O}(v^3)$ and $z=\mathcal{O}(v)$. This implies, in particular, that if we are interested in the gravitational waves emitted to infinity, we may simply impose the regularity at z = 0 as the inner boundary condition of the outer solution for the calculation up to 6PN order beyond the quadrupole formula.

The above fact that a non-trivial boundary condition due to the presence of the black hole horizon appears at $\mathcal{O}(\epsilon^{2\ell+2})$ in the post-Minkowski expansion can be more easily seen as follows. Since $j_\ell=\mathcal{O}(z^\ell)$ as z → 0, we have ${X_\ell } \to O({\epsilon^{\ell + 1}}){e^{ - iz*}}$, or $A_\ell ^{\text{trans}} = \mathcal{O}\left( {{\epsilon^{\ell + 1}}} \right)$, where z* = z + ε ln(z - ε). On the other hand, from the asymptotic behavior of j_ℓ at z = ∞, the coefficients A ^inc_ℓ and A ^ref_ℓ must be of order unity. Then, using the Wronskian argument, we find

$$|A_\ell ^{inc}| - |A_\ell ^{ref}| = \frac{{|A_\ell ^{trans}{|^2}}}{{|A_\ell ^{inc}| + |A_\ell ^{ref}|}} = \mathcal{O}({ \in ^{2\ell + 2}}).$$

((93))

Thus, we immediately see that a non-trivial boundary condition appears at $\mathcal{O}\left( {{\epsilon^{2\ell + 2}}} \right)$.

It is also useful to keep in mind the above fact when we solve for ξ_ℓ under the post-Minkowski expansion. It implies that we may choose a phase such that A ^inc_ℓ and A ^ref_ℓ are complex conjugate to each other, to $\mathcal{O}\left( {{\epsilon^{2\ell + 1}}} \right)$. With this choice, the imaginary part of X_ℓ, which reflects the boundary condition at the horizon, does not appear until $\mathcal{O}(\epsilon^{2\ell+2})$ because the Regge-Wheeler equation is real. Then, recalling the relation of ξ_ℓ to X_ℓ, Equation (71), Im (ξ ⁽ⁿ⁾_ℓ ) for a given n ≤ 2ℓ+ 1 is completely determined in terms of Re (ξ ^(r)_ℓ ) for r ≤ n - 1. That is, we may focus on solving only the real part of Equation (84).

3.5 Structure of the ingoing wave function to $\mathcal{O}(\epsilon^2)$

With the boundary condition discussed in Section 2, we can integrate the ingoing wave Regge-Wheeler function iteratively to higher orders of ε in the post-Minkowskian expansion, ε ≪ 1. This was carried out in [91] to $\mathcal{O}(\epsilon^2)$ and in [105] to $\mathcal{O}(\epsilon^3)$ (See [71] for details). Here, we do not recapitulate the details of the calculation since it is already quite involved at $\mathcal{O}(\epsilon^2)$, with much less space for physical intuition. Instead, we describe the general properties of the ingoing wave function to $\mathcal{O}(\epsilon^2)$.

As discussed in Section 2, the ingoing wave Regge-Wheeler function X_ℓ can be made real up to $\mathcal{O}(\epsilon^{2\ell+1})$, or $\mathcal{O}(\epsilon^5)$ of the post-Minkowski expansion, if we recall ℓ ≥ 2. Choosing the phase of X_ℓ in this way, let us explicitly write down the expressions of Im(ξ ⁽ⁿ⁾_ℓ ) (n = 1, 2) in terms of Re(ξ ^(m)_ℓ ) (m ≤ n - 1). We decompose the real and imaginary parts of ξ ⁽ⁿ⁾_ℓ as

$$\xi _\ell ^{(n)} = f_\ell ^{(n)} + ig_\ell ^{(n)}.$$

((94))

Inserting this expression into Equation (71), and expanding the result with respect to ε (and noting f ⁽⁰⁾_ℓ = j_ℓ and g ⁽⁰⁾_ℓ = 0), we find

$$\begin{gathered} {X_\ell } = {e^{ - i \in \ln (z - \in )}}z\left[ {{j_\ell } + \in \left( {f_\ell ^{(1)} + ig_\ell ^{(1)}} \right) + { \in ^2}\left( {f_\ell ^{(2)} + ig_\ell ^{(2)}} \right) + ...} \right] \hfill \\ \;\;\;\;\; = z\left[ {{j_\ell } + \in f_\ell ^{(1)} + { \in ^2}\left( {f_\ell ^{(2)} + g_\ell ^{(1)}\ln z - \frac{1}{2}{j_\ell }(\ln z) + ...} \right)} \right] \hfill \\ \;\;\;\;\; + iz\left[ { \in (g_\ell ^{(1)} - {j_\ell }\ln z) + { \in ^2}\left( {g_\ell ^{(2)} + \frac{1}{z}{j_\ell } - f_\ell ^{(1)}\ln z} \right) + ...} \right]. \hfill \\ \end{gathered} $$

((95))

Hence, we have

$$g_\ell ^{(1)} = {j_\ell }\ln z,\;\;\;g_\ell ^{(2)} = - \frac{1}{z}{j_\ell } + f_\ell ^{(1)}\ln z,\;\;\;...$$

((96))

We thus have the post-Minkowski expansion of X_ℓ as

$${X_\ell } = \sum\limits_{n = 0}^\infty {{ \in ^n}X_\ell ^{(n)}} ,\;\;\;with\;X_\ell ^{(0)} = z{j_\ell },\;\;\;X_\ell ^{(1)} = zf_\ell ^{(1)},\;\;\;X_\ell ^{(2)} = z\left( {f_\ell ^{(2)} + \frac{1}{2}{j_\ell }{{(\ln z)}^2}} \right),\;\;...$$

((97))

Now, let us consider the asymptotic behavior of X_ℓ at z ≪ 1. As we know that ξ ⁽¹⁾_ℓ and ξ ⁽²⁾_ℓ are regular at z = 0, it is readily obtained by simply assuming Taylor expansion forms for them (including possible In z terms), inserting them into Equation (84), and comparing the terms of the same order on both sides of the equation. We denote the right-hand side of Equation (84) by S ⁽ⁿ⁾_ℓ .

For n = 1, we have

$$\begin{gathered} \operatorname{Re} (S_\ell ^{(1)}) = \frac{1}{z}\left( {{{j''}_\ell } + \frac{1}{2}{{j'}_\ell } - \frac{{4 + {z^2}}}{{{z^2}}}{j_\ell }} \right) \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; = \left\{ \begin{gathered} \mathcal{O}(z)\;\;\;\;\;\;\;for\;\ell = 2, \hfill \\ \mathcal{O}({z^{\ell - 3}})\;\;\;for\;\ell \geqslant 3. \hfill \\ \end{gathered} \right. \hfill \\ \end{gathered} $$

((98))

Inserting this into Equation (84) with n = 1, we find

$$\operatorname{Re} (\xi _\ell ^{(1)}) = f_\ell ^{(1)} = \left\{ \begin{gathered} \mathcal{O}({z^3})\;\;\;\;\;\;for\;\ell = 2, \hfill \\ \mathcal{O}({z^{\ell - 1}})\;\;\;for\;\ell \geqslant 3. \hfill \\ \end{gathered} \right.$$

((99))

Of course, this behavior is consistent with the full post-Minkowski solution given in Equation (87).

For n = 2, we then have

$$\begin{gathered} \operatorname{Re} (S_\ell ^{(2)}) = \frac{1}{z}\left( {f{{_\ell ^{(1)}}^{\prime \prime }} = \frac{1}{z}f{{_\ell ^{(1)}}^\prime } - \frac{{4 + {z^2}}}{{{z^2}}}f_\ell ^{(1)}} \right) - \frac{1}{z}\left( {2g{{_\ell ^{(1)}}^\prime } + \frac{1}{z}g_\ell ^{(1)}} \right) \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\; = - \frac{1}{z}({j_\ell }\ln z)' - \frac{1}{{{z^2}}}{j_\ell }\ln z + \left\{ \begin{gathered} \mathcal{O}({z^{\ell - 2}}\;\;\;for\;\ell = 2,3, \hfill \\ \mathcal{O}({z^{\ell - 4}}\;\;\;for\;\ell \geqslant 4. \hfill \\ \end{gathered} \right. \hfill \\ \end{gathered} $$

((100))

This gives

$$\operatorname{Re} (\xi _\ell ^{(2)}) = f_\ell ^{(2)} = \left\{ \begin{gathered} \mathcal{O}({z^\ell }) + \mathcal{O}({z^\ell })\ln z - \tfrac{1}{2}{j_\ell }{(\ln z)^2}\;\;\;\;\;\;for\;\ell = 2,3, \hfill \\ \mathcal{O}({z^{\ell - 2}}) + \mathcal{O}({z^\ell })\ln z - \tfrac{1}{2}{j_\ell }{(\ln z)^2}\;\;\;for\;\ell \geqslant 4. \hfill \\ \end{gathered} \right.$$

((101))

Note that the In z terms in Equation (100) arising from g ⁽¹⁾_ℓ give the (In z)² term in f ⁽²⁾_ℓ that just cancels the j_ℓ(ln z)²/2 term of X ⁽²⁾_ℓ in Equation (97).

Inserting Equations (99) and (101) into the relevant expressions in Equation (97), we find

$$\begin{gathered} {X_2} = {z^3}\{ \mathcal{O}(1) + \in \mathcal{O}(z){ \in ^2}[\mathcal{O}(1) + \mathcal{O}(1)\ln z] + ...\} , \hfill \\ {X_3} = {z^3}\{ \mathcal{O}(z) + \in \mathcal{O}(1){ \in ^2}[\mathcal{O}(z) + \mathcal{O}(z)\ln z] + ...\} , \hfill \\ {X_\ell } = {z^3}\{ \mathcal{O}({z^{\ell - 2}}) + \in \mathcal{O}({z^{\ell - 3}}){ \in ^2}[\mathcal{O}({z^{\ell - 4}}) + \mathcal{O}({z^{\ell - 2}})\ln z] + ...\} \;\;\;\;\;for\ell \geqslant 4. \hfill \\ \end{gathered} $$

((102))

Note that, for ℓ = 2 and 3, the leading behavior of X ⁽ⁿ⁾_ℓ at n = ℓ - 1 is more regular than the naively expected behavior, ∼ z^ℓ+1−n, which propagates to the consecutive higher order terms in ε. This behavior seems to hold for general ℓ, but we do not know a physical explanation for it.

Given a post-Newtonian order to which we want to calculate, by setting $z=\mathcal{O}(v)$ and $\epsilon=\mathcal{O}(v^3)$, the above asymptotic behaviors tell us the highest order of X ⁽ⁿ⁾_ℓ we need. We also see the presence of In z terms in X ⁽²⁾_ℓ . The logarithmic terms appear as a consequence of the mathematical structure of the Regge-Wheeler equation at z ≪ 1. The simple power series expansion of X ⁽ⁿ⁾_ℓ in terms of z breaks down at $\mathcal{O}(\epsilon^2)$, and we have to add logarithmic terms to obtain the solution. These logarithmic terms will give rise to In v terms in the wave-form and luminosity formulae at infinity, beginning at $\mathcal{O}(v^6)$ [99, 100]. It is not easy to explain physically how these In v terms appear. But the above analysis suggests that the In v terms in the luminosity originate from some spatially local curvature effects in the near-zone.

Now we turn to the asymptotic behavior at z = ∞. For this purpose, let the asymptotic form of f ⁽ⁿ⁾_ℓ be

$$f_\ell ^{(n)} \to P_\ell ^{(n)}{j_\ell } + Q_\ell ^{(n)}{n_\ell }\;\;\;\;\;\;\;as\;z \to \infty .$$

((103))

Noting Equation (97) and the equality ε^{−iε ln(z-ε)} = e^−iz* e^iz, the asymptotic form of X_ℓ is expressed as

$${X_\ell } \to A_\ell ^{inc}e{}^{ - i(z* - \in \ln \in )} + A_\ell ^{ref}{e^{i(z* - \in \ln \in )}},$$

((104))

$$\begin{gathered} A_\ell ^{inc} = \frac{1}{2}{i^{\ell + 1}}{e^{ - i \in \ln \in }}\left\{ {1 + \in \left[ {P_\ell ^{(1)} + i\left( {Q_\ell ^{(1)} + \ln z} \right)} \right]} \right. \hfill \\ \;\;\;\;\;\;\left. { + { \in ^2}\left[ {\left( {P_\ell ^{(2)} - Q_\ell ^{(1)}\ln z\} + i\left( {Q_\ell ^{(2)} + P_\ell ^{(1)}\ln z} \right)} \right.} \right] + ...} \right\}. \hfill \\ \end{gathered} $$

((105))

Note that

$$\omega r* = \omega \left( {r + 2M\ln \frac{{r - 2M}}{{2M}}} \right) = z* - \in \ln \in ,$$

((106))

because of our definition of z*, z* = z + ε + ln(z - ε). The phase factor ε^{−iε In ε} of A ^inc_ℓ originates from this definition, but it represents a physical phase shift due to wave propagation on the curved background.

As one may immediately notice, the above expression for A ^inc_ℓ contains In z-dependent terms. Since A ^inc_ℓ should be constant, P ⁽ⁿ⁾_ℓ and Q ⁽ⁿ⁾_ℓ should contain appropriate ln z-dependent terms which exactly cancel the ln z-dependent terms in Equation (105). To be explicit, we must have

$$\begin{gathered} P_\ell ^{(1)} = p_\ell ^{(1)}, \hfill \\ Q_\ell ^{(1)} = q_\ell ^{(1)} - \ln z, \hfill \\ P_\ell ^{(2)} = p_\ell ^{(2)} + q_\ell ^{(1)}\ln z - {(\ln z)^2}, \hfill \\ Q_\ell ^{(2)} = q_\ell ^{(2)} - p_\ell ^{(1)}\ln z, \hfill \\ \end{gathered} $$

((107))

where p ⁽ⁿ⁾_ℓ and q ⁽ⁿ⁾_ℓ are constants. These relations can be used to check the consistency of the solution f⁽ⁿ⁾ obtained by integration. In terms of p ⁽ⁿ⁾_ℓ and q ⁽ⁿ⁾_ℓ , A ^inc_ℓ is expressed as

$$A_\ell ^{inc} = \frac{1}{2}{i^{\ell + 1}}{e^{ - i \in \ln \in }}\left[ {1 + \in \left( {p_\ell ^{(1)} + iq_\ell ^{(1)}} \right) + { \in ^2}\left( {p_\ell ^{(2)} + iq_\ell ^{(2)}} \right) + ...} \right].$$

((108))

Note that the above form of A ^inc_ℓ implies that the so-called tail of radiation, which is due to the curvature scattering of waves, will contain In v terms as phase shifts in the waveform, but will not give rise to such terms in the luminosity formula. This supports our previous argument on the origin of the ln v terms in the luminosity. That is, it is not due to the wave propagation effect but due to some near-zone curvature effect.

4 Analytic Solutions of the Homogeneous Teukolsky Equation by Means of the Series Expansion of Special Functions

In this section, we review a method developed by Mano, Suzuki, and Takasugi [68], who found analytic expressions of the solutions of the homogeneous Teukolsky equation. In this method, the exact solutions of the radial Teukolsky equation (14) are expressed in two kinds of series expansions. One is given by a series of hypergeometric functions and the other by a series of the Coulomb wave functions. The former is convergent at horizon and the latter at infinity. The matching of these two solutions is done exactly in the overlapping region of convergence. They also found that the series expansions are naturally related to the low frequency expansion. Properties of the analytic solutions were studied in detail in [69]. Thus, the formalism is quite powerful when dealing with the post-Newtonian expansion, especially at higher orders.

In many cases, when we study the perturbation of a Kerr black hole, it is more convenient to use the Sasaki-Nakamura equation, since it has the form of a standard wave equation, similar to the Regge-Wheeler equation. However, it is not quite suited for investigating analytic properties of the solution near the horizon. In contrast, the Mano-Suzuki-Takasugi (MST) formalism allows us to investigate analytic properties of the solution near the horizon systematically. Hence, it can be used to compute the higher order post-Newtonian terms of the gravitational waves absorbed into a rotating black hole.

We also note that this method is the only existing method that can be used to calculate the gravitational waves emitted to infinity to an arbitrarily high post-Newtonian order in principle.

4.1 Angular eigenvalue

The solutions of the angular equation (15) that reduce to the spin-weighted spherical harmonics in the limit aω → 0 are called the spin-weighted spheroidal harmonics. They are the eigenfunctions of Equation (15), with λ being the eigenvalues. The eigenvalues λ are necessary for discussions of the radial Teukolsky equation. For general spin weight s, the spin weighted spheroidal harmonics obey

$$\left\{ {\frac{1}{{\sin \theta }}\frac{d}{{d\theta }}\left[ {\sin \theta \frac{d}{{d\theta }}} \right] - {a^2}{\omega ^2}{{\sin }^2}\theta - \frac{{{{(m + s\cos \theta )}^2}}}{{{{\sin }^2}\theta }} - 2a\omega s\cos \theta + s + 2ma\omega + \lambda } \right\}{}_s{S_{\ell m}} = 0.$$

((109))

In the post-Newtonian expansion, the parameter aω is assumed to be small. Then, it is straightforward to obtain a spheroidal harmonic _sS_ℓm of spin-weight s and its eigenvalue λ perturbatively by the standard method [86, 101, 94].

It is also possible to obtain the spheroidal harmonics by expansion in terms of the Jacobi functions [35]. In this method, if we calculate numerically, we can obtain them and their eigenvalues for an arbitrary value of aω.

Here we only show an analytic formula for the eigenvalue λ accurate to $\mathcal{O}({(a\omega )^2})$, which is needed for the calculation of the radial functions. It is given by

$$\lambda = {\lambda _0} + a\omega {\lambda _1} + a{}^2{\omega ^2}{\lambda _2} + \mathcal{O}({(a\omega )^3}),$$

((110))

where

$$\begin{gathered} {\lambda _0} = \ell (\ell + 1) - s(s + 1), \hfill \\ {\lambda _1} = - 2m\left( {1\tfrac{{{s^2}}}{{\ell (\ell + 1)}}} \right), \hfill \\ {\lambda _2} = - H(\ell + 1) - H(\ell ), \hfill \\ \end{gathered} $$

((111))

$$H(\ell ) = \frac{{2({\ell ^2} - {m^2}){{({\ell ^2} - {s^2})}^2}}}{{(2\ell - 1){\ell ^3}(2\ell + 1)}}.$$

((112))

4.2 Horizon solution in series of hypergeometric functions

As in Section 3, we focus on the ingoing wave function of the radial Teukolsky equation (14). Since the analysis below is applicable to any spin, |s| = 0, 1/2, 1, 3/2, and 2, we do not specify it except when it is needed. Also, the analysis is not restricted to the case aω ≪ 1 unless so stated explicitly. For general spin weight s, the homogeneous Teukolsky equation is given by

$${\Delta ^{ - s}} = \frac{d}{{dr}}\left( {{\Delta ^{s + 1}}\frac{{d{R_{\ell m\omega }}}}{{dr}}} \right) + \left( {\frac{{{K^2} - 2is(r - M)K}}{\Delta } + 4is\omega r - \lambda } \right){R_{\ell m\omega }} = 0.$$

((113))

As before, taking account of the symmetry ${\bar R_{\ell m\omega }} = {R_{\ell - m - \omega }}$, we may assume ε = 2Mω > 0 if necessary.

The Teukolsky equation has two regular singularities at r = r_±, and one irregular singularity at r = ∞. This implies that it cannot be represented in the form of a single hypergeometric equation. However, if we focus on the solution near the horizon, it may be approximated by a hypergeometric equation. This motivates us to consider the solution expressed in terms of a series of hypergeometric functions.

We define the independent variable x in place of z (= ωr) as

$$x = \frac{{{z_ + } - z}}{{ \in \kappa }},$$

((114))

where

$${z_ \pm } = \omega r{z_ \pm },\;\;\;\;\;\;\kappa = \sqrt {1 - {q^2}} ,\;\;\;\;\;\;q = \frac{a}{M}.$$

((115))

For later convenience, we also introduce gt = (ε - mq)/κ and ε_± = (ε ± τ)/2. Taking into account the structure of the singularities at r = r_±, we put the ingoing wave Teukolsky function R ⁱⁿ_ℓmω as

$$R_{\ell m\omega }^{in} = {e^{i \in \kappa x}} = {( - x)^{ - s - i( \in + \tau )/2}}{(1 - x)^{i( \in - \tau )/2}}{p_{in}}(x).$$

((116))

Then the radial Teukolsky equation becomes

$$\begin{gathered} x(1 - x){p_{in}}^{\prime \prime } + [1 - s - i \in - i\tau - (2 - 2i\tau )x]{p_{in}}^\prime + [i\tau (1 - i\tau ) + \lambda + s(s + 1)]{p_{in}} = \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;2i \in \kappa [ - x(1 - x){p_{in}}^\prime + (1 - s + i \in - i\tau )x{p_{in}}] + [{ \in ^2} - i \in \kappa (1 - 2s)]p{}_{in}, \hfill \\ \end{gathered} $$

((117))

where a prime denotes d/dx. The left-hand side of Equation (117) is in the form of a hypergeometric equation. In the limit ε → 0, noting Equation (110), we find that a solution that is finite at x = 0 is given by

$${p_{in}}( \in \to 0) = F( - \ell - i\tau ,\ell + 1 - i\tau ,1 - s - i\tau ,x).$$

((118))

For a general value of ε, Equation (117) suggests that a solution may be expanded in a series of hypergeometric functions with ε being a kind of expansion parameter. This idea was extensively developed by Leaver [64]. Leaver obtained solutions of the Teukolsky equation expressed in a series of the Coulomb wave functions. The MST formalism is an elegant reformulation of the one by Leaver [64].

The essential point is to introduce the so-called renormalized angular momentum ν, which is a generalization of ℓ, to a non-integer value such that the Teukolsky equation admits a solution in a convergent series of hypergeometric functions. Namely, we add the term [ν(ν+1) - λ - s(s+1)] p_in to both sides of Equation (117) to rewrite it as

$$\begin{gathered} x(1 - x){p_{in}}^{\prime \prime } + [1 - s - i \in - i\tau - (2 - 2i\tau )x]{p_{in}}^\prime + [i\tau (1 - i\tau ) + v(v + 1)]{p_{in}} = \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;2i \in \kappa [ - x(1 - x){p_{in}}^\prime + (1 - s + i \in - i\tau )x{p_{in}}] \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + [v(v + 1) - \lambda - s(s + 1) + \in {}^2 - i \in \kappa (1 - 2s){p_{in}}. \hfill \\ \end{gathered} $$

((119))

Of course, no modification is done to the original equation, and ν is just an irrelevant parameter at this stage. A trick is to consider the right-hand side of the above equation as a perturbation, and look for a formal solution specified by the index ν in a series expansion form. Then, only after we obtain the formal solution, we require that the series should converge, and this requirement determines the value of ν. Note that, if we take the limit ε → 0, we must have ν → ℓ (or ν → - ℓ - 1) to assure [ν(ν + 1) - λ - s(s + 1)] → 0 and to recover the solution (118).

Let us denote the formal solution specified by a value of ν by p ^ν_in . We express it in the series form,

$$\begin{gathered} \;\;\;\;\;\;\;\;\;p_{in}^v = \sum\limits_{n = - \infty }^\infty {{a_n}{p_{n + v}}(x),} \hfill \\ {p_{n + v}}(x) = F(n + v + 1 - i\tau , - n - v - i\tau ;1 - s - i \in - i\tau ;x). \hfill \\ \end{gathered} $$

((120))

Here, the hypergeometric functions p_n+ν (x) satisfy the recurrence relations [68],

$$\begin{gathered} x{p_{n + v}} = - \frac{{(n + v + 1 - s - i \in )(n + v + 1 - i\tau )}}{{2(n + v + 1)(2n + 2v + 1)}}{p_{n + v + 1}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\; + \frac{1}{2}\left[ {1 + \frac{{i\tau (s + i \in )}}{{(n + v)(n + v + 1)}}} \right]{p_{n + v}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\; - \frac{{(n + v + s + i \in )(n + v + i\tau )}}{{2(n + v)(2n + 2v + 1)}}{p_{n + v - 1}}, \hfill \\ \end{gathered} $$

((121))

$$\begin{gathered} x(1 - x){{p'}_{n + v}} = \frac{{(n + v + i\tau )(n + v + 1 - i\tau )(n + v + 1 - s - i \in )}}{{2(n + v + 1)(2n + 2v + 1)}}{p_{n + v + 1}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\; + \frac{1}{2}(s + i \in )\left[ {1 + \frac{{i\tau (1 - i\tau )}}{{(n + v)(n + v + 1)}}} \right]{p_{n + v}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\; - \frac{{(n + v + 1 - i\tau )(n + v + i\tau )(n + v + s + i \in )}}{{2(n + v)(2n + 2v + 1)}}{p_{n + v - 1}}, \hfill \\ \end{gathered} $$

((122))

Inserting the series (120) into Equation (119) and using the above recurrence relations, we obtain a three-term recurrence relation among the expansion coefficients a_n. It is given by

$$\alpha _n^v{a_{n + 1}} + \beta _n^v{a_n} + \gamma _n^v{a_{n - 1}} = 0,$$

((123))

where

$$\begin{gathered} \alpha _n^v = \frac{{i \in \kappa (n + v + 1 + s + i \in )(n + v + 1 + s - i \in )(n + v + 1 + i\tau )}}{{(n + v + 1)(2n + 2v + 3)}}, \hfill \\ \beta _n^v = - \lambda - s(s + 1) + (n + v)(n + v + 1) + { \in ^2} + \in ( \in - mq) + \frac{{ \in ( \in - mq)({s^2} + { \in ^2})}}{{(n + v)(n + v + 1)}}, \hfill \\ \gamma _n^v = \frac{{i \in \kappa (n + v - s + i \in )(n + v - s - i \in )(n + v - i\tau )}}{{(n + v)(2n - 2v - 1)}}. \hfill \\ \end{gathered} $$

((124))

The convergence of the series (120) is determined by the asymptotic behaviors of the coefficients a ^ν_n at n → ±t8. We thus discuss properties of the three-term recurrence relation (123) and the role of the parameter ν in detail.

The general solution of the recurrence relation (123) is expressed in terms of two linearly independent solutions {f ⁽¹⁾_n } and {f ⁽²⁾_n } (n - ±1, ±2.... ). According to the theory of three-term recurrence relations (see [49], Page 31) when there exists a pair of solutions that satisfy

$$\mathop {\lim }\limits_{n \to \infty } \frac{{f_n^{(1)}}}{{f_n^{(2)}}} = 0\;\;\;\;\left( {\mathop {\lim }\limits_{n \to - \infty } \frac{{f_n^{(1)}}}{{f_n^{(2)}}} = 0} \right),$$

((125))

then the solution {f ⁽¹⁾_n } is called minimal as n → ∞ (n → - ∞). Any non-minimal solution is called dominant. The minimal solution (either as n → ∞ or as n → - ∞) is determined uniquely up to an overall normalization factor.

The three-term recurrence relation is closely related to continued fractions. We introduce

$${R_n} \equiv \frac{{{a_n}}}{{{a_{n - 1}}}},\;\;\;{L_n} \equiv \frac{{{a_n}}}{{{a_{n + 1}}}}.$$

((126))

We can express R_n and L_n in terms of continued fractions as

$${R_n} = \frac{{\gamma _n^v}}{{\beta _n^v + \alpha _n^v{R_{n + 1}}}} = \frac{{\gamma _n^v}}{{\beta _n^v - }}.\frac{{{\alpha _n}\gamma _{n + 1}^v}}{{\beta _{n + 1}^v - }}.\frac{{{\alpha _{n + 1}}\gamma _{n + 2}^v}}{{\beta _{n + 2}^v - }}....,$$

((127))

$${L_n} = \frac{{\alpha _n^v}}{{\beta _n^v + \gamma _n^v{L_{n - 1}}}} = \frac{{\alpha _n^v}}{{\beta _n^v - }}.\frac{{{\alpha _{n - 1}}\gamma _n^v}}{{\beta _{n - 1}^v - }}.\frac{{{\alpha _{n - 2}}\gamma _{n - 1}^v}}{{\beta _{n - 2}^v - }}.....$$

((128))

These expressions for R_n and L_n are valid if the respective continued fractions converge. It is proved (see [49], Page 31) that the continued fraction (127) converges if and only if the recurrence relation (123) possesses a minimal solution as n → ∞, and the same for the continued fraction (128) as n → - ∞.

Analysis of the asymptotic behavior of (123) shows that, as long as ν is finite, there exists a set of two independent solutions that behave as (see, e.g., [49], Page 35)

$$\mathop {\lim }\limits_{n \to \infty } \;n\frac{{a_n^{(1)}}}{{a_{n - 1}^{(1)}}} = \frac{{i \in \kappa }}{2},\;\;\;\;\;\mathop {\lim }\limits_{n \to \infty } \frac{{a_n^{(2)}}}{{na_{n - 1}^{(2)}}} = \frac{{2i}}{{ \in \kappa }},$$

((129))

and another set of two independent solutions that behave as

$$\mathop {\lim }\limits_{n \to - \infty } \;n\frac{{b_n^{(1)}}}{{b_{n + 1}^{(1)}}} = \frac{{i \in \kappa }}{2},\;\;\;\;\;\mathop {\lim }\limits_{n \to - \infty } \frac{{b_n^{(2)}}}{{nb_{n + 1}^{(2)}}} = \frac{{2i}}{{ \in \kappa }},$$

((130))

Thus, {a ⁽¹⁾_n } is minimal as n → ∞ and {b ⁽¹⁾_n } is minimal as n → - ∞.

Since the recurrence relation (123) possesses minimal solutions as n → ±∞, the continued fractions on the right-hand sides of Equations (127) and (128) converge for a_n = a ⁽¹⁾_n and a_n = b ⁽¹⁾_n . In general, however, a ⁽¹⁾_n and b ⁽¹⁾_n do not coincide. Here, we use the freedom of ν to obtain a consistent solution. Let {f ^ν_n } be a sequence that is minimal for both n → ±∞. We then have expressions for f ^ν_n /f ^ν_n-1 and f ^ν_n /f ^ν_n+1 in terms of continued fractions as

$${\tilde R_n} \equiv \frac{{{f_n}}}{{{f_{n - 1}}}} = - \frac{{\gamma _n^v}}{{\beta _n^v - }}.\frac{{{\alpha _n}\gamma _{n + 1}^v}}{{\beta _{n + 1}^v - }}.\frac{{{\alpha _{n + 1}}\gamma _{n + 2}^v}}{{\beta _{n + 2}^v - }}....,$$

((131))

$${\tilde L_n} \equiv \frac{{{f_n}}}{{{f_{n + 1}}}} = - \frac{{\alpha _n^v}}{{\beta _n^v - }}.\frac{{{\alpha _{n - 1}}\gamma _n^v}}{{\beta _{n - 1}^v - }}.\frac{{{\alpha _{n - 2}}\gamma _{n - 1}^v}}{{\beta _{n - 2}^v - }}.....$$

((132))

This implies

$${\tilde R_n}{\tilde L_{n - 1}} = 1.$$

((133))

Thus, if we choose ν such that it satisfies the implicit equation for ν, Equation (133), for a certain n, we obtain a unique minimal solution {f_n^ν} that is valid over the entire range of n, - ∞ < n < ∞, that is

$$\mathop {\lim }\limits_{n \to \infty } n\frac{{f_n^v}}{{f_{n - 1}^v}} = \frac{{i \in \kappa }}{2},\;\;\;\;\;\;\;\;\mathop {\lim }\limits_{n \to - \infty } n\frac{{f_n^v}}{{f_{n + 1}^v}} = \frac{{i \in \kappa }}{2}.$$

((134))

Note that if Equation (133) for a certain value of n is satisfied, it is automatically satisfied for any other value of n.

The minimal solution is also important for the convergence of the series (120). For the minimal solution {f ^ν_n }, together with the properties of the hypergeometric functions p_n+ν for large |n|, we find

$$\mathop {\lim }\limits_{n \to \infty } n\frac{{f_{n + 1}^v{p_{n + v + 1(x)}}}}{{f_n^v{p_{n + v(x)}}}} = - \mathop {\lim }\limits_{n \to \infty } n\frac{{f_{n - 1}^v{p_{n + v - 1(x)}}}}{{f_n^v{p_{n + v(x)}}}} = \frac{{i \in \kappa }}{2}[1 - 2x + {({(1 - 2x)^2} - 1)^{1/2}}].$$

((135))

Thus, the series of hypergeometric functions (120) converges for all x in the range 0 ≥ x > - ∞ (in fact, for all complex values of x except at |x| = ∞), provided that the coefficients are given by the minimal solution.

Instead of Equation (133), we may consider an equivalent but practically more convenient form of an equation that determines the value of ν. Dividing Equation (123) by a_n, we find

$$\beta _n^v + \alpha _n^v{R_{n + 1}} + \gamma _n^v{L_{n - 1}} = 0,$$

((136))

where R_n+1 and L_n+1 are those given by the continued fractions (131) and (132), respectively. Although the value of n in this equation is arbitrary, it is convenient to set n = 0 to solve for ν.

For later use, we need a series expression for Rⁱⁿ with better convergence properties at large |x|. Using analytic properties of hypergeometric functions, we have

$${R^{in}} = R_0^v + R_0^{ - v - 1},$$

((137))

where

$$\begin{gathered} R_0^v = {e^{i \in \kappa x}}{( - x)^{ - s - (i/2)( \in + \tau )}}{(1 - x)^{(i/2)( \in + \tau ) + v}} \hfill \\ \;\;\;\;\;\;\;\; \times \sum\limits_{n = - \infty }^\infty {f_n^v\frac{{\Gamma (1 - s - i \in - i\tau )\Gamma (2n + 2v + 1)}}{{\Gamma (n + v + 1 - i\tau )\Gamma (n + v + 1 - s - i \in )}}} \hfill \\ \;\;\;\;\;\;\;\; \times {(1 - x)^n}F( - n - v - i\tau , - n - v - s - i \in ; - 2n - 2v;\frac{1}{{1 - x}}). \hfill \\ \end{gathered} $$

((138))

This expression explicitly exhibits the symmetry of R_in under the interchange of ν and -ν - 1. This is a result of the fact that ν(ν+1) is invariant under the interchange ν ↔ - ν - 1. Accordingly, the recurrence relation (123) has the structure that f ^-ν-1_-n satisfies the same recurrence relation as f ^ν_n .

Finally, we note that if ν is a solution of Equation (133) or (136), ν + k with an arbitrary integer k is also a solution, since ν appears only in the combination of n+ν. Thus, Equation (133) or (136) contains an infinite number of roots. However, not all of these can be used to express a solution we want. As noted in the earlier part of this section, in order to reproduce the solution in the limit ε → 0, Equation (118), we must have ν → ℓ (or ν → -ℓ - 1 by symmetry). Thus, we impose a constraint on ν such that it must continuously approach ℓ as ε → 0.

4.3 Outer solution as a series of Coulomb wave functions

The solution as a series of hypergeometric functions discussed in Section 4.2 is convergent at any finite value of r. However, it does not converge at infinity, and hence the asymptotic amplitudes, B^inc and B^ref, cannot be determined from it. To determine the asymptotic amplitudes, it is necessary to construct a solution that is valid at infinity and to match the two solutions in a region where both solutions converge. The solution convergent at infinity was obtained by Leaver as a series of Coulomb wave functions [64]. In this section, we review Leaver’s solution based on [69].

In this section again, by noting the symmetry ${\bar R_{\ell m\omega }} = {R_{\ell - m - \omega }}$, we assume ω > 0 without loss of generality.

First, we define a variable $\hat z = \omega (r - r\_) = \epsilon\kappa (1 - x)$. Let us denote a Teukolsky function by R_C. We introduce a function $f(\hat z)$ by

$${R_C} = {\hat z^{ - 1 - s}}{\left( {1 - \frac{{ \in \kappa }}{{\hat z}}} \right)^{ - s - i( \in + \tau )/2}}f(\hat z).$$

((139))

Then the Teukolsky equation becomes

$$\begin{gathered} {{\hat z}^2}f''[{{\hat z}^2} + (2 \in + 2is)\hat z - \lambda - s(s + 1)]f = \in \kappa \hat z(f'' + f) + \in \kappa (s - 1 + 2i \in )f' \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \frac{ \in }{{\hat z}}[\kappa - i( \in - mq)](s - 1 + i \in )f \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + [ - 2{ \in ^2} + \in mq + \kappa ({ \in ^2} + i \in s)]f. \hfill \\ \end{gathered} $$

((140))

We see that the right-hand side is explicitly of $\mathcal{O}(\epsilon)$ and the left-hand side is in the form of the Coulomb wave equation. Therefore, in the limit ε → 0, we obtain a solution

$$f(\hat z) = {F_\ell }( - is - \in ,\hat z),$$

((141))

where ${F_L}(\eta ,\hat z)$ is a Coulomb wave function given by

$${F_L}(\eta ,\hat z) = {e^{ - i\hat z}}{2^L}{\hat z^{L + 1}}\frac{{\Gamma (L + 1 - i\eta )}}{{\Gamma (2L + 2)}}\Phi (L + 1 - i\eta ,2L + 2;2i\hat z),$$

((142))

and Φ is the regular confluent hypergeometric function (see [1], Section 13) which is regular at $\hat z = 0$.

In the same spirit as in Section 4.2, we introduce the renormalized angular momentum ν. That is, we add $\left[ {\lambda + s(s + 1) - v(v + 1)} \right]f(\hat z)$ to both sides of Equation (140) to rewrite it as

$$\begin{gathered} {{\hat z}^2}f'' + ({{\hat z}^2} + (2 \in + 2is)\hat z - v(v + 1)f = \in \kappa \hat z(f'' + f) + \in \kappa (s - 1 + 2i{ \in _ + })f' \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \frac{ \in }{{\hat z}}[\kappa - i( \in - mq)](s - 1 + i \in )f \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + [ - v(v + 1) + \lambda + s(s + 1) - 2{ \in ^2} + \in mq + \kappa ({ \in ^2} + i \in s)]f. \hfill \\ \end{gathered} $$

((143))

We denote the formal solution specified by the index ν by ${f_v}(\hat z)$, and expand it in terms of the Coulomb wave functions as

$${f_v} = \sum\limits_{n = - \infty }^\infty {{{( - i)}^n}} \frac{{{{(v + 1 + s - i \in )}_n}}}{{{{(v + 1 - s + i \in )}_n}}}{b_n}{F_{n + v}}( - is - \in ,\hat z),$$

((144))

where (a)_n = Γ(a + n)/Γ(a). Then, using the recurrence relations among F_n+ν,

$$\begin{gathered} \frac{1}{z}{F_{n + v}} = \frac{{(n + v + 1 + s - i \in )}}{{(n + v + 1)(2n + 2v + 1)}}{F_{n + v + 1}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + \frac{{is + \in }}{{(n + v)(n + v + 1)}}{F_{n + v}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + \frac{{(n + v - s + i \in )}}{{(n + v)(2n + 2v + 1)}}{F_{n + v - 1,}} \hfill \\ \end{gathered} $$

((145))

$$\begin{gathered} {{F'}_{n + v}} = \frac{{(n + v)(n + v + 1 + s - i \in )}}{{(n + v + 1)(2n + 2v + 1)}}{F_{n + v + 1}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + \frac{{is + \in }}{{(n + v)(n + v + 1)}}{F_{n + v}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + \frac{{(n + v + 1)(n + v - s + i \in )}}{{(n + v)(2n + 2v + 1)}}{F_{n + v - 1,}} \hfill \\ \end{gathered} $$

((146))

we can derive the recurrence relation among b_n. The result turns out to be identical to the one given by Equation (123) for a_n. We mention that the extra factor (ν+1+s-iε)_n/(ν+1-s+iε)_n in Equation (144) is introduced to make the recurrence relation exactly identical to Equation (123).

The fact that we have the same recurrence relation as Equation (123) implies that if we choose the parameter ν in Equation (144) to be the same as the one given by a solution of Equation (133) or (136), the sequence {f ^ν_n } is also the solution for {b_n}, which is minimal for both n → ± ∞. Let us set

$$g_n^v = {(i - )^n}\frac{{{{(v + 1 + s - i \in )}_n}}}{{{{(v + 1 - s + i \in )}_n}}}f_n^v.$$

((147))

By choosing ν as stated above, we have the asymptotic value for the ratio of two successive terms of g ^ν_n as

$$\mathop {\lim }\limits_{n \to \infty } \;n\frac{{g_n^v}}{{g_{n - 1}^v}} = \mathop {\lim }\limits_{n \to - \infty } \;n\frac{{g_n^v}}{{g_{n + 1}^v}} = \frac{{ \in \kappa }}{2}.$$

((148))

Using an asymptotic property of the Coulomb wave functions, we have

$$\mathop {\lim }\limits_{n \to \infty } \frac{{g_n^v{F_{n + v}}(z)}}{{g_{n - 1}^v{F_{n + v - 1}}(z)}} = \mathop {\lim }\limits_{n \to - \infty } \;n\frac{{g_n^v{F_{n + v}}(z)}}{{g_{n - 1}^v{F_{n + v + 1}}(z)}} = \frac{{ \in \kappa }}{z}.$$

((149))

We thus find that the series (144) converges at $\hat z > \epsilon\kappa$ or equivalently r > r₊.

The fact that we can use the same ν as in the case of hypergeometric functions to obtain the convergence of the series of the Coulomb wave functions is crucial to match the horizon and outer solutions.

Here, we note an analytic property of the confluent hypergeometric function (see [34], Page 259),

$$\Phi (a,c;x) = \frac{{\Gamma (c)}}{{\Gamma (c - a)}}{e^{ia\pi }}\psi (a,c;x) + \frac{{\Gamma (c)}}{{\Gamma (a)}}{e^{i\pi (a - c)}}{e^x}\psi (c - a,c; - x),$$

((150))

where Ψ is the irregular confluent hypergeometric function, and Im(x) > 0 is assumed. Using this with the identities

$$\begin{gathered} a = n + v + 1 - s + i \in , \hfill \\ c = 2(n + v + 1), \hfill \\ x = 2i\hat z, \hfill \\ \end{gathered} $$

((151))

we can rewrite R ^ν_C (for ω > 0) as

$$R_C^v = R_ + ^v + R_ - ^v,$$

((152))

where

$$\begin{array}{*{20}{c}} {R_{+}^{v}=2^{v}e^{-\pi\epsilon}e^{i\pi(v+1-s)}\frac{\Gamma(v+1-s+i\epsilon) }{\Gamma (v+1+s-i\epsilon)}e^{-i\hat{z}}\hat{z}^{v+i\epsilon+}(\hat{z}-\epsilon\kappa)^{-s-i\epsilon+}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \\ {\times \sum\limits_{n=-\infty}^{\infty}\;i^{n}f_{n}^{v}(2\hat{z})^{n}\Psi(n+v+1-s+i\epsilon,\;2n+2v+2;\;2i\hat{z})}, \\ {R_-^v=2^{v}e^{-\pi\epsilon}e^{-i\pi(v+1+s)}e^{i\hat{z}}\hat{z}^{v+i\epsilon+}(\hat{z}-\epsilon\kappa)^{-s-i\epsilon+}} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\\ {\times \sum\limits_{n=-\infty}^{\infty}\;i^{n}\frac{(v+1+s-i\epsilon)_{n}}{(v+1-s+i\epsilon)_{n}}f_{n}^{v}(2\hat{z})^{n}} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\\ {\times\Psi(n+v+1+s-i\epsilon,\;2n+2v+2;-i\hat{z})}.\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \end{array}$$

((153))

By noting an asymptotic behavior of Ψ (a, c; x) at large |x|,

$$\Psi (a,c;x) \to {x^{ - a}}\;\;\;\;\;\;\;\;as|x| \to \infty ,$$

((154))

we find

$$R_ + ^v = A_ + ^v{z^{ - 1}}{e^{ - i(z + \in \ln z)}},$$

((155))

$$R_ - ^v = A_ - ^v{z^{ - 1 - 2s}}{e^{i(z + \in \ln z)}},$$

((156))

where

$$A_ + ^v = {e^{ - (\pi /2) \in }}{e^{(\pi /2)i(v + 1 - s)}}{2^{ - 1 + s - i \in }}\frac{{\Gamma (v + 1 - s + i \in )}}{{\Gamma (v + 1 + s - i \in )}}\sum\limits_{n = - \infty }^{ + \infty } {f_n^v} ,$$

((157))

$$A_ - ^v = {2^{ - 1 + s + i \in }}{e^{ - (\pi /2)i(v + 1 + s)}}{e^{ - (\pi /2) \in }}\sum\limits_{n = - \infty }^{ + \infty } {{{( - 1)}^n}\frac{{{{(v + 1 + s - i \in )}_n}}}{{{{(v + 1 - s + i \in )}_n}}}f_n^v} ,$$

((158))

We can see that the functions R ^ν₊ and R ^ν_{_} are incoming-wave and outgoing wave solutions at infinity, respectively. In particular, we have the upgoing solution, defined for s = -2 by the asymptotic behavior (20), expressed in terms of a series of Coulomb wave functions as

$${R^{up}} = R_ - ^v.$$

((159))

4.4 Matching of horizon and outer solutions

Now, we match the two types of solutions R ^ν₊ and R ^ν_{_} . Note that both of them are convergent in a very large region of r, namely for $\epsilon\kappa < \hat z < \infty$. We see that both solutions behave as ${\hat z^v}$ multiplied by a single-valued function of $\hat z$ for large $\left| {\hat z} \right|$. Thus, the analytic properties of R ^ν₀ and R ^ν_C are the same, which implies that these two are identical up to a constant multiple. Therefore, we set

$$R_0^v = {K_v}R_C^v.$$

((160))

In the region $\epsilon\kappa < \hat z < \infty$, we may expand both solutions in powers of $\hat z$ except for analytically non-trivial factors. We have

$${R_{0}^{v}=e^{i\epsilon\kappa}e^{-i\hat{z}}(e\kappa)^{-v-i\epsilon+}\hat{z}^{v+i\epsilon+}\left(\frac{\hat{z}}{\epsilon\kappa}-1\right)^{-s-i\epsilon+}\sum\limits_{n=-\infty}^{\infty}\;\sum\limits_{j=0}^{\infty}C_{n,j}\hat{z}^{n-j}} \; \; \; \; \; \\ {=e^{i\epsilon\kappa}e^{-i\hat{z}}(\epsilon\kappa)^{-v-i\epsilon+}\hat{z}^{v+i\epsilon+}\left({\frac{\hat{z}}{\epsilon\kappa}}-1\right)^{-s-i\epsilon+}\sum\limits_{k=-\infty}^{\infty}\sum\limits_{n=k}^{\infty}C_{n,n-k}\hat{z}^{k}},$$

((161))

$${R_{C}^{v}=e^{-i\hat{z}}2^{v}(\epsilon\kappa)^{-s-i\epsilon+}\hat{z}^{v+i\epsilon+}\left(\frac{\hat{z}}{\epsilon\kappa}-1\right)^{-s-i\epsilon+}\sum\limits_{n=-\infty}^{\infty}\sum\limits_{j=0}^{\infty}D_{n,j}\hat{z}^{n+j}} \; \; \; \; \; \; \; \; \; \; \\ =e^{-i\hat{z}}2^{v}(\epsilon\kappa)^{-s-i\epsilon+}\hat{z}^{v+i\epsilon+}\left(\frac{\hat{z}}{\epsilon\kappa}-1\right)^{-s-i\epsilon+}\sum\limits_{k=-\infty}^{\infty}\sum\limits_{n=-\infty}^{k}D_{n,k-n}\hat{z}^{k},$$

((162))

where

$$\begin{gathered} {C_{n,j}} = \frac{{\Gamma (1 - s - 2i{ \in _ + })\Gamma (2n + 2v + 1)}}{{\Gamma (n + v + 1 - i\tau )\Gamma (n + v + 1 - s - i \in )}} \hfill \\ \;\;\;\;\;\;\;\;\; \times \frac{{{{( - n - v - i\tau )}_j}){{( - n - v - s - i \in )}_j}}}{{{{( - 2n - 2v)}_j}(j!)}}{( \in \kappa )^{ - n + j}}{f_n}, \hfill \\ \end{gathered} $$

((163))

$$\begin{gathered} {D_{n,j}} = {( - 1)^n}{(2i)^{n + j}}\frac{{\Gamma (n + v + 1 - s + i \in ){{(v + 1 + s - i \in )}_n}}}{{\Gamma (2n + 2v + 2)\;\;\;{{(v + 1 - s + i \in )}_n}}} \hfill \\ \;\;\;\;\;\;\;\;\; \times \frac{{{{(n + v + 1 - s + i \in )}_j}}}{{{{(2n + 2v + 2)}_j}(j!)}}{f_n}, \hfill \\ \end{gathered} $$

((164))

Then, by comparing each integer power of $\hat z$ in the summation, in the region $\epsilon\kappa \ll \hat z < \infty$, and using the formula Γ(z) Γ(1 - z) = π /sin πz, we find

$$\begin{gathered} {K_v} = {e^{i \in \kappa }}{( \in \kappa )^{s - v}}{2^{ - v}}{\left( {\sum\limits_{n = - \infty }^r {D{}_{n,r - n}} } \right)^{ - 1}}\left( {\sum\limits_{n = r}^\infty {C{}_{n,n - r}} } \right) \hfill \\ \;\;\;\;\; = \frac{{{e^{i \in \kappa }}{{(2 \in \kappa )}^{s - v - r}}{2^{ - s}}{i^r}\Gamma (1 - s - 2i{ \in _ + })\Gamma (r + 2v + 2)}}{{\Gamma (r + v + 1 - s + i \in )\Gamma (r + v + 1 + i\tau )\Gamma (r + v + 1 + s + i \in )}} \hfill \\ \;\;\;\;\;\;\;\; \times \left( {\sum\limits_{n = r}^\infty {{{( - 1)}^n}\frac{{\Gamma (n + r + 2v + 1)\Gamma (n + v + 1 + s + i \in )\Gamma (n + v + 1 + i\tau )}}{{(n - r)!\;\;\;\;\;\;\;\;\;\Gamma (n + v + 1 - s - i \in )\Gamma (n + v + 1 - i\tau )\;\;\;\;\;}}f_n^v} } \right) \hfill \\ \;\;\;\;\;\;\;\; \times {\left( {\sum\limits_{n = - \infty }^r {\frac{{{{( - 1)}^n}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{(v + 1 + s - i \in )}_n}}}{{(r - n)!\;{{(r + 2v + 2)}_n}{{(v + 1 - s + i \in )}_n}\;\;\;\;\;}}f_n^v} } \right)^{ - 1}}, \hfill \\ \end{gathered} $$

((165))

where r can be any integer, and the factor K_ν should be independent of the choice of r. Although this fact is not manifest from Equation (165), we can check it numerically, or analytically by expanding it in terms of ε.

We thus have two expressions for the ingoing wave function Rⁱⁿ. One is given by Equation (116), with p ^ν_in expressed in terms of a series of hypergeometric functions as given by Equation (120) (a series which converges everywhere except at r = ∞). The other is expressed in terms of a series of Coulomb wave functions given by

$${R^{in}} = {K_v}R_C^v + {K_{ - v - 1}}R_C^{ - v - 1},$$

((166))

which converges at r > r₊, including r = ∞. Combining these two, we have a complete analytic solution for the ingoing wave function.

Now we can obtain analytic expressions for the asymptotic amplitudes of Rⁱⁿ, B^trans, B^inc, and B^ref. By investigating the asymptotic behaviors of the solution at r → ∞ and r → r₊, they are found to be¹

$${B^{trans}} = {\left( {\frac{{ \in \kappa }}{\omega }} \right)^{2s}}{e^{i\kappa \in + (1 + \tfrac{{2\ln k}}{{1 + \kappa }})}}\sum\limits_{n = - \infty }^\infty {f_n^v,} $$

((167))

$${B^{inc}} = {\omega ^{ - 1}}\left( {{K_v} - i{e^{ - i\pi v}}\frac{{\sin \pi (v - s + i \in )}}{{\sin \pi (v + s - i \in )}}{K_{ - v - 1}}} \right)A_ + ^v{e^{ - i( \in \ln \in - \tfrac{{1 - \kappa }}{2} \in )}},$$

((168))

$${B^{ref}} = {\omega ^{ - 1 - 2s}}\left( {{K_v} + i{e^{i\pi v}}{K_{ - v - 1}}} \right)A_ - ^v{e^{i( \in \ln \in - \tfrac{{1 - \kappa }}{2} \in )}}.$$

((169))

Incidentally, since we have the upgoing solution in the outer region (159), it is straightforward to obtain the asymptotic outgoing amplitude at infinity C^trans from Equation (153). We find

$${C^{trans}} = {\omega ^{ - 1 - 2s}}A_ - ^v{e^{i( \in \ln \in - \tfrac{{1 - \kappa }}{2} \in )}}.$$

((170))

4.5 Low frequency expansion of the hypergeometric expansion

So far, we have considered exact solutions of the Teukolsky equation. Now, let us consider their low frequency approximations and determine the value of v. We solve Equation (136) with n = 0,

$$\beta _0^v + \alpha _0^v{R_1} + \gamma _0^v{L_{ - 1}} = 0,$$

((171))

with a requirement that ν → ℓ as ε → 0.

To solve Equation (174), we first note the following. Unless the value of ν is such that the denominator in the expression of α ^ν_n or γ ^ν_n happens to vanish, or β ^ν_n happens to vanish in the limit ε → 0, we have $\alpha _n^v = \mathcal{O}(\epsilon),\gamma _n^v = \mathcal{O}(\epsilon)$, and $\beta _n^v = \mathcal{O}(1)$. Also, from the asymptotic behavior of the minimal solution f ^ν_n as n → ±t8 given by Equation (134), we have ${R_n}(v) = \mathcal{O}(\epsilon)$ and ${L_{ - n}}(v) = \mathcal{O}(\epsilon)$ for sufficiently large n. Thus, except for exceptional cases mentioned above, the order of a ^ν_n in ε increases as |n| increases. That is, the series solution naturally gives the post-Minkowski expansion.

First, let us consider the case of R_n(ν) for n > 0. It is easily seen that $\alpha _n^v = \mathcal{O}(\epsilon),\gamma _n^v = \mathcal{O}(\epsilon)$, and $\beta _n^v = \mathcal{O}(1)$ for all n > 0. Therefore, we have ${R_n}(v) = \mathcal{O}(\epsilon)$ for all n > 0.

On the other hand, for n < 0, the order of L_-n(ν) behaves irregularly for certain values of n. For the moment, let us assume that ${L_{ - 1}}(v) = \mathcal{O}(\epsilon)$. We see from Equations (124) that $\alpha _0^v = \mathcal{O}(\epsilon)$, $\gamma _0^v = \mathcal{O}(\epsilon)$, since $v = \ell + \mathcal{O}(\epsilon)$. Then, Equation (174) implies $\beta _0^v = \mathcal{O}({\epsilon^2})$. Using the expansion of λ given by Equation (110), we then find $v = \ell + \mathcal{O}({\epsilon^2})$ (i.e., there is no term of $\mathcal{O}(\epsilon)$ in ν). With this estimate of ν, we see from Equation (128) that ${L_{ - 1}}(v) = \mathcal{O}(\epsilon)$ is justified if L_-2 (ν) is of order unity or smaller.

The general behavior of the order of L_-n(ν) in e for general values of s is rather complicated. However, if we assume s to be anon-integer and ℓ ≥ |s|, and $\tau = (\epsilon - mq)/\kappa = \mathcal{O}(1)$, it is relatively easily studied. With the assumption that $v = \ell + \mathcal{O}({\epsilon^2})$, we find there are three exceptional cases:

• For n = -2ℓ - 1, we have ${\alpha _n} = \mathcal{O}(\epsilon),\beta = \mathcal{O}({\epsilon^2})$, and ${\gamma _n} = \mathcal{O}(\epsilon)$.

• For n = -ℓ - 1, we have $\alpha _n^v = \mathcal{O}(1/\epsilon),\beta _n^v = \mathcal{O}(1/\epsilon)$, and ${\gamma _n} = \mathcal{O}(\epsilon)$.

• For n = -ℓ, we have $\mathcal{O}(\epsilon),{\beta _n} = \mathcal{O}(1/\epsilon)$, and ${\gamma _n} = \mathcal{O}(1/\epsilon)$.

These imply that ${L_{2\ell - 1}}(v) = \mathcal{O}(1/\epsilon),{L_{ - \ell - 1}}(v) = \mathcal{O}(1)$, and ${L_{ - \ell }}(v) = \mathcal{O}({\epsilon^2})$, respectively. To summarize, we have

$$\begin{gathered} \;\;\;\;\;\;\;\;{R_n}(v) = \frac{{f_n^v}}{{f_{n - 1}^v}} = \mathcal{O}( \in )\;\;\;\;\;\;\;\;\;\;\;\;for\;all\;n > 0, \hfill \\ \;\;\;\;\;{L_{ - \ell }}(v) = \frac{{f_{ - \ell }^v}}{{f_{ - \ell + 1}^v}} = \mathcal{O}({ \in ^2}),\;\; \hfill \\ \;\;{L_{ - \ell - 1}}(v) = \frac{{f_{ - \ell - 1}^v}}{{f_{ - \ell }^v}} = \mathcal{O}(1), \hfill \\ {L_{ - 2\ell - 1}}(v) = \frac{{f_{ - 2\ell - 1}^v}}{{f_{ - 2\ell }^v}} = \mathcal{O}(1/ \in ), \hfill \\ \;\;\;\;\;{L_n}(v) = \frac{{f_n^v}}{{f_{n + 1}^v}} = \mathcal{O}( \in )\;\;\;\;\;\;\;\;\;\;\;\;for\;all\;the\;other\;n < 0. \hfill \\ \end{gathered} $$

((172))

With these results, we can calculate the value of ν to $\mathcal{O}(\epsilon)$, which is given by

$$v = \ell + \frac{1}{{2\ell + 1}}\left( { - 2 - \frac{{{s^2}}}{{\ell (\ell + 1)}} + \frac{{{{[(\ell + 1){}^2 - {s^2})}^2}}}{{(2\ell + 1)(2\ell + 2)(2\ell + 3)}} - \frac{{{{({\ell ^2} - {s^2})}^2}}}{{(2\ell - 1)(2\ell + 1)}}} \right){ \in ^2} + \mathcal{O}({ \in ^3}).$$

((173))

Now one can take the limit of an integer value of s. In particular, the above holds also for s = 0. Interestingly, ν is found to be independent of the azimuthal eigenvalue m to $\mathcal{O}(\epsilon^2)$.

The post-Minkowski expansion of homogeneous Teukolsky functions can be obtained with arbitrary accuracy by solving Equation (123) to a desired order, and by summing up the terms to a sufficiently large |n|. The first few terms of the coefficients f ^ν_n are explicitly given in [68]. A calculation up to a much higher order in $\mathcal{O}(\epsilon)$ was performed in [98], in which the black hole absorption of gravitational waves was calculated to $\mathcal{O}({v^8})$ beyond the lowest order.

4.6 Property of ν

In this section, we discuss the property of the solution of Equation (136) which we recapitulate:

$${g_n}(v) \equiv \beta _n^v + \alpha _n^v{R_{n + 1}} + \gamma _n^v{L_{n - 1}} = 0.$$

((174))

The numerical evaluation of this equation was not done very much before. Leaver [64] briefly mentioned a numerical implementation of a code to obtain ν by solving Equation (174). In the Schwarzschild case, Tagoshi and Nakamura [99] solved Equation (174) numerically to obtain ν. They evaluated the homogeneous solution numerically based on Leaver’s method [64] by using the value of ν obtained numerically. Later, Fujita and Tagoshi [40] tried a numerical implementation of the MST method. They found that, as ω becomes large, it becomes impossible to obtain a solution of (174) if ν is restricted to a real number. In a subsequent paper [41], they found that when the real solution ceases to exist, a complex solution appears. They also found that when ν is complex, the real part of ν is always either an integer or half-integer. As an example, we show the value of ν as a function of Mω in Table 1.

Table 1 The value of ν for various value of Mω in the case s = −2, l = m = 2 and q = 0.

The fact that we only have an integer or half-integer as the real part of ν is strongly suggested from the property of g_n(ν) [37]. Let us summarize the argument. We first convert v to y as v = p/2 + iy, where p is an arbitrary integer and y is an arbitrary complex number. We note that for an arbitrary integer n,

$$\begin{gathered} \;\;\;\;\;\;\;\;\;\;(\beta _n^{p/2 + iy})* = \beta _{ - n - p - 1}^{p/2 + iy*}, \hfill \\ (\alpha _n^{p/2 + iy}\alpha _{n + 1}^{p/2 + iy})* = \alpha _{ - n - p - 2}^{p/2 + iy*}\gamma _{ - n - p - 1}^{p/2 + iy*}, \hfill \\ \end{gathered} $$

((175))

where * denotes complex conjugation. From these relations, we find that if v - p/2 + iy is a solution of g_n(v) - 0, we have

$$0 = ({g_n}(p/2 + iy))* = {g_{ - n - p - 1}}(p/2 + iy*) = {g_n}(p/2 - p - 2n - 1 + iy*),$$

((176))

where we used the relation, g_n+1 (v) - g_n (v+1). We see that in this case v′ = p/2-p-2n-1+iy* is also an solution of g_n = 0. As already discussed at the end of Section 4.2, when v is a solution of g_n - 0, v+k and -v-1+k with an arbitrary integer k are also solutions. We assume that there are no other solutions. Although we do not have a formal proof of it, numerical investigations suggest that this is so. Under this assumption, since both v - p/2 + iy and v′ - p/2 - p - 2n - 1 + iy* are solutions of g_n - 0, we have two possibilities about the property of y:

$$iy* = iy,\;\;\;\;\;\;\;\;\;or\;\;\;\;\;\;\;\;iy* = - iy.$$

((177))

In the former case, y is real. In this case, v is complex with real part p/2 (integer or half-integer). In the later case, y is pure imaginary. In this case, v is real.

It becomes possible to determine v in the wide range of ω by allowing Im(v) ≠ 0. The MST formalism is now very useful in the fully numerical evaluation of homogeneous solutions of the Teukolsky equation. As a first step, Fujita, Hikida and Tagoshi [38] considered generic bound geodesic orbits around a Kerr black hole and evaluate the energy and angular momentum flux to infinity as well as the rate of change of the Carter constant in a wide range of orbital parameters.

The critical value of ω when v becomes complex is not very small. The complex v does not appear in the analytic evaluation of v in the low frequency expansion in powers of ε - 2Mω. Thus, at the first glance, it seems impossible to express the complex v in the power series expansion of ε. However, Hikida et al. [52] pointed out that it is possible to evaluate sin² (πν) very accurately in terms of the power series expansion of ε, even if ω is larger than a critical value and v is complex. Such an analytical expression of v is very useful in the numerical root finding of Equation (174) as well as in the analytical calculation of the homogeneous solutions.

5 Gravitational Waves from a Particle Orbiting a Black Hole

Based on the ingoing wave functions discussed in Section 3 and 4, we can derive the gravitational wave energy and angular momentum flux emitted to infinity. The formula for the energy and the angular momentum luminosity to infinity are given by Equations (48) and (49). Since most of the calculations are very long, we show only the final results. In [71], some details of the calculations are summarized. We define the post-Newtonian expansion parameter by x - (MΩ_ϕ)^1/3, where M is the mass of the black hole and Ω_ϕ. is the orbital angular frequency of the particle. Since the parameter x is directly related to the observable frequency, this result can be compared with the results by another method easily.

5.1 Circular orbit around a Schwarzschild black hole

First, we present the gravitational wave luminosity for a particle in a circular orbit around a Schwarzschild black hole [100, 105]. In this case, Ω_ϕ. is given by Ω_ϕ = (M/r ³₀ )^1/2 ≡ Ω_c, where r₀ is the orbital radius, in standard Schwarzschild coordinates. The luminosity to $\mathcal{O}({x^{11}})$ is given by

$$\begin{gathered} \left\langle {\frac{{dE}}{{dt}}} \right\rangle = {\left( {\frac{{dE}}{{dt}}} \right)_N} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\; \times \left[ {1 - \frac{{1247}}{{336}}{x^2} + 4\pi {x^3} - \frac{{44711}}{{9072}}{x^4} - \frac{{8191}}{{672}}\pi {x^5}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left( {\frac{{6643739519}}{{69854400}} - \frac{{1712}}{{105}}\gamma + \frac{{16}}{3}{\pi ^2} - \frac{{3424}}{{105}}\ln 2 - \frac{{1712}}{{105}}\ln x} \right){x^6} - \frac{{16285}}{{504}}\pi {x^7} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left( { - \frac{{323105549467}}{{3178375200}}} \right. + \frac{{232597}}{{4410}}\gamma - \frac{{1369}}{{126}}{\pi ^2} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left. {\frac{{39931}}{{294}}\ln 2 - \frac{{47385}}{{1568}}\ln 3 + \frac{{232597}}{{4410}}\ln x} \right){x^8} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left( {\frac{{265978667519}}{{745113600}}\pi - \frac{{6848}}{{105}}\pi \gamma - \frac{{13696}}{{105}}\pi \ln 2 - \frac{{6848}}{{105}}\pi \ln x} \right){x^9} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left( { - \frac{{2500861660823683}}{{2831932303200}} + \frac{{916628467}}{{7858620}}\gamma - \frac{{424223}}{{6804}}{\pi ^2}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left. {\frac{{83217611}}{{1122660}}\ln 2 + \frac{{47385}}{{196}}\ln 3 + \frac{{916628467}}{{7858620}}\ln x} \right){x^{10}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left( {\frac{{8399309750401}}{{101708006400}}\pi + \frac{{177293}}{{1176}}\pi \gamma } \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left. {\frac{{8521283}}{{17640}}\pi \ln 2 - \frac{{142155}}{{784}}\pi \ln 3 + \frac{{177293}}{{1176}}\pi \ln x} \right)\left. {{x^{11}}} \right], \hfill \\ \end{gathered} $$

((178))

where (dE/dt)_N is the Newtonian quadrupole luminosity given by

$${\left( {\frac{{dE}}{{dt}}} \right)_N} = \frac{{32{\mu ^2}{M^3}}}{{5r_0^5}} = \frac{{32}}{5}{\left( {\frac{\mu }{M}} \right)^2}{x^{10}}.$$

((179))

This is the 5.5PN formula beyond the lowest, Newtonian quadrupole formula. We can find that our result agrees with the standard post-Newtonian results up to $\mathcal{O}({x^5})$ [13] in the limit μ/M ≪ 1.

5.2 Circular orbit on the equatorial plane around a Kerr black hole

Next, we consider a particle in a circular orbit on the equatorial plane around a Kerr black hole [101]. In this case, the orbital angular frequency Ω_ϕ. is given by

$${\Omega _\varphi } = {\Omega _c}[1 - q{\upsilon ^3} + {q^2}\upsilon {}^6 + \mathcal{O}({\upsilon ^9})],$$

((180))

where Ω_c, is the orbital angular frequency of the circular orbit in the Schwarzschild case, v = (M/r₀)^1/2, q = a/M, and r₀ is the orbital radius in the Boyer-Lindquist coordinate. The effect of the angular momentum of the black hole is given by the corrections depending on the parameter q. Here, q is arbitrary as long as |q| < 1. The luminosity is given up to $\mathcal{O}({x^8})$ (4PN order) by

$$\begin{gathered} \left\langle {\frac{{dE}}{{dt}}} \right\rangle = {\left( {\frac{{dE}}{{dt}}} \right)_N}\left[ {1 + (q - independent\;terms) - \frac{{11}}{4}q{x^3} + \frac{{33}}{{16}}{q^2}{x^4} - \frac{{59}}{{16}}q{x^5}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left( { - \frac{{65}}{6}\pi q + \frac{{611}}{{504}}{q^2}} \right){x^6} + \left( {\frac{{162035}}{{3888}}q + \frac{{65}}{8}\pi {q^2} - \frac{{71}}{{24}}{q^3}} \right){x^7} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left( { - \frac{{359}}{{14}}\pi q + \frac{{22667}}{{4536}}{q^2} + \frac{{17}}{{16}}{q^4}} \right){x^8}} \right]. \hfill \\ \end{gathered} $$

((181))

5.3 Waveforms in the case of circular orbit

In the previous two subsections, we only considered the luminosity formulas for the energy and the angular momentum. Here, focusing on circular orbits, we review the previous calculation of the gravitational waveforms.

On the other hand, the gravitational waveforms have also been calculated. In the case of circular orbit around a Schwarzschild black hole, Poisson [83] derived the 1.5PN waveform and Tagoshi and Sasaki [100] derived the 4PN waveform. These were done by using the post-Newtonian expansion of the Regge-Wheeler equation discussed in Section 3. Recently, Fujita and Iyer [39] derived the 5.5PN waveform by using the MST formalism. They also discussed factorized re-summed waveforms which is useful to obtain better agreement with accurate numerical data.

In the case of circular orbit around a Kerr black hole, Poisson [83] derived the 1.5PN waveform under the assumption of slow rotation of the black hole. In [94] and [101], although the luminosity was derived up to 2.5PN and 4PN order respectively, the waveform was not derived up to the same order. Recently, Tagoshi and Fujita [97] computed the all multipolar modes ${\tilde Z_{lm\omega }}$ necessary to derive the waveform up to 4PN order, and the results were used to derive the factorized, resummed, multipolar waveform in [78].

From Equations (46) and (47), we have

$$\begin{gathered} {h_ + } - i{h_ \times } = - \frac{2}{r}\sum\limits_{\ell mn} {{Z_{\ell m{\omega _n}}}\frac{{{}_{ - 2}S_{\ell m}^{a{\omega _n}}}}{{\sqrt {2\pi } }}{e^{i{\omega _n}(r* - t) + im\varphi }}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\; \equiv \sum\limits_{\ell m} {(h_{\ell m}^ + - ih_{\ell m}^ \times )} \hfill \\ \end{gathered} $$

((182))

Since the formulas for the waveform are very complicated, we only show the mode for ℓ = m = 2 up to 4PN order in the Schwarzschild case. We define ζ ^+,×_ℓ-m as [83]

$$h_{\ell m}^{ + , \times } + h_{\ell , - m}^{ + , \times } = - \frac{\mu }{r}{(M{\Omega _\varphi })^{2/3}}\varsigma _{\ell ,m}^{ + , \times }.$$

((183))

$$\begin{gathered} \varsigma _{2,2}^ + = (3 + \cos (2\theta ))\left( {\cos (2\psi ) - \frac{{107{\upsilon ^2}\cos (2\psi )}}{{42}}} \right) \hfill \\ \;\;\; + {\upsilon ^3}\left( {2\pi \cos (2\psi ) + \left( { - \frac{{17}}{3} + 4\ln 2} \right)\sin (2\psi )} \right) - \frac{{2173{\upsilon ^4}\cos (2\psi )}}{{1512}} \hfill \\ \;\;\; + {\upsilon ^5}\left( {\frac{{ - 107\pi \cos (2\psi )}}{{21}} + \left( {\frac{{1819}}{{126}} - \frac{{214\ln 2}}{{21}}} \right)\sin (2\psi )} \right) \hfill \\ \;\;\; + {\upsilon ^6}\left( {\cos (2\psi )\left( {\frac{{49928027}}{{1940400}} - \frac{{856\gamma }}{{105}} + \frac{{2{\pi ^2}}}{3} + \frac{{668\ln 2}}{{105}} - 8{{(\ln 2)}^2} - \frac{{856\ln \upsilon }}{{105}}} \right)} \right. \hfill \\ \;\;\;\;\;\;\;\left. { + \left( {\frac{{ - 254\pi }}{{35}} + 8\pi \ln 2} \right)\sin (2\psi )} \right) \hfill \\ \;\;\; + {\upsilon ^7}\left( {\frac{{ - 2173\pi \cos (2\psi )}}{{756}} + \left( {\frac{{36941}}{{4536}} - \frac{{2173\ln 2}}{{378}}} \right)\sin (2\psi )} \right) \hfill \\ \;\;\; + {\upsilon ^8}\left( {\cos (2\psi )\left( { - \frac{{326531600453}}{{12713500800}} + \frac{{45796\gamma }}{{2205}} - \frac{{107{\pi ^2}}}{{63}}} \right.} \right. \hfill \\ \;\;\;\;\;\; - \left. {\frac{{35738\ln 2}}{{2205}} + \frac{{428{{(\ln 2)}^2}}}{{21}} + \frac{{45796\ln \upsilon }}{{2205}}} \right) \hfill \\ \;\;\;\;\;\;\left. {\left. { + \left( {\frac{{13589\pi }}{{735}} - \frac{{428\pi \ln 2}}{{21}}} \right)\sin (2\psi )} \right)} \right), \hfill \\ \end{gathered} $$

((184))

$$\begin{gathered} \varsigma _{2,2}^ \times = 4\cos (\theta )\left( {in(2\psi ) - \frac{{107{\upsilon ^2}\sin (2\psi )}}{{42}} + {\upsilon ^3}\left( {\cos (2\psi )\left( {\frac{{17}}{3} - 4\log (2)} \right) + 2\pi \sin (2\psi )} \right)} \right. \hfill \\ \;\; - \frac{{2173{\upsilon ^4}\sin (2\psi )}}{{1512}} + {\upsilon ^5}\left( {\cos (2\psi )\left( { - \frac{{1819}}{{126}} + \frac{{214\log (2)}}{{21}}} \right) - \frac{{107\pi \sin (2)}}{{21}}} \right) \hfill \\ \;\; + {\upsilon ^6}\left( {\cos (2\psi )\left( {\frac{{254\pi }}{{35}} - 8\pi \log (2)} \right)} \right. \hfill \\ \;\left. {\; + \left( {\frac{{49928027}}{{1940400}} - \frac{{856\lambda }}{{105}} + \frac{{2{\pi ^2}}}{3} + \frac{{668\log (2)}}{{105}} - 8\log {{(2)}^2} - \frac{{856\log (\upsilon )}}{{105}}} \right)\sin 2(\psi )} \right) \hfill \\ \;\; + \;{\upsilon ^7}\left( {\cos (2\psi )\left( { - \frac{{36941}}{{4536}} + \frac{{2173\log (2)}}{{378}} - \frac{{2173\pi \sin (2\psi )}}{{756}}} \right)} \right) \hfill \\ \;\; + \;{\upsilon ^8}\left( {\cos (2\psi )\left( { - \frac{{13589\pi }}{{735}} + \frac{{428\pi (2)}}{{21}}} \right)\left( { + \frac{{326531600453}}{{12713500800}} + \frac{{45796\gamma }}{{2205}}} \right.} \right. \hfill \\ \;\;\;\left. {\left. {\left. { - \frac{{107{\pi ^2}}}{{63}} - \frac{{35738\log (2)}}{{2205}} + \frac{{428\log {{(2)}^2}}}{{21}} + \frac{{45796\log (\upsilon )}}{{2205}}} \right)\sin (2\psi )} \right)} \right), \hfill \\ \end{gathered} $$

((185))

where

$$\psi = \Omega (t - r*) - \varphi - 2{\upsilon ^3}(\gamma + 2\ln 2 + 3\ln \upsilon ).$$

((186))

Other modes are given in [100] up to 4PN order. (Note however the following errors which were pointed out in [61, 17, 4, 5, 39]. The signs of ζ ⁺_ℓm are opposite. The sign of ζ ^×_7,3 is also opposite. ζ ^×_8,7 , and ζ ⁺_10,6 have errors and the corrected formulas are given in Equations(6.6) and (6.7) in [39].) In the literature on the post-Newtonian approximation [17, 4, 5], the post-Newtonian waveforms are express in terms of H ^(n/2)_+,× ) defined as

$${h_{ + , \times }} = \frac{{2\mu x}}{r}\sum\limits_n {{x^{n/2}}H_{ + , \times }^{(n/2)}} ,$$

((187))

where x - (MΩ_ϕ)^2/3. The expression of H ^(n/2)_+,× up to 5.5PN order are given in [39].

5.4 Slightly eccentric orbit around a Schwarzschild black hole

Next, we consider a particle in slightly eccentric orbit on the equatorial plane around a Schwarzschild black hole (see [71], Section 7). We define r₀ as the minimum of the radial potential R(r)/r⁴. We also define an eccentricity parameter ε from the maximum radius of the orbit r_max, which is given by r_max = r₀(1 + e). These conditions are explicitly given by

$$\frac{{\partial (R/{r^4})}}{{\partial r}}\left| {_{r = {r_0}}} \right. = 0,\;\;\;and\;R({r_0}(1 + e)) = 0.$$

((188))

We assume e ≪ 1. In this case, Ω_ϕ is given to $\mathcal{O}({\epsilon^2})$ by

$${\Omega _\varphi } = {\Omega _c}[1 - f(\upsilon ){e^2}],\;\;\;f(\upsilon ) = \frac{{3(1 - 3{\upsilon ^2})(1 - 8{\upsilon ^2})}}{{2(1 - 2{\upsilon ^2})(1 - 6{\upsilon ^2})}},$$

((189))

where Ω_c = (M/r ³₀ )^1/2 is the orbital angular frequency in the circular orbit case. We now present the energy and angular momentum luminosity, accurate to $\mathcal{O}({\epsilon^2})$ and to $\mathcal{O}({x^8})$ beyond Newtonian order. They are given by

$$\begin{gathered} \left\langle {\frac{{dE}}{{dt}}} \right\rangle = {\left( {\frac{{dE}}{{dt}}} \right)_N}\left\{ {1 + (e - independent\;terms)} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + {e^2}\left[ {\frac{{157}}{{24}} - \frac{{6781}}{{168}}{x^2} + \frac{{2335}}{{48}}\pi {x^3} - \frac{{14929}}{{189}}{x^4} - \frac{{773}}{3}\pi {x^5}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left( {\frac{{156066596771}}{{69854400}} - \frac{{106144}}{{315}}\gamma } \right. + \frac{{992}}{9}{\pi ^2} - \frac{{80464}}{{315}}\ln 2 \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. { - \frac{{234009}}{{560}}\ln 3 - \frac{{106144}}{{315}}\ln x} \right){x^6} - \frac{{32443727}}{{48384}}\pi {x^7} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left( { - \frac{{3045355111074427}}{{671272842240}} + \frac{{507208}}{{245}}\gamma } \right. - \frac{{31271}}{{63}}{\pi ^2} - \frac{{151336}}{{441}}\ln 2 \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left. {\left. {\frac{{12887991}}{{3920}}\ln 3 + \frac{{507208}}{{245}}\ln x} \right){x^8}} \right]} \right\}, \hfill \\ \end{gathered} $$

((190))

and

$$\begin{gathered} \left\langle {\frac{{d{J_z}}}{{dt}}} \right\rangle = {\left( {\frac{{dJ}}{{dt}}} \right)_N}\left\{ {1 + (e - independent\;terms)} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + {e^2}\left[ {\frac{{23}}{8} - \frac{{3259}}{{168}}{x^2} + \frac{{209}}{8}\pi {x^3} - \frac{{1041349}}{{18144}}{x^4} - \frac{{785}}{6}} \right.\pi {x^5} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left( {\frac{{91721955203}}{{69854400}} - \frac{{41623}}{{210}}\gamma + \frac{{389}}{6}{\pi ^2} - \frac{{24503}}{{210}}\ln 2} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. { - \frac{{78003}}{{280}}\ln 3 - \frac{{41623}}{{210}}\ln x} \right){x^6} - \frac{{91565}}{{168}}\pi {x^7} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left( { - \frac{{105114325363}}{{72648576}} + \frac{{696923}}{{630}}\gamma - \frac{{4387}}{{18}}{\pi ^2} - \frac{{7051}}{{10}}\ln 2} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\left. {\;\left. { + \frac{{3986901}}{{1960}}\ln 3 + \frac{{696923}}{{630}}\ln x} \right){x^8}} \right]} \right\}, \hfill \\ \end{gathered} $$

((191))

where (dJ/dt)_N is the Newtonian angular momentum flux expressed in terms of x,

$${\left( {\frac{{d{J_z}}}{{dt}}} \right)_N} = \frac{{32}}{5}{\left( {\frac{\mu }{M}} \right)^2}M{x^7},$$

((192))

and the e-independent terms in both (dE/dt) and <dJ/dt> are the same and are given by the terms in the case of circular orbit, Equation (178).

5.5 Slightly eccentric orbit around a Kerr black hole

Next, we consider a particle in a slightly eccentric orbit on the equatorial plane around a Kerr black hole [95, 96]. We define the orbital radius ro and the eccentricity in the same way as in the Schwarzschild case by

$$\frac{{\partial (R/{r^4})}}{{\partial r}}\left| {_{r = {r_0}} = 0,\;\;\;and\;R({r_0}(1 + e)) = 0.} \right.$$

((193))

We also assume e ≪ 1. In this case, Ω_ϕ. is given to $\mathcal{O}({\epsilon^2})$ by

$${\Omega _\varphi } = {\Omega _c}\left[ {1 - q{\upsilon ^3} + {e^2}\left( { - \frac{3}{2} + \frac{9}{2}{\upsilon ^2} - \frac{9}{2}q{\upsilon ^3} + 3(6 + q{}^2){\upsilon ^4} - 60q{\upsilon ^5}} \right) + \mathcal{O}({\upsilon ^6})} \right].$$

((194))

We now give the energy and angular momentum luminosity that are accurate to $\mathcal{O}({\epsilon^2})$ and to $\mathcal{O}({x^5})$ beyond Newtonian order:

$$\begin{gathered} \left\langle {\frac{{dE}}{{dt}}} \right\rangle = {\left( {\frac{{dE}}{{dt}}} \right)_N}\left\{ {1 - \frac{{1247}}{{336}}{x^2}} \right. - \left( {\frac{{11}}{4}q + 4\pi } \right){x^3} - \left( {\frac{{44711}}{{9072}} + \frac{{33}}{{16}}{q^2}} \right){x^4} - \left( {\frac{{59}}{{16}}q - \frac{{8191}}{{672}}\pi } \right){x^5} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + {e^2}\left[ {\frac{{157}}{{24}} - \frac{{6781}}{{168}}{x^2} + \left( { - \frac{{2009}}{{72}}q + \frac{{2335}}{{48}}\pi } \right){x^3} + \left( { - \frac{{14929}}{{189}} + \frac{{281}}{{16}}{q^2}} \right){x^4}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\left. { + \left( { + \frac{{3223}}{{168}}q - \frac{{773}}{3}\pi } \right){x^5}} \right]} \right\}, \hfill \\ \end{gathered} $$

((195))

$$\begin{gathered} \left\langle {\frac{{d{J_z}}}{{dt}}} \right\rangle = {\left( {\frac{{d{J_z}}}{{dt}}} \right)_N}\left\{ {1 - \frac{{1247}}{{336}}{x^2}} \right. - \left( {\frac{{11}}{4}q + 4\pi } \right){x^3} - \left( {\frac{{44711}}{{9072}} + \frac{{33}}{{16}}{q^2}} \right){x^4} - \left( {\frac{{59}}{{16}}q + \frac{{8191}}{{672}}\pi } \right){x^5} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + {e^2}\left[ {\frac{{23}}{8} - \frac{{3259}}{{168}}{x^2} + \left( { - \frac{{371}}{{24}}q + \frac{{209}}{8}\pi } \right){x^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\left. { + \left( { - \frac{{1041349}}{{18144}} + \frac{{171}}{{16}}{q^2} + \frac{{949}}{{56}}q - \frac{{785}}{6}\pi } \right){x^5}} \right]} \right\}. \hfill \\ \end{gathered} $$

((196))

5.6 Circular orbit with a small inclination from the equatorial plane around a Kerr black hole

Next, we consider a particle in a circular orbit with small inclination from the equatorial plane around a Kerr black hole [94]. In this case, apart from the energy ɛ and z-component of the angular momentum l_z, the particle motion has another constant of motion, the Carter constant C. The orbital plane of the particle precesses around the symmetry axis of the black hole, and the degree of precession is determined by the value of the Carter constant. We introduce a dimensionless parameter y defined by

$$y = \frac{{\hat C}}{{\hat l_z^2 + {a^2}(1 - {{\hat {\mathcal{E}}}^2})}}.$$

((197))

Given the Carter constant and the orbital radius r₀, the energy and angular momentum is uniquely determined by R(r) = 0, and ∂R(r)/∂r = 0. By solving the geodesic equation with the assumption γ ≪ 1, we find that y^1/2 is equal to the inclination angle from the equatorial plane. The angular frequency Ω_ϕ is determined to $\mathcal{O}(y)$ and $\mathcal{O}({v^4})$ as

$${\Omega _\varphi } = {\Omega _c}\left[ {1 - q{\upsilon ^3} + \frac{3}{2}y(q{\upsilon ^3} - {q^2}{\upsilon ^4}) + \mathcal{O}({\upsilon ^6})} \right].$$

((198))

We now present the energy and the z-component angular flux to $\mathcal{O}({v^5})$

$$\begin{gathered} \left\langle {\frac{{dE}}{{dt}}} \right\rangle = {\left( {\frac{{dE}}{{dt}}} \right)_N}\left[ {1 - \frac{{1247}}{{336}}{\upsilon ^2} + \left( {4\pi - \frac{{73}}{{12}}\left( {1 - \frac{y}{2}} \right)q} \right){\upsilon ^3} + \left( { - \frac{{44711}}{{9072}} + \frac{{33}}{{16}}{q^2} - \frac{{527}}{{96}}{q^2}y} \right){\upsilon ^4}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\; + \left( { - \frac{{8191}}{{672}}\pi + \frac{{3749}}{{336}}q\left( {1 - \frac{y}{2}} \right)} \right){\upsilon ^5}} \right]. \hfill \\ \end{gathered} $$

((199))

$$\begin{gathered} \left\langle {\frac{{d{J_z}}}{{dt}}} \right\rangle = {\left( {\frac{{d{J_z}}}{{dt}}} \right)_N}\left\{ {\left( {1 - \frac{y}{2}} \right) - \frac{{1247}}{{336}}\left( {1 - \frac{y}{2}} \right){\upsilon ^2} + \left[ {4\pi \left( {1 - \frac{y}{2}} \right) - \frac{{61}}{{12}}\left( {1 - \frac{y}{2}} \right)q} \right]} \right.{\upsilon ^3} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left[ { - \frac{{44711}}{{9072}}\left( {1 - \frac{y}{2}} \right) + \left( {\frac{{33}}{{16}} - \frac{{229}}{{32}}y} \right){q^2}} \right]\upsilon {}^4 \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. { + \left[ { - \frac{{8191}}{{672}}\left( {1 - \frac{y}{2}} \right)\pi + \left( {\frac{{417}}{{56}} - \frac{{4301}}{{224}}y} \right)q} \right]{\upsilon ^5}} \right\}. \hfill \\ \end{gathered} $$

((200))

Using Equation (198), we can express v in terms of x = (MΩ_ϕ)^1/3 as

$$\upsilon = x\left( {1 + \frac{q}{3}{x^3} + \frac{1}{2}y\left( { - q{x^3} + {q^2}{x^4}} \right)} \right).$$

((201))

We then express Equations (199) and (200) in terms of x as

$$\begin{gathered} \left\langle {\frac{{dE}}{{dt}}} \right\rangle = {\left( {\frac{{dE}}{{dt}}} \right)_N}\left\{ {1 - \frac{{1247}}{{336}}{x^2} + \left( {4\pi - \frac{{11}}{4}q - \frac{{47}}{{24}}qy} \right){x^3} + \left[ { - \frac{{44711}}{{9072}} + \left( {\frac{{33}}{{16}} - \frac{{47}}{{96}}y} \right)q{}^2} \right]{x^4}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. { + \left[ { - \frac{{8191}}{{672}}\pi + \left( { - \frac{{59}}{{16}} + \frac{{11215}}{{672}}y} \right)q} \right]{x^5}} \right\}, \hfill \\ \end{gathered} $$

((202))

$$\begin{gathered} \left\langle {\frac{{d{J_z}}}{{dt}}} \right\rangle = {\left( {\frac{{d{J_z}}}{{dt}}} \right)_N}\left\{ {\left( {1 - \frac{y}{2}} \right) - \frac{{1247}}{{336}}\left( {1 - \frac{y}{2}} \right){x^2} + \left[ {4\lambda \left( {1 - \frac{y}{2}} \right) - \left( {\frac{7}{2}y + \frac{{11}}{4}\left( {1 - \frac{y}{2}} \right)} \right)q} \right]{x^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left[ { - \frac{{44711}}{{9072}}\left( {1 - \frac{y}{2}} \right) + \left( {\frac{{33}}{{16}} - \frac{{117}}{{32}}y} \right){q^2}} \right]{x^4} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\;\; + \left[ { - \frac{{8191}}{{672}}\left( {1 - \frac{y}{2}} \right)\pi + \left( { - \frac{{59}}{{16}} + \frac{{687}}{{224}}y} \right)q} \right]{x^5}} \right\}. \hfill \\ \end{gathered} $$

((203))

5.7 Absorption of gravitational waves by a black hole

In this section, we evaluate the energy absorption rate by a black hole. The energy flux formula is given by [107]

$$\left( {\frac{{d{E_{hole}}}}{{dt\;d\Omega }}} \right) = \sum\limits_{\ell m} {\int {d\omega } } \frac{{{}_2S_{\ell m}^2}}{{2\pi }}\frac{{128\omega k({k^2} + 4{{\tilde \in }^2})({k^2} + 16{{\tilde \in }^2}){{(2M{r_ + })}^5}}}{{|C{|^2}}}|\tilde Z_{\ell m\omega }^H{|^2},$$

((204))

where $\tilde \epsilon = \kappa /(4{r_ + })$, and

$$\begin{gathered} |C{|^2} = ({(\lambda + 2)^2} + 4a\omega m - 4{a^2}{\omega ^2})({\lambda ^2} + 36a\omega m - 36{a^2}{\omega ^2}) \hfill \\ \;\;\;\;\;\;\;\;\; + (2\lambda + 3)(96{a^2}{\omega ^2} - 48a\omega m) + 144{\omega ^2}({M^2} - {a^2}). \hfill \\ \end{gathered} $$

((205))

In calculating $\tilde Z_{\ell m\omega }^\text{H}$, we need to evaluate the upgoing solution R^up, and the asymptotic amplitude of ingoing and upgoing solutions, B^inc, B^trans, and C^trans in Equations (19) and (20). Evaluation of the incident amplitude B^trans of the ingoing solution is essential in the calculation. Poisson and Sasaki [85] evaluated them, in the case of a circular orbit around the Schwarzschild black hole, up to $\mathcal{O}(\epsilon)$ beyond the lowest order, and obtained the energy flux at the lowest order, using the method we have described in Section 3. Later, Tagoshi, Mano, and Takasugi [98] evaluated the energy absorption rate in the Kerr case to $\mathcal{O}({v^8})$ beyond the lowest order using the method in Section 4. Since the resulting formula is very long and complicated, we show it here only to $\mathcal{O}({v^3})$ beyond the lowest order. The energy absorption rate is given by

$$\begin{gathered} {\left( {\frac{{dE}}{{dt}}} \right)_H} = \frac{{32}}{5}{\left( {\frac{\mu }{M}} \right)^2}{x^{10}}{x^5}\left[ { - \frac{1}{4}q - \frac{3}{4}{q^3} + \left( { - q - \frac{{33}}{{16}}{q^3}} \right)} \right.{x^2} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left( {2q{B_2} + \frac{1}{2} + \frac{{13}}{2}\kappa {q^2} + \frac{{35}}{6}{q^2} - \frac{1}{4}{q^4} + \frac{1}{2}\kappa + 3{q^4} + 6{q^3}{B_2}} \right){x^3}} \right]. \hfill \\ \end{gathered} $$

((206))

where

$${B_n} = \frac{1}{{2i}}\left[ {{\psi ^{(0)}}\left( {3 + \frac{{niq}}{{\sqrt {1 - q{}^2} }}} \right) - {\psi ^{(0)}}\left( {3 - \frac{{niq}}{{\sqrt {1 - q{}^2} }}} \right)} \right],$$

((207))

and ψ⁽ⁿ⁾(z) is the polygamma function. We see that the absorption effect begins at $\mathcal{O}({v^5})$ beyond the quadrupole formula in the case q ≠ 0. If we set q = 0 in the above formula, we have

$${\left( {\frac{{dE}}{{dt}}} \right)_H} = {\left( {\frac{{dE}}{{dt}}} \right)_N}({x^8} + \mathcal{O}({\upsilon ^{10}})),$$

((208))

which was obtained by Poisson and Sasaki [85].

We note that the leading terms in (dE/dt)_H are negative for q > 0, i.e., the black hole loses energy if the particle is co-rotating. This is because of the superradiance for modes with k < 0.

5.8 Adiabatic evolution of Carter constant for orbit with small eccentricity and small inclination angle around a Kerr black hole

In the Schwarzschild case, the particle’s trajectory is characterized by the energy E and the z-component angular momentum L_z. In the adiabatic approximation, the rates of change of E and L_z are equated with those radiated to infinity as gravitational waves or with those absorbed into the black hole horizon, in accordance with the conservation of E and L_z. On the other hand, in the Kerr case, the Carter constant, Q, is also necessary to specify the particle’s trajectory. In this case, the rate of change of Q is not directly related to the asymptotic gravitational waves. Mino [70] proposed a new method for evaluating the average rate of change of the Carter constant by using the radiative field in the adiabatic approximation. He showed that the average rate of change of the Carter constant as well as the energy and angular momentum can be obtained by the radiative field of the metric perturbation. The radiative field is defined as half the retarded field minus half the advanced field. Mino’s work gave a proof of an earlier work by Gal’tsov [46] in which the radiative field is used to evaluate the average rate of change of the energy and the angular momentum without proof. Inspired by this new development, it was demonstrated in [53] and [31] that the time-averaged rates of change of the energy and the angular momentum can be computed numerically for generic orbits. A first step toward explicit calculation of the rate of change of the Carter constant was done in the case of a scalar charged particle in [30]. After that a simpler formula for the average rate of change of the Carter constant was found in [90, 89]. This new formula relates the rate of change of the Carter constant to the flux evaluated at infinity and on the horizon. Based on the new formula, in [89], the rate of change of the Carter constant for orbits with small eccentricities and inclinations is derived analytically up to O(v⁵) by using the MK method discussed in Section 4. In Ref. [48], the method was extended to the case of the orbits with small eccentricity but arbitrary inclination angle.

First we show the results for the small inclination case [89]. We define the radius r₀ and the eccentricity e as

$$\frac{{dR}}{{dr}}\left| {_{r = {r_0}} = 0,\;\;\;and\;\;\;R({r_0}(1 + e)) = 0.} \right.$$

((209))

In Ref. [89], the parameter which expresses a small inclination from the equatorial plane is defined as

$$y' = \hat C/\hat l_z^2.$$

((210))

Note that this definition of y′ is different from y in Equation (197). By solving Equation (209) with respect to E and ${\hat l_z}$, we obtain

$$\hat {\mathcal{E}} = 1 - \frac{{{\upsilon ^2}}}{2} + \frac{{3{\upsilon ^4}}}{8} - q{\upsilon ^5} - {e^2}\left( {\frac{{{\upsilon ^2}}}{2} - \frac{{\upsilon 4}}{4} + 2q{\upsilon ^5}} \right) + \frac{1}{2} + q{\upsilon ^5}y' + q{\upsilon ^5}{e^2}y',$$

((211))

$$\begin{gathered} {{\hat l}_z} - {r_0}\upsilon \left[ {1 + \frac{{3{\upsilon ^2}}}{2} - 3q{\upsilon ^3} + \frac{{27{\upsilon ^4}}}{8} + {q^2}{\upsilon ^4} - \frac{{15q{\upsilon ^5}}}{2}} \right. \hfill \\ \;\;\;\;\;\; + {e^2}\left( { - 1 + \frac{{3{\upsilon ^2}}}{2} - 6q{\upsilon ^3} + \frac{{81{\upsilon ^4}}}{8} + \frac{{7{q^2}{\upsilon ^4}}}{2} - \frac{{63q{\upsilon ^5}}}{2}} \right) \hfill \\ \;\;\;\;\;\; + y'\left( {\frac{1}{2} - \frac{{3{\upsilon ^2}}}{4} + 3q{\upsilon ^3} - \frac{{27{\upsilon ^4}}}{{16}} - \frac{{3{q^2}{\upsilon ^4}}}{2} + \frac{{15q{\upsilon ^5}}}{2}} \right) \hfill \\ \;\;\;\;\;\; + \left. {{e^2}y'\left( {\frac{1}{2} - \frac{{3{\upsilon ^2}}}{4} + 6q{\upsilon ^3} - \frac{{81{\upsilon ^4}}}{{16}} - \frac{{19{q^2}{\upsilon ^4}}}{4} + \frac{{63q{\upsilon ^5}}}{2}} \right)} \right]. \hfill \\ \end{gathered} $$

((212))

By using Equations (4.9) and (4.12) in [89], we find the azimuthal angular frequency Ω_ϕ as

$$\begin{gathered} {\Omega _\varphi } = {\Omega _c}\left[ {1 - q{\upsilon ^3} + {e^2}\left( {\frac{3}{2} + \frac{{9{\upsilon ^2}}}{2} - \frac{{21q{\upsilon ^3}}}{2} + 18{\upsilon ^4} + 3{q^2}{\upsilon ^4} - 60q{\upsilon ^5}} \right)} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\left. { + \left( {\frac{{3q{\upsilon ^3}}}{2} - \frac{{3{q^2}{\upsilon ^4}}}{2}} \right)y' + \left( {\frac{{39q{\upsilon ^3}}}{4} - \frac{{51{q^2}{\upsilon ^4}}}{4} + \frac{{75q{\upsilon ^5}}}{2}} \right){e^2}y'} \right], \hfill \\ \end{gathered} $$

((213))

where v = (M/r₀)^1/2 and Ω_c = (M/r₀)^1/2. From the definition of x ≡ (MΩ_ϕ)^1/3, the parameter v is expressed with x as

$$\begin{gathered} \upsilon = x\left[ {1 + \frac{{q{x^3}}}{3} + {e^2}\left( { - \frac{1}{2} - \frac{{3{x^2}}}{2} + \frac{{7q{x^3}}}{3} - 6{x^4} - {q^2}{x^4} + \frac{{31q{x^5}}}{2}} \right)} \right. \hfill \\ \;\;\;\;\;\; + \left. {\left( { - \frac{{q{x^3}}}{2} + \frac{{{q^2}{x^4}}}{2}} \right)y' + {e^2}y'\left( { - \frac{{3q{x^3}}}{2} + \frac{{9{q^2}{x^4}}}{4} - \frac{{23q{x^5}}}{4}} \right)} \right]. \hfill \\ \end{gathered} $$

((214))

By using the above formula, we can express y defined in Equation (197) in terms of y′ above as

$$\begin{gathered} y = \frac{{\hat C}}{{\hat l_z^2 + {a^2}(1 - {{\hat {\mathcal{E}}}^2})}} \hfill \\ \;\;\; = \frac{{\hat l_z^2y'}}{{\hat l_z^2 + {a^2}(1 - {{\hat {\mathcal{E}}}^2})}} \hfill \\ \;\;\; = (1 - {q^2}{x^4} + 4{q^2}{x^6})y' + {e^2}( - {q^2}{x^4} + 12{q^2}{x^6})y'. \hfill \\ \end{gathered} $$

((215))

The average rate of change of ɛ, l_z and Q become

$$\begin{gathered} \left\langle {\frac{{d\mathcal{E}}}{{dt}}} \right\rangle = - \frac{{32}}{5}{\left( {\frac{\mu }{M}} \right)^2}{\upsilon ^{10}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; \times \left[ {1 - \frac{{1247}}{{336}}{\upsilon ^2} - \left( {\frac{{73}}{{12}}q - 4\pi } \right){\upsilon ^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left( {\frac{{44711}}{{9072}} - \frac{{33}}{{16}}{q^2}} \right){\upsilon ^4} + \left( {\frac{{3749}}{{336}}q - \frac{{8191}}{{672}}\pi } \right){\upsilon ^5} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left\{ {\frac{{277}}{{24}} - \frac{{4001}}{{84}}{\upsilon ^2} + \left( {\frac{{3583}}{{48}} - \frac{{457}}{4}q} \right){\upsilon ^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left( {42{q^2} - \frac{{1091291}}{{9072}}} \right){\upsilon ^4} + \left( {\frac{{58487}}{{672}}q - \frac{{364337}}{{1344}}\pi } \right){\upsilon ^5}} \right\}{e^2} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left( {\frac{{73}}{{24}}q{\upsilon ^3} - \frac{{527}}{{96}}{q^2}{\upsilon ^4} - \frac{{3749}}{{672}}q{\upsilon ^5}} \right)y' \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left( {\frac{{457}}{8}q{\upsilon ^3} - \frac{{5407}}{{48}}{q^2}{\upsilon ^4} - \frac{{58487}}{{1344}}q{\upsilon ^5}} \right){e^2}y'} \right], \hfill \\ \end{gathered} $$

((216))

$$\begin{gathered} \left\langle {\frac{{d{l_z}}}{{dt}}} \right\rangle = - \frac{{32}}{5}{\left( {\frac{\mu }{M}} \right)^2}M{\upsilon ^7} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; \times \left[ {1 - \frac{{1247}}{{336}}{\upsilon ^2} - \left( {\frac{{61}}{{12}}q - 4\pi } \right){\upsilon ^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left( {\frac{{44711}}{{9072}} - \frac{{33}}{{16}}{q^2}} \right){\upsilon ^4} + \left( {\frac{{417}}{{56}}q - \frac{{8191}}{{672}}\pi } \right){\upsilon ^5} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left\{ {\frac{{51}}{8} - \frac{{17203}}{{672}}{\upsilon ^2} + \left( { - \frac{{781}}{{12}}q + \frac{{369}}{8}\pi } \right){\upsilon ^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left( {\frac{{929}}{{32}}{q^2} - \frac{{1680185}}{{18144}}} \right){\upsilon ^4} + \left( {\frac{{1809}}{{224}}q - \frac{{48373}}{{336}}\pi } \right){\upsilon ^5}} \right\}{e^2} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left\{ { - \frac{1}{2} + \frac{{1247}}{{672}}{\upsilon ^2} + \left( {\frac{{61}}{8}q - 2\pi } \right){\upsilon ^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left. {\left( {\frac{{213}}{{32}}{q^2} - \frac{{44711}}{{18144}}} \right){\upsilon ^4} - \left( {\frac{{4301}}{{224}}q - \frac{{8191}}{{1344}}\pi } \right){\upsilon ^5}} \right\}y' \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left\{ { - \frac{{51}}{{16}} + \frac{{17203}}{{1344}}{\upsilon ^2} + \left( {\frac{{1513}}{{16}}q - \frac{{369}}{{16}}\pi } \right){\upsilon ^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left. {\left( {\frac{{1680185}}{{36288}} - \frac{{5981}}{{64}}{q^2}} \right){\upsilon ^4} - \left( {168q - \frac{{48373}}{{672}}\pi } \right){\upsilon ^5}} \right\}{e^2}y'} \right], \hfill \\ \end{gathered} $$

((217))

$$\begin{gathered} \left\langle {\frac{{dQ}}{{dt}}} \right\rangle = - \frac{{64}}{5}{\left( {\frac{\mu }{M}} \right)^3}{M^3}{\upsilon ^6} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; \times \left[ {1 - q\upsilon - \frac{{743}}{{336}}{\upsilon ^2} - \left( {\frac{{1637}}{{336}}q - 4\pi } \right){\upsilon ^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left( {\frac{{439}}{{48}}{q^2} - \frac{{129193}}{{18144}} - 4\pi q} \right){\upsilon ^4} + \left( {\frac{{151765}}{{18144}}q - \frac{{4159}}{{672}}\pi - \frac{{33}}{{16}}{q^3}} \right){\upsilon ^5} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left\{ {\frac{{43}}{8} - \frac{{51}}{8}q\upsilon - \frac{{2425}}{{224}}{\upsilon ^2} - \left( {\frac{{14869}}{{224}}q - \frac{{337}}{8}\pi } \right){\upsilon ^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left( {\frac{{453601}}{{4536}} - \frac{{3631}}{{324}}q{}^2 + \frac{{369}}{8}\pi q} \right){\upsilon ^4} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left( {\frac{{141049}}{{9072}}q - \frac{{38029}}{{672}}\pi - \frac{{929}}{{32}}{q^3}} \right){\upsilon ^5}} \right\}{e^2} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left\{ {\frac{1}{2}q\upsilon + \frac{{1637}}{{672}}q{\upsilon ^3} - \left( {\frac{{1355}}{{96}}{q^2} - 2\pi q} \right){\upsilon ^4}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left. {\left( {\frac{{151765}}{{36288}}q - \frac{{213}}{{32}}{q^3}} \right){\upsilon ^5}} \right\}y' \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left\{ {\frac{{51}}{{16}}q\upsilon + \frac{{14869}}{{448}}q{\upsilon ^3} + \left( {\frac{{369}}{{16}}\pi q - \frac{{33257}}{{192}}{q^2}} \right){\upsilon ^4}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left. {\left( { - \frac{{141049}}{{18144}}q + \frac{{5981}}{{64}}{q^3}} \right){\upsilon ^5}} \right\}{e^2}y'} \right]. \hfill \\ \end{gathered} $$

((218))

The rate of change of C becomes

$$\begin{gathered} \left\langle {\frac{{dC}}{{dt}}} \right\rangle = \left\langle {\frac{{dQ}}{{dt}}} \right\rangle - 2(a\mathcal{E} - {l_z})\left( {a\left\langle {\frac{{d\mathcal{E}}}{{dt}}} \right\rangle - \left\langle {\frac{{d{l_z}}}{{dt}}} \right\rangle } \right) \hfill \\ \;\;\;\;\;\;\;\;\; = - \frac{{64}}{5}{\left( {\frac{\mu }{M}} \right)^3}{M^3}{\upsilon ^6}y'\left[ {1 - \frac{{743}}{{336}}{\upsilon ^2} - \left( {\frac{{85}}{8}q - 4\pi } \right)\upsilon {}^3} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left( {\frac{{129193}}{{18144}} - \frac{{307}}{{96}}{q^2}} \right){\upsilon ^4} + \left( {\frac{{2553}}{{224}}q - \frac{{4159}}{{672}}\pi } \right){\upsilon ^5} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left\{ {\frac{{43}}{8} - \frac{{2425}}{{224}}{\upsilon ^2} + \left( {\frac{{337}}{8}\pi - \frac{{1793}}{{16}}q} \right){\upsilon ^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left( {\frac{{453601}}{{4536}} - \frac{{7849}}{{192}}{q^2}} \right){\upsilon ^4} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left. {\left( {\frac{{3421}}{{224}}q - \frac{{38029}}{{672}}\pi } \right){\upsilon ^5}} \right\}{e^2}} \right]. \hfill \\ \end{gathered} $$

((219))

We can use Equation (214) to express these equations in terms of x. We obtain

$$\begin{gathered} \left\langle {\frac{{d\mathcal{E}}}{{dt}}} \right\rangle = - \frac{{32}}{5}{\left( {\frac{\mu }{M}} \right)^2}{x^{10}}\left[ {1 - \frac{{1247{x^2}}}{{336}} + \left( {4\pi - \frac{{11q}}{4}} \right){x^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + \left( { - \frac{{44711}}{{9072}} + \frac{{33{q^2}}}{{16}}} \right){x^4} + \left( {\frac{{8191\pi }}{{672}} - \frac{{59q}}{{16}}} \right){x^5} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + {e^2}\left\{ {\frac{{157}}{{24}} - \frac{{6781x{}^2}}{{168}} + \left( {\frac{{2335\pi }}{{48}} - \frac{{2009q}}{{72}}} \right){x^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left( { - \frac{{14929}}{{189}} + \frac{{281{q^2}}}{{16}}} \right){x^4} + \left( { - \frac{{773\pi }}{3} + \frac{{3223q}}{{168}}} \right){x^5}} \right\} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + \left( { - \frac{{47q{x^3}}}{{24}} - \frac{{47{q^2}{x^4}}}{{96}} + \frac{{11215q{x^5}}}{{672}}} \right)y' \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + \left. {{e^2}y'\left( { - \frac{{617q{x^3}}}{{48}} - \frac{{1585{q^2}{x^4}}}{{96}} + \frac{{60187q{x^5}}}{{336}}} \right)} \right], \hfill \\ \end{gathered} $$

((220))

$$\begin{gathered} \left\langle {\frac{{d{l_z}}}{{dt}}} \right\rangle = - \frac{{32}}{5}{\left( {\frac{\mu }{M}} \right)^2}M{x^7}\left[ {1 - \frac{{1247{x^2}}}{{336}} + \left( {4\pi - \frac{{11q}}{4}} \right){x^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + \left( { - \frac{{44711}}{{9072}} + \frac{{33{q^2}}}{{16}}} \right){x^4} + \left( {\frac{{8191\pi }}{{672}} - \frac{{59q}}{{16}}} \right){x^5} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + {e^2}\left\{ {\frac{{23}}{8} - \frac{{3259x{}^2}}{{168}} + \left( {\frac{{209\pi }}{8} - \frac{{371q}}{{24}}} \right){x^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left( { - \frac{{1041349}}{{18144}} + \frac{{171{q^2}}}{{16}}} \right){x^4} + \left( { - \frac{{785\pi }}{6} + \frac{{949q}}{{56}}} \right){x^5}} \right\} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + {e^2}y'\left\{ { - \frac{{23}}{{16}} + \frac{{3259{x^2}}}{{336}} + } \right.\left( { - \frac{{209\pi }}{{16}} + \frac{{1057q}}{{48}}} \right){x^3} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left( {\frac{{1041349}}{{36288}} - \frac{{825{q^2}}}{{32}}} \right){x^4} + \left( {\frac{{785\pi }}{{12}} - \frac{{925q}}{{14}}} \right){x^5}} \right\} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + \left\{ { - \frac{1}{2} + \frac{{1247{x^2}}}{{672}} + \left( { - 2\pi + \frac{{71q}}{{24}}} \right){x^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left. {\left( {\frac{{44711}}{{18144}} - \frac{{101{q^2}}}{{32}}} \right){x^4} + \left( {\frac{{8191\pi }}{{1344}} + \frac{{687q}}{{224}}} \right){x^5}} \right\}y'} \right], \hfill \\ \end{gathered} $$

((221))

$$\begin{gathered} \left\langle {\frac{{dC}}{{dt}}} \right\rangle = - \frac{{64}}{5}{\left( {\frac{\mu }{M}} \right)^3}{M^3}{x^6}y'\left[ {1 - \frac{{743{x^2}}}{{336}} + \left( {4\pi - \frac{{69q}}{8}} \right){x^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + \left( { - \frac{{129193}}{{18144}} + \frac{{307{q^2}}}{{96}}} \right){x^4} + \left( {\frac{{4159\pi }}{{672}} - \frac{{11089q}}{{16}}} \right){x^5} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + {e^2}\left\{ {\frac{{19}}{8} - \frac{{7379x{}^2}}{{672}} + \left( {\frac{{193\pi }}{8} - \frac{{89q}}{2}} \right){x^3}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left. {\left( { - \frac{{1340159}}{{18144}} + \frac{{1209{q^2}}}{{64}}} \right){x^4} + \left( { - \frac{{34295\pi }}{{448}} + \frac{{502051q}}{{4032}}} \right){x^5}} \right\}} \right]. \hfill \\ \end{gathered} $$

((222))

We note that if we set y′ = 0, <dɛ/dt> and <dl_z/dt> agree, respectively, with the rate of emission of the energy and the azimuthal angular momentum radiated to infinity, Equations (195) and (196) in Section 5.5. We also note that by using the transformation from y to y′, Equation (215), we can directly show that Equations (202) and (203) in Section 5.6 agree respectively with <dɛ/dt> and <dl_z/dt> in Equations (220) and (221) with e = 0.

5.9 Adiabatic evolution of constants of motion for orbits with generic inclination angle and with small eccentricity around a Kerr black hole

The calculation in Section 5.8 was extended to orbits with generic inclination angle by Ganz et al. [48]. We specify the geodesics by the semi-latus rectum p and the eccentricity e and a dimensionless inclination parameter y′. The outer and inner turning point of the radial motion is here define as

$$R(p/1 - e)) = 0,\;\;\;\;\;R(p/1 + e)) = 0.$$

((223))

The inclination parameter is defined by

$$y' = \hat C/\hat l_z^2,$$

((224))

which is the same as Equation (210). We define $v = \sqrt {M/p}$. By solving these equations with respect to $\hat \varepsilon$ and ${\hat l_z}$, we obtain

$$\begin{gathered} \hat {\mathcal{E}} - 1 - \frac{1}{2}{\upsilon ^2} + \frac{8}{8}{\upsilon ^4} - qY{\upsilon ^5} + {e^2}\left\{ {\frac{1}{2}{\upsilon ^2} - \frac{3}{4}{\upsilon ^4} + 2qY{\upsilon ^5}} \right\}, \hfill \\ \frac{{{{\hat l}_z}}}{p} = \upsilon Y + \frac{{3Y}}{2}{\upsilon ^3} - 3q{Y^2}{\upsilon ^4} + \left( {{q^2}{Y^3} + \frac{{9EY}}{4}} \right){\upsilon ^5} \hfill \\ \;\;\;\;\;\; + {e^2}\left\{ {\frac{Y}{2}{\upsilon ^3} - q{Y^2}{\upsilon ^4} + \left( {{q^2}{Y^3} + \frac{{9Y}}{4}} \right){\upsilon ^5}} \right\}, \hfill \\ \end{gathered} $$

((225))

where

$$Y \equiv \frac{1}{{\sqrt {1 + y'} }} = \frac{{{{\hat l}_z}}}{{\hat C + \hat l_z^2}}.$$

((226))

The average rate of change of ɛ, l_z and C become up to O(e², v⁵),

$$\begin{gathered} \left\langle {\frac{{d\mathcal{E}}}{{dt}}} \right\rangle = - \frac{{32}}{5}{\left( {\frac{\mu }{M}} \right)^2}{\upsilon ^{10}}{(1 - {e^2})^{3/2}}\left[ {\left( {1 + \frac{{73}}{{24}}{e^2}} \right) - \left( {\frac{{1247}}{{336}} + \frac{{9181}}{{672}}{e^2}} \right){\upsilon ^2}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left( {\frac{{73Y}}{{12}} + \frac{{823Y}}{{24}}{e^2}} \right)q{\upsilon ^3} + \left( {4 + \frac{{1375}}{{48}}{e^2}} \right)\pi {\upsilon ^3} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left( {\frac{{44711}}{{9072}} + \frac{{172157}}{{2592}}{e^2}} \right){\upsilon ^4} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left( {\frac{{329}}{{96}} - \frac{{527{Y^2}}}{{96}} + \left\{ {\frac{{4379}}{{192}} - \frac{{6533{Y^2}}}{{192}}} \right\}{e^2}} \right){q^2}{\upsilon ^4} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left. {\left( {\frac{{8191}}{{672}} + \frac{{44531}}{{336}}{e^2}} \right)\pi {\upsilon ^5} + \left( {\frac{{3749Y}}{{336}} + \frac{{1759Y}}{{56}}{e^2}} \right)q{\upsilon ^5}} \right], \hfill \\ \left\langle {\frac{{d{l_z}}}{{dt}}} \right\rangle = - \frac{{32}}{5}{\left( {\frac{\mu }{M}} \right)^2}M{\upsilon ^7}{(1 - {e^2})^{3/2}}\left[ {\left( {Y + \frac{{7Y}}{8}{e^2}} \right) - \left( {\frac{{1247Y}}{{336}} + \frac{{425Y}}{{336}}{e^2}} \right){\upsilon ^2}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left( {\frac{{61}}{{24}} - \frac{{61{Y^2}}}{8} + \left\{ {\frac{{63}}{8} - \frac{{91{Y^2}}}{4}} \right\}{e^2}} \right)q{\upsilon ^3} + \left( {4Y + \frac{{97Y}}{8}{e^2}} \right)\pi {\upsilon ^3} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left( {\frac{{44711Y}}{{9072}} + \frac{{302893Y}}{{6048}}{e^2}} \right){\upsilon ^4} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left( {\frac{{57Y}}{{16}} - \frac{{45{Y^3}}}{8} + \left\{ {\frac{{201Y}}{{16}} - \frac{{37{Y^3}}}{2}} \right\}{e^2}} \right){q^2}{\upsilon ^4} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left( {\frac{{8191Y}}{{672}} + \frac{{48361Y}}{{1344}}{e^2}} \right)\pi {\upsilon ^5} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left. {\left( {\frac{{2633}}{{224}} - \frac{{4301{Y^2}}}{{224}} + \left\{ {\frac{{66139}}{{1344}} - \frac{{18419{Y^2}}}{{448}}} \right\}{e^2}} \right)q{\upsilon ^5}} \right], \hfill \\ \left\langle {\frac{{dC}}{{dt}}} \right\rangle = - \frac{{64}}{5}{\left( {\frac{\mu }{M}} \right)^3}{M^3}{\upsilon ^6}{(1 - {e^2})^{3/2}}(1 - {Y^2})\left[ {\left( {1 + \frac{7}{8}{e^2}} \right)} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left( {\frac{{743}}{{336}} - \frac{{23}}{{42}}{e^2}} \right){\upsilon ^2} - \left( {\frac{{85Y}}{8} + \frac{{211Y}}{8}{e^2}} \right)q{\upsilon ^3} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left( {4 + \frac{{97}}{8}{e^2}} \right)\pi {\upsilon ^3} - \left( {\frac{{129193}}{{18144}} + \frac{{84035}}{{1728}}{e^2}} \right){\upsilon ^4} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left( {\frac{{329}}{{96}} - \frac{{53{Y^2}}}{8} + \left\{ {\frac{{929}}{{96}} - \frac{{163{Y^2}}}{8}} \right\}{e^2}} \right){q^2}{\upsilon ^4} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \left. {\left( {\frac{{2553Y}}{{224}} - \frac{{553Y}}{{192}}{e^2}} \right)q{\upsilon ^5} - \left( {\frac{{4159}}{{672}} + \frac{{21229}}{{1344}}{e^2}} \right)\pi {\upsilon ^5}} \right]. \hfill \\ \end{gathered} $$

((227))

Here, a term (1 - e²)^3/2 is factored out. We can express v in terms of x ≡ (MΩ_ϕ)^1/3 by using Equation (3.15) in [48] as

$$\begin{gathered} \upsilon = x\left[ {1 + q{x^3}\left( { - \frac{2}{3} + Y} \right) + {q^2}{x^4}\left( {\frac{1}{4} + \frac{Y}{2} - \frac{{3{Y^2}}}{4}} \right)} \right. \hfill \\ \;\;\; + {e^2}\left( {\frac{1}{2} - {x^2} + q{x^3}\left( { - \frac{4}{3} + 3Y} \right) + {x^4}\left( { - 3 + {q^2}\left( {\frac{5}{8} + \frac{{5Y}}{4} - \frac{{19{Y^2}}}{8}} \right)} \right)} \right. \hfill \\ \;\;\; + \left. {q{x^5}(3 + 3Y))} \right]. \hfill \\ \end{gathered} $$

((228))

The average rate of change of ɛ, l_z and C are rewritten as

$$\begin{gathered} \left\langle {\frac{{d\mathcal{E}}}{{dt}}} \right\rangle = - \frac{{32}}{5}{\left( {\frac{\mu }{M}} \right)^2}{(1 - {e^2})^{3/2}}{x^{10}}\left[ {1 + \frac{{193{e^2}}}{{24}} + \left( { - \frac{{1247}}{{336}} - \frac{{30865{e^2}}}{{672}}} \right){x^2}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + {x^3}\left( {4\pi - \frac{{20q}}{3} + \frac{{47qY}}{{12}} + {e^2}\left( {\frac{{2623\pi }}{{48}} - \frac{{1145q}}{{18}} + \frac{{379qY}}{{12}}} \right)} \right) \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + {x^4}\left. {\left( { - \frac{{44711}}{{9072}} - \frac{{89{q^2}}}{{96}} + 5{q^2}Y - \frac{{193{q^2}{Y^2}}}{{96}}} \right)} \right) \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + {e^2}\left. {\left( { - \frac{{522439}}{{6048}} - \frac{{4165{q^2}}}{{192}} + \frac{{1205{q^2}Y}}{{24}} - \frac{{503{q^2}{Y^2}}}{{64}}} \right)} \right) \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\; + {x^5}\left( { - \frac{{8191\pi }}{{672}} + \frac{{1247q}}{{42}} - \frac{{11215qY}}{{336}}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left. {{e^2}\left( { - \frac{{370877\pi }}{{1344}} + \frac{{17723q}}{{42}} - \frac{{39199qY}}{{96}}} \right)} \right)} \right]. \hfill \\ \end{gathered} $$

((229))

$$\begin{gathered} \left\langle {\frac{{d{l_z}}}{{dt}}} \right\rangle = - \frac{{32}}{5}{\left( {\frac{\mu }{M}} \right)^2}M{x^7}{(1 - {e^2})^{3/2}}\left[ {Y + \frac{{35{e^2}Y}}{8} + {x^2}\left( { - \frac{{1247Y}}{{336}} - \frac{{16777{e^2}Y}}{{672}}} \right)} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + {x^3}\left( {\frac{{61q}}{{24}} + 4\pi Y - \frac{{14qY}}{3} - \frac{{5q{Y^2}}}{8} + {e^2}\left( {\frac{{247q}}{{12}} + \frac{{257\pi Y}}{8} - \frac{{329qY}}{{12}} - \frac{{51q{Y^2}}}{4}} \right)} \right) \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + {x^4}\left( { - \frac{{44711Y}}{{9072}} - \frac{{29{q^2}Y}}{{16}} + \frac{{7{q^2}{Y^2}}}{2} + \frac{{3{q^2}{Y^3}}}{8}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + {e^2}\left. {\left( { - \frac{{83963Y}}{{1296}} - 21{q^2}Y + \frac{{357{q^2}{Y^2}}}{{16}} + \frac{{399{q^2}{Y^3}}}{{32}}} \right)} \right) \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\; + {x^5}\left( { - \frac{{2633q}}{{224}} - \frac{{8191\pi Y}}{{672}} + \frac{{1247qY}}{{56}} - \frac{{3181q{Y^2}}}{{224}}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left. {{e^2}\left( { - \frac{{65029q}}{{448}} - \frac{{200413\pi Y}}{{1344}} + \frac{{10651qY}}{{56}} - \frac{{15065q{Y^2}}}{{448}}} \right)} \right)} \right]. \hfill \\ \end{gathered} $$

((230))

$$\begin{gathered} \left\langle {\frac{{dC}}{{dt}}} \right\rangle = - \frac{{64}}{5}{\left( {\frac{\mu }{M}} \right)^3}{M^3}{x^6}{(1 - {e^2})^{3/2}}(1 - {Y^2})\left[ {1 + \frac{{31{e^2}}}{8} + \left( { - \frac{{743}}{{336}} - \frac{{1201{e^2}}}{{84}}} \right){x^2}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + {x^3}\left( {4\pi - 4q - \frac{{37qY}}{8} + {e^2}\left( {\frac{{241\pi }}{8} - \frac{{43q}}{2} - \frac{{575qY}}{{16}}} \right)} \right) \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\; + {x^4}\left( { - \frac{{129193}}{{18144}} - \frac{{185{q^2}}}{{96}} + 3{q^2}Y\frac{{17{q^2}{Y^2}}}{8}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + {e^2}\left. {\left( { - \frac{{438271}}{{5184}} - 18{q^2} + \frac{{141{q^2}Y}}{8} + \frac{{385{q^2}{Y^2}}}{{16}}} \right)} \right) \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\; + {x^5}\left( { - \frac{{4159\pi }}{{672}} + \frac{{743q}}{{63}} - \frac{{4229qY}}{{672}}} \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left. {\left. {{e^2}\left( { - \frac{{19227\pi }}{{224}} + \frac{{1799q}}{{18}} + \frac{{3151qY}}{{96}}} \right)} \right)} \right]. \hfill \\ \end{gathered} $$

((231))

If we assume that the inclination angle is small and Y = 1 - y′/2 + O(y′²), we find that Equations (229) – (231) reduce respectively to (220) – (222) in Section 5.8. As discussed in [48], in the case of largely inclined orbits, the fundamental frequency of gravitational waves is expressed not only with Ω_ϕ but also the frequency of θ-ocillation, Ω_θ.

6 Conclusion

In this article, we described analytical approaches to calculate gravitational radiation from a particle of mass μ orbiting a black hole of mass M with M ≫ μ, based upon the perturbation formalism developed by Teukolsky. A review of this formalism was given in Section 2. The Teukolsky equation, which governs the gravitational perturbation of a black hole, is too complicated to be solved analytically. Therefore, one has to adopt a certain approximation scheme. The scheme we employed is the post-Minkowski expansion, in which all the quantities are expanded in terms of a parameter ε = 2Mω where ω is the Fourier frequency of the gravitational waves. For the source term given by a particle in bound orbit, this naturally gives the post-Newtonian expansion.

In Section 3, we considered the case of a Schwarzschild background. For a Schwarzschild black hole, one can transform the Teukolsky equation into the Regge-Wheeler equation. The advantage of the Regge-Wheeler equation is that it reduces to the standard Klein-Gordon equation in the flat-space limit, and hence it is easier to understand the post-Minkowskian or post-Newtonian effects. Therefore, we adopted this method in the case of a Schwarzschild background. However, the post-Minkowski expansion of the Regge-Wheeler equation is not quite systematic, and as one goes to higher orders, the equations to be solved become increasingly complicated. Furthermore, for a Kerr background, although one can perform a transformation similar to the Chandrasekhar transformation, it can be done only at the expense of losing the reality of the equation. Thus, the resulting equation is not quite suited for analytical treatments.

In Section 4, we described a different method, developed by Mano, Suzuki, and Takasugi [69, 68] that directly deals with the Teukolsky equation, and we considered the case of a Kerr background with this method. Although the method is mathematically rather complicated and it is hard to obtain physical insights into relativistic effects, it has great advantage in that it allows a systematic post-Minkowski expansion of the Teukolsky equation, even on the Kerr background. We gave a thorough review on how this method works and how it gives a systematic post-Minkowski expansion.

Finally, in Section 5, we recapitulated the results of calculations of the gravitational waves for various orbits that had been obtained by various authors using the methods described in Sections 3 and 4. These results are useful not only by themselves for the actual case of a compact star orbiting a supermassive black hole, but also because they give us useful insights into higher order post-Newtonian effects even for a system of equal-mass binaries.