New Kludge Scheme for the Construction of Approximate Waveforms for Extreme-Mass-Ratio Inspirals

New Kludge Scheme for the Construction of Approximate Waveforms for Extreme-Mass-Ratio Inspirals

Carlos F. Sopuerta Institut de Ciències de l’Espai (CSIC-IEEC), Facultat de Ciències, Campus UAB, Torre C5 parells, Bellaterra, 08193 Barcelona, Spain.    Nicolás Yunes Department of Physics, Montana State University, Bozeman, MT 59717, USA. Department of Physics and MIT Kavli Institute, 77 Massachusetts Avenue, Cambridge, MA 02139, USA. Princeton University, Physics Department, Princeton, NJ 08544, USA.
July 12, 2019

We introduce a new kludge scheme to model the dynamics of generic extreme mass-ratio inspirals (stellar compact objects spiraling into a spinning supermassive black hole) and to produce the gravitational waveforms that describe the gravitational-wave emission of these systems. This scheme combines tools from different techniques in General Relativity: It uses a a multipolar, post-Minkowskian expansion for the far-zone metric perturbation (which provides the gravitational waveforms, here taken up to mass hexadecapole and current octopole order) and for the local prescription of the self-force (since we are lacking a general prescription for it); a post-Newtonian expansion for the computation of the multipole moments in terms of the trajectories; and a BH perturbation theory expansion when treating the trajectories as a sequence of self-adjusting Kerr geodesics. The orbital evolution is thus equivalent to solving the geodesic equations with time-dependent orbital elements, as dictated by the multipolar post-Minkowskian radiation-reaction prescription. To complete the scheme, both the orbital evolution and wave generation require to map the Boyer-Lindquist coordinates of the orbits to the harmonic coordinates in which the different multipolar post-Minkowskian quantities have been derived, a mapping that we provide explicitly in this paper. This new kludge scheme is thus a combination of approximations that can be used to model generic inspirals of systems with extreme mass ratios to systems with more moderate mass ratios, and hence can provide valuable information for future space-based gravitational-wave observatories like the Laser Interferometer Space Antenna and even for advanced ground detectors. Finally, due to the local character in time of our multipolar post-Minkowskian self-force, this scheme can be used to perform studies of the possible appearance of transient resonances in generic inspirals.


I Introduction

Gravitational waves (GWs) hold the promise to provide detailed information about astrophysical bodies that are obscure in the electromagnetic spectrum, such as binary black hole (BH) systems. Moreover, such waves will allow for the first studies of the nature of the gravitational interaction and of the validity of General Relativity (GR) in the strongest regimes Will (2006); Yunes and Pretorius (2009). An accurate modeling of such GWs is essential for the extraction and characterization of weak signals buried in detector noise. This is because waveform templates act as an optimal linear filter that maximizes the signal-to-noise ratio (SNR) in the presence of stochastic noise. The absence of such templates for certain GW sources renders suboptimal any GW search strategy. Therefore, the construction and modeling of GWs to construct accurate templates for data analysis is of paramount importance in the blossoming of GW astrophysics.

One of the staple GW sources of the planned space-based observatory Laser Interferometer Space Antenna (LISA) Danzmann and Rudiger (2003); *Danzmann:2003ad; *Prince:2003aa; lis () that require accurate templates for their detection and analysis are extreme-mass-ratio inspirals (EMRIs) Amaro-Seoane et al. (2007). These events consist of a small compact object (SCO), such as a stellar-mass BH or neutron star (NS), spiraling in a generic orbit into a spinning, (super)massive black hole (MBH), whose evolution is GW dominated. In such inspirals, the SCO spends up to millions of cycles in close orbits around the MBH, possibly with large pericenter velocities and eccentricities, sampling the strong gravitational field of the MBH.

Many astrophysical scenarios predict the existence of such EMRIs. One such scenario postulates that the SCO exchanges energy and angular momentum with other stars in a stellar core/cusp near a MBH at a galactic center, via two-body relaxation and dynamical friction Amaro-Seoane et al. (2007). If so, the SCO can be swung sufficiently close to the MBH to be gravitationally captured Hils and Bender (1995); *1997MNRAS.284..318S; Amaro-Seoane et al. (2007), at which point it would slowly inspiral until being swallowed by the MBH. Of course, such a scenario is complicated by mass segregation Freitag et al. (2006); *2006ApJ...645L.133H; *2006RPPh...69.2513M, triaxial density profiles Holley-Bockelmann et al. (2002); *2005ApJ...629..362H, resonant relaxation Rauch and Tremaine (1996); *2006ApJ...645.1152H, etc. Other channels of EMRI formation include binary tidal separation Miller et al. (2005) (where a binary is disrupted with one component captured by the MBH) and massive star capture or production in accretion discs Levin (2003); *2007MNRAS.374..515L, where the stellar-mass BH is directly formed in the accretion disk of the MBH. It turns out that the EMRI GWs produced by EMRIs in each astrophysical mechanism have distinct characteristic that may be use to distinguish them from EMRI observations Yunes et al. (2011); *2011PhRvD..84b4032K.

The inspiral is by far the dominant GW phase for EMRI data analysis purposes. This can be understood rather easily (see, e.g. Finn and Thorne (2000); Yunes et al. (2011)). The number of cycles accumulated in the inspiral scales with the inverse of the mass ratio , while in the plunge and merger it scales with the MBH mass only. Since the mass ratio for EMRIs is in the range -, the characteristic duration and cycle accumulation during the inspiral phase is several orders of magnitude larger than during plunge and merger. In turn, the SNR is in the range - for EMRIs at realistic distances and it scales linearly with the total number of cycles. Thus, the contribution of the plunge and merger to the SNR is reduced by a factor of relative to the inspiral contribution. In addition, there is no detectable ringdown for EMRIs, as the SCO barely perturbs the background geometry as it crosses the MBH’s event horizon. Since the SCO is not disrupted in EMRI plunges, the SCOs internal structure is erased or effaced almost completely, without affecting the inspiral signal.

