A Mapping between Spheroidal and Spherical Harmonics

Extreme Mass-Ratio Inspirals in the Effective-One-Body Approach: Quasi-Circular, Equatorial Orbits around a Spinning Black Hole


We construct effective-one-body waveform models suitable for data analysis with LISA for extreme-mass ratio inspirals in quasi-circular, equatorial orbits about a spinning supermassive black hole. The accuracy of our model is established through comparisons against frequency-domain, Teukolsky-based waveforms in the radiative approximation. The calibration of eight high-order post-Newtonian parameters in the energy flux suffices to obtain a phase and fractional amplitude agreement of better than radian and respectively over a period between and months depending on the system considered. This agreement translates into matches higher than over a period between and months, depending on the system. Better agreements can be obtained if a larger number of calibration parameters are included. Higher-order mass ratio terms in the effective-one-body Hamiltonian and radiation-reaction introduce phase corrections of at most radians in a one year evolution. These corrections are usually one order of magnitude larger than those introduced by the spin of the small object in a one year evolution. These results suggest that the effective-one-body approach for extreme mass ratio inspirals is a good compromise between accuracy and computational price for LISA data analysis purposes.

I Introduction

Extreme mass-ratio inspirals (EMRIs) are one of the most promising sources of gravitational waves (GWs) expected to be detected with the proposed Laser Interferometer Space Antenna (LISA)  Bender et al. (1998); Danzmann and Rüdiger (2003); Prince (2003); Gair et al. (2004). These sources consist of a small compact object, such as a neutron star or stellar-mass black hole (BH), in a close orbit around a spinning, supermassive BH Amaro-Seoane et al. (2007). Gravitational radiation losses cause the small object to spiral closer to the supermassive BH and eventually merge with it. Hence, the GW signal from such events encodes information about strong gravity, allowing tests of general relativity Sopuerta and Yunes (2009) and of the Kerr metric Collins and Hughes (2004); Vigeland and Hughes (2010); Hughes (2006); Apostolatos et al. (2009); Lukes-Gerakopoulos et al. (2010); Glampedakis and Babak (2006); Gair et al. (2008); Barausse and Rezzolla (2008); Barausse et al. (2007); Kesden et al. (2005); Barack and Cutler (2007), as well as measurements of the spins and masses of massive BHs Barack and Cutler (2004).

Unfortunately, EMRIs are very weak sources of GWs at their expected distances from us, and thus, they must be observed over many cycles to be detectable Amaro-Seoane et al. (2007). For example, a typical EMRI at a distance of would produce GWs with signal-to-noise ratios (SNRs) on the order of - depending on the observation time. Therefore, matched filtering is essential to extract EMRIs from LISA noise and the foreground of unresolved GWs from white dwarf binaries in our galaxy.

Matched filtering consists of cross-correlating the data stream with a certain noise-weighted waveform template Andrzej Kr—lak (2005). If the latter is similar to a GW event hidden in the data, then this cross-correlation filters it out of the noise. Of course, for matched filtering to be effective, one must construct accurate template filters. Otherwise, real events can be missed, or if an event is detected, parameter estimation can be strongly biased Cutler and Vallisneri (2007). The construction of accurate EMRI waveforms is extremely difficult due to the long duration of the signal and the strong-field nature of the orbits. A one-year EMRI signal contains millions of radians in phase information. To avoid significant dephasing, its waveform modeling must be accurate to at least one part in  Lindblom et al. (2008).

Such an exquisite accuracy requirement is complicated further by the strong-field nature of the orbit. An EMRI can reach orbital velocities of two-thirds the speed of light and orbital separations as small as a few times the mass of the supermassive companion. This automatically implies that standard, post-Newtonian (PN) Taylor-expanded waveforms fail to model such EMRI orbits Mandel and Gair (2009). PN theory relies on the assumptions that all orbital velocities are much smaller than the speed of light and that all objects are at separations much larger than the total mass of the system Blanchet (2006). A better approximation scheme to model EMRIs is BH perturbation theory, where one only assumes that the mass ratio of the system is much less than unity Mino et al. (1997). This is clearly the case for EMRIs, as the mass ratio is in the range . Perturbation theory, however, is computationally and analytically expensive. Only recently have generic orbits been computed around a non-spinning BH to linear order in the mass ratio Barack and Sago (2009); Sago (2009), and it is unlikely that these will be directly used for EMRI data analysis Gair et al. (2004).

EMRIs involve complicated inspiral analysis, but unlike comparable-mass coalescences, the merger and ringdown phase can be completely neglected. To see this, note that the instantaneous amplitude of the waves from a binary scales as , where is the reduced mass, is the total mass and are the component masses. The inspiral lasts for a time and releases an energy flux . In contrast, the merger and ringdown last for a time (independent of ), and thus, release an energy flux . For an EMRI, and the inspiral clearly dominates the signal. Based on this argument, we neglect the merger and ringdown, focusing on the inspiral for our analysis.

i.1 Summary of Previous Work

The modeling of EMRIs has been attempted in the past with various degrees of success. One approach is to compute the self-field of the test particle to understand how it modifies the orbital trajectory. This task, however, is quite involved, both theoretically and computationally, as the self-field contains a divergent piece that is difficult to regularize (see, e.g. Ref. Barack (2009) for a recent review). Recently, a breakthrough was achieved, with the full calculation of the self-force for generic EMRIs about non-spinning supermassive BHs Barack and Sago (2009); Sago (2009). Such calculations, however, are computationally prohibitive if the goal is to populate a waveform template space.

Another approach is to use more approximate methods to model the EMRI trajectories. One such approach was developed by Hughes Hughes (2000, 2001a), following the pioneering work of Poisson Poisson (1993). In this radiative-adiabatic scheme, the inspiral is treated as a sequence of adiabatically shrinking geodesics. The degree of shrinkage is determined by solving the Teukolsky equation on each individual geodesic. Its solution encodes how the constants of the motion (the energy, angular momentum and Carter constant) change due to GW emission. By interpolating across such sequence of geodesics, one then obtains a continuous inspiral and waveform. The calculation of a single waveform, however, is rather computationally expensive, as it requires the mapping of the entire orbital phase space, which for generic orbits is likely to be prohibitive. It is also worth noting that the radiative approximation neglects the impact of conservative effects which, especially for eccentric orbits, are likely to be important Pound et al. (2005).

Other, perhaps more rough approximations can also be used to model EMRIs. The templates obtained through these methods are sometimes called kludge waveforms to emphasize their approximate nature. The goal of their construction was never to provide sufficiently accurate templates for real data analysis. Instead, kludge waveforms were built to carry out descoping or parameter estimation studies to determine roughly the accuracy to which parameters could be extracted, given an EMRI detection with LISA.

The first kludge waveforms were constructed by Barack and Cutler Barack and Cutler (2004). These waveforms employ the quadrupole formula to build templates as a function of the orbital trajectories. The latter are simply Keplerian ellipses with varying orbital elements. The variation of these is determined by low-order PN expressions, constructed from the GW energy and angular momentum fluxes. An improvement of these fluxes was developed by Gair and Glampedakis Gair and Glampedakis (2006), who fitted these low-order PN expression to more accurate fluxes constructed from solutions to the Teukolsky equation. A further improvement was developed by Babak et al. Babak et al. (2007), who modeled the waveforms via a quadrupole-octopole formula and the orbital trajectories via solutions to the geodesic equations, augmented with PN–orbit-averaged evolution equations for the orbital elements.

