text:
# Zoom & Whirl: Eccentric equatorial orbits around spinning black holes and their evolution under gravitational radiation reaction.

## I Introduction

### A Background

Binary star systems, consisting of compact objects such as black holes and neutron stars, are relatively strong sources of gravitational radiation and are expected to be prime sources for the terrestrial network of kilometer-sized interferometric gravitational wave detectors, which will soon be fully operational, or for space-based detectors such as the proposed LISA mission . In order to detect gravitational radiation and subsequently study the physics of these sources it is absolutely necessary to have a prior theoretical knowledge of their dynamics. This is especially true because of the method of matched filtering (see  for a recent review) that is likely to be employed in order to identify true gravitational wave signals “buried” inside the detector’s noisy output. The success of this method depends on the use of an accurate template of the incoming waveform.

This paper will focus on the case of extreme mass ratio systems, modelling a massive central object which is a spinning (Kerr) black hole while the orbiting body is “light” and compact enough to be considered as a test-particle moving in the gravitational field of its companion. There are two important reasons for studying such a model.

Firstly, due to the extreme mass ratio, the motion of the small mass can be accurately approximated by a geodesic trajectory (which is well known ) and the system’s gravitational radiation is well described by first-order black hole perturbation theory techniques. The celebrated Teukolsky formalism  has proven particularly successful for this task. One thus has the opportunity to make a detailed study of a fully relativistic celestial system. For this reason, black hole perturbative studies can be used as a test for numerical relativity simulations of two-body systems (and vice versa!) .

Secondly, in recent years there has been an accumulation of evidence of the existence of supermassive black holes (of mass range $`10^610^9M_{}`$) in galactic nuclei (including our own Milky Way) . It is expected that scattered stellar-mass $`110M_{}`$ compact objects from the surrounding stellar population will be captured by the central black hole as a result of two-body encounters and interactions with the inhomogeneities of the background gravitational potential. The same scenario can of course work equally well for normal stars; however they will soon be tidally disrupted as they approach the black hole ,  .

Once in a bound orbit, the compact object will slowly inspiral towards the central black hole due to the emission of gravitational radiation. As the frequency of the emitted waves scales as $`1/M`$ (where $`M`$ denotes the central black hole’s mass), they will potentially lie in the low-frequency band ($`10^510^1Hz`$) where LISA will have its peak sensitivity. The (still uncertain) estimate for the number of such events is around $`1/\text{year}`$, or better, out to a distance of $`1`$Gpc and they should be detectable by LISA, with typical signal to noise ratios of $`10100`$, , assuming the use of some optimal filtering technique, such as matched filtering .

A huge payoff from direct observations of such events is to be expected, provided we have an accurate a priori description of the emitted waveform. In principle, for instance, the black hole parameters (masses,spins) can be measured to a high accuracy. Similarly, information on the mass-function of compact stellar populations in galactic nuclei could be provided. Because the total luminosity of the source depends only on its mass it may be possible to work out the distance to the source, which would be very useful to cosmologists. Moreover, one might be able to identify the massive object as a Kerr black hole, as opposed to some other, more speculative object (for example a boson star ). This was demonstrated by Ryan  who showed in detail how the massive body’s multipole moments are encrypted in the waveform emitted by an orbiting particle. In the near term, precise numerical results in the low-mass-ratio limit will be useful for testing the accuracy of post-Newtonian (PN) derived templates aimed at ground based detectors like LIGO , .

LISA will monitor the last year of inspiral of a compact body into a massive black hole by tracking the phase of the emitted waveform. It has been suggested that for astrophysically likely scenarios, drag forces operating on the orbiting body due to gas accreting onto the black hole, will operate on a timescale much longer than the radiation reaction timescale of the particle , . Based on requirements that the initial highly eccentric orbit in which the particle finds itself as the result of some scattering event should have a small enough periastron so that the radiation reaction timescale is shorter than the timescale for a second scattering event at apastron, we expect that the initial periastron should be rather close, so that $`r_p<20M`$ while the apastron will extend to a distance $`10^410^6M`$ , ,. Newtonian order estimates suggest that although radiation reaction will considerably reduce this enormous initial eccentricity during the course of the inspiral, the eccentricity will remain finite and non-negligible when the particle enters the strong-field region of interest to this paper (see section VB below). Exactly how much eccentricity remains will depend critically on the initial periastron distance and is largely insensitive to the initial apastron distance (and thus to the initial eccentricity).

We can thus argue that for a sufficiently bound orbit the system of the massive black hole and the orbiting compact object will evolve under its own spacetime dynamics. This tends to justify our “black hole plus particle” model. However, even in this simplified picture there are problems. The particle, in general, will move along a non-equatorial eccentric orbit (as the galactic central stellar population is almost spherically symmetric, capture orbits of arbitrary inclination are to be expected). The Teukolsky formalism cannot, at present, deal with such orbits, for reasons discussed in , in particular the problem of determining the rate of change of the “Carter constant” of the motion due to the emission of gravitational waves (much effort, towards this goal, is being focused on building a framework for calculating the gravitational self-force acting on the orbiting particle ). For this reason we restrict our attention to equatorial orbits around the central body. In such a case the rate of change of the orbital parameters can be deduced by reading the gravitational wave fluxes for the energy and angular momentum at infinity and the black hole horizon.

There is, however, a factor that cannot be accounted for by the previous flux-balance argument which lies at the heart of our approach. As has been observed recently ,  the gravitational self-force contains a conservative piece which is not associated with any radiation emission. Although the effect of this conservative force is negligible (scaling as $`\mu ^2`$) over short timescales (say, one orbital period), it is conceivable that the same will not be true for the accumulated effect after $`10^410^5`$ orbits  (this is, roughly, the number of orbits that LISA will record).

Another issue that has risen recently concerns the possible difficulty in defining the notion of adiabaticity and averaged flux for generic, i.e. eccentric and non-equatorial, orbits in Kerr spacetime. This is related to the belief that generic orbits have no well-defined orbital periods as they show an apparently non-periodic behaviour. For that reason, it has been suggested  that an “ergodicity” criterion would be more appropriate. However, recent work suggests  that it is possible, after all, to rigorously define (by means of Hamilton-Jacobi theory) a triplet of fundamental frequencies for generic Kerr orbits. Consequently, one may still be able to define adiabaticity for these orbits too.

Serious complications can also arise at the level of free motion, where radiation reaction is neglected. In general, the small body will have its own intrinsic spin. In such a case, due to the coupling of the particle’s spin with the background gravitational field, the motion is no longer geodesic. Although the (specific) spin magnitude is small, i.e. $`S𝒪(\mu )`$, spin-induced effects could become important over timescales much longer than, say, one orbital period. A particularly dramatic possibility is that when the test-particle is allowed to have spin, “chaotic” features may appear in the orbital motion . Presently it is unclear whether chaotic behaviour will be important for extreme mass ratio systems likely to be observed by LISA.

When radiation reaction is “switched-on” in the spinning particle case, one finds, not surprisingly, that the radiative fluxes at infinity and the horizon are inadequate for determining the evolution of the orbit. This is, in part, due to the fact that there is no known analog of the Carter constant (so there is one less constant of motion available), and also due to the existence of additional spin-degrees of freedom. A Newtonian order, weak-field estimation for the radiative change of the spin has been worked out by Apostolatos et.al. . Some speculations of what could happen to circular orbits under strong field conditions can be found in ,. For generic orbits, most likely only a self-force calculation will be able to describe the full orbital evolution.

As we are still far away from dealing with all of these challenges we make two major simplifications for this paper, that the orbiting particle has no spin and that it always remains in an equatorial orbit. Although Ryan has shown , for orbits in the weak field region, that non-equatorial orbits are forced by radiation reaction towards becoming retrograde equatorial orbits, the effect is small. The effect remains small even in the strong field region, as was recently shown by Hughes . Precisely equatorial orbits, pro- or retrograde, will remain equatorial under radiation reaction. Therefore it is reasonable to expect that detectors such as LISA will actually observe signals from particles in near equatorial orbits.

Previous studies have shown that slightly eccentric orbits of particles around Schwarzschild  and Kerr black holes  and arbitrarily eccentric orbits around Schwarzschild black holes  decrease in their eccentricity until shortly before the innermost stable circular orbit (ISCO) when a point is reached after which the eccentricity begins to increase. The present work comes as an additional piece to this series of papers. Specifically, we consider equatorial eccentric orbits of particles around a Kerr black hole and study their evolution under gravitational radiation reaction. This class of orbits is not exactly what we would expect in reality, but it is an important step towards a more realistic view of gravitational waves from this type of low-frequency source, because it includes two very important features which we know will be present in all or most sources, black hole spin and orbital eccentricity. Our study is, in this respect, a useful companion piece to Hughes’ discussion of non-equatorial (but circular) orbits . Eccentric equatorial orbits were first investigated by Shibata  who calculated fluxes and waveforms, without, however, discussing the impact of radiation reaction on the orbital motion. Our approach is similar to previous papers investigating eccentric orbits around non-spinning black holes , and nearly circular orbits around spinning black holes  and the results are qualitatively similar to those of both papers. In addition, we compute gravitational waveforms produced by moderate/high eccentricity, strong-field orbits (not discussed in Shibata’s study ) which we call “zoom-whirl” orbits. We find that these waveforms are a very characteristic, though complex, signal that might be important from an observational point of view for the planned LISA space antenna.

In this paper we focus on the final part of the inspiral, when the particle is at small radii, relatively close to the last stable bound orbit. In consequence we deal with orbits with moderate eccentricities, between 0.1 and 0.7. In a future paper  we intend to study the full inspiral, thus expanding our scope to cover orbits with large radii and larger eccentricities, on the order of 1. In that paper we plan to present wavetrains and spectra associated with a long stretch of the inspiral, covering many orbital periods, along the lines of .

Our results in this paper can be summarised as follows. Moderate eccentricities will be a feature of the signals from many inspiralling compact binaries right up to the final plunge. Immediately before plunge there will be an eccentricity increasing phase in all cases, particularly noticeable for retrograde orbits. The total amount of eccentricity gained in this phase will generally be small, on the order of 10% or less for low-eccentricity ($`e<0.1`$) prograde orbits, but perhaps as much as 50% for low-eccentricity retrograde orbits. Orbits with moderate eccentricities will gain much less in eccentricity. Where $`e>0.3`$ and the orbit is prograde, zoom-whirl features will be prominent in the waveform in the very last stages of the inspiral. Where these orbits are observed from a position away from the polar axis of the source there will relatively strong high-frequency component to the signal due to beaming of higher multipoles in the radiation in the direction of the orbiting particle’s motion. One expects that these signals will present particular problems for signal analysis, a situation which may be ameliorated when a positive detection of the source has been made during the earlier part of the inspiral when the waveform, though highly eccentric, will be less complex.

### B Organisation of the paper

The remainder of this paper is organised as follows. In Part II we discuss the geodesic motion of eccentric equatorial orbits (Sections IIA & IIB), paying particular attention to the so-called “zoom-whirl” orbits (Section IIC). In those sections we define useful orbital parameters such as the semi-latus rectum $`p`$ and the eccentricity $`e`$. Some analytic approximations on the orbital periods and number of revolutions for a particle in an orbit close to becoming dynamically unstable, are presented in Section IID. Sections IIIA & IIIB of Part III contain a review of the Teukolsky-Sasaki-Nakamura formalism for the calculation of gravitational waveforms and fluxes. In Section IIIC, we give a preliminary discussion on the orbital evolution under radiation reaction (definition of adiabaticity, general formulae for the rate of change of orbital parameters). Part IV is entirely devoted to analytic results. Section IVA discusses the weak-field limit for the orbital parameter’s rates of change.In Section IVB we derive an approximate formula relating the energy flux to the angular momentum flux, emitted by orbits close to becoming unstable. We subsequently use this formula to find strong-field approximate expressions for the rate of change of $`p`$ and $`e`$. In Section IVC we study the particularly interesting family of (equatorial) horizon-skimming orbits that can exist around a rapidly rotating black hole. The main (numerical) results of this paper are contained in Part V. In Section VA, we sketch the methods used in our numerical code and, moreover, give estimates for the various introduced errors. In Section VB we give results on the averaged rate of change of the parameters $`p,e`$ (which determine the evolution of any given orbit). This allows us to draw conclusions for the “global” behaviour of bound equatorial orbits under the influence of radiation reaction. Section VC contains calculations of waveforms generated from some zoom-whirl orbits. Part V ends with a presentation of the spectral content of the radiation emitted at infinity and at the black hole horizon (Section VD). Our conclusions are summarised in Part VI, where we also discuss prospects for future work. Tables with samples of our numerical data can be found at the end of the paper’s main body. Three Appendices are devoted to some technical details. Throughout this paper we have adopted geometrised units ($`c=G=1`$).

## II Geodesic motion

### A Equations of motion

We start by considering a test body moving in a Kerr gravitational field. For the moment, we neglect any radiation reaction effects and focus on purely geodesic motion. Working in the usual Boyer-Lindquist coordinate frame, the equations of motion, specialised for an equatorial orbit, are given by ,

$`r^2{\displaystyle \frac{dr}{d\tau }}`$ $`=`$ $`\pm (V_r)^{1/2},`$ (1)
$`r^2{\displaystyle \frac{d\varphi }{d\tau }}`$ $`=`$ $`V_\varphi (aEL)+{\displaystyle \frac{aT}{\mathrm{\Delta }}},`$ (3)
$`r^2{\displaystyle \frac{dt}{d\tau }}`$ $`=`$ $`V_\mathrm{t}a(aEL)+{\displaystyle \frac{(r^2+a^2)T}{\mathrm{\Delta }}},`$ (5)
$`\theta (\tau )`$ $`=`$ $`\pi /2,`$ (7)
where $`T=E(r^2+a^2)La,V_r=T^2\mathrm{\Delta }(r^2+(LaE)^2),\mathrm{\Delta }=r^22Mr+a^2.`$ The two constants of motion $`E`$, $`L`$ denote the orbit’s specific energy and $`z`$-component of angular momentum (for notational simplicity we drop the subscript $`z`$ for the angular momentum). We have prograde (retrograde) orbits according to whether $`L>0(<0)`$ (note that at certain points, where there is no danger of confusion, we shall label retrograde orbits by a negative value for the spin parameter $`a`$). Moreover, since we shall be discussing bound orbits, $`0<E<1`$. A general bound equatorial orbit can be equivalently described  either by the constants $`E`$ and $`L`$ or by a semi-latus rectum $`p`$ and an eccentricity $`e`$ (with $`0e<1`$). The restriction on the values of $`p`$ is discussed below. We define these parameters in terms of the two turning points of the orbit ($`r_p`$ is the periastron and $`r_a`$ the apastron, see Fig. 1 for a typical illustration),

$$r_p=\frac{p}{1+e},r_a=\frac{p}{1e}.$$
(8)
A turning point $`r_\mathrm{o}`$ by definition satisfies $`V_r(r_\mathrm{o})=0`$, or explicitly,

$$(E^21)r^3+2Mr^2(x^2+a^2+2aEx)r+2Mx^2=0,$$
(9)
where we have further defined $`x=LaE`$. Writing this polynomial in the form $`(E^21)(rr_p)(rr_a)(rr_3)`$ we can immediately write an expression for the energy,

$$E=\left[1\left(\frac{M}{p}\right)(1e^2)\left\{1\frac{x^2}{p^2}(1e^2)\right\}\right]^{1/2}.$$
(10)
Similarly, the third root $`r_3`$ of (9) is found to be

$$r_3=\frac{2M(1e^2)x^2}{p^2(1E^2)}.$$
(11)
It then follows that,

$$x^2=\frac{N(p,e)\mathrm{\Delta }_x^{1/2}(p,e)}{2F(p,e)}.$$
(12)
The explicit forms of the functions $`N,F`$ and $`\mathrm{\Delta }_x`$ are given in Appendix A. In this expression, the upper (lower) sign corresponds to a prograde (retrograde) orbit. The same convention will be followed throughout the paper.

The radial coordinate can be parameterised as

$$r(\chi )=\frac{p}{1+e\mathrm{cos}\chi },$$
(13)
where $`\chi `$ is a monotonically varying parameter, running from $`\chi =0`$ (at $`r=r_p`$) to $`\chi =\pi `$ (at $`r=r_a`$) and finally up to $`\chi =2\pi `$ (back to $`r=r_p`$ ). The radial motion can be separated into two distinct branches, namely, the motion from $`r_p`$ to $`r_a`$, and the “inverse” motion from $`r_a`$ back to $`r_p`$ again. Integration of (5) gives