A positive consequence of the large number of EMRI GW cycles is that these waves carry a detailed map of the MBH geometry, so that we expect to determine the EMRI physical parameters with high precision Barack and Cutler (2004). This information will be very useful, in particular to test the spacetime geometry of MBHs Collins and Hughes (2004); *Glampedakis:2005cf; *Barack:2006pq; *Yunes:2008gb; *Brink:2008xx; *Vigeland:2009pr; *Hughes:2010xf; Sopuerta (2010) and even alternative theories of gravity (see, e.g. Hughes (2006); Schutz (2009); Sopuerta (2010); Babak et al. (2011); *Vigeland:2011ji; *Gair:2011ym). Moreover, given that the expected even rate is in the range EMRIsyr Gair et al. (2004); Hopman and Alexander (2006), EMRI observations will allow us to understand better the stellar dynamics near galactic nuclei, populations of stellar BHs, etc. (for a review see Amaro-Seoane et al. (2007)), and it also possible that they will tell us about cosmology MacLeod and Hogan (2008); *Gair:2008bx.

Although one is left to model the inspiral phase of EMRIs, this task remains a gargantuan endeavor for several reasons. First, the accuracy requirements for EMRI templates are much more stringent than for comparable-mass binaries. For detection and parameter estimation one usually demands an absolute accuracy of better than and SNR radians in the GW phase respectively over the entire time of observation. In a yr inspiral, a typical EMRI can have cycles in-band, which then translates into a relative radian accuracy of and for detection and parameter estimation respectively. It is important to note that these precision requirements are just simple estimates that do not take into account data analysis strategies that may relax them (see, e.g. Gair et al. (2004)). In contrast, the above relative measure becomes of for ground-based data analysis of comparable-mass plunge-merger-ringdowns with current numerical relativity simulations, because these accumulate only GW cycles in these phases. In between EMRIs and comparable-mass inspirals, there are Intermediate-Mass-Ratio Inspirals (IMRIs), which can be potential sources for both space-based detectors (where an Intermediate-Mass BH (IMBH) inspirals into a MBH) and advanced ground-based detectors (where a SCO inspirals into an IMBH).

The extreme mass ratios involved in the problem also lead to the appearance of two different spatial and time scales. The two different spatial scales are represented by the very different sizes of the MBH and the SCO, , whereas the different time scales are the orbital one and the one associated with the radiation-reaction effects, . An illustration of how this complicates the EMRI problem is the recent work of Lousto and Zlochower Lousto and Zlochower (2011); *Nakano:2011pb, who evolved the first mass-ratio binary over the last two orbits before merger and plunge; this simulation took approximately days of computational time using full numerical relativity. This means that with present numerical relativity techniques, full numerical simulations are out of the question for EMRI modeling.

Another key reason for the difficulty of EMRI modeling is their intrinsic strong-relativistic nature. In the interesting part of the EMRI dynamics, the SCO is moving around the strong-field region of the MBH, acquiring large pericenter velocities and even sampling regions inside the ergosphere, leading to large relativistic -factors over tens of thousands of GW cycles. Approximate techniques employed in the comparable-mass regime, such as low-velocity expansions in the post-Newtonian (PN) approximation, are then ill-suited for EMRIs. So neither numerical relativity nor PN theory are, for very different reasons, suitable schemes to model EMRIs.

A better-suited framework that exploits the extreme mass ratios involved is BH perturbation theory, where one treats the SCO as a small perturbation of the MBH background geometry. In this context, the inspiral can be described as the action of a self-force. This local vector force is made out of the regularized metric perturbations generated by the SCO, after eliminating divergences due to the particle description of the SCO. The SCO’s motion is then governed by the MiSaTaQuWa equation of motion equation, derived in Mino et al. (1997); Quinn and Wald (1997) (see also Sec. II). The MiSaTaQuWa equation is considered the foundation of a self-consistent scheme to describe EMRIs in an adiabatic way by coupling it to the partial differential equations that describe the perturbations produced by the SCO. For recent discussions on these issues see Gralla and Wald (2008); Pound (2010, 2010) and for general reviews see Glampedakis (2005); Poisson (2004); Barack (2009).

At present, the gravitational self-force has been computed for the case of a nonrotating MBH using time-domain techniques Barack and Sago (2009); *Barack:2010tm (see also Barack et al. (2010); Barack and Sago (2011); Tiec et al. (2011) for the study of the physical consequences of the self-force) and progress is being made towards calculations for the more astrophysically relevant case of a spinning MBH Shah et al. (2011). In the meantime, a number of new techniques in the frequency and time domains are being developed to produce accurate and efficient self-force calculations Sopuerta et al. (2006); *Barack:2007jh; *Barack:2007we; *Vega:2007mc; *Lousto:2008mb; *Canizares:2008dp; *Canizares:2009ay; *Field:2009kk; *Canizares:2010yx; *Thornburg:2010tq; *Canizares:2011ft; *Warburton:2011hp; *Dolan:2011dx. In any case, given the amount of cycles required for EMRI GWs and the present complexity of self-force calculations, we cannot expect to generate complete waveform template banks by means of full self-force calculations. Instead, the goal of these studies should be to understand all the details of the structure of the self-force so that we can formulate efficient and precise algorithms to create the waveforms needed for LISA data analysis, perhaps complementing some of the existent approximating schemes that we review below.

i.1 Existent EMRI Waveform Models

In parallel to the efforts to make progress in the self-force program, there have also been some efforts to build certain less-reliable approximation schemes to model EMRIs. These are very useful, for instance, for parameter estimation studies. We compare and contrast these in Table 1, ordered by level of complexity from simplest (top) to most complex (bottom), where we also include the new Kludge scheme at the bottom for comparison.

Scheme Name Orbital Motion Radiation-Reaction Waveform Generation
Peters & Mathews Peters and Mathews (1963); Peters (1964) Keplerian Ellipses Newtonian Order Multipolar decomposition ()
Analytic Kludge Barack and Cutler (2004) Keplerian Ellipses Low-Order PN Multipolar decomposition ()
Numerical Kludge Babak et al. (2007) Kerr Geodesics Calibrated Low-Order PN Multipolar decomposition ()
EOB Yunes et al. (2010); Yunes (2010); Yunes et al. (2011) Kerr Geodesics Calibrated, Resummed 5.5PN Resummed 5.5PN ()
Teukolsky-based Hughes (2000, 2001) Kerr Geodesics Averaged Teukolsky Adiabatic Teukolsky ()
This work Kerr Geodesics Local Multipolar Post-Minkowskian Multipolar decomposition ()
Table 1: Comparison of existent modeling schemes.