All of these improvements, however, do not mean that kludge waveforms would be effectual or faithful for realistic data analysis with LISA. One cannot exactly quantify this statement because exact EMRI waveforms are not available and will not be in the near future. One can nonetheless predict that these approaches will be insufficient because critical components of the fluxes are not being taken into account. For example, GWs do not only escape to infinity, but they are also absorbed by the supermassive BH, contributing to the overall fluxes of energy and angular momentum. This contribution is non-negligible if one considers sufficiently long waveforms (longer than a few weeks). In fact, as we shall show in this paper, even the inclusion of such terms and very high order PN expressions in the fluxes is still insufficient for accurate waveform models that last more than a couple of months.

i.2 The Effective-One-Body Approach

The effective-one-body (EOB) formalism was introduced in Refs. Buonanno and Damour (1999, 2000) to model the inspiral, merger, and ringdown of comparable-mass BH binaries. This scheme was then extended to higher PN orders Damour et al. (2000), spinning BHs Damour (2001); Buonanno et al. (2006); Damour et al. (2008); Barausse and Buonanno (2010), small mass-ratio mergers Nagar et al. (2007); Damour and Nagar (2007); Bernuzzi and Nagar (2010), and improved by resumming the radiation reaction-force and waveforms  Damour et al. (1998); Damour and Nagar (2007); Damour et al. (2009); Pan et al. (2010a); Fujita and Iyer (2010). In the comparable mass case, phase and amplitude agreement was achieved between EOB and numerical-relativity waveforms, after calibrating a few parameters Damour and Nagar (2009); Buonanno et al. (2009); Pan et al. (2010b). By calibrating the EOB model to the comparable mass case, one can also improve the agreement of the model with the self-force predictions  Barack and Sago (2009); Damour (2010). The combination of EOB and BH perturbation theory tools for LISA data-analysis purposes was first carried out in Refs. Yunes et al. (2010); Amaro-Seoane et al. (2010). In these papers, the EOB scheme was found successful for the coherent modeling of EMRIs about a non-spinning background for a year period. Here we extend these results to non-precessing EMRIs about a spinning background.

As a first step toward the construction of accurate EMRI waveforms, we concentrate on quasi-circular, equatorial EMRIs about a spinning, supermassive Kerr BHs. The modeling of such EMRIs is simpler than that of inclined and eccentric ones, as only a single component of the radiation-reaction force is non-vanishing and entirely controlled by the GW energy flux (the Carter constant vanishes by symmetry). Moreover, such EMRIs are expected in at least one astrophysical scenario Levin (2007). In this setup, stellar-mass compact objects are either created in the accretion disk surrounding the supermassive BH or are captured by the disk, and hence move with the disk. The accretion disk is expected to be in the spin equatorial plane within a few hundred gravitational radii of the supermassive BH Bardeen and Petterson (1975).

We first calibrate the EOB energy flux to the leading-order energy flux computed in BH perturbation theory through the solution to the Teukolsky equation. This calibration is more complicated than for non-spinning systems because it must now be performed globally, i.e., as a function of both spin and velocity. This increases the computational cost of the calibration and the number of calibration parameters, as a bivariate series generically contains more terms than a monovariate one. After calibrating parameters, we find that the fluxes agree to within one part in for all spins [] and velocities [] considered, where is the velocity at the innermost circular orbit (ISCO).

Once the energy flux has been calibrated, we evolve the Hamilton equations in the adiabatic approximation and compare the amplitude and phase evolution to that obtained with an approximate BH perturbation theory, numerical result. For the latter, we employ the so-called radiative approximation Hughes (2000, 2001a), where one models the EMRI as an adiabatic sequence of geodesics with varying orbital elements, as prescribed by the solution to the Teukolsky equation. We find that the EOB and Teukolsky-based waveforms agree in phase and relative amplitude to better than radian and respectively after or months of evolution, depending on the system considered. Better agreements can be obtained if a larger number of calibration parameters were included.

Our EOB waveforms differ from previous kludge models on several fronts. First, the radiation-reaction force is here computed differently than in the kludge approach. In the latter this force is calculated from PN, Taylor-expanded fluxes that encode the GW that escape to infinity only. These fluxes were then improved by fitting a very large number of parameters to more accurate Teukolsky-fluxes with a log-independent, power-series expansion for the fitting functions Gair and Glampedakis (2006). In the EOB approach, the radiation-reaction force is computed directly from the factorized resummed waveforms Damour et al. (2009); Pan et al. (2010a). These are enhanced through the addition of BH absorption terms and then the calibration of eight high PN-order parameters to Teukolsky fluxes with a log-dependent, power-series expansion for the fitting functions. Second, the conservative dynamics are also treated here differently than in the kludge approach. In the latter, the Hamiltonian is either a two-body, Newtonian one Barack and Cutler (2004) or the full test-particle limit one, i.e., Schwarzschild or Kerr Babak et al. (2007). In the EOB approach, the conservative dynamics not only encodes the exact test-particle limit Hamiltonian, but they also allow for the inclusion of finite mass-ratio terms and of the spin of the small body.

i.3 Data Analysis Implications

The waveforms computed here are thus suitable for coherent data analysis over periods of several months. This can be established by computing the overlap between the EOB and Teukolsky-based waveforms, after maximizing over extrinsic parameters (an overall phase and time shift). We find that, when eight calibration parameters are used, the overlap remains higher than over to months of evolution, depending on the system considered. This is to be compared with numerical kludge waveforms Babak et al. (2007) whose overlap drops to and after and months respectively, even when forty-five calibration parameters are used to fit the flux Gair and Glampedakis (2006) . Of course, one could obtain higher overlaps by maximizing over intrinsic parameters, such as the chirp mass or the spin of the background, but this would naturally bias parameter estimation. Also, when integrating over only two weeks, the overlap increases, remaining higher than at Gpc regardless of the model used.

The benefit of coherently integrating over longer periods of time is that the recovered SNR naturally increases, thus allowing us to detect signals farther out and improving parameter estimation. One can see this by simply noting that the SNR scales with the square root of the time of observation. For example, coherent integration over or months instead of two weeks increases the SNR at Gpc from to and from to for two prototypical EMRIs. Such a large increase in SNR by coherently integrating over long observation times brings EMRIs not only to a confidently detectable range, but would also allow interesting tests of GR.

We conclude the paper by studying the error introduced in these waveforms due to neglecting second-order mass-ratio terms in the radiation-reaction force (dissipative PN self-force) and first-order in the mass-ratio terms in the Hamiltonian (conservative PN self-force). Such mass-ratio dependent effects can easily be included in the EOB prescription, as they are known in the PN/EOB framework. Of course, since these are known to finite PN order, we cannot include full second-order effects. These effects should be considered estimates, since the complete result may differ from the PN prediction. We find that such PN radiation-reaction effects modify the phase of the waveform by radians in a one year evolution, provided the EMRI samples the strong-gravity regime close to the ISCO. In a two-month period, however, the inclusion of finite mass-ratio effects increases the mismatch from to at 1 Gpc. This implies that such effects will only be seen if one coherently integrates over a sufficiently long time of observation. We find a somewhat smaller final dephasing when we allow the second body to be spinning and neglect any self-force corrections. The relative importance of the conservative or dissipative PN self-force terms and that of the spin of the second body depends on somewhat on the EMRI considered. Generically, we find all such corrections to be larger than one radian after a full year of integration, while they are negligible over a two month period.