$$t(r)=\{\begin{array}{cc}\widehat{t}(r)\hfill & \text{first branch}\hfill \\ T_r\widehat{t}(r)\hfill & \text{second branch}\hfill \end{array}$$
(14)
where

$$\widehat{t}(r)=_{r_1}^r\frac{1}{r^2}\left(\frac{dr}{d\tau }\right)^1\left[ax+\frac{r^2+a^2}{\mathrm{\Delta }}(Er^2ax)\right]𝑑r$$
(15)
We have also denoted as $`T_r`$ the period of the radial motion. For the $`\varphi `$-motion we similarly write,

$$\varphi (r)=\{\begin{array}{cc}\widehat{\varphi }(r)\hfill & \text{first branch}\hfill \\ \mathrm{\Delta }\varphi \widehat{\varphi }(r)\hfill & \text{second branch},\hfill \end{array}$$
(16)
where

$$\varphi (r)=_{r_1}^r\frac{1}{r^2}\left(\frac{dr}{d\tau }\right)^1\left[x+\frac{a}{\mathrm{\Delta }}(Er^2ax)\right]𝑑r.$$
(17)
and $`\mathrm{\Delta }\varphi `$ is the change of $`\varphi `$ during an interval $`T_r`$. Both integrands in (14), (16) are (unphysically) divergent at the turning points, an undesirable feature in a numerical calculation. This difficulty can be avoided by choosing $`\chi `$ as the integration parameter. Using

$$\frac{dr}{d\tau }=\frac{e\mathrm{sin}\chi }{p}[x^2+a^2+2xaE\frac{2Mx^2}{p}(3+e\mathrm{cos}\chi )]^{1/2},$$
(18)
we get

$$\varphi (\chi )=_0^\chi 𝑑\chi ^{}\frac{\stackrel{~}{V}_\varphi (\chi ^{},p,e)}{J(\chi ^{},p,e)\stackrel{~}{V}_\mathrm{r}^{1/2}(\chi ^{},p,e)},$$
(19)
$$t(\chi )=_0^\chi 𝑑\chi ^{}\frac{\stackrel{~}{V}_\mathrm{t}(\chi ^{},p,e)}{J(\chi ^{},p,e)\stackrel{~}{V}_\mathrm{r}^{1/2}(\chi ^{},p,e)},$$
(20)
where

$`\stackrel{~}{V}_\mathrm{r}(\chi ,p,e)`$ $`=`$ $`x^2+a^2+2axE{\displaystyle \frac{2Mx^2}{p}}(3+e\mathrm{cos}\chi ),`$ (21)
$`\stackrel{~}{V}_\varphi (\chi ,p,e)`$ $`=`$ $`x+aE{\displaystyle \frac{2Mx}{p}}(1+e\mathrm{cos}\chi ),`$ (22)
$`\stackrel{~}{V}_\mathrm{t}(\chi ,p,e)`$ $`=`$ $`a^2E{\displaystyle \frac{2aMx}{p}}(1+e\mathrm{cos}\chi )+{\displaystyle \frac{Ep^2}{(1+e\mathrm{cos}\chi )^2}},`$ (23)
$`J(\chi ,p,e)`$ $`=`$ $`1{\displaystyle \frac{2M}{p}}(1+e\mathrm{cos}\chi )+{\displaystyle \frac{a^2}{p^2}}(1+e\mathrm{cos}\chi )^2.`$ (24)
The integrand quantities in (19), (20) are well behaved and, moreover, these equations are valid for both branches of the radial motion. The radial period is simply given by $`T_r=t(2\pi )=2t(\pi )`$, and similarly, $`\mathrm{\Delta }\varphi \varphi (2\pi )=2\varphi (\pi )`$.

A general bound equatorial orbit is the combination of two separable motions: the radial motion which is, strictly speaking, periodic (in the sense that the radial coordinate returns to its original value after a certain time interval $`T_r`$ has elapsed) and the azimuthal motion which is not purely periodic (in the sense that the $`\varphi `$-coordinate monotonically increases but, nevertheless, the orbit returns to the same configuration after $`\varphi `$ has increased by some value $`\mathrm{\Delta }\varphi `$). The former motion is known in classical mechanics  as “libration” while the latter motion is called “rotation”. For such a combination of motions, it is generally known that there is a fundamental period (the period of libration) which fully describes the motion (see Appendix B for further details). We shall, therefore, call $`T_r`$ the orbital period. The fact that the orbit is periodic in a strict sense will enable us to rigorously define adiabaticity when radiative effects are to be included.

In line with the foregoing discussion, we define the orbital frequency to be $`\mathrm{\Omega }_r=2\pi /T_r`$. We can similarly refer to the frequency of the $`\varphi `$-motion as $`\mathrm{\Omega }_\varphi =\mathrm{\Delta }\varphi /T_r`$. The gravitational waves emitted by our systems will have frequencies which depend on these orbital frequencies. Below we shall see how they form a spectrum of discrete frequencies parameterised by the following wave numbers: $`\mathrm{}`$, which identifies the multipole of the emitted waves ($`\mathrm{}=2`$ for quadrupole, for instance), $`m`$ which runs from $`\mathrm{}`$ to $`+\mathrm{}`$ and $`k`$ which counts the harmonics created by the linear composition of the two orbital frequencies. The frequency of the waves emitted by a given harmonic $`k`$ of a given multipolar contribution $`m`$ is

$$\omega =k\mathrm{\Omega }_r+\frac{m\mathrm{\Delta }\varphi }{2\pi }\mathrm{\Omega }_r.$$
(25)
In the calculational scheme to be outlined below we will evaluate the fluxes of energy and angular momentum which are carried by waves of a given frequency (that is, a given multipole and harmonic of the frequency spectrum) and sum the fluxes for all frequencies of the discreet spectrum to get the total radiated fluxes of these quantities.

### B Separatrix curve

In general eqn. (9) has three distinct real roots. The case with $`r_p=r_3`$ corresponds to a marginally stable orbit: once at the periastron, the particle will enter into a circular orbit of radius $`r_{\mathrm{isbo}}=r_p=r_3`$ (ISBO stands for Innermost Stable Bound Orbit). At this stage the orbit has become unstable, so that a slight inwards “push” will drive the particle to catastrophically plunge into the black hole. Therefore, stable bound orbits should satisfy $`r_3<r_p`$. This translates to the inequality,

$$x^2(1+e)(3e)<p^2.$$
(26)
We can imagine a division of the $`(p,e)`$ plane into regions of stable and unstable orbits. The boundary curve $`p_s(e)`$ satisfying the equality in (26), defines the separatrix of bound orbits. In Fig. 2 we illustrate separatrices for a variety of black hole spins. A sample of numerical data used to generate this figure can be found in Table I. As one might have anticipated, spinning up the black hole will cause the separatrix curve for prograde (retrograde) orbits to move to the left (right) with respect to the Schwarzschild curve $`p_s(a=0)=(6+2e)M`$ . This behaviour can be seen most easily by a slow rotation approximation to (26). At leading order we find,

$$p_s=(6+2e)M8a\left[\frac{1+e}{6+2e}\right]^{1/2}+𝒪(a^2).$$
(27)
On the other hand, as can be verified by direct substitution in (26), for extreme rotation ($`a=M`$) the prograde separatrix becomes $`p_s(e)=M(1+e)`$, i.e. for all eccentricities, the periastron “descends” into the black hole “throat” at $`r=M`$, but is still separated by a finite proper distance from the horizon itself .

### C Zoom-Whirl orbits

From the short discussion in the previous Section one can imagine that as the orbit gradually approaches the separatrix, the particle will spend a considerable amount of its orbital “life” close to the periastron (see Fig. 3). An approximation for $`T_\mathrm{r}`$ as $`pp_s`$, derived in the following Section, gives

$$T_\mathrm{r}\mathrm{ln}(pp_s),$$
(28)
which shows that the period will grow (and eventually diverge) as the separatrix is approached. In that region, the particle will trace a quasi-circular path before being reflected back to the apastron. Such behaviour will be particularly prominent for high eccentricity orbits: the particle will “zoom in” from its apastron position, and perform a certain number of quasi-circular revolutions (“whirls”) reaching the periastron (which should have a value close to $`r_{\mathrm{isbo}}(e)=p_s(e)/(1+e)`$). Finally, the particle will be reflected and “zoom out” towards the apastron again. We shall heuristically (but quite descriptively) name these orbits “Zoom-Whirl” orbits. They resemble a set of orbits known in the literature as homoclinic orbits . Zoom-whirl orbits can exist in both Kerr and Schwarzschild geometries, and their potential significance for the detection of gravitational waves by space-based instruments was first pointed out some years ago by Curt Cutler and Eric Poisson<sup>*</sup><sup>*</sup>*The name “Zoom-Whirl” originated with the work of these two at Caltech. It may have been suggested by Kip Thorne., who concluded that the small number of whirls in the Schwarzschild case made the phenomenon less interesting for spinless central bodies. But as we shall shortly see, they are more pronounced in the case of near-extreme Kerr black holes, for prograde orbits. A typical example of such an orbit is illustrated in Fig. 3, for the case of a rapidly spinning ($`a=0.99M`$) black hole.

It is straightforward to calculate the total number of azimuthal revolutions $`N_\mathrm{r}=\mathrm{\Delta }\varphi /2\pi `$ during one orbital period, by numerically integrating (19). Results obtained by such a calculation are presented in Fig. 4. In this figure we have considered orbits of a given eccentricity ($`e=0.9`$ and $`e=0.3`$) and for a variety of black hole spins. For all depicted cases, the smallest value of $`p`$ resides at the same distance $`\delta p`$ from the corresponding separatrix value $`p_s(e)`$. As can be seen, the number of revolutions increases as the separatrix is approached, in agreement with our intuitive expectations. In fact, an approximate formula (valid for $`pp_s`$) derived in Section IID shows that,

$$N_\mathrm{r}\mathrm{ln}(pp_s).$$
(29)
We can furthermore deduce that the “whirling” of the particle near the separatrix becomes more pronounced as the black hole spin increases. Although for small and moderate spins $`N_r`$ stays close to the corresponding Schwarzschild value, it grows rapidly as $`aM`$, basically due to the intense “frame-dragging” induced by the black hole’s rotation in the very strong field region close to the horizon which can be reached by particles in prograde orbits. The overall behaviour can be understood as an extreme example of perihelion advance (as in the celebrated case of the planet Mercury).

In principle, as (29) suggests, the number of revolutions can be made arbitrarily large irrespective of the black hole spin, provided the particle approaches the separatrix sufficiently closely. However, as we discuss in Section IIIC the adiabatic assumption upon which our formalism relies breaks down in this regime. Sufficiently close to the separatrix, radiation reaction makes a significant correction to the particle’s motion in each orbital period. Before long this causes the particle to cross the separatrix and plunge into the black hole. These transition/plunging regimes have been studied recently by Ori and Thorne  for the case of circular equatorial orbits in the Kerr geometry. More relevant to the present discussion is the work of O’Shaughnessy and Thorne  which concerns the transition regime of zoom-whirl orbits. They show that for the case of an extreme Kerr black hole and eccentricity close to unity, the particle may experience more than 20 whirls per orbit before plunging. These have to be added to the number of whirls performed during the adiabatic phase of the orbit.

In a realistic scenario, we should not expect to find (apart from chance cases where the particle enters a near-separatrix orbit as a result of its initial scattering) very high eccentricity zoom-whirl orbits, as it is well known that the orbit has a general tendency to circularise . However, despite the decrease in eccentricity over the greater part of the inspiral, a substantial amount of eccentricity will survive, in many cases, up to the point where the orbit is about to plunge. These orbits will probably become zoom-whirl orbits, especially when a rapidly spinning black hole is involved and especially for prograde orbits. Keep in mind that many scattered particles will be in highly non-equatorial orbits. Zoom-whirl behaviour should also be seen in these cases as there is still a separatrix present, close to which the particle can spend a considerable amount of time.

A compact body in a zoom-whirl orbit will spend a considerable fraction of the orbital period in strong field regions (it can even travel close to the event horizon if the central black hole is spinning rapidly enough) and hence will radiate strongly. Our numerical results together with analytic approximations, reveal that a good fraction of the averaged flux is radiated during the motion near the periastron. As the orbit approaches the separatrix it tends to radiate as if it was a circular orbit of angular frequency $`\mathrm{\Omega }_\varphi `$ (see also  for a similar statement in the Schwarzschild case). This is clear evidence that most of the radiation is coming from the whirl part of the orbit, during which the radius hardly changes and there is a single dominant frequency characterised by the azimuthal ($`\varphi `$-dependent) orbital period. However, the most important feature of a zoom-whirl orbit is the characteristic form of the gravitational wave it emits, which is a series of rapid “quasi-circular” oscillations separated by relatively “quiet” intervals. In Section VC below we calculate some waveforms of this type.

### D Approximations near the separatrix (I)

Orbits that reside near the separatrix of the $`(p,e)`$ plane are amenable to analytic approximation, basically due to the fact that the turning point $`r_p`$ is close to a local minimum of the radial potential $`V_r`$. In this Section we derive approximate expressions for $`T_\mathrm{r}`$ and $`\mathrm{\Delta }\varphi `$. We already know that (see eqn. (20)),

$`T_\mathrm{r}`$ $`=`$ $`2{\displaystyle _0^\pi }𝑑\chi {\displaystyle \frac{\stackrel{~}{V}_\mathrm{t}(\chi ,p,e)}{J(\chi ,p,e)\stackrel{~}{V}_\mathrm{r}^{1/2}(\chi ,p,e)}},`$ (30)
$`\mathrm{\Delta }\varphi `$ $`=`$ $`2{\displaystyle _0^\pi }𝑑\chi {\displaystyle \frac{\stackrel{~}{V}_\varphi (\chi ,p,e)}{J(\chi ,p,e)\stackrel{~}{V}_\mathrm{r}^{1/2}(\chi ,p,e)}}.`$ (33)
We take the “distance” $`ϵ=pp_s`$ from the separatrix to be small, i.e. $`ϵ/M1`$. Also, we shall exclude small eccentricity orbits (more precisely, orbits with $`eϵ/M`$), or marginally bound orbits ($`e1`$). The former case of nearly circular orbits has already been discussed in . In what follows, quantities with an “s” subscript are to be evaluated exactly at the separatrix. From (21) we get,

$$\stackrel{~}{V}_\mathrm{r}(\chi ,p,e)=\frac{M}{p_s}\{1+𝒪(ϵ)\}[ϵS+2ex_\mathrm{s}^2(1\mathrm{cos}\chi )+𝒪(ϵ^2,ϵ(1\mathrm{cos}\chi )],$$
(34)
where

$$S=2p_s(1+e)(3e)\left(\frac{x^2}{p}\right)_{p=p_s}.$$
(35)
We see that $`\stackrel{~}{V}_\mathrm{r}0`$ as $`ϵ0`$ and $`\chi =0`$ (i.e. the periastron “touches” the separatrix). At the same limit, $`\stackrel{~}{V}_\mathrm{t}`$ and $`J`$ remain nonzero. We can write then, at leading order in $`ϵ`$,

$`T_\mathrm{r}`$ $``$ $`\left({\displaystyle \frac{p_s}{M}}\right)^{1/2}{\displaystyle _0^\pi }𝑑\chi {\displaystyle \frac{A_\mathrm{t}(1\mathrm{cos}\chi )}{\left[ϵS+2ex_s^2(1\mathrm{cos}\chi )\right]^{1/2}}},`$ (36)
$`\mathrm{\Delta }\varphi `$ $``$ $`\left({\displaystyle \frac{p_s}{M}}\right)^{1/2}{\displaystyle _0^\pi }𝑑\chi {\displaystyle \frac{A_\varphi (1\mathrm{cos}\chi )}{\left[ϵS+2ex_s^2(1\mathrm{cos}\chi )\right]^{1/2}}}.`$ (39)
We have defined the functions

$`A_\mathrm{t}(y)`$ $`=`$ $`{\displaystyle \frac{\left[a^2E_\mathrm{s}(1+eey)^22aMx_\mathrm{s}(1+eey)^3/p_s+E_\mathrm{s}p_s^2\right]}{(1+eey)^2\left[12M(1+eey)/p_s+a^2(1+eey)^2/p_s^2\right]}},`$ (40)
$`A_\varphi (y)`$ $`=`$ $`{\displaystyle \frac{\left[x_\mathrm{s}+aE_\mathrm{s}2Mx_\mathrm{s}(1+eey)/p_s\right]}{\left[12M(1+eey)/p_s+a^2(1+eey)^2/p_s^2\right]}},`$ (43)
with argument $`y=1\mathrm{cos}\chi `$. In order to isolate the divergent pieces in the integrals (36), (39) we split the functions (40),(43)

$$A_{\mathrm{t},\varphi }(y)=A_{\mathrm{t},\varphi }(0)+B_{\mathrm{t},\varphi }(y).$$
(44)
These expressions are just Taylor expansions around the regular point $`y=0`$ (with $`B_{t,\varphi }`$ containing the first and all higher derivatives of $`A_{t,\varphi }`$). Not surprisingly, both functions $`B_{\mathrm{t},\varphi }(y)`$ take the form

$$B_{\mathrm{t},\varphi }(y)=ey\stackrel{~}{B}_{\mathrm{t},\varphi }(y).$$
(45)
Although we don’t write the functions $`\stackrel{~}{B}_{\mathrm{t},\varphi }(y)`$ explicitly here (as they don’t take a simple form and they are not needed in what follows) we have verified that $`\stackrel{~}{B}_{\mathrm{t},\varphi }(0)0`$. It follows that the contribution to the integrals from $`B_{\mathrm{t},\varphi }(1\mathrm{cos}\chi )`$ is finite when $`ϵ,\chi 0`$. On the other hand, the contribution from $`A_{\mathrm{t},\varphi }(0)`$ is found to be divergent at the same limit,

$$_0^\pi \frac{d\chi }{\left[ϵS+2ex_\mathrm{s}^2(1\mathrm{cos}\chi )\right]^{1/2}}=\frac{1}{2}(ex_\mathrm{s}^2)^{1/2}\mathrm{ln}\left[\frac{64ex_\mathrm{s}^2}{ϵS}\right]+𝒪(\frac{ϵ}{e}\mathrm{ln}\left[\frac{e}{ϵ}\right].)$$
(46)
Hence at leading order in $`ϵ`$ (therefore close to the separatrix),

$`T_\mathrm{r}`$ $``$ $`A_\mathrm{t}(0)\left[{\displaystyle \frac{(1+e)(3e)}{eMp_s}}\right]^{1/2}\mathrm{ln}\left[{\displaystyle \frac{64ep_s^2}{ϵS(1+e)(3e)}}\right],`$ (47)
$`\mathrm{\Delta }\varphi `$ $``$ $`A_\varphi (0)\left[{\displaystyle \frac{(1+e)(3e)}{eMp_s}}\right]^{1/2}\mathrm{ln}\left[{\displaystyle \frac{64ep_s^2}{ϵS(1+e)(3e)}}\right].`$ (50)
The divergence of $`T_\mathrm{r}`$ and $`\mathrm{\Delta }\varphi `$ at the separatrix is the result of the particle being trapped in an unstable circular orbit at the location of the minimum of the radial potential $`V_\mathrm{r}`$.

## III RADIATION REACTION: formulation of the problem

### A The Teukolsky formalism

In this paper, we shall employ Teukolsky’s formalism  for the calculation of gravitational fluxes and waveforms. His eponymous equation describes the evolution of linearised radiative perturbative fields in a Kerr geometry background. In particular, instead of dealing directly with metric perturbations, the Teukolsky formalism considers perturbations on the Weyl curvature scalar $`\psi _4`$. This quantity is a result of the projection of the Weyl tensor on the null vectors $`n^\alpha `$, $`\overline{m}^\beta `$ which are members of the Newman-Penrose null tetrad , that is $`\psi _4=C_{\alpha \beta \gamma \delta }n^\alpha \overline{m}^\beta n^\gamma \overline{m}^\delta `$. The feature that makes this formalism attractive to our problem is that the radiative fluxes (at infinity and at the horizon) as well as the two wave polarisations $`h_+`$,$`h_\mathrm{x}`$ can all be extracted from $`\psi _4`$. The “master” perturbation equation is separable in the Fourier domain by means of a decomposition

$$\psi _4(t,r,\theta ,\varphi )=\rho ^4\underset{\mathrm{}m}{}𝑑\omega e_2^{i\omega t+im\phi }S_\mathrm{}m^{a\omega }(\theta )R_{\mathrm{}m\omega }(r),$$
(51)
where $`\rho =ria\mathrm{cos}\theta `$. The radial function $`R_{\mathrm{}m\omega }(r)`$ satisfies the Teukolsky equation

$$\mathrm{\Delta }^2\frac{d}{dr}\left(\frac{1}{\mathrm{\Delta }}\frac{dR_{\mathrm{}m\omega }}{dr}\right)V(r)R_{\mathrm{}m\omega }=T_{\mathrm{}m\omega }.$$
(52)
The potential $`V(r)`$ is given by

$$V(r)=\frac{K^2+4i(rM)K}{\mathrm{\Delta }}+8i\omega r+\lambda ,$$
(53)
where $`K=(r^2+a^2)\omega ma`$ and $`\lambda =E_{lm}+a^2\omega ^22am\omega `$. The angular functions $`{}_{2}{}^{}S_{lm}^{a\omega }(\theta )`$ are $`s=2`$ spin-weighted spheroidal harmonics  which satisfy the following eigenvalue equation,

$`[{\displaystyle \frac{1}{\mathrm{sin}\theta }}{\displaystyle \frac{d}{d\theta }}\left\{\mathrm{sin}\theta {\displaystyle \frac{d}{d\theta }}\right\}+a^2\omega ^2\mathrm{cos}^2\theta {\displaystyle \frac{m^2}{\mathrm{sin}^2\theta }}+4a\omega \mathrm{cos}\theta `$ (54)
$`+{\displaystyle \frac{4m\mathrm{cos}\theta }{\mathrm{sin}^2\theta }}4\mathrm{cot}^2\theta 2+E_\mathrm{}m]_2S_\mathrm{}m^{a\omega }=0.`$ (55)
We have adopted the following normalisation for the spheroidal harmonics (hereafter we drop the subscript $`2`$ for notational simplicity),

$$_0^\pi |S_\mathrm{}m^{a\omega }|^2\mathrm{sin}\theta d\theta =1.$$
(56)
The source term $`T_{\mathrm{}m\omega }`$ present in (52) is constructed directly from the particle’s energy-momentum tensor and this is the point where the particle’s motion enters explicitly in the perturbation equation. Its explicit form is given below. Let us now return to the radial equation (52). A particular solution of this equation can be found in terms of two independent solutions $`R_{\mathrm{}m\omega }^{\mathrm{in}}`$, $`R_{\mathrm{}m\omega }^{\mathrm{up}}`$ of the homogeneous equation,

$$R_{\mathrm{}m\omega }(r)=\frac{R_{\mathrm{}m\omega }^{\mathrm{up}}(r)}{W}_{r_+}^r𝑑r^{}\frac{T_{\mathrm{}m\omega }(r^{})R_{\mathrm{}m\omega }^{\mathrm{in}}(r^{})}{\mathrm{\Delta }^2(r^{})}+\frac{R_{\mathrm{}m\omega }^{\mathrm{in}}(r)}{W}_r^+\mathrm{}𝑑r^{}\frac{T_{\mathrm{}m\omega }(r^{})R_{\mathrm{}m\omega }^{\mathrm{up}}(r^{})}{\mathrm{\Delta }^2(r^{})},$$
(57)
where $`W`$ the (constant) Wronskian $`W[\mathrm{\Delta }^{1/2}R_{\mathrm{}m\omega }^{\mathrm{in}},\mathrm{\Delta }^{1/2}R_{\mathrm{}m\omega }^{\mathrm{up}}]`$. The solutions $`R_{\mathrm{}m\omega }^{\mathrm{in}}`$, $`R_{\mathrm{}m\omega }^{\mathrm{up}}`$ are chosen such as to have, respectively, purely ingoing behaviour at the horizon, and purely outgoing behaviour at infinity. Explicitly,

$$R_{\mathrm{}m\omega }^{\mathrm{in}}\{\begin{array}{cc}\mathrm{\Delta }^2e^{ikr^{}}\hfill & \text{for }rr_+\hfill \\ r^3B^{\mathrm{out}}e^{i\omega r^{}}+r^1B^{\mathrm{in}}e^{i\omega r^{}}\hfill & \text{for }r+\mathrm{},\hfill \end{array}$$
(58)
$$R_{\mathrm{}m\omega }^{\mathrm{up}}\{\begin{array}{cc}C^{\mathrm{out}}e^{ikr^{}}+\mathrm{\Delta }^2C^{\mathrm{in}}e^{ikr^{}}\hfill & \text{for }rr_+,\hfill \\ r^3e^{i\omega r^{}}\hfill & \text{for }r+\mathrm{}\hfill \end{array},$$
(59)
where $`k=\omega ma/2Mr_+`$, $`r_+=M+(M^2a^2)^{1/2}`$ is the outer event horizon, and $`r_{}`$ is the usual tortoise coordinate defined by $`dr_{}/dr=(r^2+a^2)/\mathrm{\Delta }`$. From these expressions we have that $`W=2i\omega B^{\mathrm{in}}`$. The solution (57) describes ingoing waves at the horizon and outgoing waves at infinity as it should be required on physical grounds. That is,

$`R_{\mathrm{}m\omega }(rr_+)`$ $``$ $`{\displaystyle \frac{\mathrm{\Delta }^2e^{ikr^{}}}{2i\omega B^{\mathrm{in}}}}{\displaystyle _{r_+}^{\mathrm{}}}𝑑r^{}{\displaystyle \frac{T_{\mathrm{}m\omega }(r^{})R_{\mathrm{}m\omega }^{\mathrm{up}}(r^{})}{\mathrm{\Delta }^2(r^{})}}Z_{\mathrm{}m\omega }^{\mathrm{}}\mathrm{\Delta }^2(r)e^{ikr^{}}`$ (60)
$`R_{\mathrm{}m\omega }(r\mathrm{})`$ $``$ $`{\displaystyle \frac{r^3e^{i\omega r^{}}}{2i\omega B^{\mathrm{in}}}}{\displaystyle _{r_+}^{\mathrm{}}}𝑑r^{}{\displaystyle \frac{T_{\mathrm{}m\omega }(r^{}),R_{\mathrm{}m\omega }^{\mathrm{in}}(r^{})}{\mathrm{\Delta }^2(r^{})}}Z_{\mathrm{}m\omega }^\mathrm{H}r^3e^{i\omega r^{}}.`$ (63)
The source term $`T_{\mathrm{}m\omega }`$ is given by

$$T_{\mathrm{}m\omega }=4𝑑\mathrm{\Omega }𝑑t\rho ^5\overline{\rho }^1(B_2^{}+B_2^{})e^{im\phi +i\omega t}\frac{{}_{2}{}^{}S_{\mathrm{}m}^{a\omega }}{\sqrt{2\pi }},$$
(64)
where

$`B_2^{}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\rho ^8\overline{\rho }L_1[\rho ^4L_0(\rho ^2\overline{\rho }^1T_{nn})]`$ (65)
$``$ $`{\displaystyle \frac{1}{2\sqrt{2}}}\rho ^8\overline{\rho }\mathrm{\Delta }^2L_1[\rho ^4\overline{\rho }^2J_+(\rho ^2\overline{\rho }^2\mathrm{\Delta }^1T_{\overline{m}n})],`$ (66)
$`B_2^{}`$ $`=`$ $`{\displaystyle \frac{1}{4}}\rho ^8\overline{\rho }\mathrm{\Delta }^2J_+[\rho ^4J_+(\rho ^2\overline{\rho }T_{\overline{m}\overline{m}})]`$ (69)
$``$ $`{\displaystyle \frac{1}{2\sqrt{2}}}\rho ^8\overline{\rho }\mathrm{\Delta }^2J_+[\rho ^4\overline{\rho }^2\mathrm{\Delta }^1L_1(\rho ^2\overline{\rho }^2T_{\overline{m}n})].`$ (70)
We have defined the operators

$`L_s`$ $`=`$ $`_\theta +{\displaystyle \frac{m}{\mathrm{sin}\theta }}a\omega \mathrm{sin}\theta +s\mathrm{cot}\theta ,`$ (72)
$`J_+`$ $`=`$ $`_r+iK/\mathrm{\Delta }.`$ (75)
The quantities $`T_{nn}`$, $`T_{\overline{m}n},T_{\overline{m}\overline{m}}`$ are the result of the projection of the particle’s energy-momentum tensor $`T^{\mu \nu }`$ on the tetrad vectors, i.e. $`T_{nn}=T_{\mu \nu }n^\mu n^\nu `$ etc. The energy-momentum tensor for a particle in an arbitrary orbit $`(t,r(t),\theta (t),\varphi (t))`$ is given by

$$T^{\mu \nu }=\mu \frac{u^\mu u^\nu }{\mathrm{\Sigma }\mathrm{sin}\theta u^t}\delta (rr(t))\delta (\theta \theta (t))\delta (\varphi \varphi (t)),$$
(77)
where $`u^\mu =dx^\mu /d\tau `$ and $`\mathrm{\Sigma }=r^2+a^2\mathrm{cos}\theta `$. We obtain for the individual projections ,

$`T_{nn}`$ $`=`$ $`\mu {\displaystyle \frac{C_{nn}}{\mathrm{sin}\theta }}\delta (rr(t))\delta (\theta \theta (t))\delta (\varphi \varphi (t)),`$ (78)
$`T_{\overline{m}n}`$ $`=`$ $`\mu {\displaystyle \frac{C_{\overline{m}n}}{\mathrm{sin}\theta }}\delta (rr(t))\delta (\theta \theta (t))\delta (\varphi \varphi (t)),`$ (81)
$`T_{\overline{m}\overline{m}}`$ $`=`$ $`\mu {\displaystyle \frac{C_{\overline{m}\overline{m}}}{\mathrm{sin}\theta }}\delta (rr(t))\delta (\theta \theta (t))\delta (\varphi \varphi (t)),`$ (84)
with

$`C_{nn}`$ $`=`$ $`{\displaystyle \frac{1}{4\mathrm{\Sigma }^3}}(u^\mathrm{t})^1\left[E(r^2+a^2)aL+\mathrm{\Sigma }u^\mathrm{r}\right]^2,`$ (85)
$`C_{\overline{m}n}`$ $`=`$ $`{\displaystyle \frac{\rho }{2\sqrt{2}\mathrm{\Sigma }^2}}(u^\mathrm{t})^1\left[E(r^2+a^2)aL+\mathrm{\Sigma }u^\mathrm{r}\right]\left[i\mathrm{sin}\theta \left(aE{\displaystyle \frac{L}{\mathrm{sin}^2\theta }}\right)+\mathrm{\Theta }(\theta )\right],`$ (86)
$`C_{\overline{m}\overline{m}}`$ $`=`$ $`{\displaystyle \frac{\rho ^2}{2\mathrm{\Sigma }}}(u^\mathrm{t})^1\left[i\mathrm{sin}\theta \left(aE{\displaystyle \frac{L_z}{\mathrm{sin}^2\theta }}\right)+\mathrm{\Theta }(\theta )\right]^2.`$ (87)
The quantity $`\mathrm{\Theta }(\theta )`$ represents the effective latitudinal potential, i.e, $`(\mathrm{\Sigma }u^\theta )^2=\mathrm{\Theta }(\theta )`$. Substitution in (64) yields

$`T_{\mathrm{}m\omega }`$ $`=`$ $`{\displaystyle \frac{4\mu }{\sqrt{2\pi }}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑t{\displaystyle _0^\pi }𝑑\theta e^{i\omega tim\phi (t)}`$ (88)
$`\times `$ $`[{\displaystyle \frac{1}{2}}L_1^{}\left\{\rho ^4L_2^{}(\rho ^3S_\mathrm{}m^{a\omega })\right\}C_{nn}\rho ^2\overline{\rho }^1\delta (rr(t))\delta (\theta \theta (t))`$ (89)
$`+`$ $`{\displaystyle \frac{\mathrm{\Delta }^2\overline{\rho }^2}{\sqrt{2}\rho }}\left(L_2^{}S_\mathrm{}m^{a\omega }+ia(\overline{\rho }\rho )\mathrm{sin}\theta S_\mathrm{}m^{a\omega }\right)J_+\left\{C_{\overline{m}n}\rho ^2\overline{\rho }^2\mathrm{\Delta }^1\delta (rr(t))\delta (\theta \theta (t))\right\}`$ (90)
$`+`$ $`{\displaystyle \frac{1}{2\sqrt{2}}}L_2^{}\left\{\rho ^3S_\mathrm{}m^{a\omega }(\overline{\rho }^2\rho ^4)_{,r}\right\}C_{\overline{m}n}\mathrm{\Delta }\rho ^2\overline{\rho }^2\delta (rr(t))\delta (\theta \theta (t))`$ (91)
$``$ $`{\displaystyle \frac{1}{4}}\rho ^3\mathrm{\Delta }^2S_\mathrm{}m^{a\omega }J_+\left\{\rho ^4J_+\left(\overline{\rho }\rho ^2C_{\overline{m}\overline{m}}\delta (rr(t))\delta (\theta \theta (t))\right)\right\}],`$ (92)
where

$$L_s^{}=_\theta \frac{m}{\mathrm{sin}\theta }+a\omega \mathrm{sin}\theta +s\mathrm{cot}\theta .$$
(94)
The $`\theta `$-integration can be performed directly to give

$`T_{\mathrm{}m\omega }`$ $`=`$ $`\mu {\displaystyle _{\mathrm{}}^{\mathrm{}}}dte^{i\omega tim\phi (t)}\mathrm{\Delta }^2[(A_{nn0}+A_{\overline{m}n0}+A_{\overline{m}\overline{m}0})\delta (rr(t))`$ (95)
$`+`$ $`\left\{(A_{\overline{m}n1}+A_{\overline{m}\overline{m}1})\delta (rr(t))\right\}_{,r}+\left\{A_{\overline{m}\overline{m}2}\delta (rr(t))\right\}_{,rr}]_{\theta =\theta (t)},`$ (96)
where

$`A_{nn0}`$ $`=`$ $`{\displaystyle \frac{2}{\sqrt{2\pi }\mathrm{\Delta }^2}}C_{nn}\rho ^2\overline{\rho }^1L_1^+\{\rho ^4L_2^+(\rho ^3S_\mathrm{}m^{a\omega })\},`$ (97)
$`A_{\overline{m}n0}`$ $`=`$ $`{\displaystyle \frac{2}{\sqrt{\pi }\mathrm{\Delta }}}C_{\overline{m}n}\rho ^3[\left(L_2^+S_\mathrm{}m^{a\omega }\right)({\displaystyle \frac{iK}{\mathrm{\Delta }}}+\rho +\overline{\rho })`$ (101)
$`a\mathrm{sin}\theta (t)S_\mathrm{}m^{a\omega }{\displaystyle \frac{K}{\mathrm{\Delta }}}(\overline{\rho }\rho )],`$
$`A_{\overline{m}\overline{m}0}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}\rho ^3\overline{\rho }C_{\overline{m}\overline{m}}S\left[i\left({\displaystyle \frac{K}{\mathrm{\Delta }}}\right)_{,r}{\displaystyle \frac{K^2}{\mathrm{\Delta }^2}}+2i\rho {\displaystyle \frac{K}{\mathrm{\Delta }}}\right],`$ (102)
$`A_{\overline{m}n\mathrm{\hspace{0.17em}1}}`$ $`=`$ $`{\displaystyle \frac{2}{\sqrt{\pi }\mathrm{\Delta }}}\rho ^3C_{\overline{m}n}[L_2^+S_\mathrm{}m^{e\omega }+ia\mathrm{sin}\theta (t)(\overline{\rho }\rho )S_\mathrm{}m^{a\omega }],`$ (103)
$`A_{\overline{m}\overline{m}1}`$ $`=`$ $`{\displaystyle \frac{2}{\sqrt{2\pi }}}\rho ^3\overline{\rho }C_{\overline{m}\overline{m}}S_\mathrm{}m^{a\omega }\left(i{\displaystyle \frac{K}{\mathrm{\Delta }}}+\rho \right),`$ (104)
$`A_{\overline{m}\overline{m}2}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}\rho ^3\overline{\rho }C_{\overline{m}\overline{m}}S_\mathrm{}m^{a\omega }.`$ (105)
Note that all functions of $`\theta `$ are evaluated at $`\theta =\theta (t)`$. The amplitudes $`Z_{\mathrm{}m\omega }^{\mathrm{},H}`$ defined in (63) can be written as

$`Z_{\mathrm{}m\omega }^H`$ $`=`$ $`{\displaystyle \frac{\mu }{2i\omega B^{\mathrm{in}}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑te^{i\omega tim\phi (t)}_{\mathrm{}m\omega }^H(r(t),\theta (t)),`$ (106)
$`Z_{\mathrm{}m\omega }^{\mathrm{}}`$ $`=`$ $`{\displaystyle \frac{\mu }{2i\omega B^{\mathrm{in}}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑te^{i\omega tim\phi (t)}_{\mathrm{}m\omega }^{\mathrm{}}(r(t),\theta (t)),`$ (109)
where

$`_{\mathrm{}m\omega }^\mathrm{H}`$ $`=`$ $`[R_{\mathrm{}m\omega }^{\mathrm{in}}\{A_{nn0}+A_{\overline{m}n0}+A_{\overline{m}\overline{m}\mathrm{\hspace{0.17em}0}}\}`$ (111)
$`{\displaystyle \frac{dR_{\mathrm{}m\omega }^{\mathrm{in}}}{dr}}\{A_{\overline{m}n1}+A_{\overline{m}\overline{m}1}\}+{\displaystyle \frac{d^2R_{\mathrm{}m\omega }^{\mathrm{in}}}{dr^2}}A_{\overline{m}\overline{m}2}]_{r=r(t),\theta =\theta (t)},`$
$`_{\mathrm{}m\omega }^{\mathrm{}}`$ $`=`$ $`[R_{\mathrm{}m\omega }^{\mathrm{up}}\{A_{nn0}+A_{\overline{m}n0}+A_{\overline{m}\overline{m}0}\}`$ (115)
$`{\displaystyle \frac{dR_{\mathrm{}m\omega }^{\mathrm{up}}}{dr}}\{A_{\overline{m}n1}+A_{\overline{m}\overline{m}1}\}+{\displaystyle \frac{d^2R_{\mathrm{}m\omega }^{\mathrm{up}}}{dr^2}}A_{\overline{m}\overline{m}2}]_{r=r(t),\theta =\theta (t)}.`$
Up to this point, all expressions listed in this Section are valid for an arbitrary orbit. We now specialise our discussion to equatorial orbits by setting $`\theta (t)=\pi /2`$. In this case, $`_{\mathrm{}m\omega }^{\mathrm{},H}`$ are functions of $`r(t)`$ only. As discussed in detail in Appendix B, the quantities

$$\alpha ^{\mathrm{},H}(t)=^{\mathrm{},H}(r(t))e^{im[\varphi (t)\mathrm{\Omega }_\varphi t]},$$
(116)
are periodic functions of time (with period equal to $`T_r`$). Consequently, they can be expanded in a Fourier series

$$\alpha ^{\mathrm{},H}(t)=\underset{k=\mathrm{}}{\overset{+\mathrm{}}{}}\alpha _k^{\mathrm{},H}e^{ik\mathrm{\Omega }_rt},$$
(117)
with $`\mathrm{\Omega }_\mathrm{r}=2\pi /T_\mathrm{r}`$. Inverting, we obtain the Fourier coefficients

$$\alpha _k^{\mathrm{},H}=\frac{1}{T_r}_0^{T_r}𝑑t\alpha ^{\mathrm{},H}(t)e^{ik\mathrm{\Omega }_rt}.$$
(118)
Using the Fourier series (117) in (106),(109) we arrive at

$$Z_{\mathrm{}m\omega }^{\mathrm{},H}=\underset{k=\mathrm{}}{\overset{+\mathrm{}}{}}Z_{\mathrm{}mk}^{\mathrm{},H}\delta (\omega \omega _{mk}),$$
(119)
where $`\omega _{mk}=m\mathrm{\Omega }_\varphi +k\mathrm{\Omega }_r`$ and

$$Z_{\mathrm{}mk}^{\mathrm{},H}=\frac{\mu \mathrm{\Omega }_r}{2i\omega _{mk}B^{in}}_0^{T_r}𝑑t^{\mathrm{},H}(r(t))e^{i\omega _{mk}tim\varphi (t)}.$$
(120)
Due to symmetries of the Teukolsky equation (52) we have the following conjugation relation,

$`Z_{\mathrm{},m,\omega }^{\mathrm{},H}`$ $`=`$ $`(1)^{\mathrm{}}\overline{Z}_{\mathrm{},m,\omega }^{\mathrm{},H},`$ (121)
$`Z_{\mathrm{},m,k}^{\mathrm{},H}`$ $`=`$ $`(1)^{\mathrm{}}\overline{Z}_{\mathrm{},m,k}^{\mathrm{},H},`$ (122)
where an overbar denotes complex conjugate.

We proceed by writing (120) as an integral over $`\chi `$,

$$Z_{\mathrm{}mk}^{\mathrm{},H}=\frac{\mu \mathrm{\Omega }_r}{2i\omega _{mk}B^{in}}_0^{2\pi }𝑑\chi \frac{\stackrel{~}{V}_\mathrm{t}(\chi )}{J(\chi )\stackrel{~}{V}_\mathrm{r}^{1/2}(\chi )}_{\mathrm{}m\omega }^{\mathrm{},H}(r(\chi ))e^{i\omega _{mk}t(\chi )im\varphi (\chi )}.$$
(123)
As in the case of eqns. (19), (20), this expression is well-behaved at the orbital turning points. Moreover, noting that the $`\chi `$-dependence of the integrand in (123) appears in terms of the form $`\mathrm{cos}\chi `$ (in terms with $`r(\chi )`$)) and $`\mathrm{sin}\chi `$ (in terms with $`u^\mathrm{r}`$) we can write,

$`Z_{\mathrm{}mk}^{\mathrm{},H}`$ $`=`$ $`{\displaystyle \frac{\mu \mathrm{\Omega }_r}{2i\omega _{mk}B^{in}}}{\displaystyle _0^\pi }d\chi {\displaystyle \frac{\stackrel{~}{V}_\mathrm{t}(\chi )}{J(\chi )\stackrel{~}{V}_\mathrm{r}^{1/2}(\chi )}}[_{\mathrm{}m\omega (+)}^{\mathrm{},H}(r(\chi ))e^{i\omega _{mk}t(\chi )im\varphi (\chi )}`$ (124)
$`+`$ $`_{\mathrm{}m\omega ()}^{\mathrm{},H}(r(\chi ))e^{i\omega _{mk}t(\chi )+im\varphi (\chi )}].`$ (125)
The subscripts $`(\pm )`$ mean “$`\mathrm{sin}\chi \pm \mathrm{sin}\chi `$” in the functions $`_{\mathrm{}m\omega }^{\mathrm{},H}`$. The numerical calculation of the amplitudes $`Z_{\mathrm{}mk}^{\mathrm{},H}`$ is the “backbone” of our radiation reaction code (see Section V). Finally, we express the $`A`$ and $`C`$ amplitudes (87),(105), in terms of $`\chi ,p,e`$,

$`C_{nn}(\chi ,p,e)`$ $`=`$ $`{\displaystyle \frac{J(\chi ,p,e)}{4p^4\stackrel{~}{V}_\mathrm{t}(\chi ,p,e)}}\left[p^2Eax(1+e\mathrm{cos}\chi )^2+ep\mathrm{sin}\chi \stackrel{~}{V}_\mathrm{r}^{1/2}(\chi ,p,e)\right]^2,`$ (126)
$`C_{\overline{m}n}(\chi ,p,e)`$ $`=`$ $`{\displaystyle \frac{ixJ(\chi ,p,e)}{2\sqrt{2}p^3\stackrel{~}{V}_\mathrm{t}(\chi ,p,e)}}(1+e\mathrm{cos}\chi )\left[p^2Eax(1+e\mathrm{cos}\chi )^2+ep\mathrm{sin}\chi \stackrel{~}{V}_\mathrm{r}^{1/2}(\chi ,p,e)\right],`$ (128)
$`C_{\overline{m}\overline{m}}(\chi ,p,e)`$ $`=`$ $`{\displaystyle \frac{x^2J(\chi ,p,e)}{2p^2\stackrel{~}{V}_\mathrm{t}(\chi ,p,e)}}(1+e\mathrm{cos}\chi )^2,`$ (130)
and

$`A_{\overline{m}n0}(u)`$ $`=`$ $`{\displaystyle \frac{2}{\sqrt{\pi }}}{\displaystyle \frac{C_{\overline{m}n}}{u(12Mu+a^2u^2)^2}}\left[2a^2u^3+\{ia(a\omega m)4M\}u^2+2u+i\omega \right]`$ (132)
$`\times \left[{\displaystyle \frac{S_\mathrm{}m^{a\omega }}{\theta }}(\pi /2)+(a\omega m)S_\mathrm{}m^{a\omega }(\pi /2)\right],`$
$`A_{\overline{m}\overline{m}0}(u)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \frac{C_{\overline{m}\overline{m}}S_\mathrm{}m^{a\omega }(\pi /2)}{u^2(12Mu+a^2u^2)^2}}[2ia^3(a\omega m)u^5+a(a\omega m)\{6iM+a(a\omega m)\}u^4`$ (136)
$`4ia(a\omega m)u^3+2\omega \{iM+a(a\omega m)\}u^22i\omega u+\omega ^2],`$
$`A_{\overline{m}n1}(u)`$ $`=`$ $`{\displaystyle \frac{2}{\sqrt{\pi }}}{\displaystyle \frac{C_{\overline{m}n}}{u(12Mu+a^2u^2)}}\left[{\displaystyle \frac{S_\mathrm{}m^{a\omega }}{\theta }}(\pi /2)+(a\omega m)S_\mathrm{}m^{a\omega }(\pi /2)\right],`$ (139)
$`A_{\overline{m}\overline{m}1}(u)`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{C_{\overline{m}\overline{m}}S_\mathrm{}m^{a\omega }(\pi /2)}{u^2(12Mu+a^2u^2)}}\left[a^2u^3+\{ia(a\omega m)2M\}u^2+u+i\omega \right],`$ (142)
$`A_{\overline{m}\overline{m}2}(u)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \frac{C_{\overline{m}\overline{m}}S_\mathrm{}m^{a\omega }(\pi /2)}{u^2}},`$ (145)
$`A_{nn0}(u)`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{C_{nn}}{(12Mu+a^2u^2)^2}}[2ia({\displaystyle \frac{S_\mathrm{}m^{a\omega }}{\theta }}(\pi /2)+(a\omega m)S_\mathrm{}m^{a\omega }(\pi /2))u`$ (150)
$`+{\displaystyle \frac{^2S_\mathrm{}m^{a\omega }}{\theta ^2}}(\pi /2)+2(a\omega m){\displaystyle \frac{S_\mathrm{}m^{a\omega }}{\theta }}(\pi /2)+\{(a\omega m)^22\}S_\mathrm{}m^{a\omega }(\pi /2)],`$
where $`u(\chi ,p,e)=(1+e\mathrm{cos}\chi )/p`$. By means of (119) one can obtain the following expressions for $`\psi _4`$ at infinity and on the horizon,

$$\psi _4(t,r,\theta ,\varphi )\{\begin{array}{cc}\rho (r_+)^4_{\mathrm{}mk}\psi _{\mathrm{}mk}^H\hfill & \text{for }rr_+\hfill \\ & \\ r^1_{\mathrm{}mk}\psi _{\mathrm{}mk}^{\mathrm{}}\hfill & \text{for }r+\mathrm{},\hfill \end{array}$$
(151)
where

$$\psi _{\mathrm{}mk}^{H,\mathrm{}}=\frac{1}{\sqrt{2\pi }}Z_{\mathrm{}mk}^{\mathrm{},H}S_\mathrm{}m^{a\omega _{mk}}(\theta )e^{i\omega _{mk}(tr_{})+im\varphi }.$$
(152)
Once the Weyl scalar $`\psi _4`$ is known, we can immediately relate it to the two polarisation components $`h_+,h_\mathrm{x}`$ of the transverse-traceless metric perturbation as $`r\mathrm{}`$ ,

$$\psi _4\frac{1}{2}\left(\frac{^2h_+}{t^2}i\frac{^2h_\mathrm{x}}{t^2}\right).$$
(153)
It follows from (153), (152) that $`h_{+,\mathrm{x}}(t,r,\theta ,\varphi )`$ are given by (here the coordinates $`t,r,\theta ,\varphi `$ are referred to the observation point),

$$h_+ih_\mathrm{x}=\frac{2}{r}\underset{\mathrm{}mk}{}\frac{Z_{\mathrm{}mk}^\mathrm{H}}{\omega _{mk}^2}\frac{S_\mathrm{}m^{a\omega _{mk}}(\theta )}{\sqrt{2\pi }}e^{i\omega _{mk}(tr_{})+im\varphi }.$$
(154)
Note that the gravitational waveform is exclusively radiated at harmonics of the two orbital frequencies $`\mathrm{\Omega }_r`$, $`\mathrm{\Omega }_\varphi `$. The gravitational wave energy and angular momentum flux at infinity can be found in terms of the Landau-Lifschitz pseudotensor ,

$`\left({\displaystyle \frac{dE}{dt}}\right)_{GW}^{\mathrm{}}`$ $`=`$ $`{\displaystyle \frac{1}{16\pi }}{\displaystyle \left\{\left(\frac{h_+}{t}\right)^2+\left(\frac{h_\mathrm{x}}{t}\right)^2\right\}r^2𝑑\mathrm{\Omega }},`$ (155)
$`\left({\displaystyle \frac{dL}{dt}}\right)_{GW}^{\mathrm{}}`$ $`=`$ $`{\displaystyle \frac{1}{16\pi }}{\displaystyle \left\{\frac{h_+}{t}\frac{h_+}{\varphi }+\frac{h_\mathrm{x}}{t}\frac{h_\mathrm{x}}{\varphi }\right\}r^2𝑑\mathrm{\Omega }}.`$ (157)
We define the averaged (over one orbital period) fluxes to be $`(C=E,L)`$,

$$\dot{C}_{GW}\frac{1}{T_r}_0^{T_r}𝑑t\left(\frac{dC}{dt}\right)_{GW}^{\mathrm{}}.$$
(158)
With the help of (155), (157), (154) we arrive at

$`\dot{E}_{GW}^{\mathrm{}}`$ $`=`$ $`{\displaystyle \underset{\mathrm{},m,k}{}}{\displaystyle \frac{|Z_{\mathrm{}mk}^\mathrm{H}|^2}{4\pi \omega _{mk}^2}},`$ (159)
$`\dot{L}_{GW}^{\mathrm{}}`$ $`=`$ $`{\displaystyle \underset{\mathrm{},m,k}{}}{\displaystyle \frac{m|Z_{\mathrm{}mk}^\mathrm{H}|^2}{4\pi \omega _{mk}^3}}.`$ (160)
The calculation of the respective fluxes at the black hole horizon is a more complicated issue as it is not possible to use expressions such as (155), (157). Despite this difficulty, Teukolsky and Press  were able to derive formulae for the horizon fluxes using the approach of Hawking and Hartle  who studied the deformation of the hole’s event horizon under the influence of infalling radiation. These formulae are,

$`\dot{E}_{GW}^H`$ $`=`$ $`{\displaystyle \underset{\mathrm{},m,k}{}}\alpha _{\mathrm{}mk}{\displaystyle \frac{|Z_{\mathrm{}mk}^{\mathrm{}}|^2}{4\pi \omega _{mk}^2}},`$ (161)
$`\dot{L}_{GW}^H`$ $`=`$ $`{\displaystyle \underset{\mathrm{},m,k}{}}\alpha _{\mathrm{}mk}{\displaystyle \frac{m|Z_{\mathrm{}mk}^{\mathrm{}}|^2}{4\pi \omega _{mk}^3}},`$ (162)
where

$$\alpha _{\mathrm{}mk}=\frac{256(2Mr_+)^5p_{mk}(p_{mk}^2+4ϵ^2)(p_{mk}^2+16ϵ^2)\omega _{mk}^3}{C_{\mathrm{}mk}},$$
(163)
with $`ϵ=\sqrt{M^2a^2}/4Mr_+`$ and

$`C_{\mathrm{}mk}`$ $`=`$ $`[(\lambda +2)^2+4am\omega _{mk}4a^2\omega _{mk}^2](\lambda ^2+36am\omega _{mk}36a^2\omega _{mk}^2)`$ (165)
$`+(2\lambda +3)(96a^2\omega _{mk}^248am\omega _{mk})+144\omega _{mk}^2(M^2a^2),`$
is the so-called Starobinsky constant. Note that eqns. (159)-(162) have to be divided by $`\mu `$ in order to convert them to fluxes of specific energy and angular momentum. Moreover, we can exploit the conjugation relations (122) in the numerical calculation of the amplitudes $`Z_{\mathrm{}mk}`$ and reduce by one half the required computational time.

### B The Sasaki-Nakamura equation

From eqn. (125) it is obvious that in order to calculate the amplitudes $`Z_{\mathrm{}mk}^{\mathrm{},H}`$, which will give us the gravitational waveform and fluxes (154), (159)- (162), we need to evaluate the quantity $`B^{\mathrm{in}}`$. In principle, one could numerically integrate the Teukolsky equation (52) from the horizon out to “infinity” and extract the amplitudes $`B^{\mathrm{in},\mathrm{out}}`$. But this is a poor strategy, since the effective potential $`V(r)`$ is long-ranged and the $`B^{\mathrm{in}}`$ term drops off towards infinity much faster than the $`B^{\mathrm{out}}`$ term and can only be extracted with very low accuracy . A way to circumvent this difficulty is to integrate, instead, the Sasaki-Nakamura equation ,

$$\frac{d^2X}{dr_{}^2}F(r)\frac{dX}{dr_{}}U(r)X=0.$$
(166)
The “potentials” $`F(r),U(r)`$ are given in Appendix C. The solutions of this equation are related to the solutions of the Teukolsky equation via the transformation,

$$R_{\mathrm{}m\omega }(r)=\frac{1}{\eta }\left[\left(\alpha +\frac{\beta _{,r}}{\mathrm{\Delta }}\right)\frac{\mathrm{\Delta }X_{\mathrm{}m\omega }}{(r^2+a^2)^{1/2}}\frac{\beta }{\mathrm{\Delta }}\frac{d}{dr}\left(\frac{\mathrm{\Delta }X_{\mathrm{}m\omega }}{(r^2+a^2)^{1/2}}\right)\right].$$
(167)
The functions $`\alpha (r),\beta (r)`$ are also given in Appendix C. The key property of (166) is that it encompasses a short-range potential. This can be demonstrated more easily if we shift to the function,

$$Y(r)=\eta ^{1/2}(r)X(r).$$
(168)
Then, eqn. (166) transforms into the Schrödinger-type equation,

$$\frac{d^2Y}{dr_{}^2}+QY=0,$$
(169)
with the effective potential,

$$Q=U\frac{1}{4}F^2+\frac{\mathrm{\Delta }}{2\eta (r^2+a^2)^2}\left\{\mathrm{\Delta }\eta _{,rr}\frac{\mathrm{\Delta }}{\eta }(\eta _{,r})^2+2M\eta _{,r}\left(\frac{r^2a^2}{r^2+a^2}\right)\right\}.$$
(170)
The functions $`F(r)`$, $`U(r)`$ have the following behaviour at infinity and at the horizon,

$$F(r)\{\begin{array}{cc}0+𝒪(rr_+)\hfill & \text{for }rr_+\hfill \\ r^2c_1/c_0+𝒪(r^3)\hfill & \text{for }r+\mathrm{},\hfill \end{array}$$
(171)
$$U(r)\{\begin{array}{cc}k^2+𝒪(rr_+)\hfill & \text{for }rr_+\hfill \\ \omega ^2+r^2\left[\lambda +2(1+am\omega a^2\omega ^2)i\omega c_1/c_0\right]+𝒪(r^3),\hfill & \text{for }r+\mathrm{}.\hfill \end{array}$$
(172)
It follows that,

$$Q(r_{})\{\begin{array}{cc}\omega ^2r_{}^2[\lambda +2(1+am\omega a^2\omega ^2i\omega c_1/c_0]+𝒪(r_{}^3\mathrm{ln}r_{})\hfill & \text{for }r_{}+\mathrm{}\hfill \\ k^2+𝒪(e^{cr_{}})\hfill & \text{for }r_{}\mathrm{},\hfill \end{array}$$
(173)
where $`c=(r_+r_{})/2M`$ is a positive constant. From this expression is obvious that $`Q`$ is short-ranged. Consequently, equation (166) admits a solution (“in” mode) which is purely ingoing at the horizon and a mixture of ingoing/outgoing waves at infinity:

$$X^{\mathrm{in}}\{\begin{array}{cc}A^{\mathrm{down}}e^{ikr_{}}\hfill & \text{for }rr_+\hfill \\ A^{\mathrm{in}}e^{i\omega r_{}}+A^{\mathrm{out}}e^{i\omega r_{}}\hfill & \text{for }r+\mathrm{},\hfill \end{array}$$
(174)
Another useful independent solution to (166) is the “up” mode,

$$X^{\mathrm{up}}\{\begin{array}{cc}D^{\mathrm{in}}e^{ikr_{}}+D^{\mathrm{out}}e^{ikr_{}}\hfill & \text{for }rr_+\hfill \\ D^{\mathrm{up}}e^{i\omega r_{}}\hfill & \text{for }r+\mathrm{}.\hfill \end{array}$$
(175)
The relation between the asymptotic amplitudes appearing in (58) and (174) can be deduced from (167),

$`B^{\mathrm{in}}`$ $`=`$ $`{\displaystyle \frac{1}{4\omega ^2}}A^{\mathrm{in}},`$ (176)
$`B^{\mathrm{out}}`$ $`=`$ $`{\displaystyle \frac{4\omega ^2}{c_0}}A^{\mathrm{out}},`$ (177)
where the constant $`c_0`$ is given in Appendix C. Hence, we can simply integrate equation (166) instead of (52) and easily identify the ingoing and outgoing waves and evaluate their respective amplitudes. We can then simply find the desired amplitudes $`B^{\mathrm{in}/\mathrm{out}}`$ from (177). Similarly, knowledge of the wavefunction $`X(r)`$ and its derivative at a given point immediately leads to the Teukolsky radial function $`R(r)`$ and its derivative via the rule (167). In conclusion, all the quantities (apart from the spheroidal harmonics) required for the calculation of the gravitational flux and waveform, can be obtained by numerical integration of the Sasaki-Nakamura equation (166).

### C Orbital evolution: adiabaticity and flux balance

Due to the emission of gravitational radiation the orbit of a particle around a black hole will slowly evolve in time and the orbital constants $`E,L`$ (or equivalently $`p,e`$) will no longer be conserved. Radiation reaction effects become noticeable on a timescale that scales as $`M^2/\mu `$, i.e. they are always tiny in a timescale $`𝒪(M)`$, provided the system’s mass ratio is sufficiently small. We can define as the radiation reaction timescale,

$$T_{\mathrm{RR}}=\text{min}[T_p,T_e],$$
(178)
where $`T_e=e/|\dot{e}|`$ and $`T_p=p/|\dot{p}|`$ are the radiative timescales for $`p`$ and $`e`$ respectively (approximate expressions for these timescales are given in the following Section). We will then say that an orbit evolves adiabatically if

$$T_{\mathrm{RR}}T_\mathrm{r}.$$
(179)
In other words, it is a good approximation to assume the motion of the particle to be geodesic, as long as we are interested in timescales much shorter than $`T_{\mathrm{RR}}`$. On the other hand, by making this simplification we “freeze” the evolution of the orbit, as if there was no radiation reaction. That is, within the adiabatic approximation, we cannot know the exact evolution of the functions $`E(t),L(t)`$ (or of $`p(t),e(t)`$). It is still possible, however, to calculate an averaged rate of change of such quantities. This can be done by assuming the following “flux-balance” relation

$$\dot{C}=\dot{C}_{\mathrm{GW}}=(\dot{C}_{\mathrm{GW}}^{\mathrm{}}+\dot{C}_{\mathrm{GW}}^H),$$
(180)
where $`C=E,L`$. We have separately denoted the gravitational wave fluxes at infinity and down to the horizon by $`\dot{C}_{\mathrm{GW}}^{\mathrm{}},\dot{C}_{\mathrm{GW}}^H`$ respectively. The overdot symbol stands for the averaged (over one orbital period) rate, see Section IIIB. We can equally well describe an orbit by means of the parameters $`(p,e)`$ and calculate the relevant averaged rates of change of those quantities. Since $`E=E(p,e)`$ and $`L=L(p,e)`$ we have that (commas denote partial derivatives),

$`\dot{E}`$ $`=`$ $`E_{,p}\dot{p}+E_{,e}\dot{e},`$ (181)
$`\dot{L}`$ $`=`$ $`L_{,p}\dot{p}+L_{,e}\dot{e}.`$ (182)
These relations can be inverted to obtain,

$`\dot{p}`$ $`=`$ $`H^1(E_{,e}\dot{L}+L_{,e}\dot{E}),`$ (183)
$`\dot{e}`$ $`=`$ $`H^1(E_{,p}\dot{L}L_{,p}\dot{E}),`$ (184)
with $`H=E_{,p}L_{,e}E_{,e}L_{,p}`$. Eventually, all partial derivatives of $`E`$ and $`L`$ can be found in terms of the corresponding partial derivatives of the functions $`F,N`$ and $`\mathrm{\Delta }_\mathrm{x}`$ which are given explicitly in Appendix A. However, the resulting formulae are quite messy so we do not present them here.

Note that although the formalism adopted in our analysis offers only a “local” information on the radiative orbital evolution, it can be further manipulated in order to obtain additional information. A recipe for “evolving” orbits under radiation reaction, using the known averaged rates of change of the relevant orbital constants, was given recently by Hughes  in the context of circular non-equatorial orbits. In effect, one is able to construct a series of “snapshots” of the radiation-induced inspiral, and make predictions of the evolution of the emitted waveform close to the point where the orbit becomes unstable (that is until the adiabatic condition no longer holds).

As we have already mentioned, adiabaticity will eventually break down near the separatrix, no matter how small the mass ratio is. This can be immediately seen from (179) and recalling that $`T_r\mathrm{}`$ at the separatrix, as predicted by (47). For an order-of-magnitude estimation, we can use the quadrupole approximation for the fluxes (see next Section) and translate (179) into a constraint on the mass ratio,

$$\frac{\mu }{M}\frac{5}{128\pi }\left(\frac{p}{M}\right)^{5/2}f_3^1(e)\left[1\pm \frac{a}{M}\left(\frac{M}{p}\right)^{3/2}f_3^1(e)\left(\frac{169}{12}+\frac{185}{12}e^2+\frac{223}{96}e^4\right)\right].$$
(185)
The function $`f_3(e)`$ is defined in the following Section. Equation (185) is accurate to leading order in $`M/p`$ and $`a/M`$, and to derive it we have used the corresponding order expression for the orbital period

$$T_r=2\pi M(1e^2)^{3/2}\left(\frac{p}{M}\right)^{3/2}\left[1\frac{3a}{M}\left(\frac{M}{p}\right)^{3/2}(1e^2)\right].$$
(186)
The mass-ratio constraint (185) is automatically satisfied as the black hole perturbation scheme we employ requires $`\mu /M1`$.

On the other hand, in the strong-field regime near the separatrix we find (using results derived in Section IVB),

$$\frac{\mu }{M}\delta H\left(\mathrm{ln}\left[\frac{64ep_s^2}{ϵS(1+e)(3e)}\right]\right)^1,$$
(187)
where $`\delta `$ is a combination of $`A_t(0),A_\varphi (0),E_{,p/e},L_{,p/e}`$ and is of order unity. As we discuss in Section IVB, the quantity $`H`$ also becomes zero when $`ϵ0`$ (unless $`a=M`$, in which case it remains finite). This is clearly the most severe restriction for the mass ratio. Fortunately, the real astrophysical systems we are trying to model are typically characterised by a mass ratio $`\mu /M10^6`$. Therefore we can approach the separatrix closely, probably to the point where the physical body would begin its plunge into the black hole, in the cases which interest us.

## IV Analytical results

### A Weak-field approximations

Orbits with $`pM`$ are well described by weak-field approximate results. In particular, the energy and angular momentum fluxes should be given with sufficient accuracy by the quadrupole-order formulae as given in ,. However, these authors make use of a different set of orbital parameters. For example, Ryan’s semi-major axis $`\overline{a}`$ and eccentricity $`\overline{e}`$  are related to our parameters via the transformation,

$`1e^2`$ $`=`$ $`(1\overline{e}^2)\left[1{\displaystyle \frac{4a}{M}}\left({\displaystyle \frac{M}{p}}\right)^{3/2}e^2\mathrm{cos}\iota \right],`$ (188)
$`p`$ $`=`$ $`\overline{a}(1e^2)\left[1{\displaystyle \frac{2a}{M}}\left({\displaystyle \frac{M}{p}}\right)^{3/2}e^2\mathrm{cos}\iota \right].`$ (189)
Note that at a Newtonian level the two sets are consistent to each other. Rewriting Ryan’s fluxes in terms of our parameters we obtain :

$`\dot{E}_{quad}`$ $`=`$ $`{\displaystyle \frac{32}{5}}{\displaystyle \frac{\mu ^2}{M^2}}\left({\displaystyle \frac{M}{p}}\right)^5(1e^2)^{3/2}\left[f_1(e){\displaystyle \frac{a}{M}}\left({\displaystyle \frac{M}{p}}\right)^{3/2}f_2(e)\right],`$ (190)
$`\dot{L}_{quad}`$ $`=`$ $`{\displaystyle \frac{32}{5}}{\displaystyle \frac{\mu ^2}{M}}\left({\displaystyle \frac{M}{p}}\right)^{7/2}(1e^2)^{3/2}\left[f_3(e)+{\displaystyle \frac{a}{M}}\left({\displaystyle \frac{M}{p}}\right)^{3/2}(f_4(e)f_5(e))\right],`$ (191)
where

$`f_1(e)`$ $`=`$ $`1+{\displaystyle \frac{73}{24}}e^2+{\displaystyle \frac{37}{96}}e^4,`$ (192)
$`f_2(e)`$ $`=`$ $`{\displaystyle \frac{73}{12}}+{\displaystyle \frac{823}{24}}e^2+{\displaystyle \frac{949}{32}}e^4+{\displaystyle \frac{491}{192}}e^6,`$ (193)
$`f_3(e)`$ $`=`$ $`1+{\displaystyle \frac{7}{8}}e^2`$ (194)
$`f_4(e)`$ $`=`$ $`{\displaystyle \frac{61}{24}}+{\displaystyle \frac{63}{8}}e^2+{\displaystyle \frac{95}{64}}e^4,`$ (195)
$`f_5(e)`$ $`=`$ $`{\displaystyle \frac{61}{8}}+{\displaystyle \frac{91}{4}}e^2+{\displaystyle \frac{461}{64}}e^4.`$ (196)
The formulae (190),(191) can be utilised for order of magnitude estimations even in the strong field regime though becoming increasingly inaccurate with decreasing $`p/M`$ (this has been verified by comparing them to the fully numerical results ). We can now estimate the timescales of radiative evolution for $`p`$,$`e`$. For $`pM`$ the energy and angular momentum, at leading order in $`M/p`$ and $`a/M`$, are given by,

$`E`$ $``$ $`1{\displaystyle \frac{M}{2p}}(1e^2){\displaystyle \frac{a}{M}}(1e^2)^2\left({\displaystyle \frac{M}{p}}\right)^{5/2},`$ (197)
$`L`$ $``$ $`\pm \sqrt{Mp}{\displaystyle \frac{aM}{p}}(3+e).`$ (198)
Accordingly, eqns. (183), (184) become,

$`\dot{p}`$ $`=`$ $`{\displaystyle \frac{64}{5}}{\displaystyle \frac{\mu }{M}}(1e^2)^{3/2}\left({\displaystyle \frac{M}{p}}\right)^3\left[f_3(e){\displaystyle \frac{a}{4M}}\left({\displaystyle \frac{M}{p}}\right)^{3/2}f_6(e)\right],`$ (199)
$`\dot{e}`$ $`=`$ $`{\displaystyle \frac{304}{15}}{\displaystyle \frac{\mu }{M^2}}e(1e^2)^{3/2}\left({\displaystyle \frac{M}{p}}\right)^4\left[f_7(e){\displaystyle \frac{a}{M}}\left({\displaystyle \frac{M}{p}}\right)^{3/2}f_8(e)\right],`$ (201)
where

$`f_6(e)`$ $`=`$ $`{\displaystyle \frac{133}{12}}+{\displaystyle \frac{379}{24}}e^2+{\displaystyle \frac{475}{96}}e^4,`$ (202)
$`f_7(e)`$ $`=`$ $`1+{\displaystyle \frac{121}{304}}e^2,`$ (203)
$`f_8(e)`$ $`=`$ $`{\displaystyle \frac{879}{76}}+{\displaystyle \frac{699}{76}}e^2+{\displaystyle \frac{1313}{608}}e^4.`$ (204)
The equations above demonstrate the well known fact  that, in the weak-field regime, eccentric orbits tend to circularise under radiation reaction (while they slowly shrink towards the central body). For the associated timescales one finds,

$$\frac{T_p}{T_e}=\frac{19}{12}\left(1+\frac{7}{8}e^2\right)^1\left(1+\frac{121}{304}e^2\right)\left[1\frac{a}{M}\left(\frac{M}{p}\right)^{3/2}f_3^1(e)\left(\frac{55}{114}+\frac{6431}{1824}e^2+\frac{9593}{1824}e^4+\frac{9191}{4864}e^6\right)\right].$$
(205)
According to this equation, in the weak-field regime the eccentricity decays faster than the size of the orbit . The leading order spin term furthermore implies that this behaviour is more pronounced for retrograde orbits.

### B Approximations near the separatrix (II)

The previous Section discussed results which are already familiar from the existing literature . We now present new results regarding strong-field orbits which reside near the separatrix.

The analysis of Section IIIA has shown that the gravitational wave spectrum will essentially contain harmonics of $`\mathrm{\Omega }_\mathrm{r},\mathrm{\Omega }_\varphi `$. We can use the approximate expressions (47),(50) to deduce that for orbits near the separatrix, i.e. $`pp_s(e)=ϵM`$ (and as long as $`e`$ is not close to zero and is not unity),

$`\mathrm{\Omega }_\mathrm{r}`$ $``$ $`{\displaystyle \frac{2\pi }{A_\mathrm{t}(0)}}\left[{\displaystyle \frac{eMp_s}{(1+e)(3e)}}\right]^{1/2}\left(\mathrm{ln}\left[{\displaystyle \frac{64ep_s^2}{ϵS(1+e)(3e)}}\right]\right)^1,`$ (206)
$`\mathrm{\Omega }_\varphi `$ $``$ $`{\displaystyle \frac{A_\varphi (0)}{A_\mathrm{t}(0)}}.`$ (209)
Hence, for $`ϵ0`$ we have $`\mathrm{\Omega }_r0`$ and in effect, the spectrum becomes almost “circular”:

$$\omega _{\mathrm{mk}}m\mathrm{\Omega }_\varphi .$$
(210)
Furthermore, by substitution in (159), (160) we get

$$\dot{E}_{GW}^{\mathrm{},H}\mathrm{\Omega }_\varphi \dot{L}_{GW}^{\mathrm{},H}.$$
(211)
We conclude that orbits near the separatrix radiate energy and angular momentum at rates so that the ratio $`\dot{E}_{GW}/\dot{L}_{\mathrm{GW}}`$ is almost equal to the respective ratio of a circular orbit with the same $`\mathrm{\Omega }_\varphi `$. The effective radius of this fiducial orbit is

$$r_{\mathrm{eff}}=M^{1/3}\left(\frac{A_t(0)}{A_\varphi (0)}a\right)^{2/3}.$$
(212)
For example, for the prograde orbit $`p=2.11M,e=0.7`$ we find $`r_{\mathrm{eff}}=1.88M>r_p=1.24M`$ while for the orbit $`p=2.35M,e=0.9`$ we find $`r_{\mathrm{eff}}=3.90M>r_p=1.24M`$ (for both cases we have taken $`a=0.99M`$ ).

Equation (211) suggests that particles in zoom-whirl orbits lose most of their energy and angular momentum while they revolve near the periastron, which is what we would intuitively expect.

We next discuss approximations for $`\dot{p}`$ and $`\dot{e}`$ near the separatrix. Unfortunately, the lack of a simple analytic expression for $`p_s(e)`$ makes such a task difficult, and the resulting formulae are quite cumbersome with little analytic transparency. Nevertheless, we can follow a much simpler path and still gain some significant insight. For $`pp_s`$ and using eq. (211), eqns. (183), (184) become,

$`\dot{p}`$ $``$ $`\left[H^1(L_{,e}\mathrm{\Omega }_\varphi E_{,e})\dot{L}\right]_{pp_s},`$ (213)
$`\dot{e}`$ $``$ $`\left[H^1(L_{,p}\mathrm{\Omega }_\varphi +E_{,p})\dot{L}\right]_{pp_s}.`$ (215)
By direct substitution of (26), it can be shown that the function $`H(p,e)`$ becomes exactly zero at the separatrix. On the other hand, and as long as $`aM`$, one can verify numerically that the numerators in (213), (215) remain finite near and at the separatrix. It follows that for non-extreme Kerr holes both $`\dot{p}`$ and $`\dot{e}`$ diverge at the location of the separatrix. This pathological behaviour signals the breakdown of the adiabaticity assumption upon which our method stands. A proper discussion of this transition regime should take into account the rapid radiative evolution of the orbit (which now varies in a timescale comparable to the orbital period).

Moving on, we divide (215) into (213) to get,

$$\frac{\dot{e}}{\dot{p}}\left[\frac{(L_{,p}\mathrm{\Omega }_\varphi +E_{,p})}{(L_{,e}\mathrm{\Omega }_\varphi E_{,e})}\right]_{pp_s}.$$
(216)
Exploring the numerical value of this quantity for numerous very near-separatrix orbits and black hole spins $`a<M`$, we have found it to be always negative and finite. This means that $`\dot{e}`$ and $`\dot{p}`$ have opposite signs near the separatrix. Since the latter is always negative (the orbit always shrinks) we conclude that very close to the separatrix $`\dot{e}>0`$, i.e. the orbit gains eccentricity, (as previously found, in less general cases, in , , , ). Since weak field orbits always lose eccentricity, there must be a critical curve $`p_{crit}(e)`$ on the $`(p,e)`$ plane at which $`\dot{e}=0`$. As eqn. (216) is formally accurate (within the constraints imposed by the adiabaticity condition) not only at the separatrix but also in its vicinity, we can actually study the behaviour of $`\dot{e}`$ in a thin zone near the separatrix. For a given small or moderate black hole spin, we find that the ratio (216) is again negative. However, for high eccentricities $`e1`$ we initially get a positive value which gradually passes from zero and becomes negative as $`pp_\mathrm{e}`$. With increasing $`a/M`$ we observe the same behaviour at even lower eccentricities, provided we are considering prograde orbits. The opposite behaviour is observed for retrograde orbits. For $`aM`$, (216) becomes negative only very close to the separatrix for all eccentricities. These results suggest that, at least for $`e1`$, the critical curve $`p_{crit}(e)`$ is located close to $`p_s(e)`$ (this has been shown to be true for $`a=0`$ ), and that (for prograde orbits) $`p_{crit}(e)p_s(e)0`$ as $`aM`$ (which resembles the situation for nearly circular equatorial Kerr orbits ). All of our (semi)analytic predictions are fully supported by the numerical results that are presented in Section VB.

Prograde orbits near the separatrix of an extreme Kerr hole are discussed separately, and fully analytically, in the following Section. Here we should emphasise once again that all the approximations presented in this Section are valid provided $`eϵ/M`$. This excludes nearly circular orbits, which have been explored in detail in .

We can now write approximate expressions for the timescales $`T_e,T_p`$ for an orbit close to the separatrix.

$$\frac{T_e}{T_p}\left[\frac{e}{p}\frac{|\mathrm{\Omega }_\varphi L_{,e}E_{,e}|}{|E_{,p}\mathrm{\Omega }_\varphi L_{,p}|}\right]_{pp_s}.$$
(217)
For example, for a $`a=0.99M`$ Kerr hole, we have $`T_e/T_p=0.81`$ for $`p=1.7M,e=0.3`$ while for $`p=2.11M,e=0.7`$ we get $`T_e/T_p=4.5`$ (for both orbits, eqn. (217) is a good approximation). This situation is typical for non-extreme holes. As we move along the separatrix keeping a fixed distance from it, the ratio $`T_e/T_p`$ tends to increase (and become larger than unity) with increasing eccentricity. In comparison, the corresponding weak-field timescales ratio, eqn. (205), is relatively insensitive to variations of eccentricity.

### C Horizon-skimming orbits

A particularly interesting class of prograde strong-field orbits are those that potentially “graze” the black hole horizon. These orbits can only exist provided the black hole is near extremally rotating, $`aM`$ (this can be deduced from Fig. 2). Circular, non-equatorial horizon-skimming orbits were first studied by Wilkins  and more recently by Hughes . Here, on the other hand, we discuss equatorial horizon-skimming orbits of arbitrary (but not equal to unity or close to zero) eccentricity around an extreme Kerr black hole.

As the separatrix for these orbits takes the very simple form $`p_s(e)=M(1+e)`$ we can duplicate the analysis of the previous section following a purely analytical path. Expanding (12) around $`p=p_s=M(1+e)`$ we find that

$$x^2=M^2\left(\frac{1+e}{3e}\right)+𝒪(ϵ).$$
(218)
We then get for the energy and angular momentum,

$`E`$ $`=`$ $`\sqrt{{\displaystyle \frac{1+e}{3e}}}+g(e){\displaystyle \frac{ϵ}{M}}+𝒪(ϵ^2),`$ (219)
$`L`$ $`=`$ $`2M\sqrt{{\displaystyle \frac{1+e}{3e}}}+f(e){\displaystyle \frac{ϵ}{M}}+𝒪(ϵ^2).`$ (221)
The explicit form of the functions $`g(e),f(e)`$ is not required for the following analysis. We use these equations (together with (211), noting that $`\mathrm{\Omega }_\varphi =1/2M`$ on the separatrix for the orbits under discussion) to obtain,

$`H`$ $``$ $`{\displaystyle \frac{2[2Mg(e)f(e)]}{(1+e)^{1/2}(3e)^{3/2}}},`$ (222)
$`E_{,e}\dot{L}+L_{,e}\dot{E}`$ $``$ $`M[2Mg(e)f(e)]\dot{E},`$ (224)
$`E_{,p}\dot{L}L_{,p}\dot{E}`$ $``$ $`[2Mg(e)f(e)]\dot{E}.`$ (226)
Here, unlike the non-extreme case, the function $`H`$ remains finite as $`ϵ0`$. The above formulae, as well as the following ones, have a fractional error $`𝒪(ϵ/eM)`$. Hence for horizon skimming orbits,

$`\dot{p}`$ $``$ $`{\displaystyle \frac{1}{2}}M(1+e)^{1/2}(3e)^{3/2}\dot{E},`$ (227)
$`\dot{e}`$ $``$ $`{\displaystyle \frac{1}{2}}(1+e)^{1/2}(3e)^{3/2}\dot{E}.`$ (228)
We see that both rates are finite all the way down to the separatrix, unlike the $`a<M`$ case. However, the adiabaticity condition (179) is still invalidated at $`p=p_s`$.

More interesting is the behaviour of the ratio of the rates (227), (228),

$$\frac{\dot{e}}{\dot{p}}=\frac{1}{M}+𝒪(ϵ/eM).$$
(229)
This is always positive, which means that the $`\dot{e}>0`$ region that exists for $`a<M`$ shrinks to zero for extreme Kerr black holes. In other words, the critical curve $`p_{crit}(e)`$ has the same value of the Boyer-Lindquist coordinate as the separatrix itself. This conclusion completes the discussion of the previous subsection. For the ratio of the respective timescales we get,

$$\frac{T_e}{T_p}=\frac{e}{1+e}+𝒪(ϵ/eM).$$
(230)
As (230) predicts, the eccentricity always evolves more rapidly than the semi-latus rectum, which again is contrary to the situation with weak-field orbits.

## V Numerical results

### A Method and error estimates

Numerical solution of the Teukolsky equation or related equations has been a minor industry for nearly thirty years, since the pioneering work of Press and Teukolsky , Detweiler  and Sasaki and Nakamura . Our method is based on a numerical algorithm outlined in  and employs subroutines found in . It involves the use of Bulirsch-Stoer integration to solve the Sasaki-Nakamura equation. Our code is a direct descendent of the codes used in  and , since we deal with both arbitrarily eccentric orbits and black holes of non-zero spin. Romberg integration is used both to integrate Eqns. (19) and (20) and to integrate Eqn. (125). To calculate the spheroidal harmonic functions $`S_\mathrm{}m^{a\omega }`$ we use the “spectral decomposition” method described in . The reliability of these methods in general is well known.

We were able to check our numerical results for circular equatorial orbits with those of , where our agreement was good to 5 or 6 significant digits, and with the codes on which this code was based , , for eccentric orbits in Schwarzschild and for nearly circular orbits in Kerr, and again our agreement was good to 4 to 6 significant digits. A similarly good agreement was achieved by comparison with the published results of Tanaka et.al.  for equatorial eccentric orbits in Schwarzschild and with those of Shibata  for circular equatorial orbits in Kerr. We were also able to compare our results with those given by Shibata  for equatorial eccentric orbits in Kerr. In this case, however, we did find some disagreement of about $`1\%`$. The cause of this disagreement is not apparent, while it stays roughly at the same level for moderate and high eccentricities. The disagreement does not seem to be due to the problems of maintaining accuracy with the long runtimes and large number of harmonics required for moderate/large eccentricities. The error introduced by the truncation of the $`\mathrm{},k`$ sums in the flux calculation does not seem to be the source of the disagreement. We cannot say at present which code might be at fault.

Finally, we have also been able to compare our code with some results for circular orbits in Kerr from , and again our agreement is good to several significant figures. Similarly, comparison with post Newtonian results for eccentric Kerr orbits (as found, for example, in  ) reveals good agreement in the weak-field regime. In Table III we compare some sample results.

In view of the lack of any check for strong field orbits with high/moderate eccentricities and high spins, it is of obvious importance that we present some estimate of the likely error in our numerical results. The main sources of numerical error in our code are as follows:

1. Inaccuracy in the Bulirsch-Stoer integration routines, from . We set the relative accuracy parameter $`eps`$, which governs the convergence of the final result, at $`10^6`$.

2. Inaccuracies in the Romberg integrator, also from . We set $`eps=10^6`$ for the routines which integrated Eqs. (20) and (19). This parameter governs the level of convergence which the routine demands in the final result, before it stops iterating. However, in the case of the Romberg routine which governed the main program loop, i.e. the integration of Eq. (125), we typically set $`eps=10^5`$ in many cases in order to achieve large savings in computing times. In a few runs designed to produce data for illustration of waveforms only (not numerical data on flux quantities), we used $`eps=10^4`$.

3. Our method requires that we calculate the quantity $`B^{\mathrm{in}}`$ (see Eqn. (125) above). To do this we integrated the Sasaki-Nakamura equation (166) out to $`r=100/\omega _{mk}`$ and then successively doubled the limit of integration, until our Richardson extrapolator told us we had achieved convergence to a relative accuracy of $`10^5`$.

4. Our method for calculating spheroidal harmonic functions $`S_\mathrm{}m^{a\omega }`$ involved writing them as an expansion of an infinite set of spherical harmonic functions. Fortunately this expansion can be truncated at 30 terms and remain very precise in most cases, but for high black hole spins, $`a`$ and high angular frequencies, $`\omega _{mk}`$ we were obliged to use 40 terms to avoid truncation errors causing small high frequency ripples in the wave forms. However, in our numerical results of flux rate and orbital evolution this source of error appears to be considerably less than $`10^6`$.

5. In principle our calculation of fluxes must be summed over an infinite number of harmonics in each of the integers $`\mathrm{}`$,$`m`$ and $`k`$. In practice truncating these sums for the $`\mathrm{}`$ and $`m`$ harmonics was not difficult. Fluxes for a sequence of these harmonics usually monotonically decrease after a peak at some value of $`\mathrm{}`$ and $`m`$ and so we demanded that the loop through these variables halt once fluxes went below a factor of $`10^5`$ times the peak contribution. However the spectrum of fluxes in the $`k`$ harmonic was much more complicated, typically involving several peaks (see figures in Section VD) before finally monotonically decreasing after a number of peaks which increased for increasingly eccentric orbits. We examined spectra in $`k`$ by hand to confirm the machine’s results and experimented widely to convince ourselves that we had caught all significant contributions to the total flux from different frequencies.

Clearly there are several significant independent sources of error, so that we can only offer our best judgement of the total relative error in our code in those cases where we have no independent check on our results. As we have addressed every systematic source of error that we encountered, and as we are confident that the code is running correctly in all of the cases dealt with in this paper, we estimate the relative error for numerical results quoted in this paper as no greater than $`10^310^4`$, in the case of fluxes, $`\dot{E}`$ and $`\dot{L}`$, and no greater than $`10^210^3`$ for quantities such as $`\dot{e}`$ and $`\dot{p}`$ because cancellations between terms when converting flux numbers into orbital evolution quantities tends to increase the size of the relative errors. This is especially true near the critical point where the rate of change of the eccentricity becomes zero, due to the complete cancellation of these terms. As an illustration of this, we will note in passing that the mysterious ”bump” seen in Fig. 2 of Ref.  turns out to be due to a rare case where flux errors which appear insignificant themselves are greatly magnified when the numerical flux data is combined to produce $`\dot{e}`$ and $`\dot{p}`$.

It is useful, in this context, to mention that in comparing the results from our code to the code in , the flux quantities for radiation emitted toward infinity agree to about 1 part in $`10^6`$, the flux quantities for radiation towards the horizon agree to about 1 part in $`10^4`$ and the position of the critical curve, as calculated by the two codes, can disagree by about 1%. This suggests that the only way in which our numerical error is large enough to make a visible difference in our figures would be as a slight change in the position of the critical curve in Fig. 5 (the retrograde case).

### B Backreaction on the orbit

In this Section we present numerical results on the evolution of bound equatorial orbits in terms of $`p`$ and $`e`$. This pair of parameters is preferable to the equivalent set $`E,L`$, because of their clearer geometrical meaning. We have calculated the averaged rates $`\dot{p}`$,$`\dot{e}`$ for a number of prograde and retrograde orbits and for two different black hole spins: $`a=0.5M`$, and $`a=0.99M`$. A representative part of our numerical results is presented in Fig. 5. Each individual orbit is represented as a point on the $`(p,e)`$ plane. At each one of those points we have attached a vector with components $`(M/\mu )(\dot{p},M\dot{e})`$ that indicates the direction at which the orbit adiabatically evolves under radiation reaction. Moreover, all orbits shown are chosen so as to be strongly adiabatic for the typical mass ratio $`10^6`$ (the most severe constraint is $`\mu /M10^2`$ ). These figures (together with the analytic approximations of Sections IVA, IVB and IVC) lead to the following conclusions :

* The semi-latus rectum $`p`$ always decreases (the orbit is shrinking). In fact, $`|\dot{p}|`$ grows monotonically, and finally diverges, as the separatrix is reached. This divergence is an artificial feature of our formalism, associated to the breakdown of the adiabatic approximation.
* The eccentricity $`e`$ shows a more complicated behaviour. For sufficiently large $`p`$, we always find $`\dot{e}<0`$. However, as the orbit approaches the separatrix, $`\dot{e}`$ changes sign and becomes positive, i.e near the separatrix, eccentricity increases. As in the case of $`\dot{p}`$, $`\dot{e}`$ will also diverge at the separatrix, due to the failure of adiabaticity.
* As the black hole spins up, in the case of prograde orbits the critical radius after which $`\dot{e}>0`$ moves closer to the separatrix (in coordinate terms) for a given $`e`$. The same is true for a fixed black hole spin, but for increasing $`e`$. For retrograde orbits, spinning up the black hole tends to move the $`\dot{e}=0`$ curve away from the separatrix.
* In a sense, the increasing eccentricity regime is a precursor of orbital instability and plunging. This is hinted by the proximity of the critical curve $`p_{crit}`$, where $`\dot{e}`$ flips sign and becomes positive, to the separatrix curve $`p_s(e)`$ which is the boundary between stable and unstable bound orbits. Qualitively speaking, at this stage of the inspiral, the radial potential $`V_r`$ is quite “flat” and as a consequence the particle has more room to move radially, even as it continues to “sink” towards the bottom of the potential well (the behaviour which is responsible for the characteristic “circularizing” tendency).

These results agree with, and at the same time generalise previous results concerning bound orbits (of arbitrary $`e`$) around Schwarzschild black holes  and slightly eccentric equatorial orbits around Kerr black holes . As we have discussed in Section VC, some of the conclusions above must be modified when the black hole is extreme ($`a=M`$).

In Table VII we give a sample of our numerical data, for the energy and angular momentum fluxes as well as for $`\dot{p},\dot{e}`$, for some of the orbits presented in Fig. 5. As we have discussed, we believe that these numbers have fractional accuracy at least $`10^3`$.

Another important result concerns the significance of the horizon fluxes on the evolution of orbits with relatively small periastrii. Specifically, we have encountered very-strong field orbits for which $`|\dot{E}_{GW}^\mathrm{H}|0.1|\dot{E}_{GW}^{\mathrm{}}|`$. However, the most intriguing property of the horizon fluxes is that they assist, in a sense, the orbiting body. This is most easily illustrated by plotting the evolution of the set of orbits of the top graph of Fig. 5, without including the horizon flux (represented by dashed arrows in Fig. 5; the solid arrows represent the total rates). For very strong field orbits, when the horizon fluxes are taken into account, the shrinking of the orbit is noticeably stalled. Typically, when $`\dot{E}_{GW}^\mathrm{H}`$ is non-negligible it also happens that it represents energy gain instead of energy loss (in other words the fluxes $`\dot{E}_{GW}^\mathrm{H},\dot{E}_{GW}^{\mathrm{}}`$ have opposite signs). The orbiting particle is effectively draining energy from the black hole itself. This is just a manifestation of the so-called superradiance phenomenon, well known in black hole physics ,: waves scattered within the black hole’s ergoregion and having frequencies (as measured at infinity) that lie in the interval $`0<\omega <m\omega _+`$, effectively appear (for a distant observer) as emerging from the horizon, and amplified at the expense of the hole’s rotational energy. The outgoing reflected waves “push” the particle outwards, and this interaction is manifested as a gain of orbital energy and angular momentum. Our result can be easily understood if we recall that the ergoregion is growing for increasing black hole spin. At the same time, because the boundary of instability moves in to lower radii with increasing spin, a particle can enter regions with much stronger fields and therefore emit a substantial amount of radiation towards the ergoregion (see Fig. 2). Hence we find a significant negative (superradiant) horizon flux. An alternative way of viewing this phenomenon is as an exchange of energy and angular momentum via tidal coupling analogous to tidal friction in the Earth-Moon system (and elsewhere). For an exposition of this intuitively instructive viewpoint see , and references therein.

It was recognised long ago  that if the superradiance effect ever became large enough a floating orbit would result when it balanced the energy loss due to radiation emitted towards infinity. Our results confirm earlier work in collaboration with Scott Hughes, suggesting that even for orbits very close to the horizon of very rapidly rotating black holes, the gain in energy from superradiance is only 10% of the energy lost by the system as a whole. For further details, see .

The results presented so far in this Section, although very insightful, are still incomplete in the sense that they do not describe the radiative evolution of a single given orbit on the ($`p,e`$) plane. Instead, they provide local information on the evolution of an orbit at a given point. An effort to “paste” together a sequence of such points, in order to follow the full orbit is currently underway . Meanwhile, we can use certain approximations to foresee what the full inspiral trajectory will look like. As a starting point we use the leading quadrupole-order expressions for $`\dot{p},\dot{e}`$, i.e. setting $`a=0`$ in eqns (199),(201) and derive,

$$p(e)=p_i\left(\frac{e}{e_i}\right)^{12/19}\left[\frac{1+\frac{121}{304}e^2}{1+\frac{121}{304}e_i^2}\right]^{870/2299}.$$
(231)
Given some initial values $`p_i,e_i`$ this relation describes, in the weak-field limit, the trajectory of the orbit on the ($`p,e`$) plane. Such curves for astrophysically relevant initial parameters, are shown in Fig. 6 (these curves remain essentially unchanged when the spin terms are retained in eqns. (199), (201) ). One feature that is immediately seen in the Newtonian-order inspiral is the absence of the critical curve $`\dot{e}=0`$ and the subsequent $`\dot{e}>0`$ behaviour. This should not come as a a surprise, as expression (231) is not formally valid unless $`pM`$. It is not safe to use weak-field approximations in strong field regimes.

A simple way to make better predictions is by using the exact expressions (183), (184) for $`\dot{p},\dot{e}`$, but still employing the weak-field formulae (190), (191) for the fluxes. The outcome of this trick is also shown in Fig. 6, for a set of orbits with initial parameters $`r_p=5,10,20M`$ and $`r_a=10^6M`$ (this translates to $`e=0.99999,0.99998,0.99996`$ and $`p=10,20,40M`$ respectively) for $`a=0,0.5M`$. We also considered a set of retrograde orbits with initial $`r_p=7,10,20M`$ and the same apastron as before, and $`a=0.99M`$ spin. Note that these curves, like the ones given by eqn. (231), are shape-invariant with respect to the mass ratio as long as $`\mu /M1`$. It is rewarding to see that these new trajectories do show the existence of the $`\dot{e}>0`$ region, and additionally are in good qualitive and quantitive agreement with the accurate numerical results . This is true as long as we do not attempt to evolve prograde orbits around rapidly spinning black holes, because the agreement quickly degrades when $`p`$ becomes small. Essentially this approach takes into account the correct form of the potential, which is the main cause behind the change in sign of $`\dot{e}`$, but the fact that the PN fluxes are increasingly inaccurate in strong-field regions precludes precise numerical agreement with the real trajectories.

The new curves clearly predict a higher residual eccentricity as compared to the pure Newtonian curves. This approximate result strongly suggests that many astrophysically relevant inspiralling orbits will have a significant amount of eccentricity left when they are close to plunging and will therefore be likely to exhibit zoom-whirl behaviour. Our full numerical results cannot at present be used to reproduce complete trajectories, but Fig. 5, which displays arrows which are tangential to these trajectories at individual points, certainly shows that if significant eccentricity remains at $`p5M`$, that this eccentricity will not disappear in the last part of the inspiral before plunge.

The work of Freitag  suggests that the initial periastra of scattered compact bodies which will eventually plunge into the black hole due to radiation reaction will be generally less than $`40M`$. However below that point the distribution of their periastra will be fairly flat, so that small initial periastra will be just as likely as large ones. As our figure shows, one expects, in the case of a Schwarzschild black hole, that bodies with initial $`r_p>20M`$ will have $`e<0.1`$ by the time of plunge. But for initial $`r_p<20M`$ the final eccentricity will be $`e>0.1`$ and can easily be as great as $`e0.7`$ (see Fig. 6) or higher. For instance, for $`r_p=10M`$, the final eccentricity will be $`e0.3`$ (see Fig. 6). We know that retrograde orbits will have less time to circularize and a longer ”de-circularizing” time, so eccentricities in this case should be greater. This is clearly seen in Fig. 6. In the case of prograde orbits we should generally expect smaller eccentricities before plunging (as compared to orbits around Schwarzschild black holes) but still at a significant level. Again Fig. 5 suggests that the change in eccentricity will not be great, despite the longer circularizing and shorter de-circularizing times. Moreover, near extreme holes allow a wider range of initial periastra (for prograde orbits), even down to $`r_p2M`$ and in these cases the residual eccentricity will be quite large. One concludes that eccentricity will play an important role in signal analysis for LISA.

The technique outlined here has been recently applied in constructing approximate, radiation reaction-driven, inspirals of test-bodies in generic Kerr orbits (for which case only weak-field results are currently available ). For more details see .

Motivated by the results discussed above, we would like to obtain simple estimates of the total amount of eccentricity gain and the number of orbits the particle will spend in the $`\dot{e}>0`$ region. Such estimates may provide a useful guideline for assessing the observational importance of this phase.

The numerical radiation reaction evolution-arrows in Fig. 5 show that the gradient $`de/dp`$ grows (in the negative sense) monotonically as soon as the critical curve is crossed. Hence, the maximum (negative) value is attained exactly at the separatrix. The approximate expression (216) for $`de/dp`$ is expected to be very accurate there. We can use this fixed gradient and extrapolate out to the critical curve, at some eccentricity $`e_i`$, for a given eccentricity $`e_f`$ at the separatrix. Then we have that $`\delta e=e_fe_i(p_s(e_f)p_c(e_f))(de/dp)_{p_s(e_f)}`$. This number should set an upper limit of the total increase in eccentricity. Results for some representative cases are given in Table V. From these numbers we deduce that, at best, there is a fractional increase of 5-50 % in eccentricity, the most favourable case being low-eccentricity retrograde orbits around rapidly spinning black holes. This gain decreases as we move upwards to larger final eccentricities (basically due to the shrinkage of the $`\dot{e}>0`$ region). We therefore conclude that we should not expect any dramatic increase in eccentricity when the orbit is about to become unstable.

A crude estimate on the number of orbits can be made by integrating (213) to find the time required to cross the $`\dot{e}>0`$ region,

$$t_c\frac{1}{\dot{L}}_{p_{crit}}^{p_s}\frac{Hdp}{(L_{,e}\mathrm{\Omega }_\varphi E_{,e})},$$
(232)
where we have factored out the angular momentum flux as it can be taken as constant within the integration interval (and recovered by our numerical data). Similarly we have assumed a fixed eccentricity. The number of orbits is then calculated by dividing $`t_c`$ with a typical period $`T_r`$ (or $`T_\varphi `$). We give some representative results in Table VI. We need to emphasise that these numbers should be viewed only as order-of-magnitude estimates, as eqn. (232) is a rough approximation. Nevertheless, we can still draw some reliable conclusions. For a small eccentricity our numbers are in agreement with existing results for nearly circular orbits around a Kerr black hole . For $`e=0.1`$ we should typically have a few thousands revolutions in the $`\dot{e}>0`$ region around a $`a=0.99M`$ hole and for a mass ratio $`\mu /M10^6`$. For the same parameters, but for retrograde orbits, there is an order of magnitude increase in the number of orbits. On the other hand, for all cases, there is an order of magnitude (or more) decrease as we move to eccentricity $`e=0.5`$. Note that it is possible to have a small number of full orbits, but yet a significant number of azimuthal revolutions or ”whirls” (see for example the $`a=0.99,e=0.5`$ case in Table VI).

### C Waveforms and fluxes from zoom-whirl orbits

Let us now focus on the class of orbits we have named zoom-whirl. As we have shown, orbits located near the separatrix should radiate in accordance with (211), as though they were nearly circular. In Table IV we list numerical fluxes for zoom-whirl orbits of various eccentricities. It is clear from these results that as the separatrix is approached, the analytic prediction (211) is indeed confirmed. However, one has to be very cautious when applying (211) to the study of a real astrophysical, extreme mass ratio, binary system. As our data reveal, in the region where this relation is fractionally accurate at the level of $`10^2`$, the adiabaticity constraint (179) on the mass ratio is quite severe, typically $`\mu /M10^210^3`$.

The zoom-whirl orbits are of interest for future detection efforts because of the characteristic waveform they generate. In Figs. 7 to 9 we show such waveforms (in particular we plot the quantity $`(\mu /r)h_+`$ as a function of retarded time $`tr_{}`$) for a couple of strong-field zoom-whirl orbits ( $`p=2.11M`$ and $`e=0.7`$, $`p=1.7M`$ and $`e=0.3`$) and for a rapidly spinning black hole with $`a=0.99M`$. Both orbits would evolve adiabatically if we consider a typical mass ratio $`\mu /M10^6`$. The corresponding gravitational wave flux data can be found in Table VII.

First we shall discuss the waveform as seen by an observer located on the black hole’s equatorial plane, see Fig. 7 and top panel of Fig. 9. Clearly, these waveforms have a very distinct appearance. A rapidly oscillating, high amplitude signal is radiated during the whirling of the particle near periastron. In between these bursts one observes low-amplitude signals produced during the particle’s zoom in and out from apastron. This contrast in amplitudes is greatest for larger eccentricities (compare Fig. 7 and Fig. 9). An interesting feature of the equatorial waveform is the prominent high frequency ripples superimposed on the waveforms, associated with the higher multipoles components ($`\mathrm{}=3`$ and higher. The illustrated waveform includes all multipoles up to $`\mathrm{}=18`$) of the wave. It is noteworthy that the high frequency features are prominent in both the zoom and the whirl parts of the orbit, although if they are solely the result of beaming we might expect that they would be purely a whirl feature, as the motion is fastest near periastron. However preliminary results from a time domain code written by one of us (KG) suggest that the high frequency features may be associated with quasi-normal modes of the black hole, which are excited by the high frequency emissions from the orbit. The quasi-normal mode ringing results in a continuous and time-delayed emission at these frequencies.

We also show waveforms seen by an observer on the polar axis of the black hole, in the top panel of Fig. 8 and the bottom panel of Fig. 9. In this case both “plus” and “cross” polarisations are present (only $`h_+`$ is nonzero for an equatorial observer) but we illustrate only $`h_+`$ because the “cross” waveform is the same except for a phase lag. The polar waveforms have the characteristic features of a high-amplitude, multi-cycle whirl part and a low-amplitude two cycle zoom part, but the high frequency features are absent. This suggests that the high frequency features are associated with beaming due to the rapid motion of the particle in the equatorial plane in the very strong field region, although the high-frequency features are not associated only with the strongest-field whirl part but are distributed throughout the whole cycle, including the weaker field (larger radius) zoom part. The polar waveforms are completely dominated by the quadrupole ($`\mathrm{}=2`$) emission. The $`\mathrm{}=3`$ and higher multipoles do not contribute significantly. In the equatorial waveforms the quadrupolar contribution is not much greater than that of the $`\mathrm{}=3`$ multipole and the fall off, in terms of the amplitude of $`h_+`$, for each subsequent multipolar waveform is slow. The transition between the polar and equatorial waveforms can by understood by looking at the waveform depicted in the bottom panel of Fig. 8 which corresponds to observation at angle $`\theta =\pi /4`$.

One can get an idea of what is going on by looking at the $`\theta `$-dependence of the energy flux from the system. In Fig. 12 one sees that the $`m=2`$ flux (dominated by the $`\mathrm{}=m=2`$ contribution) is concentrated somewhat towards the pole ($`\theta =0`$), but there is a strong shift towards the equator ($`\theta =\pi /2`$) in the $`m=3`$ flux, and further concentrations in that direction for each successively higher value of $`m`$. Therefore one sees the sort of beaming of the higher multipoles observed in the waveforms, although because of the dominance of the quadrupole the amplitude of the polar and equatorial waveforms is similar in this case.

Next, in Fig. 10, we show the waveform from the retrograde zoom-whirl orbit $`p=10.5M,e=0.5`$, retaining $`a=0.99M`$ for the black hole spin. The familiar zoom-whirling pattern is clear also in this case. However, we do not see any prominent high-frequency structure in these waveforms as the contribution coming from higher multipoles is small (because the orbit does not reside in a very strong field regime). We are not surprised, in this case, that the waveform seen when observing from a point along the hole’s polar axis (see bottom panel of Fig. 10) is not very different.

Although these zoom-whirl waveforms appear rather complex, and could present a problem for matched filtering data analysis techniques, one can take comfort in the fact that the number of harmonics of the spectrum which contribute significantly to the waveform and overall flux is not that great. Scott Hughes  has suggested that analysing individual ”voices,” which are monochromatic, of a complex wave may be an effective way to proceed in data analysis. Take the case of the zoom-whirl orbit with $`p=2.11M`$, $`e=0.7`$ and $`a=0.99M`$. If we look at Fig. 14, which shows the entire spectrum (up to $`\mathrm{}=18`$) of this signal, we see that the number of individual harmonics, or voices, each corresponding to an instance of the three numbers $`\mathrm{}`$,$`m`$ and $`k`$, which stand out are on the order of a dozen. Each of these voices has a rather simple waveform, as can be seen from Fig. 11. It is only their superposition which is complicated. Therefore following the chirp of individual voices as the orbit evolves may not be such a formidable computational task.

### D Spectra

The final part of our numerical results concerns the harmonic decomposition of the gravitational radiation fluxes. Specifically, we examine the $`k`$ distribution of the energy flux, at infinity and at the horizon, for a given multipole channel $`\mathrm{},m`$ (the angular momentum flux spectrum exhibits a similar behaviour). In all the figures in this Section, expect Fig. 14, we plot $`\dot{E}_{\mathrm{}mk}`$ versus $`k`$.

To begin with, in Fig. 13 we present the energy flux spectrum, at infinity, of the $`\mathrm{}=m=2,3,4`$ multipoles for an orbit with parameters $`p=2.11M,e=0.7`$, and for spin $`a=0.99M`$. As we have already mentioned, for such a strong-field orbit the higher multipoles give significant contributions to the total flux. The spectrum itself is composed of a series of “humps” which grow in height as $`k`$ increases, up to the point where the maximum harmonic is reached at $`k=k_{max}`$, after which the spectrum rapidly fades away. This behaviour closely resembles the one found in the Schwarzschild case , although in the present case the “humps” show a somewhat less regular structure. As we are not dealing with a very high eccentricity orbit, we find that $`k_{max}`$ is not very large. Moreover, the spectrum peak shifts to higher $`k`$ values as $`\mathrm{}`$ increases.

In Fig. 15 we graph the corresponding horizon fluxes for the same orbit and same multipoles as in Fig. 13. Both infinity and horizon fluxes peak at a common $`k_{max}`$ value. Since the selected orbit resides in a strong field regime, the horizon flux is a respectable fraction ($`10\%`$) of the flux at infinity, and moreover, is predominantly negative i.e. superradiant, as we should expect from the discussion in the previous Section. Note that the relative contribution of the quadrupole channel is much larger in the case of horizon fluxes as compared to the infinity fluxes.

The full energy flux spectrum, up to $`\mathrm{}=18`$, for the $`p=2.11`$,$`e=0.7`$ orbit is shown in Fig. 14 in terms of frequency rather that $`k`$.

Spectra belonging to retrograde orbits appear similar to the prograde spectra with the difference that they peak at negative values of $`k`$. For a typical situation see Fig. 16 for the retrograde orbit of $`p=10.4M,e=0.5`$ and black hole spin $`a=0.99M`$.

Before closing this Section, we should point out that only the harmonics around $`k_{max}`$ (typically 10-30 of them) give a significant contribution to the total flux, for given $`\mathrm{},m`$. This statement applies for all eccentricities we examined, and could have an important impact on the computation of waveforms and fluxes from highly eccentric orbits in which case $`k_{max}`$ attains very large values. For example, if one could find by other means (for example a time-domain code, see discussion in the next Section) the location of the main peak, then the calculation of the surrounding harmonics should give an sufficiently accurate result for the flux. This strategy is far more economic, from a computational point of view, than calculating all the harmonics between $`k=0`$ and $`k_{max}`$.

## VI Concluding discussion

We have examined gravitational waves from and radiation reaction of equatorial orbits of particles in the last stages of inspiral around a central black hole. We expect that such orbits will often have moderate eccentricities, not circular but significantly less than $`e=1`$. As these orbits approach the point at which they will plunge into the black hole two things will happen. The first is that the orbital eccentricity will begin to increase, having been decreasing throughout the preceding inspiral. For retrograde orbits this eccentricity increasing phase will last for many cycles, while for prograde orbits, especially ones which retain fairly high eccentricities, this phase will be fairly brief (in the case of rapidly rotating black holes). We have found a $`10\%`$ fractional increase in eccentricity, for the most favourable situations. Secondly, as the orbit draws closer to the unstable region it will tend to spend a longer and longer portion of each orbital period near periastron. For reasonably eccentric orbits this will promote the “zoom-whirl” behaviour we have described, with its characteristic waveform. Although extreme cases involving high eccentricities will presumably be rare, it is still likely that LISA will be seeing signals with $`510`$ whirls, and therefore $`1020`$ gravitational wave cycles between each apastron. Therefore templates of such signals will be very important. It is worth noting that because the orbital radius changes very little during these whirls, the whirl part of the waveform looks, in many respects, like a near circular waveform with harmonics of a single dominant frequency. This means, for instance, that for orbits with many whirls, the ratio of the fluxes of energy and angular momentum is close to $`\mathrm{\Omega }_\varphi `$ (recall the celebrated relation for circular orbits $`\dot{E}/\dot{L}=\mathrm{\Omega }_\varphi `$).

Considering inclined orbits (not confined to the equatorial plane) for a moment, we expect that zoom-whirl orbits will be found in cases with small inclination angles, because higher inclination angles feature greater plunge radii and the particle cannot ever have a very small periastron radius (as we see when the orbit is retrograde, which corresponds to the largest possible inclination angle). When the test-particle is in an inclined and eccentric orbit, close to crossing some separatrix of bound stable orbits, it will spend a considerable amount of time moving in a quasi-circular non-equatorial trajectory close to the periastron. It is plausible, by extrapolating the results of the present work, that during the “whirling” stage, the orbital energy, angular momentum and Carter constant will approximately evolve in such a way as if the orbit was circular and inclined. Assuming that most of the radiated energy, angular momentum and Carter constant are generated near the periastron (which in the present case is located in a strong field regime) then one might be able to estimate the, otherwise elusive, rate of change of the Carter constant for a near-to-plunge generic Kerr orbit. Such information could provide a very useful test for the recently adopted assumption of obtaining the rate of change of the Carter constant by keeping the orbital inclination angle fixed during the inspiral .

Waveforms from prograde orbits residing close to the horizon of the black hole will feature significant high frequency components when seen from a position on the equatorial plane of the system (i.e. when the orbit is observed “edge on”). This seems to be due to beaming, resulting from the rapid motion of the orbiting particle along the line of sight of the observer. When observed from on or near the polar axis the waveform is largely quadrupolar, dominated by a single nearly circular frequency. A glance at the equatorial zoom-whirl waveforms presented here suggests that they can be very complex and not necessarily amenable to matched filtering methods following the full wavetrain. But recall that a successful source identification may already have been made during the previous year of observation by LISA and from the source parameters deduced during this period it may be possible to search for the whirl parts of the waveform individually. On the other hand, it is worth noting that this high frequency structure could make it possible for LISA to detect late inspiral signals from very large black holes, with masses above $`10^7M_{}`$, which would otherwise be too low frequency (in the low multipole parts of the waveform) for detection. By contrast, if we are looking down on the system from the pole then the signal is much “cleaner”, without much contribution beyond the $`\mathrm{}=m=2`$ multipole. Obviously non-equatorial motion will introduce further harmonics and, as suggested in  it may prove more useful to examine different harmonics or “voices” of the signal separately, rather than trying to model the entire complex signal as one template.

It has been recently realised  that orbits during the long inspiral phase, before the radial orbital frequency is within the LISA waveband, will, in principle, emit detectable radiation. Recall that the periastron is rather close, $`r_p<20M`$, so that the frequency of the cycles in the whirl part of the waveform does fall within the LISA waveband. As the time between successive bursts (equal to the radial orbital period) will be very long (typically up to a century or even more) these bursts would not be expected to be detectable in practice. However, as the number of objects in this long inspiral stage will be rather great, one does have to take them into account as a background noise which will tend to hinder LISA’s efforts to detect signals from those objects in the last stages of their inspiral. In fact, one could in principle expect to have proper zoom-whirl orbits even at this stage, provided the central hole is rapidly spinning and the orbit is retrograde. Then, as we have shown, the corresponding separatrix translates to a minimum periastron value, $`r_{\mathrm{isbo}}6M`$, which falls within the expected periastron distribution.

In a follow-up paper  we will study orbits with larger radii and larger eccentricities. In particular we are planning to produce (using the method of  which was applied for inspiralling circular inclined orbits) the full inspiral trajectory and the resulting waveform for bound equatorial orbits that could be of importance for LISA. It is of great interest to (a) produce full inspiral waveforms for the kind of orbits described above and (b) give a good estimate of the total inspiral time, and the residual eccentricity (taking under consideration effects like the sign reversal of $`\dot{e}`$) just before plunging. However, the calculation of fluxes and waveforms produced by bodies in $`e1`$ orbits, is currently ranked as a difficult task. When pursued in the frequency domain (as it was the case in the present work) one has to calculate an enormous number of individual harmonics, as it is obvious from the spectra we presented. Moreover, the numerical computation of the integrals (125) is poorly convergent, a manifestation of the fact that the source term in the Teukolsky equation (52) diverges as $`r\mathrm{}`$ (that is, when the orbit tends to become parabolic $`e1`$). One way to cure this pathology would be to work with the inhomogeneous Sasaki-Nakamura equation  which has a well behaved source term (in the sense that it decays at spatial infinity). However, the first difficulty outlined above will still be present. A possible way to overcome it could be the calculation of the waveform/fluxes directly in the time domain by evolving the time-dependent Teukolsky equation without resorting to any separation of variables apart from $`\varphi `$. Conceivably, the required numerical code could be based on the Teukolsky codes used to study the dynamics of scalar and gravitational perturbations in a Kerr background metric , see  for a report on such an attempt. Looking further ahead, such time-domain codes could be the only practical tool for computing the waveform/fluxes generated by bodies orbiting non-black hole massive compact objects (in which case there is no known Teukolsky-like separable wave equation). We are currently working along both of these directions.

## VII Acknowledgements

The authors wish to thank Scott Hughes, Sterl Phinney, B.S. Sathyaprakash, Wolfram Schmidt and Lee Lindblom for very useful discussions and comments. In addition KG thanks the State Scholarships Foundation of Greece for financial support. He would also like to acknowledge support from PPARC grant PPA/G/0/1999/0214. DK thanks the Astronomy department of Oxford University for the generous use of their facilities and their unfailingly warm hospitality which greatly aided this research. In particular he would like to thank Julia Kennefick and Gavin Dalton. He would also like to acknowledge support from NSF grant PHY-0099568, and the help of the Physics and Astronomy department at Cardiff University.

## VIII Tables

## A Functions that appear in the solution for $`x^2`$.

The quantity $`x=L_zaE`$ satisfies the quartic equation,

$$F(p,e)x^4+N(p,e)x^2+C(p,e)=0$$
(A1)
where

$`F(p,e)`$ $`=`$ $`{\displaystyle \frac{1}{p^3}}[p^32M(3+e^2)p^2+M^2(3+e^2)^2p4Ma^2(1e^2)^2]`$ (A2)
$`N(p,e)`$ $`=`$ $`{\displaystyle \frac{2}{p}}[Mp^2+(M^2(3+e^2)a^2)pMa^2(1+3e^2)]`$ (A4)
$`C(p,e)`$ $`=`$ $`(a^2Mp)^2`$ (A6)
Defining the discriminant

$$\mathrm{\Delta }_x(p,e)=N^24FC=\frac{16a^2M}{p^3}[p^44Mp^3+2\{2M^2(1e^2)+a^2(1+e^2)\}p^24Ma^2(1e^2)p+a^4(1e^2)^2]$$
(A7)
the solution for $`x^2`$ is,

$$x^2=\frac{N\mathrm{\Delta }_x^{1/2}}{2F}$$
(A8)
where the upper (lower) sign corresponds to prograde (retrograde) motion.

## B Fundamental frequencies for bound equatorial orbits in Kerr geometry.

The fact that the function $`r(t)`$ is periodic in time (with period $`T_\mathrm{r}`$) implies that the function

$$\frac{d\varphi }{dt}=\frac{aT+\mathrm{\Delta }x}{(r^2+a^2)T+ax\mathrm{\Delta }}$$
(B1)
is also periodic and with the same period. Hence, it can be expanded in a Fourier series,

$$\frac{d\varphi }{dt}=\underset{k=\mathrm{}}{\overset{+\mathrm{}}{}}\beta _\mathrm{k}e^{ik\mathrm{\Omega }_\mathrm{r}t}$$
(B2)
By integrating this relation we get,

$$\varphi (t)=\beta _0t+\underset{k0}{}c_\mathrm{k}e^{ik\mathrm{\Omega }_\mathrm{r}t}+\text{(const)}$$
(B3)
where $`c_\mathrm{k}=i\beta _\mathrm{k}/k\mathrm{\Omega }_\mathrm{r}`$. This expression clearly shows that the $`\varphi `$-motion is of “rotation” type . Making use of the integration constant to define the term $`c_0`$, we have

$$\varphi (t)\beta _0t=\underset{\mathrm{k}}{}c_\mathrm{k}e^{ik\mathrm{\Omega }_\mathrm{k}t}$$
(B4)
That is, the function $`F(t)=\varphi (t)\beta _0t`$ is periodic with period equal to $`T_\mathrm{r}`$. As we have defined $`\mathrm{\Delta }\varphi =\varphi (t+T_\mathrm{r})\varphi (t)`$ we find that

$$\beta _0=\mathrm{\Omega }_\varphi \frac{\mathrm{\Delta }\varphi }{T_\mathrm{r}}$$
(B5)
## C Potentials of the Sasaki-Nakamura equation

We give an explicit listing of the potentials $`F(r)`$,$`U(r)`$ that appear in the Sasaki-Nakamura equation (166):

$`F(r)`$ $`=`$ $`{\displaystyle \frac{\eta _{,r}}{\eta }}{\displaystyle \frac{\mathrm{\Delta }}{r^2+a^2}}`$ (C1)
$`U(r)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }U_1(r)}{(r^2+a^2)^2}}+G(r)^2+{\displaystyle \frac{\mathrm{\Delta }G_{,r}}{r^2+a^2}}F(r)G(r)`$ (C2)
where the function $`\eta (r)`$ is given by

$$\eta (r)=c_0+c_1/r+c_2/r^2+c_3/r^3+c_4/r^4$$
(C3)
with the following coefficients

$`c_0`$ $`=`$ $`12i\omega M+\lambda (\lambda +2)12a\omega (a\omega m)`$ (C4)
$`c_1`$ $`=`$ $`8ia[3a\omega \lambda (a\omega m)]`$ (C5)
$`c_2`$ $`=`$ $`24iaM(a\omega m)+12a^2[12(a\omega m)^2]`$ (C6)
$`c_3`$ $`=`$ $`24ia^3(a\omega m)24Ma^2`$ (C7)
$`c_4`$ $`=`$ $`12a^4.`$ (C8)
In addition, the functions $`G(r)`$ and $`U_1(r)`$ are

$`G(r)`$ $`=`$ $`{\displaystyle \frac{2(rM)}{r^2+a^2}}+{\displaystyle \frac{r\mathrm{\Delta }}{(r^2+a^2)^2}}`$ (C9)
$`U_1(r)`$ $`=`$ $`V(r)+{\displaystyle \frac{\mathrm{\Delta }^2}{\beta }}\left[\left(2\alpha +{\displaystyle \frac{\beta _{,r}}{\mathrm{\Delta }}}\right)_{,r}{\displaystyle \frac{\eta _{,r}}{\eta }}\left(\alpha +{\displaystyle \frac{\beta _{,r}}{\mathrm{\Delta }}}\right)\right]`$ (C10)
with

$`\alpha `$ $`=`$ $`{\displaystyle \frac{i\beta K}{\mathrm{\Delta }^2}}+3iK_{,r}+\lambda +{\displaystyle \frac{6\mathrm{\Delta }}{r^2}}`$ (C12)
$`\beta `$ $`=`$ $`2\mathrm{\Delta }(iK+rM2\mathrm{\Delta }/r)`$ (C15)
The functions $`K(r),V(r)`$ appear in the Teukolsky equation (52).