The simplest model is that of Peters and Mathews Peters and Mathews (1963); Peters (1964), in which the SCO is assumed to be moving on Keplerian ellipses. The orbital elements of this ellipse evolve according to leading-order (Newtonian), dissipative radiation-reaction, i.e. the quadrupole formula for the loss of energy and angular momentum. Waveforms are then computed also to leading order via the quadrupole formula Misner et al. (1973).

A better model was introduced by Barack and Cutler Barack and Cutler (2004), the so-called Analytical Kludge waveform model. This model is based on the Peters & Mathews Peters and Mathews (1963); *Peters:1964zz model, but it is enhanced via different PN formulas in order to account for all the relativistic effects, both dissipative and conservative, present in a generic EMRI event. It has the advantage that the orbital evolution is decoupled from the evolution of the additional relativistic effects (which evolve effectively in the radiation-reaction time scale), and hence EMRI waveforms, also prescribed via the quadrupole formula Misner et al. (1973), can be computed very fast. For this reason, it has become the method of choice for many LISA parameter estimation and data analysis studies (see, e.g. Arnaud et al. (2007)).

A more sophisticated approach was proposed by Babak, et. al. Babak et al. (2007) and is sometimes referred to as the Numerical Kludge waveform model. In this setup, the orbital motion is given by geodesics around a Kerr BH and the radiative effects are prescribed via PN evolution equations for orbital elements (from 2PN expressions for the fluxes of energy and angular momentum) calibrated to more accurate Teukolsky fluxes with fitting parameters Gair and Glampedakis (2006). The waveforms are then modeled again via a multipolar expansion Thorne (1980), but this time taken to next-to-leading-order (quadrupole plus octopole).

Recently, a new hybrid scheme has been proposed by Yunes, et. al. Yunes et al. (2010); Yunes (2010); Yunes et al. (2011) based on effective-one-body (EOB) techniques Buonanno and Damour (1999); *Buonanno:2000ef; *Damour:2009sm. In this approach, the SCO-MBH, two-body system is mapped to an effective one-body system: a Kerr BH perturbed by a small effective object. The orbital motion is obtained by solving the Hamilton equations for the Hamiltonian of the effective system. When neglecting conservative self-force corrections, this reduces to solving the geodesic equations in the Kerr background. The radiation-reaction comes from the short-wavelength approximation of Isaacson’s Isaacson (1968); *Isaacson:1968gw, where the waveform is constructed from an orbit-averaged, but resummed PN expression Damour et al. (1998). This radiation-reaction force is then enhanced through the addition of very high PN order point-particle results Tanaka et al. (1996); Mino et al. (1997) and expressions that account for the flux of energy and angular momentum into the MBH’s horizon Shibata et al. (1995); Mino et al. (1997). Although shown to be accurate for equatorial, circular orbits Yunes et al. (2010) relative to Teukolsky waveforms, the EOB scheme has not yet been tested for eccentric or inclined orbits.

Another approach is the Teukolsky-based scheme of Hughes and others Poisson (1993); Hughes (2000, 2001); Glampedakis and Kennefick (2002); Drasco and Hughes (2006); Drasco et al. (2005); Hughes et al. (2005); Lopez-Aleman et al. (2003); Khanna (2004); Burko and Khanna (2007); Sundararajan et al. (2007, 2008) have developed. In this model, one prescribes the inspiral as a sequence of slowly-changing geodesics. The mapping between them is given by the orbit-averaged evolution of orbital elements, which in turn is obtained by a balance law relating averaged fluxes at the boundaries of spacetime and at the location of the SCO. These averaged fluxes at a given geodesic or point in orbital phase space are computed by solving the Teukolsky equation at that point. Therefore, the construction of any single waveform requires the mapping of the entire orbital phase space, which in turn is computationally prohibitive for truly generic EMRIs. Moreover, this scheme has numerical (and conceptual) difficulties when modeling EMRIs in regimes of spacetime where the evolution deviates from adiabaticity, such as in or close to plunge, around rapidly spinning MBHs (), or highly eccentric inspirals.

i.2 The New Kludge Scheme

In this paper, we devise a new kludge approximation scheme (relative to numerical kludge) that combines seemingly disparate ingredients from BH perturbation theory and the multipolar post-Minkowskian formalism of Blanchet and Damour Blanchet and Damour (1984). We shall here combine two different approximation schemes, BH perturbation theory (which assumes only that the mass-ratio is small) and the multipolar post-Minkowskian formalism (which assumes the gravitational field strength is small), together with other ingredients related to the choice of coordinate system and waveform construction method.

In the new kludge scheme, the orbital motion is prescribed as a spacetime trajectory that is piecewise geodesic (with respect to the MBH geometry) and such that the different geodesic intervals are connected via the SCO’s local, self-acceleration (due to its own gravity in the presence of the MBH). In this sense, each geodesic interval can be chosen arbitrarily small. This is in contrast to the Teukolsky approach (see, e.g. Hughes (2000)) where the mapping between sequences is given by the averaged (over several orbits) GW energy-momentum fluxes and balance laws. This fact, that the mapping is given purely in terms of quantities local to the SCO’s worldline and not by nonlocal balance laws, is a distinctive feature of the new kludge scheme. The implementation of this idea uses the evolution of geodesics with varying orbital elements; the energy, angular momentum in the spin direction, and Carter constants are then functions of time governed by the self-force.

Such an evolution scheme is a direct implementation of the osculating orbits method, proposed by Pound and Poisson Pound and Poisson (2008) and Pound Pound (2010). This method professes that at each point of the SCO’s nongeodesic worldline there is a unique geodesic that lies tangent to it. Therefore, the worldline is simply an interpolation between these tangent geodesics. Such a scheme hinges on a fundamental assumption of EMRI modeling: the adiabatic approximation, which assumes the deviation vector between adjacent tangent geodesics is small, i.e. the radiation-reaction time scale is much longer than all other timescales, particularly the orbital one. As Gralla and Wald Gralla and Wald (2008) explained, the adiabatic approximation only holds quasilocally around some small neighborhood of proper time at each point of the SCO’s worldline. This problem can be circumvented, however, if at each point on the worldline the deviation vector is recomputed, which is exactly the basis of the method of osculating orbits Pound and Poisson (2008); Pound (2010).