This paper is organized as follows. Section II describes how we model EMRIs analytically and numerically. Section III discusses how the analytical EOB model is calibrated to the Teukolsky energy flux. Section IV compares EOB evolutions to Teukolsky ones, while Sec. V discusses the data analysis implications of such a comparison. Section VI estimates the effect of mass-ratio dependent effects and Sec. VII concludes and points to future research. Appendix A presents details on the transformation between spheroidal and spherical tensor harmonics. Appendix B contains expressions for the GW energy flux absorbed by BHs. Finally, in Appendix C we write the EOB Hamiltonian derived in Ref. Barausse and Buonanno (2010) when BHs carry spins aligned or anti-aligned with the orbital angular momentum. We use geometric units, with , unless otherwise noted.

Ii EMRI modeling

ii.1 Analytical modeling: EOB-based waveforms

Consider a BH binary system with masses and , total mass , reduced mass and symmetric mass ratio . We assume that the orbital angular momentum is co-aligned or counter-aligned with the individual BH spins , where denotes the th BH’s spin parameter and denotes the dimensionless spin parameter.

We first discuss the case of a non-spinning BH () with mass orbiting a spinning BH with spin parameter and mass , to leading order in the mass ratio . Subleading terms in the mass ratio and terms proportional to introduce conservative corrections that are not included in the Teukolsky waveforms, which we shall use to calibrate our model, and we therefore neglect them during the calibration. Eventually, however, we shall turn these conservative terms on and estimate their effect using the spin EOB Hamiltonian of Ref. Barausse and Buonanno (2010) (see Appendix C).

In the EOB framework, the orbital trajectories are obtained by solving Hamilton’s equations, supplemented by a radiation-reaction force describing the backreaction of GW emission on the orbital dynamics. Neglecting conservative corrections of order and , the spin EOB Hamiltonian reduces to the Hamiltonian of a non-spinning test-particle in Kerr, :




being the Kerr metric. In Boyer-Lindquist coordinates and restricting ourselves to the equatorial plane , the relevant metric components read


where , and the metric potentials are


Although the EOB formalism includes possible non-adiabaticities in the last stages of the inspiral and plunge, it is not necessary to include non-adiabatic effects here. Generically, for the systems that we consider, we find that the inclusion of non-adiabatic corrections leads to small phase corrections (of after one year of evolution) Yunes et al. (2010); Amaro-Seoane et al. (2010). The assumption of adiabaticity allows us to simplify the evolution equations that are solved numerically. This in turn reduces the computational cost of producing EOB waveforms: an adiabatic EOB evolution requires a few CPU seconds, while a non-adiabatic one would require CPU days or weeks. The non-adiabatic model is computationally more expensive because one needs to solve all of Hamilton’s equations, with radiation reaction source terms that are expensive to evaluate.

The Hamiltonian of Eq. (2) simplifies drastically when we consider circular, equatorial orbits () with co-aligned or counter-aligned with the orbital angular momentum (see eg. Bardeen et al. (1972)). Imposing , which is a necessary condition for circular orbits, and inserting Eqs. (6a)–(6d) in Eq. (2), a straightforward calculation returns




and is the conjugate momentum to the Boyer-Lindquist coordinate or simply the orbital angular momentum. Imposing the condition , which is also satisfied by circular orbits, we can solve for as a function of and  Bardeen et al. (1972)


where corresponds to prograde or retrograde orbits, respectively. Inserting the above equation in Eq. (9) yields an expression for the energy of circular orbits in Kerr Bardeen et al. (1972)


The above quantities can also be expressed in terms of the orbital velocity , once is derived for circular orbits. Computing , using Eq. (11), we obtain


We also define the parameter .

In the adiabatic approximation, the orbital evolution is fully determined by the frequency evolution through Eq. (13). Assuming the motion follows an adiabatic sequence of quasi-circular orbits, we can use the balance equation to derive


where is the GW energy flux (see e.g. Ref. Buonanno and Damour (2000)). The multipolar factorized form of this flux, proposed in the non-spinning case in Refs. Damour and Nagar (2007); Damour et al. (2009) and extended to the spin case in Ref. Pan et al. (2010a), is given by


which under the assumption of adiabaticity reduces to




where denotes the parity of the multipolar waveform (i.e., if is even, if is odd), and


When spin effects are present, the expressions for all the terms in Eq. (17), namely , , and can be read in Ref. Pan et al. (2010a) [see Eqs. (24), (25), (26) and (29) therein]. The functions are the standard spherical harmonics, while and are numerical coefficients that depend on the mass ratio (see Eqs. - in Ref. Damour et al. (2009)). As before, we work to leading order in initially, and later study how the terms of higher-order in affect the GW phase evolution.

The solution to Eq. (14) requires that we prescribe initial data. We here choose post-circular initial conditions, as described in Ref. Buonanno and Damour (2000), to set-up a mock evolution that starts at a separation of and ends at either the ISCO or whenever the GW frequency reaches . This mock evolution is then used to read initial data one-year before the end of the mock evolution. This approach leads to an accurate initial data prescription, without any eccentricity contamination. For example, the error in the initial frequency induced by starting the mock evolution at , instead of , is on the order of , which leads to a difference in accumulated GW cycles of rads after a one year evolution.

Finite mass-ratio corrections can be incorporated into the EOB model by including subleading terms of and in the Hamiltonian, angular momentum and relation. We shall first ignore such terms to compare against Teukolsky-based waveforms. In Sec. VI, we shall study how our results change when we include such terms. To do so, we shall still assume circular, equatorial orbits and an adiabatic evolution, but employ the spin EOB Hamiltonian of Refs. Barausse et al. (2009); Barausse and Buonanno (2010) (reviewed in Appendix C), instead of the Kerr Hamiltonian of Eq. (9).

Except for this change, the EOB waveform modeling with finite mass-ratio corrections follows closely the derivation presented above. First, we compute the angular momentum associated with the Hamiltonian of Eq. (101) for circular, equatorial orbits, imposing and solving for . Then, we derive the orbital frequency to relate to , and to express in terms of . When mass-ratio corrections are present, however, the Hamiltonian becomes much more involved, so solutions for as a function of must be searched numerically. We have checked that the discretization and interpolation used to solve these equations numerically do not introduce an error larger than in the Hamiltonian (101).

ii.2 Numerical modeling: Teukolsky-based waveforms