In the new kludge scheme, the evolution of orbital elements is prescribed by the SCO’s self-acceleration, but this quantity is not known exactly (numerically or otherwise) for generic orbits around a spinning MBH. In view of this, we model the self-force via a radiative approximation, i.e., through the time-asymmetric part of the radiation field, given by the “half-retarded minus half-advanced” Green function Mino (2003). This can be implemented with a high-order multipolar, post-Minkowskian expansion, which assumes gravitational radiation is small relative to the gravitational field of the background, i.e. an expansion in powers of Newton’s gravitational constant .

The new kludge, local self-force prescription is completed once we say how the multipole moments depend on the orbital trajectories. This can be achieved by asymptotically matching a PN and a post-Minkowskian solution Burke (1971); Blanchet and Damour (1988, 1989); Damour and Iyer (1991); Blanchet (1993); Poujade and Blanchet (2002). In this paper, we use only the leading-order expressions for these multipoles, although in the future it would be trivial to including higher PN corrections. As such, we are not consistently keeping a given PN order, but instead using leading-order expressions for the multipole moments, while keeping several such moments in the expansions.

The use of the radiative, multipolar post-Minkowskian expansion is only because of the lack of a more precise self-force. Clearly, this radiative approximation neglects the conservative part of the self-force, which could be important in GW modeling Pound and Poisson (2008). Once the full self-force becomes available, however, one could easily employ it instead of its post-Minkowskian expansion. Our set-up is general and easily adaptable to other, more precise expressions for the self-force.

Once the orbital evolution has been prescribed, one can construct the GWs again through multipolar, post-Minkowskian expressions in terms of a sum over multipole moments. Since the mapping between Boyer-Lindquist and harmonic coordinates is known, there are no coordinate issues to relate the trajectories obtained from the orbital evolution to the trajectories that enter the definition of the multipole moments. We here employ an expansion to second-order in the multipole moments, including both the mass hexadecapole and the current octopole, thus keeping contributions one order higher than traditional kludge waveforms.

Let us emphasize that the new kludge waveforms are only approximate gravitational-wave solutions, meant to be used for qualitative descoping studies and investigations of resonant behavior. That is, the waveforms constructed could be used to determine the accuracy to which parameters could be extracted given a detection with new space-based gravitational-wave observatories as a function of the detector. Moreover, one could also study how the resonances found by Flanagan and Hinderer Flanagan and Hinderer (2010) affect such parameter estimation and detection. Having said that, improvements would have to be implemented before our kludge waveforms can be used as realistic templates in data analysis. As we will see, the new kludge scheme is amenable to such improvements, which will be the focus of future work.

i.3 Comparison to Standard Kludge Waveforms

The main advantage of the new kludge scheme is the fact that the radiation-reaction effects are described in terms of a local self-force. This means that the inspiral description does not need to average certain gravitational-wave fluxes over a number of periods/orbits like is traditionally done in kludge implementations. Instead, the self-force is prescribed through a multipolar, post-Minkowskian expansion (e.g. the quantity mentioned above) that contains time-derivatives of the system multipole moments in a nonaveraged form.

When implementing such a nonaveraged and local scheme, it is critical to use an exact mapping between Boyer-Lindquist coordinates, used in the integration of the geodesic equations of motion, and harmonic coordinates, employed in the calculation of the multipolar, post-Minkowskian self-force. This eliminates gauge issues that plague kludge waveforms due to the neglect of such a mapping (i.e. kludge schemes simply use Boyer-Lindquist like Cartesian coordinates in the multipolar decomposition of the GWs).

The use of a local self-force provides the freedom to choose how often to apply radiation-reaction effects in the numerical implementation of the dynamics (trajectory and waveform construction) of the new kludge scheme. The two extremes are: (i) We can apply the self-force at every single time step, which corresponds to the case of a continuous local self-force, or (ii) we can store the information about the self-force and apply it after a certain period of time, mimicking the averaging procedure of other schemes.

These two extreme ways of using the new kludge scheme are in correspondence with two very relevant potential applications. Whereas the type of application (ii) can be used to try to generate efficiently EMRI gravitational-wave templates for parameter estimation and data analysis development purposes, the type of application (i) seems to be very well fitted for studying local phenomena in the dynamics of EMRIs. In this sense, an important application of our scheme would be to study the transient resonances that Flanagan and Hinderer Flanagan and Hinderer (2010) have recently reported in generic EMRI orbits (eccentric and inclined) when the fundamental orbital frequencies become commensurate. The new kludge scheme can be in principle used to study such behavior and shed some light on the relevance that it may have for future EMRI detection with LISA-like detectors.

Let us clearly state, however, that the new kludges presented here are not intended to model comparable-mass systems, for which the effective-one-body framework has proven the most successful Buonanno and Damour (1999, 2000). Instead, we here focus on EMRIs only and we are interested in studying local phenomena that no other kludge scheme can currently handle. Perhaps in the future, one could also use effective-one-body tools for such studies, but this would require a tested framework capable of handling completely generic orbits. Work along this lines is currently underway Yunes et al. (2010); Yunes (2010); Yunes et al. (2011).

i.4 Executive Summary of Main Results

In this paper we present the new kludge scheme, introducing one by one each of its ingredients, from the form of the equations of motion for the inspiral to the GW construction, including the details of the approximations that we use to construct the local self-force that drives the inspiral. From a technical point of view, one of the main challenges of the numerical implementation of our scheme is the computation of high-order time derivatives (of the mass and current multipole moments), which are crucial for the estimation of the radiation-reaction effects (the multipolar post-Minkowskian self-force involves up to eight-order time derivatives of the trajectory in harmonic coordinates) and the GW construction (since we are using up to the mass hexadecapole and current octopole multipoles in the calculation, we require up to fourth-order time derivates of the trajectory in harmonic coordinates). The computation of these time derivatives is very challenging, forcing us to implement numerical techniques adapted to the properties of EMRIs dynamics. The key point is to use the fact that geodesic orbits have, in the generic case of eccentric and inclined trajectories, three fundamental frequencies. This then allows us to fit any quantity that needs differentiating to a multiple Fourier series, using a standard least-squares technique. Numerical derivatives of such quantities can then be obtained simply by analytically differentiating the Fourier expansion. Numerical experimentation has shown that this technique works remarkably well, even for the highest-order time derivatives that the new kludge requires.

To illustrate the capabilities of our scheme we show in this paper results from evolutions for different types of orbits: circular equatorial, eccentric equatorial, circular inclined, and the generic eccentric inclined orbits. Using these evolutions we study different aspects of the new kludge scheme. First, we consider the impact of harmonic coordinates in the trajectories and waveform observables in comparison with using other coordinate systems. We find that not properly accounting for this transformation can lead to huge errors in the amplitude and phase of the waveform, eg. up to errors of order a factor of 2 in the total accumulated cycles after a yr evolution. Second, we study the impact of the different radiation-reaction potentials in the resulting waveforms. Although these corrections have a smaller impact than the proper use of coordinates, they are still large for strong-field EMRIs. For example, including higher-order terms to the Burke-Thorne potential leads to corrections of order radians in a two month evolution. Third, we investigate the use of a quadrupole waveform prescription versus a more accurate hexadecapole-octopole prescription. We find no change in the resulting waveform phases, but an amplitude correction of less than . Fourth, we consider the time-evolution of different orbital parameters when we apply radiation-reaction effects locally in time through our multipolar, post-Minkowskian self-force. Although there are some orbital parameters whose time evolution is not affected, other parameters like the inclination angle can acquire small oscillations with period equal to the orbital period. Such small oscillations are not captured if the self-force is orbit-averaged.

Finally, we perform some tests and comparisons with results in the literature to validate the new kludge numerical implementation. In particular we test the prediction that the inclination angle remains almost constant during the evolution, a test that our scheme successfully passes.

i.5 Notation and Organization of the Paper

Throughout this paper we use the metric signature and geometric units in which . The MBH geometry, whose metric we denote by , is determined by its mass and (magnitude of the) spin angular momentum , with dimensionless, Kerr spin parameter (). The SCO is parameterized only in terms of its mass since we neglect its spin and other internal properties. The binary system’s parameters are the mass ratio and the total mass . The reduced mass of the system is therefore , while the symmetric mass ratio is . The mass difference is denoted by .

The SCO orbit can be parameterized in terms of the constants of geodesic motion , which stand for the SCO’s energy normalized with respect to , the -component of the angular momentum normalized to , and the Carter constant also normalized to ( and stand for two definitions of the Carter constant that we use in this paper). Alternatively, we will also parameterize the SCO orbit in terms of the orbital elements , which are also constants of the geodesic motion, where is the orbit eccentricity, is the semilatus rectum, and and are two measures of the orbit inclination. We present the mappings between and in Appendix E. The SCO spacetime trajectory is denoted via , where is proper time and thus, its four-velocity is the unit timelike vector .

Post-Newtonian orders always refer to a relative ordering scheme (instead of an absolute one), such that the -th PN order term refers to one of the form , where is the leading order contribution. We shall commonly drop the factor of when referring to PN expansions.

Greek letters in index lists are used to denote indices on the -dimensional spacetime, while Latin letters in the middle of the alphabet denote spatial indices only. Covariant differentiation is denoted using the symbol , while partial derivatives with respect to the coordinate are denoted as or . We denote symmetrization and antisymmetrization with parenthesis and square brackets around the indices respectively, such as and .

We use two main sets of coordinate systems: Boyer-Lindquist coordinates and harmonic coordinates . Other systems of coordinates that we also consider are asymptotic-Cartesian mass-centered (ACMC) coordinates, , and approximate harmonic coordinates, . Retarded time is denoted in harmonic coordinates via , where . When we refer to the Kerr metric, we sometimes use the label , e.g. and . In some situations, it is crucial to specify in which coordinate system the metric has to be written in a certain equation, like in a coordinate transformation. In those situations, we also incorporate in the metric (or related objects) a label associated with the coordinate system, e.g. or . As for angular coordinates, we shall find it convenient to sometimes perform multipolar decompositions as in Thorne (1980), with spin-weighted spherical harmonics and symmetric and trace-free spherical harmonic tensors , where is a multi-index (see Thorne (1980) for details on the multi-index notation).

The organization of this paper is as follows. Section II introduces the self-force approach to EMRI modeling and some basic details of the method of osculating orbits. Section III describes in detail the new kludge scheme, including the MBH geometry and the properties of the geodesic orbits, the application to them of the method of osculating orbits, the multipolar post-Minkowskian approximation to the self-force we use and the radiation-reaction potentials from which it can be obtained, the mapping between Boyer-Lindquist and harmonic coordinates, and the multipolar post-Minkowskian waveform generation formalism that we use in our scheme. Section IV explains the numerical implementation of the new kludge approach and the main numerical techniques that we use. This includes the algorithms for the integration of the ODEs governing the local geodesic motion and the accurate estimation of time derivatives of several orders of the multipole moments. Section V summarizes the different ways in which our scheme can be used and presents several numerical results that illustrate its main features, using from circular-equatorial orbits to eccentric inclined orbits. Section VI concludes and discusses several avenues for future research applying the new kludge scheme.

We have attempted to present the main ingredients of our approach in the main body of the paper, relegating some details to the Appendices. Appendix A gives explicit expressions of the different pieces of the multipolar post-Minkowskian self-force in terms of the radiation-reaction and local potentials. It also gives formulas to simplify the computation of the different spatial derivatives of the Kerr local potentials. Appendix B provides complementary formulas related to the mapping between Boyer-Lindquist and harmonic coordinates. In particular, we give expressions for the components of the Jacobian and Hessian, and for the components of the covariant and contravariant Kerr metric tensor in harmonic coordinates. Appendix C constructs a system of asymptotically, mass-centered coordinates and a system of approximate harmonic coordinates that we compare with the exact harmonic coordinates of Sec. III.4. Appendix D performs a far-field expansion of the Kerr metric coefficients in the approximate harmonic coordinates of Appendix C. Appendix E describes in detail how to implement the one-to-one mapping between the orbital elements and the constants of motion . Appendix F summarizes the main formulas for the computation of the fundamental frequencies and periods with respect to the Boyer-Lindquist coordinate time. Finally, Appendix G provides expressions for the coefficients that determine the evolution of the radius and Carter constants of an inspiral through circular nonequatorial geodesics.