Teukolsky-based waveforms is the name we give to radiative models that use the Teukolsky equation to prescribe the radiation-reaction force for an inspiral. We use BH perturbation theory, considering a background spacetime with mass and spin (recall that the spin parameter . The inspiraling object is a test-body with mass and no spin (). In principle, the masses and spins used here should be the same as those introduced in the EOB model.

The radiative approximation assumes that EMRIs can be modeled as an adiabatic sequence of geodesics with slowly-varying constants of the motion. Consider the discretization of the orbital phase space, each point of which represents a certain geodesic with a given set of constants of motion (energy , angular momentum and Carter constant ). For quasi-circular, equatorial EMRIs, the Carter constant vanishes, while the variation of the energy and angular momentum are related via , being the orbital frequency. At each point in the orbital phase space, the geodesic equations can be solved to obtain the orbital trajectory of the small compact object, given any spin of the background.

Once the geodesic trajectories are known, one can use these to solve the linearized Einstein equations and obtain the gravitational metric perturbation. This is best accomplished by rewriting the linearized Einstein equations in terms of the Newman-Penrose curvature scalar to yield the Teukolsky equation Teukolsky (1973). One can decompose into spin-weight spheroidal harmonics , using Boyer-Lindquist coordinates , in the Fourier domain:


The radial functions satisfy the radial Teukolsky equation


where is given in Eq. (7) and the radial potential is


with , , and the spheroidal harmonic eigenvalue. The source function is given explicitly in Eq.  of Ref. Hughes (2000) and it depends on the stress-energy tensor for a test-particle in a geodesic trajectory.

The Teukolsky equation admits two asymptotic solutions: one outgoing as and one ingoing as one approaches the background’s event horizon. These two solutions represent outgoing radiation at future null infinity and ingoing radiation that falls into the BH through the event horizon. Both types of radiation are critical in the modeling of EMRIs; not including BH absorption can lead to errors in the waveform of order radians Hughes (2001b, a). These solutions can then be used to reconstruct both the GW radiated out to infinity, as well as the total energy flux lost in GWs. The energy flux can then be related to the temporal rate of change of the orbital elements, such as the orbital radius.

Solving the Teukolsky equation for a geodesic orbit tells us how that orbit tends to evolve due to the dissipative action of GW emission. By doing so for each point in orbital phase space, we endow this space with a set of vectors that indicate how the binary flows from one orbit to another. We compute these vectors at a large number of points, and use cubic spline interpolation to estimate the rates of change of orbital constants between these points. This allows us to compute the temporal evolution of all relevant quantities, including the orbital trajectories and gravitational waveforms.

We implemented this algorithm, discretizing the orbital phase space from an initial separation of to the Kerr ISCO


in a point grid, equally spaced in


We cannot evolve inside the Kerr ISCO with such a frequency-domain code, as stable orbits do not exists in this regime. We cannot evolve inside the Kerr ISCO with a frequency-domain code, as stable orbits do not exist in this regime (and so such orbits do not have a well-defined frequency spectrum). The code we use to construct our waves is based on Hughes (2000, 2001a); Drasco and Hughes (2006), updated to use the spectral methods introduced by Fujita and Tagoshi Fujita and Tagoshi (2004, 2005). A detailed presentation of this code and its results is in preparation Throwe ().

The dominant error in these Teukolsky-based waveforms is due to truncation of the sums over and . In particular, to compare with PN results, we must map this spheroidal decomposition to a spherical one (see Appendix A for more details). Such a mapping requires one to include a buffer region of modes about the largest mode computed. We have been careful to use a sufficiently wide buffer and total number of modes such that the energy fluxes are accurate to for all velocities and spins. In particular, this means that in the strong field region (close to the ISCO) up to modes were included. Other sources of error are due to the intrinsic double precision in the numerical solution to the Teukolsky equation, the discretization of the orbital phase space, and its cubic-spline interpolation. All of these amount to errors of order . All in all and in terms of GW phase, the Teukolsky-based waveforms are accurate to at least rads. during an entire year of evolution.

Iii Calibrating the test-particle energy flux

We consider here a calibrated EOB model that is built from the functions in Eq. (17), but in which higher-order PN terms are included in the functions and are calibrated to the numerical results. In particular, we write


where are eulerlog-independent and eulerlog-dependent calibration coefficients that enter the mode at and proportional to . As in Refs. Damour et al. (2009); Pan et al. (2010a) the euler-log function is defined as


where is Euler’s constant. Notice that we have introduced calibration parameters in the non-spinning sector of the flux and in the spinning sector. The spin parameter denotes here the spin of the background. When we neglect mass-ratio terms, we choose . However, when we switch on the mass-ratio terms we have an ambiguity on the choice of . Following Ref. Pan et al. (2010a) we choose , where is the deformed-Kerr spin parameter defined in Appendix  C. Note that since , these choices are identical in the test-particle limit, which is when the energy flux is calibrated.

The choice of calibrating function in Eq. (24) is rather special and requires further discussion. We have chosen this function so that leading-order corrections in the two dominant GW modes, and , are included. Higher modes contribute significantly less to the GW and its associated flux. In each mode, we have included the leading-order unknown terms that are both -independent and -dependent. Since -independent terms in are known to much higher order (PN) than the -dependent ones (PN) in the test-particle limit Shibata et al. (1995); Tanaka et al. (1996); Mino et al. (1997), spin-independent calibration coefficients enter at a much higher PN order. The spin-dependence of the calibration terms is inferred from known terms at lower PN orders. We have investigated many functional forms for the calibrating functions, with a varying number of degrees of freedom, and found the one above to be optimal in the class studied.

The introduction of additional calibration terms might seem like a lot. In Ref. Yunes et al. (2010), the knowledge of non-spinning terms up to PN order was found to be crucial to obtain a sufficiently good agreement in the flux. Moreover, only additional calibration terms ( at 6PN order in and at 5PN order in ) were needed to reach a phase agreement of 1 rad after two years of evolution. Similarly here, we expect that if the remaining , , and PN order terms were calculated in the test-particle limit when the central black carries a spin, then the flux would also improve, requiring a smaller number of calibration parameters. It is quite likely that those higher-order PN terms will be computed in the near future, as they involve dramatically less complicated calculations than PN terms in the comparable-mass case.

Having in hand an improved GW energy flux carried away to infinity, this must be enhanced with expressions for the GW energy flux that is absorbed by the background BH. We do so here by simply adding the Taylor-expanded form of the latter (see Appendix B) to the flux of Eq. (15). The BH absorption terms in the GW energy flux depend on polygamma functions, which are computationally expensive to evaluate. We have empirically found that expanding this function in to th order is a sufficiently good approximation for our purposes. When performing computationally expensive calculations (like the fits described below) we shall employ such expansions, but when solving for the orbital phase and when computing the waveforms we shall return to the full polygamma expressions.

The resulting EOB energy flux, including BH absorption terms, is then calibrated via a two-dimensional, least-squares minimization relative to numerical data obtained from Teukolsky-based calculations. The fitting routine is two dimensional because when considering EMRIs about spinning backgrounds, the flux depends on two independent variables: the orbital velocity (or frequency or separation) and the spin of the background. This, in turn, increases the number of points that need to be used by more than an order of magnitude to properly calibrate Eq. (24). In all fits, we have assumed a data variance of for all velocities and spins and we have required a relative accuracy of one part in . Since the data is now two-dimensional, one must search for a global minimum in space. After doing so, we find the calibration parameters


The computational cost of the calibrations performed in this paper is much larger than those carried out in Ref. Yunes et al. (2010) for the following reasons. First, we consider twice as many calibration parameters as in Ref. Yunes et al. (2010), increasing the dimensionality of the fitting space. Second, global minimization routines require non-trivial algorithms that are numerically more expensive than those employed in one-dimensional minimizations. Third, the amount of data fitted increases by at least one order of magnitude, due to the intrinsic bi-dimensionality of the problem. Combining all of this, the computational cost of performing the calibration is now more than times larger than in Ref. Yunes et al. (2010). Even then, however, these fits require CPU minutes to complete. Once they have been carried out, this calculation does not need to be repeated again in the waveform modeling.

Figure 1 plots the fractional difference between the analytical GW energy flux and that computed with Teukolsky-based waveforms as a function of velocity, from an initial value of to the velocity at the ISCO, for five different spin values: . All comparisons are here normalized to the Newtonian value of the flux .

Figure 1: Fractional difference between PN and Teukolsky-based fluxes as a function of velocity for spins (top) and (bottom). The dotted curves employ the Taylor-expanded PN flux with BH absorption terms, while the dashed and solid curves use the uncalibrated and calibrated -resummed fluxes with BH absorption terms respectively.

The different curve styles differentiate between analytical models: the dotted curves use the total, uncalibrated Taylor-expansion; the dashed curves use the uncalibrated -resummed flux with BH absorption terms; the solid curves use the calibrated -resummed flux with BH absorption terms. Notice that the calibrated model does better than the other two by at least two orders of magnitude near the ISCO for all spin-values.

Several interesting conclusions can be drawn from Fig. 1. First, as obtained in Ref. Pan et al. (2010a) the uncalibrated -resummed model is better than the Taylor-expanded version of the flux, by up to nearly an order of magnitude at the ISCO for all spins. In turn, the calibrated model is better than the uncalibrated one by one to two orders of magnitude near ISCO for all spins. One could also calibrate the Taylor-expanded flux (not shown in Fig. 1), but this would not produce such good agreement in the entire space. This is clearly because the uncalibrated -resummed model is more accurate than the Taylor one, and thus the calibration terms have to do less work to improve the agreement. For the calibrated Taylor and -resummed models to become comparable in accuracy one would have to include up to at least calibration coefficients in the Taylor model.

The inclusion of BH absorption coefficients is crucial to obtain good agreement with the full Teukolsky-based flux, a result that was not obvious for the case of non-spinning EMRIs. Figure 2 plots the fractional difference between the uncalibrated EOB GW energy flux and Teukolsky-based one as a function of velocity for five different spin values: , from an initial value of at the ISCO, for five different spin values: . For these cases, we have . The solid curves use the uncalibrated EOB model including the Taylor-expanded BH absorption contributions, while the dotted curves do not. For the non-spinning case, observe that there is a very small difference (smaller than ) between adding the BH absorption terms or not.

For the spinning cases, however, this is not the case. For rapidly spinning backgrounds, adding the BH absorption terms improves the agreement by an order of magnitude. Presumably, resumming these terms in a multipolar-factorized manner would improve the agreement even more. The BH absorption terms play a much larger role in the spinning case because spin changes the PN order at which absorption enters in the energy flux. These terms enter at 4PN order for Schwarzschild black holes, but at 2.5PN order for non-zero spin. This change of order has a very large and important impact on the system’s evolution.

Figure 2: Fractional difference between -resummed and Teukolsky-based fluxes as a function of velocity for spins . The dotted curves do not include the Taylor-expanded BH absorption contributions, while the solid lines do.

The inclusion of calibration parameters to improve the agreement of PN-inspired fluxes and Teukolsky-based ones for EMRIs is certainly not new. In Ref. Gair and Glampedakis (2006), a similar, PN-inspired calibration was carried out for circular-inclined orbits (and more generic ones). Before calibration, their fluxes were Taylor-expanded to 2PN order and included only the contribution that escapes to infinity (not the BH absorption terms discussed above). Their fit was then done with Teukolsky-data produced by an older version of the code used here, which was accurate only to one part in . Moreover, the fit was done in the range [], so the fitted function loses accuracy rapidly outside this regime, particularly close to the ISCO. Inside the fitting regime, the flux was fitted to an accuracy of using calibration coefficients for an inclined, but fixed orbit. The accuracy decreases to for orbits which get closer to the ISCO.

To fairly compare the results of Ref. Gair and Glampedakis (2006) with our results which are restricted to circular, equatorial orbits, we implemented their model and re-calibrated it considering only circular, equatorial orbits. Using coefficients, we found an accuracy similar to ours at high velocities close to the ISCO, but worse at low velocities. This is because their fluxes before calibration are not as accurate as the one employed here (by including up to PN order terms and BH absorption terms), particularly at low velocities. It is important to emphasize that by calibrating parameters instead of we here obtain better flux accuracies than in Ref. Gair and Glampedakis (2006) for circular, equatorial EMRIs. We could obtain even better accuracy if we were using a larger number of calibration coefficients, e.g. using coefficients the agreement with the Teukolsky-based flux would be of .

Iv Comparison of the GW phase and amplitude

The comparison of EOB and Teukolsky evolutions requires that we choose a specific EMRI. We shall here follow Ref. Yunes et al. (2010) and choose system parameters that define two classes of EMRIs:

  • System I explores a region between orbital separations , which spans orbital velocities and GW frequencies in the range and respectively. Such an EMRI has masses and for a mass ratio of and it inspirals for rads of orbital phase depending on the spin.

  • System II explores a region between orbital separations , which spans orbital velocities and GW frequencies in the range and respectively. Such an EMRI has masses and for a mass ratio of and it inspirals for rads of orbital phase depending on the spin.

The evolution of Sys. I is stopped around an orbital separation of , because this coincides with a GW frequency of , which is close to the end of the LISA sensitivity band. The evolution of Sys. II is usually stopped at the orbital separation corresponding to the ISCO, or whenever its GWs reach a frequency of . For each of these systems, we shall investigate different background spin parameters.

Before proceeding, notice that Sys. I and II should not be compared on a one-to-one basis. One might be tempted to do so, as Sys. I resembles a weak-field EMRI, which inspirals at a larger orbital separation and with smaller orbital velocities than Sys. II, a more strong-field EMRI. Comparisons are not straightforward, however, as these systems accumulate a different total number of GW cycles. In fact, Sys. I usually accumulates almost twice as many GW cycles as Sys. II. Therefore, even though one might expect PN models of Sys. I to agree better with Teukolsky-based evolutions, this need not be the case, as this system has more time (as measured in GW cycles) to accumulate a phase and amplitude difference than Sys. II.

We compare the EOB and the Teukolsky-based waveforms after aligning them in time and phase. Such an alignment is done by minimizing the statistic in Eq. (23) of Ref. Buonanno et al. (2009), just as was done in Ref. Yunes et al. (2010). This is equivalent to maximizing the fitting factor over time and phase of coalescence in a matched filtering calculation with white noise Buonanno et al. (2009). The alignment is done in the low-frequency regime, inside the time interval , where is the GW wavelength. This quantity depends on the spin of the background, ranging from () to () for Sys. I (Sys. II). This corresponds to aligning the initial phase and frequency inside a window of length in the range months depending on the system and spin of the background. We have checked that increasing the size of the alignment window does not affect the final phase and amplitude difference; for example, for a spin of and Sys. I, increasing the alignment window by a factor of two changes the final phase difference by rads and the relative, fractional amplitude agreement by .

Figure 3: Absolute value of the dephasing (left) and relative, fractional amplitude difference (right) computed in the calibrated EOB model and the Teukolsky-based waveforms for the dominant mode. Different curves correspond to different background spin values.

Figure 3 shows the absolute value of the dephasing and relative, fractional amplitude difference of the dominant mode for both systems and a variety of background spins. For Sys. I, the calibrated EOB model maintains a 1 radian phase accuracy over at least the first 6 months for all spin values, while for Sys. II the same phase accuracy is maintained for up to only the first 2 months. The amplitude agreement is also excellent for all spin values, with better agreement for Sys. I. As found in Ref. Yunes et al. (2010) the GW phase agreement is primarily due to the correct modeling of the orbital phase, as the former tracks the latter extremely closely; we find that the difference between the orbital and GW phase over a one year evolution is less than rads.

The agreement in the phase as a function of background spin follows closely the flux agreement shown in Fig. 1. This is hard to see in Fig. 1, which is why Fig. 4 zooms into the velocity region sampled by Sys. I and plots all background spin cases for the calibrated -resummed system. Observe that the case has the best flux agreement, which explains why the phase and amplitude agreement is so good for this case in Fig. 3. Observe also that the and cases have the worst flux agreement, which also explains why they disagree the most in phase and amplitude in Fig. 3.

Figure 4: Fractional difference between the calibrated -resummed and Teukolsky-based fluxes for spins as a function of velocity. We plot here only the range of velocities sampled by Sys. I.

The accuracy of the calibrated EOB model is excellent relative to Taylor-expanded PN models. If one were to use an uncalibrated Taylor-expanded version of the flux, instead of the calibrated -resummed flux, one would find a phase and amplitude disagreement of rads [ rads] and () for Sys. I (Sys. II) after a one year-evolution for different spin values. The above results are consistent with the arguments in Ref. Mandel and Gair (2009), who concluded that PN accurate GW phase expressions could lead to phase errors around radians over the last year of inspiral. That analysis reached those conclusions by comparing to PN accurate, analytic expressions for the GW phase. Here, we are comparing full-numerical evolutions of the PN equations of motion carried out to much higher order, and, of course, we find that such conclusions depend sensitively on the type of EMRI considered and the spin of the background.

Figure 5: Absolute value of the dephasing computed for the dominant mode in different PN models and the Teukolsky-based waveforms. Different curve colors/shades correspond to different background spin values, while different curve types correspond to different PN models.

The increase in accuracy of the calibrated -resummed model is due both to the calibration and to the factorized resummation. This fact can be appreciated in Fig. 5, where we plot the absolute value of the dephasing for the dominant mode in different PN models and the Teukolsky-based waveforms. Light curves (orange in the color version) correspond to background spins of , dark curves (red in the color version) to and black curves to non-spinning backgrounds. Dotted curves use the uncalibrated Taylor flux model, dashed curves the uncalibrated -resummed model and solid curves the calibrated version. For Sys. I, there is a large gain in accuracy by switching from the uncalibrated Taylor model to the uncalibrated -resummed model, but then the calibration itself does not appear to improve the accuracy substantially. For Sys. II, on the other hand, the calibration can increase the accuracy up to almost 2 orders of magnitude, as in the case.

The agreement in the phase and amplitude is not only present in the dominant mode, but also in higher ones, as shown in Fig. 6. We here plot the absolute value of the dephasing and the relative, fractional amplitude difference between the calibrated EOB model and Teukolsky-based waveforms for the dominant mode, as well as the and modes. We have here shifted the higher modes using the best frequency shift that maximizes the agreement for the dominant mode. This agreement is simply a manifestation of the agreement in the orbital phase. In fact, we find that , with differences that are always less than rad for all systems considered.

Figure 6: Absolute value of the dephasing (left) and relative, fractional amplitude difference (right) computed in the calibrated EOB model and the Teukolsky-based waveforms. The solid curve corresponds to the dominant mode, while the dashed curve is for the mode and the dotted curve for the mode. Different curves stand for different background spins.

Higher- modes contribute significantly less to the SNR than the dominant mode. Figure 7 plots the relative fraction between the squared of the SNR computed with only the component of the waveform and that computed by summing over all modes. This figure uses data corresponding to Systems I and II, both with spin (results for other spin values are almost identical). Clearly, the mode is dominant, followed by the and modes. Because of this feature of quasi-circular inspirals, obtaining agreement for the mode implies one can recover over of the SNR.

Figure 7: Relative fraction between the squared of the SNR computed with only the mode and that computed with all modes. The lines connecting points are only meant to group points with (dotted), (dashed) and (dot-dashed).

V Data Analysis Implications

Although the phase agreement presented in the previous section is a good indicator of the validity of the EOB model, one is really interested in computing more realistic data analysis measures. In this section we compute the mismatch between the Teukolsky and the EOB model, maximized over extrinsic parameters and as a function of observation time.

Let us first introduce some basic terminology. Given any time series and , we can define the following inner-product in signal space


where the overhead tildes stand for the Fourier transform and the star stands for complex conjugation. The quantity is the spectral noise density curve, where here we follow Barack and Cutler (2004); Berti et al. (2005). Notice that we use the sky-averaged version of this noise curve here, which is larger than the non-sky-averaged version by a factor of . In particular, this means that our SNRs are smaller than those one would obtain with a non-sky-averaged noise curve by a factor of .

Given this inner-product, we can now define some useful measures. The SNR of signal is simply


while the overlap between signals and is simply


with the mismatch . The max label here is to remind us that this statistic must be maximized over a time shift and a phase shift (see eg. Appendix B of Damour et al. (1998) for a more detailed discussion).

The data analysis measures introduced above ( and ) depend on the length of the time-series, i.e. the observation time. Figure 8 plots the mismatch between the Teukolsky-based waveforms and a variety of models for both Sys. I and II and a background spin of as a function of observation time. The vertical lines correspond to observation times of weeks, months, months, months and moths, together with their associated SNRs at 1 Gpc. The mismatches are computed with different analytical models: black crosses stand for the calibrated -model with calibration coefficients; red circles to the uncalibrated - model; blue squares to the uncalibrated Taylor model; green circles to that EOB evolution using the original flux of Ref. Gair and Glampedakis (2006) which has calibration coefficients (denoted in the figure). For comparison, we also include the amount of dephasing (numbers next to data points in Fig. 8) at weeks, months, months, months and months. Observe that the calibrated -model maintains an overlap higher than over 9 and 4 months for Sys. I and II respectively. The uncalibrated -model performs comparably to the EOB model using the flux of Ref. Gair and Glampedakis (2006) which has calibration coefficients, both of which have an overlap higher than over 6 and month for Sys. I and II respectively. In the case of Sys. II the calibrated flux of Ref. Gair and Glampedakis (2006) perform better than the uncalibrated model. Also observe that the uncalibrated Taylor model is simply inadequate to model EMRIs for any observation time.

Figure 8: Mismatch between Teukolsky and EOB waveforms as a function of waveform duration for Systems I and II and a background spin of .

We then see that the use of the calibrated EOB model allows us to integrate over longer observation times, compared to a 2-week period or to other models. In turn, this allows us to recover a higher SNR that we would otherwise. The increase in SNR scales as the square root of the observation time, as expected. For example, since the calibrated EOB model is accurate over months, one would be able to coherently recover an SNR of at 1 Gpc for Sys. I, to be compared with an SNR of obtained after months of coherent integration if using for example the fluxes of Ref. Gair and Glampedakis (2006). In general, an integration over a period longer than two weeks gains us a large increase in SNR. Such gains in SNR are important because they allow us to see EMRIs farther out. Since the SNR scales as , where is the luminosity distance, even an SNR increase in a factor of increases our accessible volume by a factor of , since the later scales as .

Vi Higher-Order Effects

Let us now discuss how finite mass-ratio higher-order effects affect the GW phase and amplitude. Those effects are encoded in the terms present in the radiation-reaction force and in the Hamiltonian. The former are second-order effects in the dissipative dynamics, while the latter are first-order effects in the conservative dynamics. We have analytic control over the PN version of such effects within the EOB formalism, but until now, we had set both of these to zero when comparing to Teukolsky evolutions, as the latter do not account for such effects. In the EOB model, however, it is straightforward to include such terms, as PN expansion are formally known for all mass ratios at some given order in . The inclusion of conservative terms is achieved by using the spin EOB Hamiltonian of Ref. Barausse and Buonanno (2010) reviewed in Appendix C. The inclusion of dissipative terms is achieved by including relative terms in the multipolarly decomposed waveform and flux .

Whether such mass-ratio effects matter depends on the EMRI considered. In Ref. Yunes et al. (2010), it was found that such effects increase the dephasing between EOB and Teukolsky models by one order of magnitude after a two year evolution for non-spinning EMRIs, a result consistent with Ref. Pound and Poisson (2008). This effect is greatly amplified when considering spinning EMRIs. Table 1 compares the effect that the inclusion of relative terms in the EOB Hamiltonian and the radiation-reaction force has on the final dephasing after a one-year evolution.

System I I I II II II
No rel.  10.04 1.60 9.36 7.63 42.21 48.39
in H 30.38 0.083 18.96 4.38 40.35 47.13
in 19.99 5.31 4.49 7.24 41.50 46.55
in H and 40.32 6.83 14.08 3.98 39.64 45.29
Table 1: Absolute value of the total dephasing after a moths of evolution. The first row includes no relative contributions in the Hamiltonian or the radiation-reaction force. The second row includes relative terms in the Hamiltonian, while the third row includes such terms in the radiation reaction force.

In order to read out the effect of such finite mass-ratio terms in the phasing, one must compare rows two, three and four to the baseline given in the first row of Table 1. For example, the effect of the terms in the Hamiltonian are such as to increase the dephasing by radians. Observe that the conservative and dissipative terms usually push the dephasing in different directions, partially canceling out when both of them are present. Even then though, the generic effect of high-order terms is to increase the rate of dephasing by several tens of radians are a one year evolution. Notice furthermore that the magnitude of the effect is here not very large because we are considering circular equatorial orbits.

Finite mass-ratio effects are clearly suppressed when dealing with circular orbits. This is because, for such orbits outside the ISCO, the effect of the conservative self force is simply to shift the waveform from one orbital frequency to another. Thus, from an observational standpoint, such an effect is unmeasurable. Even though the conservative self-force shifts the ISCO, this effect is still degenerate with a shift of the system’s mass parameters. This discussion, however, neglects radiation-reaction, which is crucial to model a true inspiral waveform. One can think of the radiation-reaction force as defining a trajectory through the sequence of orbital energies that an orbit follows. There is gauge invariant information in this sequence, in the sense that the mapping between energies and orbital frequencies depends on the details of the orbit at each energy level. For a given radiation-reaction force, the sequence of geodesic orbits (and hence the sequence of frequencies) depends on whether the conservative self-force is included or not.

As is clear from this discussion, the “real” (gauge-invariant) effect of the conservative self-force on quasi-circular inspiral waveforms can be rather subtle. A robust effect, however, does arise if the self-force acts on a more generic orbit, such as an eccentric one. In that case, this force will act separately on the radial and the azimuthal orbital frequencies, which can leave a potentially strong imprint in the waveform. In principle, even for an inclined circular orbit there could be a strong imprint. In practice, however, the azimuthal and polar orbital frequencies are quite similar, which suggests that perhaps, even in this case, the self-force effects will be small.

With all of this in mind, let us discuss the results presented in Table 1 in more detail. Our study suggests that the overall effect of terms in both and leads to only and additional radians of phase for non-spinning, Sys. I and II respectively. This is in fact consistent with the results presented in Ref. Yunes et al. (2010), except that there one considered -year long evolutions. One might wonder whether using the non-spinning Hamiltonian of Ref. Buonanno et al. (2009) (where the deformed-Schwarzschild potential are Padè-resummed instead of being given by Eqs. (C), (87)) has an effect on this dephasing. We have investigated this question and found that the additional contribution to the phase is () and () radians for Sys. I and II respectively over the entire year of inspiral using the 3PN (4PN) Pade form of the deformed potentials. [We notice Favata (2010) that the 4PN Padè potentials of Ref. Buonanno et al. (2009) reproduce very closely the ISCO-shift of Ref. Barack and Sago (2009).] This implies that the non-spinning Hamiltonian Buonanno et al. (2009) at 3PN and 4PN order is sufficiently close to the Hamiltonian presented in Appendix B for data analysis of non-spinning EMRIs.

One can also compare the results in Table 1 to the recent study of Huerta and Gair Huerta and Gair (2009), who investigated the effect of -corrections in the determination of parameters, given an EMRI signal. Their Table I presents the number of cycles accumulated for a variety of mass ratios. Their last column happens to correspond to our Sys. II with no spin, for which they get a total dephasing of rads and rads after the last year of inspiral when including only conservative and all second-order corrections. This is to be compared to our results: rads in and rads after the last year of inspiral when including only conservative and all second-order corrections. These numbers are in excellent agreement, allowing for differences in the modeling. Their analysis suggests that such small difference will not affect parameter estimation for EMRIs similar to Sys II. For Sys I, however, the dephasing is much larger as the mass ratio is less extreme by one order of magnitude; thus, parameter estimation might be affected in this case.

Another high-order effect that one can study is the inclusion of the small object’s spin in the evolution of the binary. Since the spin angular momentum of the small body scales with its mass, its contribution to the orbital evolution is one order in suppressed. We can model this effect by allowing in the EOB Hamiltonian of Ref. Barausse and Buonanno (2010) and letting . Doing so for Sys. I (Sys. II) with and , we find that the total dephasing now becomes ( rads), as shown in Table 2. This is to be compared to the case when , which returns a dephasing of rads ( rads). Thus, the effect of the second spin contributes roughly rads in this case.

System I I I II II II
36.34 3.96 16.03 4.41 40.41 47.04
40.32 6.83 14.08 3.98 39.64 45.29
44.31 9.70 12.14 3.56 38.88 43.54
Table 2: Absolute value of the total dephasing in rads after a moths of evolution. The first row sets the second BH’s spin to zero, while the second row sets it to . In both cases, the spin of the background is set to and in both the Hamiltonian and the radiation-reaction force.

We can also compare these results to those estimated by Barack and Cutler Barack and Cutler (2004). In their Appendix C, they estimate that for quasi-circular inspirals similar to our Sys. II, the spin of the second body should induce a dephasing of roughly rads. This is in good agreement with the results corresponding to Sys. II in our Table 2. Our results are also in good agreement with an upcoming and independent investigation of spin-effects in the PN phasing Favata ().

Finally, taking into account the results of including finite mass-ratio effects, we can conclude that unless these are precisely modeled, it is not worth requiring an agreement better than rads when calibrating the phase of the test-particle EOB waveforms against the Teukolsky-based waveforms.

Vii Conclusions

We have constructed an EOB model for EMRIs in quasi-circular, equatorial orbits about spinning backgrounds. In the test-particle limit, this model consists of adiabatically evolving a test-particle in the Kerr spacetime using the factorized energy flux of Refs. Damour et al. (2009); Pan et al. (2010a), augmented by 8 calibration coefficients. The latter are determined by comparing the factorized energy flux to a Teukolsky-based flux, built from solutions to the Teukolsky equation in the radiative approximation. In the adiabatic approximation, the EOB waveforms can be constructed in CPU seconds at a very low computational cost. When finite mass-ratio effects and the small object’s spin are included, we build the EOB model by numerically solving the Hamilton equations with the spin EOB Hamiltonian of Ref. Barack and Cutler (2004) and the Teukolsky-calibrated factorized energy flux augmented by finite mass-ratio effects Pan et al. (2010a).

For both EMRI systems considered, we find excellent phase and amplitude agreement, with dephasing less than one radian, over periods of months. The exact length of the agreement depends on how relativistic the EMRI system is. We also calculated the overlap between EOB and Teukolsky-based waveforms to find it higher than over to months, depending on the EMRI system considered.

The EOB waveforms built here have higher overlaps and better phase agreements that all currently known EMRI models for spinning, equatorial systems, while requiring much fewer calibration parameters. In particular, the EOB model with 8 calibration coefficients outperforms by almost an order of magnitude the numerical kludge waveforms with the calibrated fluxes of Ref. Gair and Glampedakis (2006) which use 45 calibration coefficients. This implies that EOB waveforms with 8 calibration coefficients can be used for longer coherent integrations, allowing us to obtain a increase in SNR. In turn, this implies that EOB waveforms can see EMRIs that are farther out, increasing the accessible volume by at least a factor of two relative to numerical kludge waveforms. Furthermore if we were using 16 calibration coefficients, we could improve the dephasing from rads ( rads) to rads ( rads) for System I (System II) after 6 months of evolution. In turn, this would decrease the mismatch from () to () for System I and II after months of evolution.

Another possible avenue for future research is the calculation of high PN order terms in the energy flux and waveforms. Our EOB model relies on the use of accurate fluxes, but for spinning systems, the flux to infinity is only known up to PN order in the test-particle limit. This is in contrast to the non-spinning terms that are known to PN order or the BH absorption terms that are known to PN order. The calculation of the spin-dependent terms in the flux to infinity to , and PN order terms in the test-particle limit is not quixotic and would be invaluable. Once these coefficients are known, then presumably the EOB waveforms would be more accurate and might require less adjustable parameters.

Of course, the EOB waveforms we constructed here are less powerful than kludge waveforms Barack and Cutler (2004); Babak et al. (2007); Gair and Glampedakis (2006) in their generality. Our waveforms cannot yet model inclined or eccentric inspirals. The inclusion of inclined orbit should be relatively straightforward, but the addition of eccentricity might require some revamping of the EOB framework. Future work in this direction would be definitely worthwhile.

Ultimately, one would like to obtain a waveform model that is sufficiently fast, efficient and accurate to do realistic EMRI data-analysis for LISA. Such a model would need to include the correct self-force contributions to the conservative dynamics, the appropriate second-order terms in the radiation-reaction force and the correct terms that describe the small object’s spin. The EOB model we developed here does agree with all known PN self-force calculation to date. However, since full self-force calculations are not yet completed, we do not have a way of assessing the error of including the currently known EOB finite mass-ratio effects. In fact, so far the only comparisons between the PN/EOB and self-force results have been concerned with the non-spinning case, and have been limited to the ISCO shift Barack and Sago (2009); Damour (2010) and other gauge invariant quantities Blanchet et al. (2009, 2010). Quite interestingly, the calibration of the EOB model to comparable-mass numerical-relativity simulations improves the agreement of the model to self-force results Damour (2010); Favata (2010). Thus, we hope that future calibrations of the spin EOB Hamiltonian to comparable-mass simulations of spinning BHs will allow us to build an EOB model which includes finite mass-ratio effects in a more accurate way.

All that said, we have found that the presence of PN self-force terms in the spin EOB model of Ref. Barausse and Buonanno (2010) leads to dephasing of rad over one year depending on the EMRI system and the spin of the background. By contrast, the inclusion of the small object’s spin introduces dephasing on the order of a few radians. Taking into account those results, we can conclude that unless those finite mass-ratio effects are precisely modeled, currently, it is not worth requiring an agreement better than rads when calibrating the phase of the test-particle EOB waveforms against the Teukolsky-based waveforms. This in turn implies that calibrating parameters instead of to an EOB model is overkill as other systematics (induced by neglecting self-force effects) will be dominant.

We are grateful to Leor Barack, Carlos Sopuerta and Frans Pretorius for useful suggestions and comments. NY, AB, EB and YP, and SAH acknowledge support from the NSF grants PHY-0745779, PHY-0903631, and PHY-0449884; AB, SAH and MCM also acknowledge support from NASA grants NNX09AI81G, NNX08AL42G and NNX08AH29G. NY acknowledges support from the National Aeronautics and Space Administration through Einstein Postdoctoral Fellowship Award Number PF0-110080 issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060.

Appendix A Mapping between Spheroidal and Spherical Harmonics

a.1 Quantities at Spatial Infinity

Let us first consider the GW fluxes that describe radiation escaping to infinity, and later discuss the radiation absorbed by the BH. Frequency-domain Teukolsky equation codes expand the curvature scalar as


We have here incorporated the dependence into the spin-weighted spheroidal harmonics . These harmonics depend on the value of , and the minus superscript is a reminder that the we consider here harmonics of spin-weight . Throughout this appendix, the index refers to the spheroidal harmonic index, while refers to the spherical harmonic index.

From , we compute waveforms via , and hence, for a frequency-domain application,


We have implicitly assumed that is time-independent, or at least that its time-dependence is subleading. A better expansion is to re-express things in terms of the accumulated phase, i.e., the integral of the frequency , namely


The EOB and numerical relativity (NR) communities like to project these quantities onto a basis of spin-weighted spherical harmonics. For , they define via


and the harmonically-decomposed waveforms via


The minus superscript again denotes spin-weight . Defining the inner-product


the extraction of the and is simple:


Let us now take the Schwarzschild limit to see whether these expressions simplify. In that limit


and thus, performing the necessary projections and taking advantage of the orthonormality of spherical harmonics, we find


In spinning backgrounds, however, the mapping is more complicated. We can use the fact that spheroidal harmonics can be expressed as a sum of spherical harmonics via


The dependence on enters through the expansion coefficients (see, e.g. Ref. Hughes (2001b)). Inserting this expansion into Eq. (33) and using the inner-product definition, Eq. (40), we find


In the Schwarzschild limit, , so that Kerr simply limits as it should.

From the definition of the Isaacson stress-energy tensor, one can easily show that


We have here used the orthonormality of both the spheroidal and spherical harmonics to simplify the sums, as well as the fact that the time dependence of the disappears when its modulus is computed. Performing the angular integrals leaves us with familiar formulas:


This breaks down nicely enough that it is useful and sensible to define the modal contributions .

a.2 Quantities at Event Horizons

As a general principle, computing quantities that are related to an event horizon is usually more complicated than computing the same quantities at spatial infinity. For the fluxes, for example, this is because there is no simple generalization of the Isaacson tensor on the horizon. Instead, one must examine the shear of the horizon’s generators, look at how this shear generates entropy, and then apply the area theorem to compute fluxes Teukolsky and Press (1974). The relevant quantities at the horizon depend on the Newman-Penrose scalar , instead of , the former of which is a quantity of spin-weight , rather than .

The GW energy flux per unit solid angle at the horizon is given by


where , and where is the angular velocity of observers co-rotating with the event horizon. The quantity is the shear to the horizon’s generators as found by Ref. Hawking and Hartle (1972). This quantity is fairly simply related to , so let us introduce an expansion of in spin-weight spheroidal harmonics


where is given in Eq. (7). Notice that this quantity diverges on the event horizon because the Kinnersley tetrad, which is used to define the projection for , is ill-behaved as . This can be corrected for by converting to for any given mode Teukolsky and Press (1974)


The complex number is given by , where