Ii The Self-Force Approach to EMRIs

In this section we review some of the basic and well-established concepts related to the self-force as they are related to the new kludge approach (see the review papers Poisson (2004); Poisson et al. (2011); Barack (2009) for details).

The foundations for the first-order perturbative description of EMRIs were laid down in the papers by Mino, Tanaka, and Sasaki Mino et al. (1997) and Quinn and Wald Quinn and Wald (1997). The main result of these papers was the equation of motion of a massive pointlike object (describing the SCO) in the geometry of a MBH. The stress-energy tensor of the SCO is then given by


where g denotes the determinant of . Then, the SCO generates metric perturbations, , around the MBH background geometry that in the Lorenz gauge,


satisfy the linearized Einstein equations


where is the Riemann tensor of the MBH background geometry. However, according to this equation, the metric perturbations diverge at the particle location. Then, the gravitational backreaction on the particle motion, the self-force, is provided by the regularized part of the perturbations, say , according to a Hadamard prescription given in Mino et al. (1997). Then, the equation of motion for an EMRI, the MiSaTaQuWa equation, is


where the self-force, , is given by


Therefore, the dynamics of EMRIs is determined by the coupled system of Eqs. (4),(5), and (2) with a practical regularization scheme (like the mode sum scheme Barack et al. (2002)). A remarkable point is that Eqs. (4) and (5) can be rewritten as geodesic equations of motion for a point particle in a perturbed geometry that only take into account the regularized part of the metric perturbations, i.e. geodesic in the geometry  Detweiler (2001).

But how do we evolve the trajectory of the SCO accounting for the self-force in a self-consistent way? As it was proposed in Pound and Poisson (2008) (see also Gralla et al. (2009) and Gair et al. (2011) for a recent use of this technique), one can use a relativistic extension of the well-known method of osculating orbits. The idea is to consider always the geodesic tangent to the trajectory, so that the motion transitions from one geodesic to the next. Such smooth transition is facilitated in EMRI by the fact that EMRI trajectories are very close to a (local) geodesic for a long time. This is because of the clean separation in EMRI time scales: the radiation-reaction time scale is much larger than the orbital (geodesic) one, except for the tiny fraction corresponding to the merger-plunge phase.

The way to carry out this transition is to properly account for the time-evolution of the set of orbital elements that completely characterize a geodesic orbits. Following Pound and Poisson (2008), it is important to distinguish between two sets of orbital elements: principal orbital elements, in our case either or ; and positional orbital elements that determine the initial position in the geodesic as well as the geodesic initial spatial orientation. The radiation-reaction changes in the principal elements are due to the dissipative part of the self-force, while radiation-reaction changes in the positional elements are due to the conservative part of the self-force. In this work we only consider dissipative effects and hence we are only concerned with changes in the principal orbital elements.

The implementation of the method of osculating orbits consists in the translation of the fact that at any time there will be a geodesic trajectory, , with orbital elements whose position and velocity at that time will coincide with those of the accelerated trajectory. This can be written in the following way:


where denote the positional orbital elements. Although we have used here the constants of motion as principal orbital elements, we could also have used . Combining these osculation conditions with the equations of motion (4), we can arrive at the following equations Pound and Poisson (2008):


where is the SCO self-acceleration, which is related to the self-force of Eq. (5) by . This is due to the fact that has been defined as the SCO constants of motion per unit mass. From the inversion of these equations we can obtain the evolution of the different orbital elements. In this paper, we will focus exclusively on the dissipative effects of the self-force which only affect to the principal orbital elements (either or ), i.e.  and we ignore the first term in Eqs. (8) and (9). We will refer to the principal orbital elements simply as orbital elements.

Iii The New Kludge Scheme

In this section we present all the details of the new kludge scheme. As explained in the previous sections, we separate the problem into two parts: that of constructing the SCO’s trajectory and that of building the waveform from these trajectories. We begin then with a description of the background MBH geometry and its associated geodesics. We then show how to enhance the geodesic equation system to allow for the variation of the constants of the motion. The latter require knowledge of the self-acceleration, which we calculate in a multipolar, post-Minkowskian expansion. Finally, we present an explicit transformation between Boyer-Lindquist coordinates (used to evolve the modified geodesic system) and harmonic coordinates, needed to generate waveforms in a multipolar decomposition.

iii.1 MBH Geometry and Geodesic Motion

The geometry of the MBH is modeled by the Kerr metric Kerr (1963), a vacuum stationary and axisymmetric spacetime that describes the final state of gravitational collapse, according to the BH no-hair conjecture Ruffini and Wheeler (1971) and uniqueness theorems (see, e.g. Chandrasekhar (1992)). In Boyer-Lindquist coordinates Boyer and Lindquist (1967), , the line element corresponding to the Kerr metric, , is given by


where and , with . For convenience, we also define the quantity . The function has two roots:


The root () coincides with the location of the event horizon.

The Kerr geometry is stationary, as described by the timelike Killing vector field , and axisymmetric, as described by a spacelike Killing vector field . It also well-known that the Kerr geometry has an additional symmetry described by a 2-rank Killing tensor, , given by


where and are the two null principal directions of the Kerr geometry


When the effect of the self-force is neglected (or equivalently, in the limit of zero mass for the SCO), the SCO follows geodesics orbits of the MBH geometry. The Boyer-Lindquist coordinates of a timelike geodesic can be parameterized in terms of proper time as . For a given geodesic, we can construct three geodesic constants of motion, corresponding to each of the three symmetries of the Kerr spacetime. Stationarity leads to a conserved energy, , or equivalently an energy per unit mass


Axial symmetry leads to a conserved component of the angular momentum vector (the one along the spin axis, which we choose to be the axis)


The Killing tensor symmetry (12) leads to a conserved Carter constant, which can be defined as follows


but also we will use the alternative definition


The existence of these three symmetries makes the geodesic equations for integrable and completely separable Carter (1968); *Carter:1968ks. The separation can be carried out using the definitions of these constants of motion and also the normalization condition of the four-velocity, . The separated equations of motion for the components of obey the following set of ordinary differential equations (see Bardeen et al. (1972) for a detailed analysis of the physical properties of the solutions of these equations):


We are interested here in timelike bound and stable geodesics, what we call orbits. These orbits, apart from being characterized by the three constants of motion , can also be characterized by the orbital elements , which can be defined in terms of the turning points of the radial and polar motion (see Appendix E for more details). The turning points for the radial motion are just the minimum and maximum of , also known as pericenter and apocenter, and , which can be used to introduce the concepts of semilatus rectum and eccentricity


or equivalently


The turning point for the polar motion is just the minimum of , , which determines the interval in which oscillates, i.e.  and it can be used to introduce the concept of inclination angle, , as


Another common definition of orbital inclination angle uses the constants of motion


Alternatively, the orbits can also be characterized in terms of three fundamental frequencies (see e.g. Schmidt (2002); Drasco and Hughes (2004); Fujita and Hikida (2009)) with respect to the Boyer-Lindquist coordinate time (they can be also constructed using proper time or any other time): , associated with the radial motion (from periapsis to apoapsis and back); , associated with polar motion; and , associated with azimuthal motion. These frequencies are important because precessional orbital effects are due to mismatches between them and because they can be used to decompose, among other things, the GW in a Fourier expansion Drasco and Hughes (2004).

Expressions for these frequencies in terms of quadratures have been obtained for Kerr in Schmidt (2002), and recently also in Drasco and Hughes (2004); Fujita and Hikida (2009). Appendix F provides explicit formulas for the fundamental frequencies and periods used in this paper.

A final consideration regarding geodesic motion that is going to be important in this paper is the splitting of the four-velocity into spatial and time components. This splitting is associated with the time variable used to evolve the geodesic equations, in this paper Boyer-Lindquist time . The four-velocity of the SCO, , choosing the Boyer-Lindquist time parameterization , can then be decomposed as


where and are given by


Using the normalization condition , the factor , the GR generalization of the special-relativistic Lorentz factor, can be written in terms of the metric and the components of the velocity as follows:


iii.2 New kludge Osculating Trajectories

In this work we consider a version of our scheme in which we only include the dissipative effects of the self-force, i.e. those that only affect to principal orbital elements. In Sec. VI we discuss how to introduce conservative pieces of the self-force in the new kludge scheme. The time scale of change of , the radiation-reaction time scale , is much bigger than the orbital time scales , such that the ratio of these time scales satisfies: .

Following the method of osculating orbits described in Sec. II, we can describe the orbital evolution as given locally in time by the geodesic Eqs. (18)-(21), where the constants of motion are promoted to time-dependent quantities:


Although here we parameterize the time dependence of these quantities in terms of proper time, in practice we will use coordinate time, that corresponds to time associated with distant observers. For our discussion, we denote the solution of the evolution Eqs. (18)-(21) with the modifications of Eq. (29) by , but it is clear that it would no longer be a solution of the original geodesic equations. We do not decompose this solution, which takes into account self-force effects, into a background geodesic plus a deviation Gralla and Wald (2008), as the deviations grow secularly in time and after a certain number of cycles it cannot be considered a small deviation of the background geodesic orbit. Instead, in the spirit of the osculating orbits method, we treat as a new, self-consistent trajectory that is continuously corrected away from geodesic motion during evolution. In this sense, we are constructing an orbit that is made out of geodesic patches corresponding to different constants of motion. The transition from one patch to another is given by the (multipolar, post-Minkowskian) self-force, lacking a more accurate prescription. The length (duration) of the geodesic patches, or equivalently, the frequency at which the constants of motion are updated, is a free parameter that we can choose in the new kludge numerical implementation.

The evolution equations for can be obtained from the equations of osculating orbits, Eqs. (8) and (9). In our case, the inversion of these equations to find can be done easily by using the symmetries of the Kerr geometry. By applying , , and to Eq. (9) in combination with Eqs. (14)-(17) we obtain the evolution equations for , , and respectively.


The evolution of these quantities with respect to coordinate time introduces a factor in the above equations, due to the relation: . The most important quantity here is the self-acceleration , which as we shall see scales [see Eq. (57)] and is given by the metric perturbations . This is probably the main ingredient of self-force descriptions of EMRIs, which is described in more detail in the next section.

For certain special orbits, the rate of change of the orbital elements satisfies special relations due to the symmetries. One of these are equatorial orbits, i.e. orbits with . We can see that from Eq. (20) this implies . Therefore, equatorial orbits are characterized by and . Equation (33) allows us to write , where


The vector vanishes for equatorial geodesics, which implies that the Carter constant is always zero and can be obtained directly from a combination of and . The self-force then maps equatorial geodesics to equatorial geodesics, i.e. the equatorial character of geodesics is preserved upon self-force evolution.

Another type of orbits with extra symmetries are equatorial and circular orbits. These orbits have a helical symmetry described by an approximate Killing vector (which is exact on the geodesic intervals). This symmetry can be derived from the fact that for these special orbits [see Eq. (21)]. The angular velocity can be written as (see Bardeen et al. (1972)):


where the upper/lower sign corresponds to prograde/retrograde circular orbits (i.e. that corotate/counterrotate with the MBH spin and have /), , and is the Boyer-Lindquist radial coordinate of the circular orbit. Then, the helical Killing vector is . The associated constant of motion is then . The evolution of is and it can be shown that for locally circular-equatorial orbits. Hence, the evolution of the angular momentum in the spin direction is related to the evolution of the energy by .

Finally, let us consider circular but nonequatorial orbits, i.e. orbits with but . There is no helical symmetry for these orbits due to the fact that the MBH spin is not aligned with the orbital angular momentum. Nevertheless, as shown in Kennefick and Ori (1996); Ryan (1996) the radiation-reaction evolution of these orbits preserves their circular character (see also Hughes (2000)). These orbits are then characterized by the vanishing of the right-hand side of Eq (19), (the radial coordinate does not change), and its first radial derivative, (the orbit is always at a turning point of the radial motion). In addition, the condition has to be satisfied for the circular orbit to be stable. Following Hughes (2000), the preservation of the circularity of the orbit along the inspiral translates into the following two conditions: and . From these two conditions we can express the time evolution of the radius of the orbit and the evolution of the Carter constant in terms of the evolution of the energy and angular momentum of the orbit:


where the coefficients () and are functions of and are given in Appendix G. Therefore, the evolution of and can be obtained by evaluating and from the multipolar, post-Minkowskian self-force.

To finish this section, and as a consistency check, we take the Newtonian limit of the rate of change of the orbital elements in Eqs. (30)-(32). In this limit, we find that to leading order


which agrees with the Newtonian expressions of Flanagan and Hinderer (2007) after transforming to Cartesian coordinates. Notice, however, that Eqs. (30)-(32) contain many more terms than in Flanagan and Hinderer (2007), due to the fact that in the latter the factor is expanded in the low-speed and weak-field limit.

iii.3 Multipolar Post-Minkowskian Self-Acceleration

Here we discuss a method to obtain an approximation for the self-force (5). The most rigorous approach would be to compute this force within the framework of BH perturbation theory, but as already argued, such a task is still under development and in computational terms would be very costly. A different approach is to extract the self-force from the PN equations of motion. Such a path was taken by Pound and Poisson Pound and Poisson (2008) for a Schwarzschild background and Gair, et al. for a Kerr background Gair et al. (2011). This approach, however, is built under a PN approximation, which is not ideal for EMRI modeling as the SCO can reach relativistic velocities, with significant relativistic factor, as it orbits close to the MBH.

Instead of either of these approaches, we here approximate the self-force via a multipolar, post-Minkowskian expansion (see e.g. Blanchet and Damour (1984); Iyer and Will (1993, 1995)):


where and are time-symmetric potentials, in the sense that they are made out of the half-sum of retarded and advanced integrals of the stress-energy tensor over the source, i.e. the half-sum of the retarded and advanced Green functions associated with the perturbative equations. As a consequence, these potentials are invariant under time inversion. The quantities and are time-asymmetric radiation-reaction potentials, constructed from the half-difference of retarded and advanced waves, and hence, odd under time inversion. These potentials are given by Blanchet and Damour (1984); Blanchet (1997)


where we recall that are spacetime harmonic coordinates, is the antisymmetric Levi-Civita symbol, and


is the symmetric trace-free (STF) projection of the multi-index quantity . The first term in [Eq. (44)] corresponds to the well-known Burke-Thorne radiation-reaction potential Burke (1971):


The quantities , , and are the th-time-derivative of the STF mass quadrupole, mass octopole and current quadrupole moments. Formally, the radiation-reaction potentials depend on the integral of certain derivatives of the asymmetric sum (half-difference of retarded and advanced waves) of multipole moments (see, e.g. Eqs.  in Blanchet (1997)). Equation (45), however, is obtained by expanding these integrals in a slow-velocity approximation, after which the arguments of the radiation-reaction potentials depend on time only. This is consistent with the fact that the radiation-reaction force is to be evaluated in the source zone of the SCO, and not in the wave zone.

Let us provide explicit expressions for these multipole moments. To lowest order, the mass moments are given by


and the current moment is


where angle-brackets are STF projections. To higher order, these moments become more complicated as there are now nonlinear contributions from the nonreactive potentials (i.e. nonlinear contributions from the background) as well as tail and memory contributions. These expressions can be found for example in Eq.  of Arun et al. (2008) and Eq.  of Blanchet et al. (2006). In BH perturbation theory language, these higher-order terms would contribute conservative and dissipative corrections to the dissipative equations of motion. We neglect these contributions in the current version of the new kludge approach, but they can be easily incorporated in future improvements of the scheme.

The metric given in Eqs. (41)-(42) is expanded in the far-field limit, a resummation of which is necessary in order to use it for self-force calculations. All terms that are independent of the radiation-reaction potentials can be identified with MBH background geometry terms, corrected perhaps by the presence of the SCO. Neglecting the latter, we can exactly resum the metric in Eqs. (41)-(43) so that it can written as follows


where is the Kerr metric in harmonic coordinates (we shall discuss the issue of coordinates in more detail in Sec. III.4) and where we have introduced the metric perturbation, , which is given by


These are the metric perturbations that are going to describe the radiation-reaction effects in our new kludge scheme, and hence they are our approximation to the regularized metric perturbations in Eqs. (4) and (5).

Regarding the MBH background geometry, it is useful for the purposes of this work to introduce some scalar, vector, and tensor potentials in harmonic coordinates, which contain information equivalent to the potentials and in Eqs. (41)-(43):


or more compactly and . The expressions from these local potentials can be derived directly from the expressions of the components of the Kerr metric in harmonic coordinates given in Appendix B.2.

We are now in a position to compute the self-acceleration in terms of the radiation-reaction potentials of Eqs. (44) and (45) through Eq. (53) and Eq. (5). In what follows we develop the expression for the acceleration in order to describe its structure and the role of the different terms that appear on it. The first step is to use the decomposition of the SCO four-velocity given in Eqs. (26) and (27) in such a way as to rewrite Eq. (5) in the following form:


where we have factored out the dependence on the relativistic factor [see Eq. (105) for the expression of in terms of the potentials and the spatial velocity ] is the projector orthogonal to the SCO four-velocity, (in terms of the potentials it is given in Eq. (104)). The two pieces of the self-acceleration in Eq. (57) contain different terms. The first piece, , only contains gradients of the radiation-reaction potentials and the spatial velocity , whereas the second piece, , contains explicit couplings between the MBH potentials (actually their derivatives) and the radiation-reaction potentials. The form of these two terms is the following:






where the connection is to be computed with the Kerr background metric in harmonic coordinates. Putting all these different ingredients together and separating space and time components, we arrive at the following structure for the piece :


where the expressions for the quantities and , which depend only on the three-velocity and spacetime derivatives of the radiation-reaction potentials and , are given in Eqs. (106) and (107) of Appendix A. The second piece of the self-acceleration, , has the following structure:


where the quantities , , ,