Adiabatic waveforms for extreme mass-ratio inspirals via multivoice decomposition in time and frequency
Scott A. Hughes, Niels Warburton, Gaurav Khanna, Alvin J. K. Chua, Michael L. Katz
AAdiabatic waveforms for extreme mass-ratio inspirals viamultivoice decomposition in time and frequency
Scott A. Hughes, Niels Warburton, Gaurav Khanna,
3, 4
Alvin J. K. Chua, and Michael L. Katz Department of Physics and MIT Kavli Institute, Cambridge, MA 02139, United States School of Mathematics and Statistics, University College Dublin, Belfield, Dublin 4, Ireland Physics Department, University of Massachusetts, Dartmouth, MA 02747, United States Department of Physics, University of Rhode Island, Kingston, RI 02881, United States Theoretical Astrophysics Group, California Institute of Technology, Pasadena, CA 91125, United States Max-Planck-Institut f¨ur Gravitationsphysik, Albert-Einstein-Institut,Am M¨uhlenberg 1, 14476 Potsdam-Golm, Germany
We compute adiabatic waveforms for extreme mass-ratio inspirals (EMRIs) by “stitching” togethera long inspiral waveform from a sequence of waveform snapshots, each of which corresponds to aparticular geodesic orbit. We show that the complicated total waveform can be regarded as a sumof “voices.” Each voice evolves in a simple way on long timescales, which can be exploited toefficiently produce waveform models that faithfully encode the properties of EMRI systems. Welook at examples for a range of different orbital geometries: spherical orbits, equatorial eccentricorbits, and one example of generic (inclined and eccentric) orbits. To our knowledge, this is thefirst calculation of a generic EMRI waveform that uses strong-field radiation reaction. We examinewaveforms in both the time and frequency domains. Although EMRIs evolve slowly enough thatthe stationary phase approximation (SPA) to the Fourier transform is valid, the SPA calculationmust be done to higher order for some voices, since their instantaneous frequency can change fromchirping forward ( ˙ f >
0) to chirping backward ( ˙ f <
I. INTRODUCTIONA. Extreme mass-ratio inspirals and self forces
The large mass ratio limit of the two-body problem isa laboratory for studying strong-field motion in generalrelativity. By treating the spacetime of such a binaryas an exact black hole solution plus a perturbation dueto the less massive orbiting body, it is possible to ana-lyze the binary’s dynamics and the gravitational waves itgenerates using the tools of black hole perturbation the-ory. Because the equations of perturbation theory can, inmany circumstances, be solved very precisely, large-massratio serves as a high-precision limit for understandingthe two-body problem in general (see, e.g., [1]).This limit is also significant thanks to the importanceof extreme mass-ratio inspirals (EMRIs) as gravitationalwave sources. Binaries formed by the capture of stellarmass ( µ ∼ − M (cid:12) ) compact bodies onto relativis-tic orbits of black holes with M ∼ M (cid:12) in the coresof galaxies will generate low-frequency ( f ∼ .
01 Hz)gravitational waves, right in the sensitive band of space-based detectors like LISA. A typical EMRI will execute ∼ − orbits in LISA’s band as gravitational-wavebackreaction shrinks its binary separation. Because thesmall body orbits in the larger black hole’s strong field,EMRI waves are highly sensitive to the nature of thelarge black hole’s spacetime. EMRI events will make itpossible to precisely map black hole spacetimes, weigh-ing black holes’ masses and spins with exquisite accuracy,and testing the hypothesis that astrophysical black holesare described by the Kerr spacetime [2]. To achieve these ambitious science goals for EMRImeasurements, we will need waveform models, or tem-plates, that accurately match signals in data over theirfull duration. Such templates will provide guidance toalgorithms for finding EMRI signals in detector noise,and will be necessary for characterizing astrophysicalsources. Developing such models is one of the goals ofthe self force program, which seeks to develop equationsdescribing the motion of objects in specified backgroundspacetimes, including the interaction of that object withits own perturbation to that spacetime — i.e., includingthe small body’s “self interaction” [3]. Taking the back-ground spacetime to be that of a black hole, self forcescan be developed using tools from black hole perturba-tion theory, with the mass ratio of the system ε ≡ µ/M (where µ is the mass of the orbiting body, and M themass of the black hole) serving as a perturbative ordercounting parameter.To make this more quantitative, we sketch the gen-eral form of such equations of motion in the action-angleformulation used by Hinderer and Flanagan [4]. Let q α . = ( q t , q r , q θ , q φ ) be a set of angle variables which de-scribe the motion of the smaller body in a convenientcoordinate system, let J β . = ( J t , J r , J θ , J φ ) be a set of ac-tions associated with those motions, and let λ be a con-venient time variable. The motion of the smaller body isthen governed by a set of equations with the form dq α dλ = ω α ( J ) + ε g (1) α ( q r , q θ ; J ) + ε g (2) α ( q r , q θ ; J ) + . . . (1.1) dJ α dλ = ε G (1) α ( q r , q θ ; J ) + ε G (2) α ( q r , q θ ; J ) + . . . (1.2) a r X i v : . [ g r- q c ] F e b The frequency ω α describes the rate at which the anglesaccumulate per unit λ , neglecting the self interaction.The ω α ( J ) thus characterize the geodesic motion of thesmall body in the black hole background. The terms g ( n ) α describe how the small body’s trajectory is modifiedby the self interaction at O ( ε n ); the terms G ( n ) α describehow the actions (which are constant along geodesics) aremodified. Note that the self-interaction terms only de-pend on the angles q r and q θ — because black hole space-times are axisymmetric and stationary, these terms areindependent of q t and q φ . Many aspects of this problemare now under control at O ( ε ) (see, e.g., [5]), and workis rapidly proceeding on the problem at O ( ε ) [6].As the issue of the smaller body’s motion is brought un-der control, more attention is now being paid to the grav-itational waveforms that arise from this motion. That isthe focus of this paper. B. The adiabatic approximation and its use
The forcing terms in Eqs. (1.1) and (1.2) can be furtherdecomposed by splitting them into averages and oscilla-tions about their average. Consider the first-order forcingterm G (1) α . We put G (1) α ( q r , q θ ; J ) = (cid:104) G (1) α ( J ) (cid:105) + δG (1) α ( q r , q θ ; J ) , (1.3)where the averaged contribution is (cid:104) G (1) α ( J ) (cid:105) = 1(2 π ) (cid:90) π dq r (cid:90) π dq θ G (1) α ( q r , q θ ; J ) , (1.4)and we define the oscillations about this average as δG (1) α ( q r , q θ ; J ) ≡ G (1) α ( q r , q θ ; J ) − (cid:104) G (1) α ( J ) (cid:105) . The os-cillations vary about zero on a rapid orbital timescale T o ∼ M ; the average evolves on a much slower dissipa-tive inspiral timescale T i ∼ M /µ , or T i ∼ T o /ε . Becauseof the large separation of these two timescales for EMRIs,the oscillations nearly average away during an inspiral.Neglecting the impact of the oscillations introduces errorsof O ( T o /T i ) = O ( ε ).The simplest model for inspiral which captures thestrong-field dynamics of black hole orbits is known asthe adiabatic approximation . It amounts to solving thefollowing variants of Eqs. (1.1) and (1.2): dq α dλ = ω α ( J ) , dJ α dλ = ε (cid:104) G (1) α ( J ) (cid:105) . (1.5)In words, we treat the short-timescale orbital dynamicsas geodesic, but use the orbit-averaged impact of the selfforce on the actions J α . This amounts to including the“dissipative” part of the orbit-averaged self force, sinceto O ( ε ), the action of (cid:104) G (1) α (cid:105) is equivalent to computingthe rates at which an orbit’s energy E , axial angular mo-mentum L z , and Carter constant Q change due to thebackreaction of gravitational-wave emission. The adia-batic approximation treats inspiral as a “flow” through a sequence of geodesics, with the rate of flow determinedby the rates of change of E , L z , and Q [7].The computational cost of even this simplified model isquite high. As we outline in Sec. III, computing adiabaticbackreaction involves solving the Teukolsky equation formany multipoles of the radiation field and harmonics ofthe orbital motion; tens of thousands of multipoles andharmonics may need to be calculated for tens of thou-sands of orbits. To produce EMRI models which capturethe most important elements of the mapping betweensource physics and waveform properties, the communityhas developed several kludges as a stop-gap for data anal-ysis and science return studies. The analytic kludge ofBarack and Cutler [8] essentially pushes post-Newtonianmodels beyond their domain of validity. Their matchwith fully relativistic models is not good, but they cap-ture the key qualitative features of EMRI physics. Al-most as important, the analytic kludge is very fast andeasy to implement; as such, it has been heavily used formany LISA measurement studies.The numerical kludge of Babak et al. [9] is closer to thespirit of the adiabatic approximation we describe here, inthat it treats the small body’s motion as a Kerr geodesic,but uses semi-analytic fits to strong-field radiation emis-sion to describe inspiral. Wave emission is treated witha crude multipolar approximation based on the smallbody’s coordinate motion. Despite the crudeness of someof the underlying approximations, the numerical kludgefits relativistic waveform models quite well. It is howeverslower and harder to use than the analytic kludge, andas such has not been used very much. More recently,Chua and Gair [10] showed that one can significantly im-prove matches to relativistic models by using an analytickludge augmented with knowledge of the exact frequencyspectrum of Kerr black hole orbits.The shortcomings of the kludges illustrate that, ul-timately, one needs waveform models that capture thestrong-field dynamics of Kerr orbits and that accuratelydescribe strong-field radiation generation and propaga-tion through black hole spacetimes. Waveforms based onthe adiabatic approximation are the simplest ones thataccurately include both of these effects. Though the adi-abatic approximation misses important aspects from ne-glected pieces of the self force, they get enough of thestrong-field physics correct that they will be effective anduseful tools for understanding the scope of EMRI dataanalysis challenges, and for accurately assessing the sci-ence return that EMRI measurements will enable. C. This paper
Although the computational cost of making adiabaticinspiral waveforms is high, the most expensive step ofthis calculation — computing a set of complex numberswhich encode rates of change of (
E, L z , Q ), as well as thegravitational waveform’s amplitude — need only be doneonce. These numbers can be computed in advance for arange of astrophysically relevant EMRI orbits, and thenstored and used to assemble the waveform as needed. Thegoal of this paper is to lay out what quantities must becomputed, and to describe how to use such precomputeddata to build adiabatic waveforms.Our particular goal is to show how to organize andstore the most important and useful data needed to as-semble the waveforms in a computationally effective way.A key element of our approach is to view the compli-cated EMRI waveform as a sum of simple “voices.” Eachvoice corresponds to a mode ( l, m, k, n ) representing aparticular multipole of the radiation and harmonic of thefundamental orbital frequencies. The voice-by-voice de-composition was suggested long ago to one of us by L.S. Finn, and was first presented in Ref. [11] for the spe-cial case of spherical Kerr orbits (i.e., orbits of constantBoyer-Lindquist radius, including those inclined from theequatorial plane). This paper corrects an important er-ror in Ref. [11] (which left out a crucial phase that is setby initial conditions), and extends this analysis to fullygeneric configurations.The approach that we present has several importantfeatures. First, we find that the data which must bestored to describe the waveform voice-by-voice evolvessmoothly over an inspiral. This suggests that these datacan be sampled at a relatively low cadence, and we canthen build a high-quality waveform using interpolation.Such behavior had been seen in earlier work [11]; thisanalysis demonstrates that this behavior is not unique tospherical orbits, but is generic. We describe the meth-ods we have used for our initial exploration, and thatwere used in a related study [12] focusing on the rapidevaluation of waveforms for data analysis. Reference [12]showed that the chasm between the computational de-mands of accurate modeling and efficient analysis canbe bridged, but much work remains to “optimize” thesemethods — for example, in determining the best way forthe numerical data to be sampled and interpolated.Second, we show that this framework works well inboth the time and frequency domains. The frequency-domain construction is particularly interesting. Thanksto their extreme mass ratios, EMRIs evolve slowlyenough that the stationary-phase approximation (SPA)to the Fourier transform should provide an excellent rep-resentation of the frequency-domain form. However, forsome EMRI voices, the instantaneous frequency grows toa maximum and then decreases. For such voices, the firsttime derivative of the frequency vanishes at some pointalong the inspiral. The standard SPA is singular at thesepoints. We show that including information about thesecond derivative fixes this behavior, allowing us to com-pute frequency-domain waveforms for all EMRI voices.Third, we show that several of the waveform’s param-eters can be included in a simple way which effectivelyreduces the dimensionality of the waveform parameterspace. These parameters are angles which control theinitial position of the orbit in its eccentric radial motion,and its initial polar angle in the range of motion allowed by its orbital inclination. They are initial conditions onthe relativistic analog of the “true anomaly angles” usedin Newtonian celestial mechanics. Associated with theseinitial anomaly angles are phases, originally introducedin Ref. [13], which correct the complex amplitudes asso-ciated with the waveform’s voices. By comparing withoutput from a time-domain Teukolsky equation solver[14, 15], we show that these phases allow one to matchany allowed initial condition with very little computa-tional cost. This means that we can generate a suiteof data using only a single initial choice of the anomalyangles, and then transform to initial conditions corre-sponding to any other choice. This significantly reducesthe computational cost associated with covering the fullEMRI parameter space.Finally, it should not be difficult to extend this frame-work to include at least some important effects beyondthe adiabatic approximation, many of which are dis-cussed in detail in Refs. [16, 17]. For example, bothconservative self forces and spin-curvature coupling haveorbit-averaged effects that are small, but that secularlyaccumulate over many orbits [18–21]. Such effects canbe included in this framework by allowing the relativisticanomaly angles (which in the adiabatic limit are con-stant) to evolve over the inspiral; the phases associatedwith these angles will evolve as well. We also expectthat the impact of small oscillations (such as arise fromboth self forces and spin-curvature coupling) can likewisebe incorporated, perhaps very efficiently using a near-identity transformation [22].As we were completing the bulk of the calculationswhich appear here, a similar analysis of adiabatic EMRIwaveforms was presented by Fujita and Shibata [23].Their analysis focuses to a large extent on the measura-bility of EMRI waves by LISA, confining their discussionto eccentric and equatorial sources. Like us, they alsotake advantage of the fact that one can pre-compute themost expensive data on a grid of orbits, and then assem-ble the waveform by interpolation. We are encouraged bythe fact that their independently developed framework issubstantially similar to what we present here. D. Organization of this paper
The remainder of this paper is organized as follows.Since inspiral in the adiabatic approximation is treatedas a sequence of geodesic orbits, we begin by reviewingthe properties of Kerr geodesics in Sec. II. Nothing inthis section is new; it is included primarily to keep themanuscript self contained, and to allow us to carefullydefine our notation and the meaning of important quan-tities which are used elsewhere. Section III reviews howwe solve the Teukolsky equation and use its solutionsin order to calculate adiabatic backreaction on an orbit.This allows us to compute how a system evolves fromorbit to orbit, as well as the gravitational-wave ampli-tude produced by that orbit. This material is likewisenot new, but is included for completeness, as well as tointroduce and explain all relevant notation.In Sec. IV, we lay out how one “stitches” together datadescribing radiation from geodesics to construct an adi-abatic waveform. This construction essentially amountsto taking the solution to the Teukolsky equation for ageodesic orbit and promoting various factors which areconstants on geodesics into factors which vary slowlyalong an inspiral. This construction introduces a two-timescale expansion: some quantities vary on the “fast”timescale associated with orbital motions, T o ∼ M ; oth-ers vary on the “slow” timescale associated with the in-spiral, T i ∼ M /µ = T o /ε . The adiabatic waveform isonly accurate up to corrections of order the system’s massratio, essentially because it assumes that all time deriva-tives only include information about the system’s “fast”time variation. One way in which a post-adiabatic anal-ysis will improve on these results will be by includinginformation about time derivatives with respect to theslow variation along the inspiral.Section V describes how to compute multivoice sig-nals in the frequency domain. Because EMRI systemsare slowly evolving, the stationary phase approximation(SPA) should accurately describe the Fourier transformof EMRI signals. However, because the evolution of cer-tain voices is not monotonic, the “standard” SPA calcu-lation can fail, introducing singular artifacts at momentswhen a voice’s frequency evolution switches sign. We re-view the standard SPA Fourier transform and show howby including an additional derivative of the frequency itis straightforward to correct this artifact. We concludethis section by showing how to combine multiple voicesto construct the full frequency-domain EMRI waveform.In Sec. VI, we present various important technical de-tails describing how we implement this framework forthe results we present in this paper. We strongly empha-size that there is a great deal of room for improvementon the techniques described here. We have not, for ex-ample, carefully assessed the most effective method forlaying out the grid of data on which we store informationabout adiabatic backreaction and waveform amplitudes,nor have we thoroughly investigated efficient methodsfor interpolating these data to off-grid points (e.g., [12]).These important points will be studied in future work, aswe begin assessing how to take this framework and useit to develop EMRI waveforms in support of LISA dataanalysis and science studies.In Secs. VII and VIII we present examples of adiabaticEMRI waveforms. In both of these sections, we show ex-amples of the complete time-domain waveform producedby summing over many voices, as well as the (much sim-pler) structure of representative voices which contributeto these waveforms. Section VII shows results for inspi-ral into Schwarzschild black holes, presenting the detailsof an inspiral with small initial eccentricity ( e init = 0 . e init = 0 . adhoc correction that corrects some of the neglected “slowtime” evolution. Though this correction is not rigorouslyjustified, its form suggests that the dephasing may arisefrom a slow-time evolution of the phases variables whichare set by the orbit’s initial conditions.Our conclusions are presented in Sec. IX. Along withsummarizing the main results of this analysis, we describeplans and directions for future work. Chief among theseplans is to investigate how to optimize the implemen-tation in order to make EMRI waveforms as rapidly aspossible, which will be very important in order for thesewaveforms to be usable for LISA data and science anal-ysis studies, as well as investigating how to include post-adiabatic effects with small modifications of the frame-work we present here. We also plan to publicly releasethe code and data used in this study, and describe thestatus of our plans as this analysis is being completed. II. KERR GEODESICS
We begin by discussing bound Kerr geodesics. Themost important aspects of this content are discussed indepth elsewhere [24–28]; we briefly review this materialfor this paper to be self contained, as well as to introducenotation and conventions that we use. Certain lengthybut important formulas are given in Appendix A.
A. Mino-time formulation of orbital motion
Consider a point-like body of mass µ orbiting a Kerrblack hole with mass M and spin parameter a = | S | /M (where S is the black hole’s spin angular momentum inunits with G = 1 = c ), and use Boyer-Lindquist coordi-nates (with the angle θ measured from the black hole’sspin axis) to describe its motion. We use Mino time asour time parameter describing these orbits. An intervalof Mino time dλ is related to an interval of proper time dτ by dλ = dτ / Σ, where Σ = r +cos θ , and where r and θ are the Boyer-Lindquist radial and polar coordinates.With this parameterization, motion in Boyer-Lindquistcoordinates is governed by the equations (cid:18) drdλ (cid:19) = [ E ( r + a ) − aL z ] − ∆[ r + ( L z − aE ) + Q ] ≡ R ( r ) , (2.1) (cid:18) dθdλ (cid:19) = Q − cot θL z − a cos θ (1 − E ) ≡ Θ( θ ) , (2.2) dφdλ = csc θL z + 2 M raE ∆ − a L z ∆ ≡ Φ r ( r ) + Φ θ ( θ ) , (2.3) dtdλ = E (cid:20) ( r + a ) ∆ − a sin θ (cid:21) − M raL z ∆ ≡ T r ( r ) + T θ ( θ ) . (2.4)We have introduced ∆ = r − M r + a . The quantities E , L z , and Q are the orbit’s energy (per unit µ ), axialangular momentum (per unit µ ), and Carter constant(per unit µ ). These quantities are conserved along anygeodesic; choosing them specifies an orbit, up to initialconditions. Writing d/dλ = Σ d/dτ puts these equationsinto more familiar forms typically found in textbooks,such as Eqs. (33.32a-d) of Ref. [24].Equations (2.1) and (2.2) tell us that bound Kerr orbitsare periodic in r and θ when parameterized using λ : r ( λ ) = r ( λ + n Λ r ) , θ ( λ ) = θ ( λ + k Λ θ ) , (2.5)where n and k are each integers. Simple formulas existfor the Mino-time periods Λ r and Λ θ [27]; we define theassociated frequencies by Υ r,θ = 2 π/ Λ r,θ .The motions in t and φ are the sum of secularly accu-mulating pieces and oscillatory functions: t ( λ ) = t + Γ λ + ∆ t r [ r ( λ )] + ∆ t θ [ θ ( λ )] , (2.6) φ ( λ ) = φ + Υ φ λ + ∆ φ r [ r ( λ )] + ∆ φ θ [ θ ( λ )] . (2.7)In these equations, t and φ describe initial conditions,Γ = (cid:104) T r ( r ) (cid:105) + (cid:104) T θ ( θ ) (cid:105) , (2.8)Υ φ = (cid:104) Φ r ( r ) (cid:105) + (cid:104) Φ θ ( θ ) (cid:105) , (2.9) and ∆ t r [ r ( λ )] = T r [ r ( λ )] − (cid:104) T r ( r ) (cid:105) ≡ ∆ t r ( λ ) , ∆ t θ [ θ ( λ )] = T θ [ θ ( λ )] − (cid:104) T θ ( θ ) (cid:105) ≡ ∆ t θ ( λ ) ; (2.10)∆ φ r [ r ( λ )] = Φ r [ r ( λ )] − (cid:104) Φ r ( r ) (cid:105) ≡ ∆ φ r ( λ ) , ∆ φ θ [ θ ( λ )] = Φ θ [ θ ( λ )] − (cid:104) Φ θ ( θ ) (cid:105) ≡ ∆ φ θ ( λ ) . (2.11)The quantity Γ describes the mean rate at which observertime t accumulates per unit λ ; the Mino-time frequencyΥ φ describes the mean rate at which φ accumulates perunit λ . The associated period is Λ φ = 2 π/ Υ φ . Simpleformulas likewise exist for Γ and Υ φ [27]. The averagesused in Eqs. (2.8)–(2.11) are given by (cid:104) f r ( r ) (cid:105) = 1Λ r (cid:90) Λ r f r [ r ( λ )] dλ , (2.12) (cid:104) f θ ( θ ) (cid:105) = 1Λ θ (cid:90) Λ θ f θ [ θ ( λ )] dλ . (2.13)The ratio of the Mino-time frequencies to Γ gives theobserver-time frequencies:Ω r,θ,φ = Υ r,θ,φ Γ . (2.14)We thus have useful closed-form expressions for all fre-quencies, conjugate to both Mino time λ and observertime t , that characterize black hole orbits. B. Orbit parameterization and initial conditions
Take the orbit to oscillate over θ min ≤ θ ≤ θ max , with θ max = π − θ min , and over r min ≤ r ≤ r max , with r min / max = p ± e . (2.15)Choosing p , e , and θ min is equivalent to choosing theintegrals of motion E , L z , and Q . We have found it par-ticularly convenient to replace θ min with an inclinationangle I , defined by I = π/ − sgn( L z ) θ min . (2.16)The angle I varies smoothly from 0 for equatorial pro-grade to π for equatorial retrograde. The definition (2.16)may seem a bit awkward thanks to the branch associatedwith the sign of L z . It is simple to show thatcos θ min = sin I . (2.17) This angle has often been labeled θ inc in previous work, such asRef. [29]. We have found this label to be potentially confusing,since the Boyer-Lindquist angle θ is measured from the blackhole’s spin axis, whereas the inclination angle is measured fromthe plane normal to this axis. We have changed notation to avoidconfusion with the coordinates. We have found that x I ≡ cos I is a particularly good pa-rameter to describe inclination: x I varies smoothly from1 to − L z having the same sign as x I .Schmidt [25] first showed how to compute ( E, L z , Q ) forgeneric Kerr orbits; a particularly clean representation isprovided by van de Meent [28]. We summarize his for-mulas in Appendix A, tweaking them slightly to use ourpreferred parameter set ( p, e, x I ).Initial conditions for t and φ are set by the parameters t and φ given in Eqs. (2.6) and (2.7). To set initialconditions on r and θ , we introduce anomaly angles χ r and χ θ to reparameterize those coordinate motions: r = p e cos( χ r + χ r ) , (2.18)cos θ = (cid:113) − x I cos( χ θ + χ θ ) . (2.19)We put χ θ = 0, χ r = 0, t = t , and φ = φ when λ = 0.The phase χ r then determines the value of r at λ = 0,and χ θ determines the corresponding value of θ . When χ θ = 0, the orbit has θ = θ min when λ = 0; when χ r = 0, it has r = r min when λ = 0.We define the fiducial geodesic to be the geodesic thathas χ θ = χ r = φ = 0 = t . We denote with a“check-mark” accent any quantity which is defined alongthe fiducial geodesic. For instance, ˇ r ( λ ) is orbital radiusalong the fiducial geodesic, ˇ θ ( λ ) is the polar angle θ alongthe fiducial geodesic. For non-fiducial geodesics, we de-fine λ = λ r to be the smallest positive value of λ atwhich r = r min ; likewise λ = λ θ is the smallest positivevalue of λ at which θ = θ min . This means that r ( λ ) = ˇ r ( λ − λ r ) , θ ( λ ) = ˇ θ ( λ − λ θ ) , (2.20)There is a one-to-one correspondence between λ θ and χ θ , and between λ r and χ r . A useful corollary is∆ t r ( λ ) = ∆ˇ t r ( λ − λ r ) − ∆ˇ t r ( − λ r ) , ∆ t θ ( λ ) = ∆ˇ t θ ( λ − λ θ ) − ∆ˇ t θ ( − λ θ ) , (2.21)with analogous formulas for ∆ φ r and ∆ φ θ . Some of our definitions differ from those used in Ref.[13]. In that reference , χ θ and χ r were not used. In-stead, χ θ = 0 corresponded to λ = λ θ , and χ r = 0corresponded to λ = λ rθ . As we discuss briefly in Secs.IV and IX, the angles χ r and χ θ will play an importantrole going beyond adiabatic waveforms. In the adiabaticapproximation, the angles χ r , χ θ , and φ are constantas we move from geodesic to geodesic. When we include,for example, orbit-averaged conservative self-force effectsor orbit-averaged spin-curvature forces, we find secularlyaccumulating phases associated with each of these mo-tions. Allowing the angles χ r , χ θ , and φ to evolveduring inspiral is a simple and robust way to “upgrade”this framework to include these next-order effects. III. ADIABATIC EVOLUTION ANDWAVEFORM AMPLITUDES VIA THETEUKOLSKY EQUATION
The next critical ingredient to constructing an adia-batic inspiral is the backreaction which arises from theorbit-averaged self interaction. The quantities which en-code the backreaction also tell us the amplitude of theinspiral’s associated gravitational waveform. In this sec-tion, we briefly summarize how these quantities are calcu-lated using the Teukolsky equation. As with Sec. II, thecontents of this section are discussed at length elsewhere,but are summarized here to introduce critical quantitiesand concepts important for later parts of this analysis.
A. Solving the Teukolsky equation
The Teukolsky equation [30] describes perturbationsto the Weyl curvature of Kerr black holes. The versionthat we will use in this analysis focuses on the Newman-Penrose curvature scalar ψ = − C αβγδ n α ¯ m β n γ ¯ m δ , (3.1)where C αβγδ is the Weyl curvature tensor, and n α and¯ m α are legs of the Newman-Penrose null tetrad [31].Teukolsky showed that ψ is governed by the equation (cid:20) ( r + a ) ∆ − a sin θ (cid:21) ∂ t Ψ − (cid:20) r + ia cos θ − M ( r − a )∆ (cid:21) ∂ t Ψ + 4
M ar ∆ ∂ φ ∂ t Ψ − ∆ ∂ r (cid:0) ∆ − ∂ r Ψ (cid:1) − θ ∂ θ (sin θ∂ θ Ψ) + (cid:20) a ∆ − θ (cid:21) ∂ φ Ψ + 4 (cid:20) a ( r − M )∆ + i cos θ sin θ (cid:21) ∂ φ Ψ + (cid:0) θ + 2 (cid:1) Ψ = 4 π Σ T . (3.2)The field Ψ = ( r − ia cos θ ) ψ , and T is a source term These definitions account for the fact that these motions are periodic: r ( λ r + n Λ r ) = r min for any integer n , and θ ( λ θ + whose precise form is not needed here. See Ref. [30] foradditional details and definitions.An important point for our analysis is that ψ = 12 d dt ( h + − ih × ) as r → ∞ , (3.3)so ψ far from the source encodes the emitted gravita-tional waves. These solutions also encode contributionsto the rates of change of E , L z , and Q from gravitational-wave backreaction. This is equivalent to the orbit-averaged self interaction arising from fields which areregular far from the source (see Ref. [32], as well asadditional discussion on this point in Ref. [33]). As r → r + ≡ M + √ M − a (the coordinate radius of theevent horizon), ψ encodes tidal interactions of the or-biting body with the black hole’s event horizon. Thesesolutions encode contributions to the rates of change of E , L z , and Q from radiation absorbed by the horizon,which is equivalent to the orbit-averaged self interactionarising from fields which are regular on the event hori-zon [32, 33]. Knowledge of ψ in the limits r → ∞ and r → r + provides all the data we need to construct adia-batic inspirals.The frequency-domain approach we use to solve theTeukolsky equation begins by writing ψ in a Fourier andmultipolar expansion. Writing ψ = 1( r − ia cos θ ) (cid:90) ∞−∞ dω ∞ (cid:88) l =2 l (cid:88) m = − l R lm ( r ; ω ) S lm ( ϑ ; aω ) e i [ mϕ − ω ( t − t )] , (3.4)Eq. (3.2) separates [30], with ordinary differential equa-tions governing R lm ( r, ω ) and S lm ( θ, aω ).The field ψ is measured at the event ( t, r, ϑ, ϕ ). (Notethe distinction between the orbit’s polar and axial an-gles, θ and φ , and the polar and axial angles at whichthe field is measured, ϑ and ϕ .) The function S lm ( ϑ ; aω )is a spheroidal harmonic (of spin-weight −
2, left out forbrevity); this function and methods for computing it arediscussed at length in Appendix A of Ref. [34]. For rea-sons we will explain below, we have also introduced theinitial time t into our expansion (3.4).The separated radial dependence has simple asymp-totic behavior: R lm ( r, ω ) → Z ∞ lmω r e iωr ∗ , r → ∞ , (3.5) → Z H lmω ∆ e − i ( ω − m Ω H ) r ∗ , r → r + . (3.6)(We have absorbed a coefficient C trans lmω into the defini-tion of Z ∞ lmω , and a coefficient B trans lmω into the definitionof Z H lmω ; see Refs. [35–37] for further discussion of thesequantities.) These asymptotic forms depend on the “tor-toise coordinate,” r ∗ ( r ) = r + M r + √ M − a ln (cid:18) rr + − (cid:19) − M r − √ M − a ln (cid:18) rr − − (cid:19) , (3.7)where r ± = M ± √ M − a . The frequencyΩ H = a M r + (3.8) k Λ θ ) = θ min for any integer k . Reference [13] also used slightly different symbols for the anomalyangles, writing χ for the polar anomaly, and ψ for radial. Here,we only use ψ for the Newman-Penrose curvature scalars. is often called the rotation frequency of the horizon. Itdescribes the frequency at which an observer held in-finitesimally outside the event horizon will orbit the blackhole as seen by distant observers.The quantities Z ∞ , H lmω are computed by integrating ho-mogeneous solutions of the radial Teukolsky equationagainst the source term of the separated radial Teukol-sky equation. Further detailed discussion can be foundin Ref. [29], with updates to notation and minor correc-tions in Ref. [38]. Of importance for this analysis is thatthese quantities are computed by evaluating integrals ofthe form Z ∞ , H lmω = (cid:90) ∞−∞ dτ e iω [ t ( τ ) − t ] e − imφ ( τ ) I ∞ , H lmω [ r ( τ ) , θ ( τ )] . (3.9)The integration variable τ is proper time along thegeodesic, and we subtract off t because it is alreadyaccounted for in Eq. (3.4). The function I ∞ , H lmω ( r, θ ) isdiscussed in Refs. [29] and [38]; schematically, it is aGreen’s function used to solve the radial Teukolsky equa-tion, multiplied by this equation’s source term. Using theproperties of Kerr geodesic orbits and using the methodsdeveloped in Refs. [35–37], it is well understood how tobuild I ∞ , H lmω [ r ( τ ) , θ ( τ )].Changing integration variable from proper time τ toMino time λ , and using Eqs. (2.6) and (2.7), this becomes Z ∞ , H lmω = e − imφ (cid:90) ∞−∞ dλ e i ( ω Γ − m Υ φ ) λ J ∞ , H lmω [ r ( λ ) , θ ( λ )] , (3.10)where we have introduced J ∞ , H lmω ( r, θ ) = ( r + a cos θ ) I ∞ , H lmω ( r, θ ) × e iω [∆ t r ( r )+∆ t θ ( θ )] e − im [∆ φ r ( r )+∆ φ θ ( θ )] . (3.11)By virtue of the periodicity of orbit’s r and θ motionswith respect to Mino time, the function J ∞ , H lmω can be written as a double Fourier series: J ∞ , H lmω = ∞ (cid:88) k = −∞ ∞ (cid:88) n = −∞ J ∞ , H lmkn e − i ( k Υ θ + n Υ r ) λ , (3.12)where J ∞ , H lmkn = 1Λ r Λ θ (cid:90) Λ r dλ r e in Υ r λ r (cid:90) Λ θ dλ θ e ik Υ θ λ θ J ∞ , H lmω [ r ( λ r ) , θ ( λ θ )] . (3.13)Combining Eqs. (3.10), (3.12), and (3.13) plus the rela-tions Ω r,θ,φ = Υ r,θ,φ / Γ, we find that Z ∞ , H lmω = ∞ (cid:88) k = −∞ ∞ (cid:88) n = −∞ Z ∞ , H lmkn δ ( ω − ω mkn ) , (3.14)where ω mkn = m Ω φ + k Ω θ + n Ω r (3.15)and Z ∞ , H lmkn = e − imφ J ∞ , H lmkn / Γ . (3.16)These coefficients have the symmetry Z ∞ , H l, − m, − k, − n = ( − ( l + k ) ¯ Z ∞ , H lmkn , (3.17)where overbar denotes complex conjugation. Our coderespects the symmetry (3.17) to double-precision accu-racy. We take advantage of this by computing Z ∞ , H lmkn forall l , all m , all k , and n ≥
0, then using Eq. (3.17) to inferresults for n <
0. This cuts the amount of computationroughly in half.As we describe in the following two subsections, thecoefficients Z ∞ lmkn set the amplitude of the gravitationalwaves from a specified geodesic orbit, and also encode thecontribution to the orbit-averaged self force from fieldsthat are regular at null infinity; the coefficients Z H lmkn encode the contribution to the orbit-averaged self forcefrom fields that are regular on the event horizon. Thesesets of coefficients are thus of crucial importance for com-puting adiabatic inspiral and its gravitational waveform. B. The gravitational waveform from an orbit
As shown in Sec. 8 of Ref. [13], the phase of Z ∞ , H lmkn de-pends on the values of λ r and λ θ , which in turn dependon the initial anomaly angles χ r and χ θ . We denote byˇ Z ∞ , H lmkn the value of Z ∞ , H lmkn for the fiducial geodesic. For ageneral geodesic, Z ∞ , H lmkn = e iξ mkn ˇ Z ∞ , H lmkn , (3.18) where the correcting phase is ξ mkn = k Υ θ λ θ + n Υ r λ r + m (cid:2) ∆ ˇ φ r ( − λ r ) + ∆ ˇ φ θ ( − λ θ ) − φ (cid:3) − ω mkn (cid:2) ∆ˇ t r ( − λ r ) + ∆ˇ t θ ( − λ θ ) (cid:3) . (3.19)The form of ξ mkn we use is slightly different from thatderived in Ref. [13]. In particular, we have separated outthe dependence on the initial axial phase φ , as shownin Eq. (3.16), and the dependence on t , which is dis-cussed further in Sec. IV. See Ref. [13] for discussion andderivation of the dependence on λ θ and λ r .The fact that the initial conditions only influence theseamplitudes via the phase factor ξ mkn means that we onlyneed to compute and store quantities on the fiducialgeodesic. Using Eq. (3.19), we can then easily convertthese results to any initial condition. This vastly cutsdown on the amount of computation that must be doneto cover the space of physically important EMRI systems.In addition, notice that ξ mkn can be written ξ mkn = mξ + kξ + nξ . (3.20)For each orbit one need only compute ( ξ , ξ , ξ ) inorder to know the phases for all ( l, m, k, n ).As we will see below, the values of these phases are ir-relevant for the orbit-averaged backreaction, but they arecritical for getting the phase of the gravitational wave-form correct given initial conditions. To write down thegravitational waves, we use Eq. (3.3) to relate h to ψ farfrom the source. Combining Eqs. (3.4), (3.5), and (3.14),we further know that ψ = 1 r (cid:88) lmkn Z ∞ lmkn S lm ( ϑ, aω mkn ) e i [ mϕ − ω mkn ( t − t )] (3.21)as r → ∞ . Here and in what follows, any sum over theset ( l, m, k, n ) will be assumed to be from 2 to ∞ for l ,from − l to l for m , and from −∞ to ∞ for k and n . Letus define h ≡ h + − ih × ≡ r (cid:88) lmkn h lmkn = 1 r (cid:88) lmkn A lmkn S lm ( ϑ ; aω mkn ) e i [ mϕ − ω mkn ( t − t )]) . (3.22)Combining Eqs. (3.3), (3.21), and (3.22), we see that A lmkn = − Z ∞ lmkn ω mkn . (3.23)As with Z ∞ lmkn , we define ˇ A lmkn as the wave amplitudefor the fiducial geodesic, and we have A lmkn = e iξ mkn ˇ A lmkn . (3.24)The data ˇ A lmkn interpolate very well and should bestored for generating inspiral waveforms. It will be con-venient for later discussion to further define H lmkn = A lmkn S lm ( ϑ ; aω mkn ) (3.25)Because spheroidal harmonics slowly change along an in-spiral as ω mkn evolves, we find it useful to examine H lmkn rather than A lmkn when examining wave amplitudes dur-ing inspiral. C. Adiabatic backreaction
Here we summarize how to use the coefficients Z ∞ , H lmkn to compute the adiabatic dissipative evolution, or back- reaction, on a geodesic. We assume that a body is on aKerr geodesic orbit, and so is characterized (up to initialconditions) by the orbital integrals E , L z , and Q . Re-sults for dE/dt and dL z /dt have been known for quitea long time [39]; these quantities each split into a con-tribution from fields that are regular at infinity, whichcan be extracted from knowledge of the distant gravi-tational radiation, and a contribution from fields thatare regular on the black hole’s event horizon. Comput-ing the down-horizon contribution is a little more tricky;one must compute how tidal stresses shear the horizon’sgenerators, increasing its surface area (or its entropy),and then apply the first law of black hole dynamics toinfer the change in the hole’s mass and spin. Resultsfor dQ/dt were first derived by Sago et al. [40], and arefound by averaging the action of the dissipative self forceon a geodesic. It also separates into pieces that arise fromfields regular at infinity and fields regular on the horizon.The results which we need for our analysis are: (cid:18) dEdt (cid:19) ∞ = (cid:88) lmkn | Z ∞ lmkn | πω mkn , (cid:18) dEdt (cid:19) H = (cid:88) lmkn α lmkn | Z H lmkn | πω mkn , (3.26) (cid:18) dL z dt (cid:19) ∞ = (cid:88) lmkn m | Z ∞ lmkn | πω mkn , (cid:18) dL z dt (cid:19) H = (cid:88) lmkn α lmkn m | Z H lmkn | πω mkn , (3.27) (cid:18) dQdt (cid:19) ∞ = (cid:88) lmkn | Z ∞ lmkn | × ( L mkn + k Υ θ )2 πω mkn , (cid:18) dQdt (cid:19) H = (cid:88) lmkn α lmkn | Z H lmkn | × ( L mkn + k Υ θ )2 πω mkn . (3.28)In these equations, L mkn = m (cid:104) cot θ (cid:105) L z − a ω mkn (cid:104) cos θ (cid:105) E , (3.29) α lmkn = 256(2 M r + ) ( ω mkn − m Ω H )[( ω mkn − m Ω H ) + 4 (cid:15) ][( ω mkn − Ω H ) + 16 (cid:15) ] ω mkn | C lmkn | . (3.30)The terms (cid:104) cot θ (cid:105) and (cid:104) cos θ (cid:105) in Eq. (3.29) mean cot θ and cos θ evaluated at the θ coordinate along the orbit,and then averaged using Eq. (2.13). The terms | C lmkn | and (cid:15) in Eq. (3.30) are given by | C lmkn | = [( λ lmkn + 2) + 4 aω mkn − a ω mkn ]( λ lmkn + 36 amω mkn − a ω mkn )+ (2 λ lmkn + 3)(96 a ω mkn − amω mkn ) + 144 ω mkn ( M − a ) , (3.31) (cid:15) = √ M − a M r + . (3.32)The quantity λ lmkn appearing in Eq. (3.31) is one formof the eigenvalue of the spheroidal harmonic; see Ref. [34]for discussion of the algorithm we use to compute it, andAppendix C of Ref. [38] for discussion of the various forms of the eigenvalue in the literature (which are simple toconvert between). Note that the factor (cid:15) which appearshere (and is not used elsewhere in this paper) is distinctfrom the similar symbol ε ≡ µ/M , the mass ratio.0In evaluating Eqs. (3.26), (3.27), and (3.28), we musttruncate the infinite sums, with cutoffs determined bythe needs of the analysis in question. For the purpose ofthis paper, we have implemented the following cutoffs: • We include all l from 2 to 10. • At each l , we include all m from − l to l . • We truncate the k sum when the fractional changeto the accumulated sum is smaller than 10 − forthree consecutive values of k . Holding all otherindices fixed, contributions to this sum tend to falloff monotonically and fairly rapidly with increasing k , so in practice this condition means that neglectedterms change the sum by less than several × − . • We truncate the n sum when the fractional changeto the accumulated sum is smaller than 10 − forthree consecutive values of n . Especially for e (cid:38) .
4, contributions to this sum do not fall off mono-tonically with n until some threshold has beenpassed (see Figs. 5 and 6 of Ref. [41], and Figs.2 and 3 of Ref. [29]). Once past that threshold,convergence is quick, and we find that neglectedterms in this case also change the sum by less thanseveral × − .We emphasize that these cutoffs have not been selectedcarefully, but are simply chosen for ease of calculationand to produce results which are “converged enough” forthe exploratory purposes of this paper. A more carefulanalysis and assessment of how to truncate these sumsshould be done before using these ideas to make “pro-duction quality” waveforms (e.g., for exploring LISA dataanalysis questions, or science return with EMRI measure-ments) in order to make sure that systematic errors inwaveform modeling are understood and do not adverselyaffect one’s analysis.As E , L z , and Q change, we require the system toevolve from one geodesic to another. To do this, we letthese orbital integrals change by enforcing a balance law: (cid:18) d C dt (cid:19) orbit = − (cid:18) d C dt (cid:19) ∞ − (cid:18) d C dt (cid:19) H , (3.33)for C ∈ [ E, L z , Q ]. As these orbital integrals evolve, theorbit’s geometry slowly changes in order to keep the sys-tem on a geodesic trajectory. Appendix B shows howto relate rates of change for p , e , and x I to the rates ofchange of E , L z , and Q . IV. ADIABATIC INSPIRAL AS DISSIPATIVEEVOLUTION ALONG A SEQUENCE OFGEODESICS
We now examine solutions to the Teukolsky equationfor a slowly evolving source. Critical to our analysis is theidea of a two-timescale expansion: the waveform phase varies on an “fast” orbital timescale T o ∼ M , and orbitcharacteristics vary on a “slow” inspiral timescale T i ∼ M /µ . The two timescales differ by the system’s massratio: T o /T i = µ/M ≡ ε . The waveforms we compute inthis way are accurate up to corrections of order ε .Suppose that we have used the rates of change dE/dt , dL z /dt , dQ/dt to compute how a system evolves fromgeodesic to geodesic. We parameterize inspiral by a book-keeper time t i which measures evolution along the inspi-ral as seen by a distant observer. We treat the inspiralas a geodesic at each moment t i , and call the geodesic atthis moment the “osculating” geodesic [42–44]. At eachsuch moment, the osculating geodesic’s energy E ( t i ), itsangular momentum L z ( t i ), and its Carter constant Q ( t i )are known. By our assumption that the inspiral is ageodesic at each moment, we can reparameterize and de-termine p ( t i ), e ( t i ), and x I ( t i ). We can also computequantities such as the frequencies Ω r,θ,φ ( t i ) and the am-plitudes Z ∞ lmkn ( t i ) for each geodesic in this sequence.To leading order in ε , the curvature scalar ψ thatarises from this sequence of geodesics is given by ψ ( t i ) = 1 r (cid:88) lmkn Z ∞ lmkn ( t i ) S lm [ ϑ ; aω mkn ( t i )] × e i [ mϕ − Φ mkn ( t i )] . (4.1)This is Eq. (3.21), but with the amplitude Z ∞ lmkn and thefrequency ω mkn now functions of t i . Notice the depen-dence on harmonics of the accumulated orbital phase:Φ mkn ( t i ) = (cid:90) t i t ω mkn ( t (cid:48) ) dt (cid:48) . (4.2)This reduces to ω mkn ( t i − t ) in the limit that the orbitdoes not inspiral. Equation (4.2) builds in the depen-dence on the initial time t , which is why we leave thisparameter out of the factor ξ mkn in Eq. (3.19).To justify this inspiraling solution for ψ , substituteEq. (4.1) into the Teukolsky equation, Eq. (3.2). Doingso, one finds that it satisfies the Teukolsky equation upto errors of order the orbital timescale over the inspiraltimescale, O ( T o /T i ) ∼ O ( ε ). These errors in turn arisefrom the fact that time derivatives have both fast-timecontributions, for which ∂ t ∼ /T o ∼ ω mkn , as well asslow-time contributions, for which ∂ t ∼ /T i ∼ ε/T o . Inthe adiabatic approximation, we neglect the slow-timederivatives, expecting that at any moment their contri-bution will be small as long as the system’s mass ratiois large. Some of the errors arising from this neglect canaccumulate secularly, leading to phase errors up to sev-eral radians after an inspiral. Post-adiabatic correctionswill change the amplitude and phase of Eq. (4.1) in such away as to correct the adiabatic approximation’s fast timeover slow time errors. See [16] for further discussion.In our applications, we are typically more interested in1the waveform h ( t i ) than in ψ ( t i ). This is given by h ( t i ) ≡ (cid:88) lmkn h lmkn ( t i ) = 1 r (cid:88) lmkn H lmkn ( t i ) e i [ mϕ − Φ mkn ( t i )] , (4.3)where H lmkn ( t i ) = A lmkn ( t i ) S lm [ ϑ ; aω mkn ( t i )] (4.4)and A lmkn ( t i ) = − Z ∞ lmkn ( t i ) ω mkn ( t i ) . (4.5)Equation (4.3) showcases the “multivoice” structure ofEMRI waveforms: each h lmkn that contributes to h ( t i )constitutes a single “voice” in this waveform. As we willsee in Secs. VII and VIII, even when the waveform iscomplicated, each voice tends to be quite simple.As discussed in Sec. III, we can compute the ampli-tudes ˇ A lmkn for the fiducial geodesic and then correctusing the phase factor ξ mkn in order to get the wave am-plitude for our system’s initial conditions. Let us imag-ine we have computed ˇ A lmkn on a dense grid of orbits.Knowing [ p ( t i , e ( t i ) , x I ( t i )], it is then simple to constructthe fiducial amplitudes ˇ A lmkn ( t i ) along an inspiral.To convert from the fiducial amplitudes to A lmkn ( t i ),we need the phase factor ξ mkn ( t i ) at each moment alongthe inspiral. Recall that ξ mkn depends on the angles χ θ and χ r which set the polar and radial initial con-ditions. An important feature of the adiabatic approx-imation is that these angles are constant : χ θ and χ r do not change as inspiral proceeds. However, the map-ping between ( χ θ , χ r ) and ξ mkn does change as inspiralproceeds, leading to a slow evolution in this phase fac-tor. When post-adiabatic physics is included, this storychanges: the angles χ θ and χ r slowly evolve underthe influence of orbit-averaged conservative self forces,and orbit-averaged spin-curvature interactions. This willchange the slow evolution of ξ mkn , and is one way inwhich conservative effects leave an observationally im-portant imprint on EMRI waveforms. V. FREQUENCY-DOMAIN DESCRIPTION
Our description so far has focused on presenting adia-batic EMRI waveforms in the time domain. We showedthat the time-domain waveform can be regarded as asuperposition of different “voices,” each of which hasits own slowly evolving amplitude. In this section, wewill exploit this multivoice structure to compute theFourier transform of an EMRI waveform, thereby de-scribing these waves in the frequency domain.Because all quantities in an EMRI evolve slowly (aslong as the two-timescale approximation is valid), ourexpectation is that the stationary phase approximation (SPA) will provide a high-quality approximation to the Fourier transform. For some voices, the frequency evolu-tion is not monotonic: some voices rises to a maximumfrequency, and then their frequency decreases. As we dis-cuss below, it is conceptually straightforward to extendthe “standard” SPA calculation in such a circumstance.We begin by reviewing the standard calculation, then dis-cuss voices whose frequency evolution is not monotonic.We conclude this section by describing how to assemblea frequency-domain waveform with many voices.
A. Standard SPA calculation
Begin by assuming a single voice signal of the form h ( t ) = H ( t ) e − i Φ( t ) . (5.1)The Fourier transform of this is given by˜ h ( f ) ≡ (cid:90) ∞−∞ h ( t ) e πift dt = (cid:90) ∞−∞ H ( t ) e i [2 πft − Φ( t )] dt . (5.2)To compute the stationary phase approximation to theFourier transform, expand the signal’s phase asΦ( t ) = Φ( t S ) + 2 πF ( t − t S ) + π ˙ F ( t − t S ) + . . . . (5.3)We will define the time t S momentarily. In Eq. (5.3), wehave introduced the signal’s instantaneous frequency andthe instantaneous frequency derivative at t = t S : F ≡ π d Φ dt (cid:12)(cid:12)(cid:12)(cid:12) t S , ˙ F = dFdt ≡ π d Φ dt (cid:12)(cid:12)(cid:12)(cid:12) t S . (5.4)We will assume that ˙ F is small, in a sense to be madeprecise below.Using these definitions, we rewrite the Fourier trans-form integral:˜ h ( f ) = e − i Φ( t S ) (cid:90) ∞−∞ H ( t ) e πi [ ft − F ( t − t S ) − (1 /
2) ˙ F ( t − t S ) ] dt = e i [2 πft S − Φ( t S )] × (cid:90) ∞−∞ H ( t (cid:48) + t S ) e πi [ ft (cid:48) − F t (cid:48) − (1 /
2) ˙ F ( t (cid:48) ) ] dt (cid:48) . (5.5)On the second line, we changed the integration variableto t (cid:48) = t − t S . The integrand of Eq. (5.5) very rapidlyoscillates unless the Fourier frequency f matches the in-stantaneous frequency F . When this condition is met, thephase is stationary : it is approximately constant, vary-ing very slowly due to the contribution from ˙ F , which isassumed to be small.We take the time t S to be the time at which the phaseis stationary. Let us rewrite it t ( f ), the time at which F ( t ) = f . Under the assumption that the largest contri-bution to the integral comes from t (cid:48) (cid:39)
0, or equivalently2from t (cid:39) t ( f ), we have˜ h ( f ) (cid:39) H [ t ( f )] e i [2 πft ( f ) − Φ( t ( f ))] (cid:90) ∞−∞ e − iπ ˙ F [ t ( f )]( t (cid:48) ) dt (cid:48) . (5.6)This integral can be evaluated with standard methods,and we finally obtain˜ h ( f ) (cid:39) H [ t ( f )] e − i [2 πft ( f ) − Φ( t ( f ))+ π/ (cid:113) ˙ F [ t ( f )] . (5.7)Note that ˙ F can be positive or negative, and it appearsunder the square root. To eliminate ambiguity about thephase of this voice in the frequency domain, we clean thisup as follows:˜ h ( f ) (cid:39) H [ t ( f )] e − i [2 πft ( f ) − Φ( t ( f )) ± π/ (cid:113) | ˙ F [ t ( f )] | . (5.8)We choose the plus sign in the phase if ˙ F >
0, and theminus sign for ˙
F <
0. This approximation works wellwhen both F ( t ) and H ( t ) change slowly,1 F dHdt (cid:28)
H , F dFdt (cid:28)
F . (5.9)It also requires that the signal frequency F evolves mono-tonically — the sign of dF/dt cannot change. B. Non-monotonic frequency
What if our signal has a frequency which does notevolve monotonically? In particular, what if F rises toa maximum and then decreases, or falls to a minimumand then increases? When this occurs, two problemsarise with the standard SPA analysis. First, in this cir-cumstance there are multiple solutions to the condition F ( t ) = f . The signal at frequency f must include contri-butions from all times from which the signal frequency F equals the Fourier frequency f . Second, ˙ F = 0 at the mo-ment that the evolution of F switches sign. The standardSPA Fourier transform is singular at that point. Theseissues affect EMRI signals, since the frequency associatedwith many voices rises to a maximum and then decreases.In particular, this occurs for EMRI voices which involve harmonics of the radial frequency: Ω r reaches a maxi-mum in the strong field, then goes to zero as the inspiralapproaches the last stable orbit.The root cause of the singularity is that the standardSPA assumes F ( t ) and ˙ F ( t ) completely describe the sig-nal’s phase. If ˙ F vanishes at frequency f , then the calcu-lation assumes all times contribute to the Fourier integralat f , and the integral (5.6) diverges. (This is consistentwith the fact that the Fourier transform of a constantfrequency signal is a delta function.) However, for realEMRI signals, F ( t ) is not constant when ˙ F = 0; thesingularity in the SPA analysis is an artifact of our as-sumption that F ( t ) and ˙ F ( t ) completely characterize thesignal’s phase. To remove this artifact, we need moreinformation about the signal’s phase evolution. Let ustherefore include the next term in the expansion:Φ( t ) = Φ( t S )+2 πF ( t − t S )+ π ˙ F ( t − t S ) + π F ( t − t S ) . . . (5.10)where ¨ F ≡ π d Φ dt (cid:12)(cid:12)(cid:12)(cid:12) t S . (5.11)Now revisit Eq. (5.2), but use (5.10) to expand thephase. We again find that t S is the stationary time, forwhich F = f . However, there may be multiple roots ofthe equation F ( t ) = f . Let us define t j ( f ) to be the j th time at which F ( t ) = f , and write ˙ F j ≡ ˙ F [ t j ( f )],¨ F j ≡ ¨ F [ t j ( f )]. For EMRI waveforms, N ≤
2. Each valueof j contributes to the Fourier transform, so we have˜ h ( f ) (cid:39) N (cid:88) j =1 H [ t j ( f )] e i [2 πft j ( f ) − Φ( t j ( f ))] × (cid:90) ∞−∞ e − iπ [ ˙ F j ( t (cid:48) ) + ¨ F j ( t (cid:48) ) / dt (cid:48) . (5.12)To perform this integral, put α = γ + 2 πi ˙ F , with γ realand positive. Define β = 2 π ¨ F , and use (cid:90) ∞−∞ e − αt / − iβt / dt = 2 √ α | β | e α / β K / ( α / β ) , (5.13)where K n ( z ) is the modified Bessel function of the 2ndkind. Taking the limit γ →
0, we find˜ h ( f ) (cid:39) √ N (cid:88) j =1 H [ t j ( f )] e i [2 πft j ( f ) − Φ( t j ( f ))] i ˙ F j | ¨ F j | e − πi ˙ F j / F j K / ( − πi ˙ F j / F j ) . (5.14)This result defines our “extended” SPA.It is useful to examine two limits of Eq. (5.14). To facilitate taking these limits, we define X j ≡ π F j ¨ F j . (5.15)3First, we take ˙ F j to be arbitrary, and expand about ¨ F j = 0. To set this up this, eliminate ¨ F j in Eq. (5.14):˜ h ( f ) (cid:39) (cid:114) π N (cid:88) j =1 H [ t j ( f )] e i [2 πft j ( f ) − Φ( t j ( f ))] ˙ F j (cid:113) ˙ F j i (cid:112) X j e − iX j K / ( − iX j ) . (5.16)Expanding about ¨ F j is equivalent to examining Eq. (5.16) for X j → ±∞ . Using˙ F j (cid:113) ˙ F j = | ˙ F j | − / ˙ F j > , = i | ˙ F j | − / ˙ F j < , (5.17)and, for X real, lim X →±∞ i √ Xe − iX K / ( − iX ) = (cid:114) π e − iπ/ (cid:18) − i
72 1 X + . . . (cid:19) , (5.18)we have ˜ h ( f ) = N (cid:88) j =1 H [ t j ( f )] e i [2 πft j ( f ) − Φ( t j ( f )) ∓ π/ | ˙ F j | / (cid:32) − i π ¨ F j ˙ F j + . . . (cid:33) . (5.19)The minus sign in the exponential is for ˙ F >
0, the plus sign for ˙
F <
0. Equation (5.19) is an accurate approximationwhen | ¨ F j | (cid:28) | ˙ F j | . Notice that we recover the standard SPA result, Eq. (5.8), when ¨ F j = 0 and N = 1.Next, allow ¨ F j to be any real value, and expand about ˙ F j = 0. To do this, we replace ˙ F j for X j in Eq. (5.14),˜ h ( f ) (cid:39) / / π / N (cid:88) j =1 H [ t j ( f )] e i [2 πft j ( f ) − Φ( t j ( f ))] ¨ F / j | ¨ F j | iX j / e − iX j K / ( − iX j ) , (5.20)and expand about X j = 0. Usinglim X → iX / e − iX K / ( − iX ) = e πi/ Γ( )2 / (cid:34) − X / e πi/ Γ( − )2 / Γ( ) (cid:35) , (5.21)we find ˜ h ( f ) (cid:39) e πi/ Γ( )3 / π / N (cid:88) j =1 H [ t j ( f )] e i [2 πft j ( f ) − Φ( t j ( f ))] ¨ F / j | ¨ F j | (cid:34) − e πi/ (cid:16) π (cid:17) / Γ( − )Γ( ) ˙ F j ¨ F / j + . . . (cid:35) . (5.22)Equation (5.22) is accurate when | ˙ F j | (cid:28) | ¨ F j | / . No-tice that this result is well behaved and finite at ˙ F j = 0,demonstrating that the extended SPA cures the singular-ity at points where the rate of change of the instantaneoussignal frequency passes through zero.Figure 1 illustrates how these elements come togetherfor a voice that showcases many of the features we havediscussed here. We show the l = m = 2, k = 0, n = 20voice for an equatorial Kerr inspiral with a = 0 . M , p init = 12 M , e init = 0 .
7, and mass ratio ε = 10 − . Thefull time-domain waveform and additional voices are dis-cussed for this case in more detail in Sec. VIII A.On the left-hand side of Fig. 1, we show this voice’stime-domain amplitude and the evolution of its fre- quency. The voice frequency increases until it reachesa maximum, then rapidly decreases, at least until the in-spiral ends several hundred M after reaching this maxi-mum. The Fourier transform, computed using Eq. (5.14)and shown in the right-hand panels, has two branches:the branch with ˙ f >
0, shown as the long-dashed curve inthe upper right-hand panel of Fig. 1; and the branch with˙ f <
0, shown as the short-dashed curve in this panel.Combining the two branches yields the solid curveshown in the lower right-hand panel of Fig. 1. We over-lay on this plot a discrete Fourier transform computedusing this time-domain voice; because of the large dy-namic range in the signal’s amplitude, we consider sep-arately a low-frequency DFT (focusing on data for t i (cid:46) FIG. 1. An example of the time-domain and frequency-domain structure of one voice: l = m = 2, k = 0, n = 20, for anequatorial Kerr inspiral with a = 0 . M , p init = 12 M , e init = 0 .
7, and mass ratio ε = 10 − . On the left, we show this voice’samplitude (top) and frequency (bottom) as a function of time over the inspiral. The voice’s frequency increases until it reaches f max (cid:39) . /M ; it then reverses and falls to f final (cid:39) . /M at the end of inspiral. The spiky features occur when thisamplitude passes through zero (note the log scale). On the right, we show the frequency domain structure computed with theextended SPA. In the top right, we examine contributions to the SPA along the two branches of F ( t ) = f . The curve with longdashes is the magnitude of the SPA for the branch with ˙ f >
0; the curve with short dashes is for the branch with ˙ f <
0. Bottomright shows the final SPA Fourier transform obtained by summing contributions from the two branches, and compares this to adiscrete Fourier transform (DFT). Because of the large dynamic range between the low-frequency and high-frequency behaviorof this voice, we separately examine the DFT for the early time, low-frequency portion (corresponding to t i (cid:46) . × M ;plotted with magenta dashes) and the DFT for the late, high-frequency portion ( t i (cid:38) . × M ; plotted with blue dots).Modulo lobes at the boundaries of the two DFTs (associated with the Tukey window used to taper the time-domain signal nearthese boundaries), the DFTs coincide perfectly with the SPA. Notice the interesting behavior in the band 0 . (cid:46) Mf (cid:46) . . × M , for which f (cid:46) . M ) and a high-frequencyDFT (focusing on data for t i (cid:38) . × M ). We use aTukey window of width 500 M to taper these segmentsof the time-domain signal. Aside from lobes near theboundaries associated with these windows, we see perfectagreement between the DFT and the SPA. Notice the in-teresting structure in the band 0 . (cid:46) M f (cid:46) . f > f <
0: at each frequency inthis band, the signal contributes at two different times,and with two different phases. Additional examples ofvoices, in both the time and frequency domains, areshown in Secs. VII and VIII.
C. Multiple voices
Generalization to a multivoice signal is straightfor-ward. Let us write our signal in the time domain h ( t ) = (cid:88) V H V ( t ) e − i Φ V ( t ) , (5.23) where V labels the voice, and is shorthand for all themode indices which describe each voice of the waveform.For generic EMRIs, V ≡ ( l, m, k, n ).The calculation proceeds as before, but now the phaseis stationary for each voice at some moment. Define F V = 12 π d Φ V dt , (5.24)and likewise define ˙ F V , ¨ F V . Assume that the condition F V ( t ) = f has N V solutions for voice V , and define t j,V ( f ) to be the j th such solution. (As stated in theprevious section, N V will be either 1 or 2 for EMRIs.)The stationary phase Fourier transform is then5˜ h ( f ) = 2 √ (cid:88) V N V (cid:88) j =1 H V [ t j,V ( f )] e i [2 πft j,V ( f ) − Φ V ( t j,V ( f ))] i ˙ F j,V | ¨ F j,V | exp (cid:34) − πi F j,V ¨ F j,V (cid:35) K / (cid:34) − πi F j,V ¨ F j,V (cid:35) ≡ (cid:88) V ˜ h lmkn ( f ) . (5.25)We have introduced ˙ F j,V ≡ ˙ F V [ t j,V ( f )] and ¨ F j,V ≡ ¨ F V [ t j,V ( f )]. Depending on the relative values of vari-ous powers of ˙ F j,V and ¨ F j,V , one can expand the V thvoice as in Eqs. (5.19) or (5.22). VI. IMPLEMENTATION
In this section, we describe various technical details bywhich we implement this formalism for computing EMRIwaveforms. To make an adiabatic inspiral and its associ-ated waveform, we lay out a grid of orbits, parameterizedby each orbit’s ( p, e, x I ). We store all the data at eachgrid point needed to construct the inspiral and the wave-form. We then interpolate to estimate the values of eachdatum at locations away from the grid points. In thissection, we describe this data grid and the data whichare stored on it, and details of how we interpolate dataoff grid. We emphasize that there is surely a great deal ofroom to improve on the techniques we present here; in-deed, we used different algorithms to design our data gridand to perform interpolation in a closely-related compan-ion analysis, Ref. [12]. For this initial study, the gridsand interpolation techniques we use are chosen for easeof use. In later work, we plan to investigate how best tooptimize the grids and interpolation methods for speedand accuracy of waveform calculation. A. EMRI data grids
We store our data on a grid that is rectangular in p − p LSO , e , and x I , where p LSO parameterizes the laststable orbit (LSO). The value of p LSO is easy to calculateas a function of e and x I [45], making it simple to set upa grid in this space. We set our innermost grid point to p min = p LSO + 0 . M , slightly outside of the LSO. Theradial frequency Ω r → p → p LSO , which means thatFourier expansions in Ω r tend to be badly behaved as theLSO is approached; this can be regarded as a precursorto the small body’s plunge into the black hole [46, 47]at the end of inspiral. Very little inspiral remains whenthe small body has reached our choice of p min , so we areconfident that the error incurred by truncating at p min (rather than closer to p LSO ) is negligible. That said, itshould be emphasized that this choice of p min has notbeen carefully evaluated. It will be useful to systemati-cally examine how to select the grid’s inner edge in future“production quality” work. Many important quantities vary rapidly near the LSO.Of particular importance is the phase of ˇ A lmkn , whichtends to rotate rapidly as p → p LSO . It is crucial toresolve this behavior in order to compute accurate wave-forms. To account for this behavior, we use a grid whosedensity increases near p LSO . For this paper, our grid isuniformly spaced in u ≡ √ p − . p LSO . (6.1)Using this spacing, we have laid down 40 points between p min = p LSO + 0 . M and p max = p min + 10 M . Differ-ent choices certainly could be used; for example, a gridreaching to larger p and with a different algorithm forincreasing density near p LSO was used in Ref. [12]. Thechoice of p spacing is an example of an issue that shouldbe more carefully investigated, and perhaps empiricallydesigned depending on what works best given computingresources and accuracy needs for one’s application.For any X stored on our grid, dX/de → e →
0. Empirically, we find it is important to have densegrid coverage for small e in order for this behavior to beaccurately captured. For this paper, we have used gridsthat run over 0 ≤ e ≤ . e = 0 .
1. Workin progress [48] suggests that this spacing may introducesmall systematic errors in computing the inspiral rateas a function of initial eccentricity; using higher densityacross at least part the small e part of this range appearsto effectively address this. More detailed discussion ofthis point will be presented in later work [48].For fixed p − p LSO and fixed e , all our stored data tendsto be very smooth and indeed nearly linear as a functionof x I . Our grid covers the range − ≤ x I ≤
1, with 16points spaced by ∆ x I = 2 / (cid:39) . • The rates of change ( dE/dt ) ∞ , H , ( dL z /dt ) ∞ , H ,and ( dQ/dt ) ∞ , H obtained by summing over manymodes until convergence has been reached, usingthe convergence criteria described in Sec. III C. • The rates of change of the orbital elements( dp/dt ) ∞ , H , ( de/dt ) ∞ , H , and ( dx I /dt ) ∞ , H obtainedby using the Jacobian described in App. B withthese fluxes. • The fiducial amplitude ˇ A lmkn for all modes used tocompute ( dE/dt ) ∞ , ( dL z /dt ) ∞ , and ( dQ/dt ) ∞ .So far, we have constructed such data sets for spheri-cal ( e = 0) and equatorial orbits ( x I = ±
1) for a/M ∈ , . , . , . . . , . , . , . a = 0 . M covering 0 ≤ e ≤ .
4. These data areproduced using the code GREMLIN, a frequency-domainTeukolsky solver primarily developed by author Hughes,with significant input from collaborators. The core meth-ods of this code are described in Refs. [29, 34], updatedto use methods developed in Refs. [36, 37] for solvingthe homogeneous Teukolsky equation; see [49] for furtherdiscussion.Two-dimensional orbits (eccentric and equatorial, orspherical and inclined) typically require several hun-dred to several thousand modes in order to convergeas described in Sec. III C, reaching ∼ modes in thestrongest fields. The number of modes needed for conver-gence of generic modes is an order of magnitude or twolarger. Each spherical mode requires about 0.1 secondsof CPU time for small values of l , increasing to roughly0.25 seconds for modes with l = 10. Computing eccentricmodes is more time consuming, since each mode involvesan integral over the radial domain covered by its orbit.This cost also varies significantly by radial mode number n . For small l , modes for small eccentricity ( e (cid:46) .
2) takeon average 1 CPU second or less; medium eccentricitymodes ( e ≈ .
5) average about 5–10 CPU seconds; andlarge eccentricity modes ( e = 0 .
8) average 30–40 CPUseconds each. These averages are skewed significantly bylarger values of n , for which the integrand of Eq. (3.13)rapidly oscillates, and the integral tends to be small com-pared to the magnitude of the integrand. At l = 10, thesetimes increase: small eccentricity modes take on averageup to 20 CPU seconds; medium eccentricity modes av-erage about 40–50 CPU seconds; and large eccentricitymodes require on average 150 seconds.The CPU cost per mode is ameliorated by the fact thateach orbit ( p, e, x I ) and each mode ( l, m, k, n ) is indepen-dent of all others. As such, this problem is embarrass-ingly parallelizable, and data sets can be effectively gen-erated on distributed computing clusters. The data setsdescribed above are publicly available through the BlackHole Perturbation Toolkit [50]. Plans to extend thesesets, develop further examples, and release the GREM-LIN code, are described in Sec. IX. B. Interpolating and integrating across the grid
To find data away from the grid points, in this anal-ysis we use cubic spline interpolation in the three direc-tions. Because our grid is rectangular in ( p − p LSO , e, x I ),this can be implemented effectively, and is adequate fordemonstrating how to build adiabatic waveforms andillustrating the results. Cubic spline interpolation forthe individual mode amplitudes will not scale well to“production-level” code, in terms of computational ef-ficiency and memory considerations. In Ref. [12], a setof the present authors used reduced-order methods withmachine learning techniques to construct a global fit tothe set of mode amplitudes, finding outstanding efficiency gains in an initial study of Schwarzschild EMRI wave-forms. Future work will apply these techniques to themore generic conditions we examine here.Other data needed to construct the waveform (for ex-ample, the geodesic frequencies Ω r,θ,φ and the phases ξ mkn ) are calculated at each osculating geodesic as in-spiral proceeds. Substantial computing speed could begained by storing and interpolating such data; indeed,all such data tends to evolve smoothly and fairly slowlyover an inspiral, so it is likely that effective interpola-tion could be implemented. We leave an investigation ofwhat to store and interpolate versus what to compute forfuture work and future optimization.To make an adiabatic inspiral, we constructa sequence of osculating geodesics, parameter-ized by [ p ( t i ) , e ( t i ) , x I ( t i )] given initial conditions[ p init , e init , x I, init ]. To do this, we interpolate to find dp/dt , de/dt , and dx I /dt at each inspiral time t i , thenuse a fixed-stepsize 4th-order Runge–Kutta integratorto construct the inspiral. This simple method is anobvious point for improvement in future work; indeed,in the related analysis [12], we used a variable-stepsize8th-order Runge–Kutta integrator.Note that both ( dp/dt ) and ( de/dt ) are singular as theLSO is approached. This is because the Jacobian (App.B) relating these quantities has zero determinant at theLSO. This singularity can adversely affect the accuracy ofinterpolating these quantities. Multiplying both ( dp/dt )and ( de/dt ) by ( p − p LSO ) , which exactly describes thesingularity in the zero eccentricity limit and approxi-mately describes it in general, yields smooth data whichinterpolate very well. Another approach is to interpo-late the fluxes dE/dt , dL z /dt , and dQ/dt , normalized bypost-Newtonian expansions of these quantities in orderto “divide out” their most rapid variations across orbits.This would construct a sequence of osculating geodesicparameterized as [ E ( t i , L z ( t i ) , Q ( t i )]. Determining themost efficient way to parameterize these orbits in orderto optimize waveform construction for speed and accu-racy are very natural directions for future work. VII. RESULTS I: EXAMPLE SCHWARZSCHILDEMRI WAVEFORMS AND THEIR VOICES
We now present examples of EMRI waveforms andtheir voices in the time and frequency domains. Our goalis not to exhaustively catalog EMRI waveforms, but justto present examples which showcase the behavior that wefind, and how this behavior tends to correlate with sourceproperties. We begin here with results for Schwarzschild;the next section shows Kerr results.Thanks to spherical symmetry, Schwarzschild orbitsare confined to a plane, which we define as equatorial.We can thus set Q = 0 and focus on voices with k = 0(i.e., neglecting harmonics of the θ motion). We examinetwo cases: one starts at ( p init , e init ) = (12 M, .
2) and in-spirals to e final (cid:39) . p init , e init ) =7(12 M, .
7) and inspirals to e final (cid:39) . p final follows from the Schwarzschild last stable orbit: p LSO = (6 + 2 e ) M . We use mass ratio ε = 10 − in all thecases we examine, both here and in the following section.Results for other extreme mass ratios can be inferred byscaling durations and accumulated phases with 1 /ε . A. Waveforms in the time domain
The top two panels of Fig. 2 show h + for the twoSchwarzschild inspirals we examine, both with initialanomaly angle χ r = 2 π/
3. The waveform is shown inthe system’s equatorial plane, so h × = 0. For the casewith e init = 0 .
2, we plot contributions from all modeswith l ∈ [2 , , m ∈ [ − l, . . . , l ], n ∈ [0 , . . . ,
10] (as wellas modes simply related by symmetry); for the case with e init = 0 .
7, we use the same l and m range, but go over n ∈ [0 , . . . , ξ m n on the waveform, zooming in on early andlate times. The red curves include all ξ m n corrections,and the blue curves neglect them, showing inspirals madeusing the fiducial amplitudes ˇ A lm n . Both early and latein the inspiral, the phase correction has a noticeable in-fluence. This is not surprising, since the impact of ξ m n is to adjust the system’s initial conditions — differentchoices of ξ m n correspond to physically different inspi-rals. The influence on the large eccentricity case is par-ticularly strong.Figures 3 and 4 show how the phases ξ (top pan-els) and ξ (bottom panels) evolve over these inspi-rals. We show these phases for initial anomaly angles χ r = 0 (black curves), π/ π/ π/ π/ ξ and ξ are nearly flat over the inspiral,though they show significant variation in the very lastmoments. The variation is larger in the higher eccen-tricity case for ξ , changing by almost a radian overthe inspiral for χ r = π/ π/ ξ and ξ are smooth and well behaved. They are alsorelatively simple to calculate, only requiring informationabout the geodesic with parameters p and e . Since com-puting A lm n is an expensive operation, one should onlycompute the fiducial amplitudes ˇ A lm n and use the phase ξ m n = mξ + nξ to convert.To calibrate how well the phases ξ m n allow us to ac-count for initial conditions, we compare the waveform as-sembled voice-by-voice with one computed independentlyusing a time-domain Teukolsky equation solver. For thecomparison waveform, we compute the worldline followedby an inspiraling body, use it to build the source for the time-domain Teukolsky equation as described in Ref. [14],and then compute the waveform using the techniques de-veloped in Refs. [14, 15]. The time-domain solver projectsits output onto spherical ( l, m ) modes of spin-weight − l = 2, m = 2.Figures 5, 6, and 7 summarize the results that we findfor Schwarzschild inspiral with p init = 12 M , e init = 0 . χ r = 0. In this case, we find that the waveformassembled voice-by-voice and the time-domain compari-son remain in phase for the entire inspiral. In the figure,we compare the two waveforms for a stretch of duration∆ t i = 3000 M at the beginning of inspiral as well as astretch of duration ∆ t i = 1000 M in the middle (near t i = 1 . × M ). The two waveforms lie on top of eachother in both cases; for this choice of χ r we find anexcellent match all the way to the end at t i (cid:39) × M .Figure 6 shows the inspiral waveforms when we put χ r = 2 π/
3. In this case, we find a secular drift whichaccumulates as inspiral proceeds. The top panel is againa stretch ∆ t i = 3000 M from the beginning of inspiral; asin Fig. 5, the two computed waveforms lie on top of oneanother. However, by t i = 1 . × M , the two waveformsare about 2 radians out of phase, as can be seen in thelower panel of this figure. This mismatch grows to about4 or 5 radians by the end of the inspiral.As discussed in Sec. IV, our solution neglects the im-pact of “slow-time” derivatives on EMRI evolution, leav-ing out time derivatives of terms which vary on the in-spiral timescale T i . As such, our solution only solves theTeukolsky equation (3.2) up to errors of O ( ε ). The time-domain solver by contrast finds a solution which, up tonumerical discretization, solves Eq. (3.2) at all orders in T o /T i . Our hypothesis is that this phase offset is be-cause the time-domain solver captures at least some ofthe “slow-time” derivatives which are missed by the adi-abatic construction. Interestingly, the magnitude of theoffset depends strongly on χ r . By examining multiplevalues of χ r , we find that the effect varies at least ap-proximately as 1 − cos χ r . This suggests that a slow-timevariation in ξ m n may play a particularly important rolein this secular drift.To test the hypothesis that the offset is due to over-looked slow-time terms, we replaced the accumulatedphase, Eq. (4.2), with the following ad hoc modification:Φ mod m n ( t i ) = (cid:90) t i t (cid:20) m Ω φ (cid:18) − cos χ r ) 32 dp/dtp Ω φ (cid:19) + n Ω r (cid:18) − cos χ r ) 32 dp/dtp Ω r (cid:19)(cid:21) dt . (7.1)The factor of (1 − cos χ r ) in this modification accountsfor the empirical dependence on χ r that we found; thefactor of dp/dt connects this phase to the inspiral, and thefactors 1 /p and 1 / Ω r,φ provide dimensional consistency.The numerical factor 3 / FIG. 2. The waveform h + for inspiral into a Schwarzschild black hole with ( p init , e init , x I, init ) = (12 M, . ,
1) (left-hand panels)and ( p init , e init , x I, init ) = (12 M, . ,
1) (right-hand panels). The waves are observed in the plane of the orbit, so h × = 0, andthe mass ratio is ε = 10 − in all cases. The top panels show the complete inspiral waveform over this domain, with the initialradial anomaly angle χ r = 2 π/
3. The bottom panels compare h + including the phase corrections ξ m n (red curves) with h + neglecting these corrections (i.e., incorrectly using the fiducial amplitudes ˇ A lm n ; blue curves). In both cases, we compare earlytimes (the first 3000 M of inspiral) and late times (an interval of 1000 M near the end of inspiral). Neglecting ξ m n leads tonoticeable differences in the waveforms. Newtonian limit yields d Ω dt = − dpdt Ω p − e dedt Ω1 − e . (7.2)Our empirical phase modification appears to be consis-tent with a weak-field correction associated with the rateat which p changes due to inspiral.We strongly emphasize that Eq. (7.1) is completely adhoc , and has not been justified by any careful calculation.However, we find that it does surprisingly well improv-ing the match between the two calculations. Figure 7 isthe equivalent of Fig. 6, but with Eq. (7.1) used to com-pute the phase rather than Eq. (4.2). Notice that thetwo waveforms lie on top of one another, at least overthe domain shown here. As inspiral proceeds, our ad hoc fix becomes less accurate: we find a roughly 1 radian off-set between the two waveforms when t i (cid:39) . × M ,growing to several radians by the end of inspiral. We findnearly identically improved matches examining inspiralswith different values of χ r , and for different choices ofthe mass ratio ε . Interestingly, including terms in de/dt inspired by the weak-field rate of change of Ω, does nothelp, suggesting that the similarity to the weak-field for-mula may be a coincidence.This analysis indicates that the phase offset we findis consistent with a post-adiabatic effect, and thereforeis missed by construction when making adiabatic wave-forms. The surprising effectiveness of our ad hoc fix sug- gests it may not be too difficult to analytically model thisbehavior and improve these waveforms. B. Waveform voices
The waveforms shown and discussed in Sec. VII A arefairly complicated, especially for the high eccentricitycase. By contrast, the individual voices which contributeto these waveforms are very simple, evolving smoothlyand simply on the much longer inspiral timescale.Figures 8 and 9 show individual voices that contributefor the case with e init = 0 .
2. A handful of the voiceswe show look jagged in these figures due to how we havepresented the data: the amplitude passes through zero insome cases, so | H lm n | appears spiky on a log-linear plot.These zero passings correspond to moments when the in-stantaneous frequency F m n = (1 / π ) d Φ m n /dt changessign. In this low eccentricity case, the voices with l = 2, m = ± n = 0 have the largest amplitudes.Figure 10 shows some of the voices which contributefor the case with e init = 0 .
7. Again we see that thevoices’ amplitudes and phases are smooth and well be-haved. In contrast to the small eccentricity case, modeswith n = 0 do not dominate here. Of the voices we plot, n = 6 dominates at early times, though it falls belowseveral other voices as the inspiral ends. The voice with n = 3 starts out weakest, but becomes strongest roughly9 FIG. 3. The phase correction ξ (bottom panel) and ξ (top panel) for inspiral with p init = 12 M , e init = 0 . ε = 10 − . In both panels, the black line correspondsto χ r = 0, the red curves show χ r = π/
6, blue shows χ r = π/
2, magenta shows χ r = 11 π/
6, and green shows χ r = 3 π/
2. Notice that these curves show only gentle varia-tion until nearly the end, with many changing rapidly as theinspiral approaches the last stable orbit.FIG. 4. The same as Fig. 3, but for an inspiral with p init =12 M , e init = 0 .
7. These phases show more variation in thiscase than when e init = 0 .
2, although the curves still showonly gentle variations until the system approaches the end ofinspiral. FIG. 5. The waveform for inspiral into a Schwarzschild blackhole with p init = 12 M , e init = 0 . χ r = 0 at mass ra-tio ε = 10 − computed voice-by-voice using the methodsdescribed here (blue crosses), and computed using a time-domain Teukolsky equation solver for the worldline of an in-spiral with these initial conditions. Though only voices with l = 2, m = 2 are included, the multivoice data are otherwiseidentical to the χ r = 0 data shown in Fig. 2. Top panel showsa stretch ∆ t i = 3000 M near the beginning of inspiral; bottomshows a stretch ∆ t i = 1000 M near t i = 1 . × M , roughlythe middle of this inspiral. The two calculations agree per-fectly in these two snapshots; we in fact find that they holdthis agreement all the way to the end at t i (cid:39) × M . halfway through this inspiral.All of the voices we examine have this behavior: boththe amplitudes and phases of individual voices evolvesmoothly on the inspiral timescale T i . This property isshared by the voices’ frequency-domain behavior. Fol-lowing Sec. V, we compute the frequency-domain repre-sentation of the voices examined here; Figs. 11 and 12show our results. Again, we see that all the voices evolvesmoothly over the inspiral. The apparent spikiness insome cases (for example, the voice with l = m = 4, n = 4for e init = 0 .
2) is because this voice’s amplitude passesthrough zero, and we show its magnitude on a log scale.It’s also worth noting that the frequency range of differ-ent voices varies. This is because our analysis begins inthe time domain, and then transforms to the frequencydomain using the SPA. Different voices thus start at dif-ferent frequencies, and reach different frequencies at theend of inspiral.Both the time-domain and the frequency-domain rep-resentation of these voices can be computed quicklyand efficiently. Future work, particularly data-analysis-focused applications which compare EMRI waves to de-0
FIG. 6. Same as Fig. 5, but for χ r = 2 π/
3; the multivoicedata shown here are identical to the l = 2, m = 2 χ r =2 π/ t i (cid:39) . × M . Thisgrows to about 4 or 5 radians by the end of inspiral. Ourhypothesis is that the time-domain integrator includes slow-time contributions (i.e., contributions from source terms thatevolve on the longer inspiral timescale T i ) which are left outof the adiabatic waveform. tector noise, could benefit substantially by focusing uponthe waveform in voice-by-voice fashion, studying whichvoices are most relevant for detection and measurementas a function of source parameters. VIII. RESULTS II: EXAMPLE KERR EMRIWAVEFORMS AND THEIR VOICES
We next consider a few examples of inspiral into Kerrblack holes. Although certain important details differfrom the results discussed in Sec. VII, much of what wefind for Kerr inspiral waveforms is qualitatively quite sim-ilar to those we find for the Schwarzschild case. As such,we keep this discussion brief, focusing on the most impor-tant highlights of this analysis. Future work will explorethese waveforms and their properties in depth.We begin in Sec. VIII A with a discussion of two con-strained cases: one that initially has zero eccentricity, butis inclined with respect to the hole’s equatorial plane; anda second case that is equatorial, but starts with large ec-centricity. We then discuss one example of fully generic(inclined and eccentric) inspiral in Sec. VIII B.
FIG. 7. Same as Fig. 6, but using the ad hoc phase modi-fication introduced in Eq. (7.1). This correction adjusts theorbital phase with a term that depends on the inspiral rate.Though this correction has not been rigorously computed, itnonetheless substantially improves the agreement between thetwo waveforms, at least over this domain of t i . (As describedin the text, the waveforms drift from one another as inspi-ral continues.) This supports the hypothesis that this driftarises due to slow-time variations neglected in the adiabaticapproximation. A. Constrained orbital geometry: Spherical andequatorial inspirals
We begin our Kerr study with two cases of inspiralinto black holes with spin a = 0 . M . In the first case,we examine an orbit that is spherical, with large incli-nation to the black hole’s equatorial plane. We take( r init , x I, init ) = (10 M, . r = 3 . M , x I = 0 . x I de-creasing very slightly. A similarly slight decrease of x I isseen in all of the cases we have examined.The left-hand panels of Fig. 13 show the gravitationalwaveform we find in this case. For the time-domainwaveform, we plot contributions from all modes with l ∈ [2 , , m ∈ [ − l, . . . , l ], k ∈ [0 , . . . , n = 0 contribute.) Thestrong influence of spin-orbit modulation can be seen inthe lower left-hand panels, which zoom in on early andlate times. These lower panels also illustrate the role ofthe initial polar anomaly angle χ θ , contrasting the wave-form with χ θ = 2 π/ χ θ = 0. Thetwo waveforms are similar in shape but shifted, consis-tent with the fact that χ θ controls the system’s initial1 FIG. 8. Some of the voices which contribute to h + for the case ( p init , e init ) = (12 M, . | H lm n ( t i ) | is on the left; the phase Φ m n ( t i ) is on the right. Top panels show voices with l = 4 and m = 4, bottom panels show l = 2 and m = 2. Red curves are data with n = 0, blue is n = 1, green is n = 2, magentais n = 3, and black is n = 4. In all cases both the amplitude and the phase evolve smoothly and simply. As functions of time,these voices can be sampled much less densely than h + must be sampled in order to accurately track its behavior. conditions: when χ θ = 0, the orbit is at θ = θ min when t = 0; for χ θ , the orbit starts at a value of θ roughlymidway between the equator and θ max = π − θ min .The right-hand panels of Fig. 13 show several of thevoices with l = m = 2 which contribute to this wave-form. The trend we have found across all the cases weexamine is that voices with | k + m | = l are loudest, andfall off as | k + m | moves away from this peak. This trendcan be seen in the cases shown here: the strongest voicehas k = 0, followed by k = −
1. In the time domain,the voices with k = 2 and k = 1 have roughly the sameamplitude; the voice with k = − k = 2, k = 1,and k = − / (cid:112) ˙ F which enters the frequency-domain magnitude compensates for the fact that thisvoice’s time-domain amplitude is smaller by a factor of2 or 3. Such differences have important implications forthe measurability of these signals.The second case we examine is an equatorial eccentricinspiral, with ( p init , e init ) = (12 M, . a = 0 . M , these initial conditions lead toa much longer inspiral that goes very deep into the strongfield, becoming nearly circular before plunge: inspirallasts for ∆ t i (cid:39) . × M , roughly twice the durationof the high-eccentricity Schwarzschild inspiral, and endswith ( p final , e final ) = (2 . M, . l ∈ [2 , , m ∈ [ − l, . . . , l ], n ∈ [0 , . . . , k = 0 contributesince there is no θ motion. The early waveform is quali-tatively quite similar to the early large eccentricity wave-form we found for Schwarzschild (Fig. 2). The late wave-form, by contrast, is quite different, reflecting the factthat the orbit has nearly circularized as it approaches thefinal plunge. As in the Schwarzschild case, we see thatthe initial radial anomaly angle χ r has a large impacton the waveform.The right-hand panels of Fig. 14 show some of the l = m = 2 voices which contribute to this waveform.An interesting trend we see is that the importance of dif-ferent voices changes dramatically during inspiral. Theevolution of the n = 0 voice is especially dramatic: itis fairly weak at early times (and in fact passes throughzero during the inspiral), but dominates the waveform atlate times. This makes sense — at late times the system’sgeometry is nearly circular, and voices corresponding toradial harmonics play a substantially less important role.2 FIG. 9. Some voices with l = 2 and m = − h + for the small eccentricity case shown in Fig. 2.Colors and definitions are exactly as in Fig. 8. The behav-ior of modes with n = 3 and n = 4 are interesting: Duringthe inspiral, there are moments at which the phase evolutionreverses direction, corresponding to the voice’s instantaneousfrequency changing sign. The amplitudes corresponding tothese voices pass through zero at times very close to thesemoments. The zero passage of these amplitudes leads to thesharp appearance that we see here. The real and imaginaryparts of H − and H − are perfectly smooth. B. Generic Kerr inspiral
We conclude our discussion of results by looking ata Kerr inspiral that is both inclined from the equato-rial plane and eccentric. As discussed in Sec. VI, wehave not yet generated dense data sets covering a widerange of such orbits. The example shown here demon-strates that the techniques we have developed to buildadiabatic EMRI waveforms have no difficulty with suchcases. Although generic EMRI waveforms have been de-veloped using “kludges,” to our knowledge this is the firstgeneric example that uses strong-field backreaction andstrong-field wave generation for the entire calculation.The case we examine begins at ( p init , e init , x I, init ) =(12 M, . , . ε = 10 − , inspirallasts for ∆ t i (cid:39) . × M , at which time the smallerbody encounters the LSO at ( p final , e final , x I, final ) =(4 . M, . , . δx I = 0 .
012 corresponds to the inclination angle I increasing by about 0 . ◦ . Figure 15 shows the trajectorythat the smaller body follows in ( p, e, x I ). The left-handpanels of Fig. 16 shows the time-domain +-polarizationof the waveform that we find in this case, including all FIG. 10. Some voices with l = 2, m = 2 which contribute tothe large eccentricity waveform shown in Fig. 2. Top panelshows the voices’ phase, lower panel their amplitudes. Redcurves are data with n = 0, blue is n = 3, green is n = 6,magenta is n = 9, and black is n = 12. In contrast to the smalleccentricity case, the voice with n = 0 does not dominateover most of the inspiral. In fact, of the voices shown, theone which starts weakest ( n = 3) evolves to be become theloudest by the end of inspiral. Nonetheless, all amplitudesand phases evolve in a smooth and simple way, just as in thelow eccentricity case. modes with l ∈ [2 , , m ∈ [ − l, . . . , l ], k ∈ [0 , . . . , n ∈ [0 , . . . ,
25] (as well as modes simply related to thismodes by symmetry). The right-hand panel of Fig. 16and both panels of Fig. 17 show some of the voices thatcontribute to this waveform ( l = m = 2, k = 0 voices inFig. 16; l = m = 2, k = ± k = 0, the time-domain structure of voices thatcontribute to this waveform is similar to what we saw inthe low-eccentricity Schwarzschild case shown in Fig. 8: n = 0 is the strongest voice, with contributions steadilydecreasing as n increases. For k = ±
2, we see more varia-tion: n = 0 often starts out relatively weak, but becomes3 FIG. 11. Frequency domain representation of the voices shown in Figures 8 and 9. Top left shows voices with l = m = 2;bottom left shows voices with l = m = 4; both panels on the right show voices with l = 2, m = −
2. In all cases, the red curvesare the frequency-domain amplitudes with n = 0; blue is for n = 1; green is n = 2; magenta is n = 3; and black is n = 4. (Weuse two panels for the voices with m = − much stronger late in the inspiral after eccentricity issignificantly decreased. The frequency domain structureof the voices is consistent with this picture; given thealready large number of plots in this paper, we do notinclude figures showing this. FIG. 12. Frequency domain representation of the voicesshown in Fig. 10. All curves show voices with l = m = 2.The red curve is the frequency-domain amplitude with n = 0;blue is for n = 3; green for n = 6; magenta for n = 9; andblack for n = 12. IX. CONCLUSION
In this paper, we have shown how to to use precom-puted frequency-domain Teukolsky equation solutions tomake EMRI waveforms. In this way, one can computeand store in advance the most computationally expen-sive aspects of EMRI analysis, and then use the meth-ods we describe to rapidly assemble waveforms using thestored data. We particularly emphasize the usefulness ofthe multivoice waveform structure, which facilitates iden-tifying particularly important waveform multipoles andharmonics, both in the time and frequency domains.The framework and techniques we describe here havenot been optimized, so there is much scope for effi-ciency gains. In a companion analysis [12], some ofthe present authors examined several algorithmic im-provements, showing that we can in fact constructanalysis-length time-domain EMRI waveforms in theSchwarzschild limit in under a second. There are twoways that the results of the present work can be incorpo-rated into the framework of [12]. The first is to extendto Kerr inspirals. The main challenges here stem fromthe higher dimensional parameter space, and the need tosum over more modes, since the waveforms depend on anadditional harmonic frequency. The inspirals also extenddeeper into the strong field, so more modes are likely tomake strong contributions to the waveform. The frame-work in [12] is designed to overcome these challenges. Itsneural network interpolation is a promising technique for4
FIG. 13. Left panels: the waveform h + for an example of spherical inspiral into a black hole with a = 0 . M . This examplehas ( r init , x I, init ) = (10 M, . ε = 10 − , and is observed in the black hole’s equatorial plane. During inspiral,the orbit’s radius steadily shrinks while x I decreases very slightly (see text for details). The polar anomaly angle χ θ = 2 π/ χ θ = 2 π/ l = m = 2 voices which contributeto this waveform, both time domain (top) and frequency domain (bottom). Each line is a different k index: red is k = 0, blueis k = 1, green is k = −
1, magenta is k = 2, and black is k = −
2. Across the examples of zero eccentricity inspiral we haveexamined, the tendency is that voices with | k + m | = l are the strongest, and fall off as | k + m | moves away from this peak. dealing with the increased dimensionality, and its use ofGPU-based hardware acceleration alleviates the compu-tational burden of summing over thousands of additional( l, m, k, n )-modes. The second extension is to incorpo-rate the calculation of frequency-domain waveforms. Forthe Schwarzschild limit, Ref. [12] already demonstrateshow to efficiently interpolate waveform amplitudes; theremaining new challenge is the efficient calculation of thestationary time as a function of frequency. This canlikely be achieved by sparsely evaluating t ( f ) at a fewkey frequencies and interpolating to compute additionalvalues. For both extensions, though the details will com-plicate the development of fast waveforms, they do notchange the fundamental message of Ref. [12]. We areconfident that adapting those methods with the data setsand framework described here will be quite effective.Because of the high cost of generating the wave-form amplitude data, once computed, the data shouldbe widely shared. The data sets we have computed,described in Sec. VI, have been released through theBlack Hole Perturbation Toolkit [50] (hereafter “theToolkit”). These data were all computed using a fairlysmall (roughly 1000 core) in-house cluster at the MITKavli Institute. This scale of cluster is adequate formaking inspiral data for 2-D orbits (spherical or equa- torial), but is too small to be useful for fully generic3-D (inclined and eccentric Kerr) data sets. We haveported our frequency-domain Teukolsky equation solver,GREMLIN, to the US National Science Foundation’sXSEDE [51] environment, and plan to develop genericdata sets there. It should be noted, however, that thereis a lot to be learned from 2-D cases about how bestto lay out orbits on the grid (for example, the work inprogress mentioned in Sec. VI that shows the importanceof dense coverage at small eccentricity). Data sets re-leased through the Toolkit are likely to be revised withsome regularity as we learn more about the best way tolay out the data grids.Plans are in place to release the GREMLIN code whichwas used to generate these data sets, and which includestools for post-processing analysis and generating EMRIwaveforms. A reduced functionality version of GREM-LIN (specialized to circular, equatorial orbits) has al-ready been released. The generic version will be releasedafter it is cleaned of proprietary libraries and fully com-ports with Open Source licensing requirements.In terms of analyses of the source physics, there areseveral natural places to extend what we have done here.One clear step is to consider how physics beyond the adi-abatic approximation affects waveforms. As described in5 FIG. 14. Left panels: the waveform h + for an example of eccentric equatorial inspiral into a black hole with a = 0 . M . Thisexample has ( p init , e init ) = (12 M, . ε = 10 − , and is observed in the black hole’s equatorial plane. Theinitial conditions are similar to those used to make the large eccentricity case shown in Fig. 2, and the two waveforms aresimilar at early times. However, at this rapid spin, inspiral goes very deep into the strong field, becoming nearly circular beforeplunge; we find ( p final , e final ) = (2 . M, . χ r = 2 / χ r = 2 π/ l = m = 2 voices which contribute to this waveform, both time domain (top) and frequency domain(bottom). Red curves show the voices with n = 0; blue is for n = 3; green is n = 6; magenta is n = 9; and black is n = 12. Sec. VII, we see evidence from comparison with wave-forms computed by a time-domain Teukolsky solver thatneglecting terms which vary on the long inspiral timescaleleads to an initial-condition-dependent secular drift inthe waveform. Our hypothesis is that this arises due toa slow evolution in the ξ mkn phases described in Sec. III.In a similar vein, we expect it will not be difficult to in-corporate the orbit-averaged first-order conservative selfforce. In brief, on average the conservative self forcechanges the rate at which the orbit precesses, and canbe modeled as a slow “anomalous” change to the rate ofperiastron advance and the advance of the line of nodes.These anomalous changes can in turn be incorporatedinto an osculating geodesic framework by allowing theparameters χ r , χ θ , and φ to slowly evolve under theinfluence of this force. See, for example, Ref. [43] for fur-ther discussion. A similar effect due to the orbit-averagedcoupling of the smaller body’s spin to the backgroundcurvature probably can also be modeled in such a way.It should also be possible to include non -orbit-averagedself forces and spin-curvature couplings [18–21, 52] byusing a near-identify transformation [22]. Such simpleextensions to this framework may make it possible to in-corporate the most important post-adiabatic effects quiteeasily, significantly improving the ability of these modelsto serve as templates for EMRI measurements. Finally, we note that recent phenomenologicalfrequency-domain gravitational wave models are beingcalibrated in the large mass-ratio limit with Teukolskywaveforms [53]. The directly constructed frequency-domain waveforms presented in this work could be usefulin this effort. ACKNOWLEDGMENTS
An early presentation of some of the results and meth-ods described here was given at the 22nd Capra Meet-ing on Radiation Reaction in General Relativity, hostedby the Centro Brasiliero de Pequisas F´ısicas in Rio deJaneiro, Brazil in June 2019; we thank the partici-pants of that meeting for useful discussions, and par-ticularly Marc Casals for organizing the meeting. Thetime-domain waveform computations were performed onthe MIT/IBM
Satori
GPU supercomputer supported bythe Massachusetts Green High Performance Comput-ing Center (MGHPCC). SAH thanks the MIT KavliInstitute for providing computing resources and sup-port; L. S. Finn for (long ago) suggesting the multivoiceframework for examining EMRI waveforms; StanislavBabak for emphasizing the importance of understand-ing EMRI waveform structure in the frequency domain;6
FIG. 15. The inclined and eccentric trajectory in ( p, e, x I ) for the generic inspiral we examine. This trajectory, illustratedby the blue curve, begins at ( p init , e init , x I, init ) = (12 M, . , .
5) and proceeds until the smaller body encounters the LSO(a section of which is illustrated by the orange plane) at ( p final , e final , x I, final ) = (4 . M, . , . p, e ) plane is illustrated by the red curve, along with the projection of the LSO at the final value of x I (black line in lower plane); a projection of this trajectory into the ( p, x I ) plane is illustrated by the green curve, along with theprojection of the LSO at the final value of e (black line on back “wall” of the box).FIG. 16. Left panel: the waveform h + for an example of generic inspiral into a black hole with a = 0 . M . This examplecorresponds to ( p init , e init , x I, init ) = (12 M, . , . ε = 10 − ; the small body inspirals until it encountersthe LSO at ( p final , e final , x I, final ) = (4 . M, . , . χ r = χ θ = 2 π/ l = m = 2, k = 0; red curve is n = 0, blue is n = 1, green is n = 2, magenta is n = 3, and black is n = 4. The behavior ofthese voices is quite simple, with n = 0 strongest, and the voices falling away as n increases. FIG. 17. Additional voices that contribute to the waveform shown in Fig. 16. Left panel shows voices with k = 2; right panelshows voices with k = −
2. As in the right-hand panel of Fig. 16, the red curve is n = 0, blue is n = 1, green is n = 2, magentais n =, and black is n = 4. The behavior is somewhat more complicated for these voices; although larger values of n tend tobe weaker, the evolution for small values of n evolves rather differently than in the k = 0 case. For both the cases shown here,the n = 0 voice in particular evolves in importance. Steve Drasco and William T. Throwe for past contribu-tions to GREMLIN development; Gustavo Velez for in-sight into the Jacobian between ( dE/dt, dL z /dt, dQ/dt )and ( dp/dt, de/dt, dx I /dt ); Talya Klinger for work port-ing GREMLIN to the XSEDE computing environment;and both Talya Klinger and Sayak Datta for discus-sions of accuracy and systematic error in EMRI datasets. SAH’s work on this problem was supported byNASA ATP Grant 80NSSC18K1091 and NSF GrantPHY-1707549. NW acknowledges support from a RoyalSociety–Science Foundation Ireland Research Fellowship.This publication has emanated from research conductedwith the financial support of Science Foundation Irelandunder Grant number 16/RS-URF/3428. GK acknowl-edges support from NSF grants PHY-2106755 and DMS-1912716. AJKC acknowledges support from the NASAgrant 18-LPS18-0027. MLK acknowledges support fromNSF grant DGE-0948017. Appendix A: Formulas for E , L z , and Q In this appendix, we list formulas for the geodesic con-stants of the motion E , L z , and Q as functions of ourpreferred parameters p , e , and x I ≡ cos I . Formulas forthese constants were first worked out by Schmidt [25],and are provided in particularly clean form for genericorbits by van de Meent [28]. We list them here using ourpreferred parameterization, and also include formulas for spherical orbits.The three geodesic constants are given by E = (cid:115) κρ + 2 (cid:36)σ − x I ) (cid:112) σ ( σ(cid:36) + ρ(cid:36)κ − ηκ ) ρ + 4 ησ , (A1) L z = − g ( r a ) E − (cid:112) g ( r a ) + h ( r a ) f ( r a ) E − h ( r a ) d ( r a ) h ( r a ) , (A2) Q = (1 − x I ) (cid:20) a (1 − E ) + L z x I (cid:21) , (A3)where r a = p/ (1 − e ) is the coordinate radius of the orbit’sapoapsis. For generic orbits, the quantities appearinghere are given by κ = d ( r a ) h ( r p ) − d ( r p ) h ( r a ) , (A4) (cid:36) = d ( r a ) g ( r p ) − d ( r p ) g ( r a ) , (A5) ρ = f ( r a ) h ( r p ) − f ( r p ) h ( r a ) , (A6) η = f ( r a ) g ( r p ) − f ( r p ) g ( r a ) , (A7) σ = g ( r a ) h ( r p ) − g ( r p ) h ( r a ) , (A8)where r p = p/ (1+ e ) is the coordinate radius of the orbit’s8periapsis, and where d ( r ) = ∆( r )[ r + a (1 − x I )] , (A9) f ( r ) = r + a [ r ( r + 2 M ) + (1 − x I )∆( r )] , (A10) g ( r ) = 2 aM r , (A11) h ( r ) = r ( r − M ) + (1 − x I )∆( r ) x I , (A12)Note that we have switched notation slightly versus Refs.[25] and [28]: we use (cid:36) in these lists rather than (cid:15) to avoidnotational collision with (cid:15) = √ M − a / M r + , as wellas with the very similar ε ≡ µ/M .For spherical orbits with e = 0, r a = r p ≡ r o . In thislimit, we use a different form of these quantities: κ = d ( r o ) h (cid:48) ( r o ) − d (cid:48) ( r o ) h ( r o ) , (A13) (cid:36) = d ( r o ) g (cid:48) ( r o ) − d (cid:48) ( r o ) g ( r o ) , (A14) ρ = f ( r o ) h (cid:48) ( r o ) − f (cid:48) ( r o ) h ( r o ) , (A15) η = f ( r o ) g (cid:48) ( r o ) − f (cid:48) ( r o ) g ( r o ) , (A16) σ = g ( r o ) h (cid:48) ( r o ) − g (cid:48) ( r o ) h ( r o ) . (A17)In these formulas, (cid:48) denotes ∂/∂r . Appendix B: Jacobian from ( dE/dt, dL z /dt, dQ/dt ) to ( dp/dt, de/dt, dx I /dt ) Once the rates of change dE/dt , dL z /dt , and dQ/dt have been computed, we use them to compute how in the adiabatic limit a system evolves from one geodesic toanother. As part of this, we would like to know how thegeodesic geometry parameters p , e , and x I change dueto this backreaction. In this appendix, we write out thedetails of this procedure.A generic Kerr orbit has turning points in its radialmotion at apoapsis, r a = p/ (1 − e ), and at periapsis, r p = p/ (1 + e ). It also has a turning point in its po-lar motion at θ m . (The second polar turning point at π − θ m yields no new information because of reflectionsymmetry about θ = π/ θ m is related toour inclination angle I by Eq. (2.16). The radial turningpoints mean that R ( r a ) = 0 and R ( r p ) = 0, where R ( r )is defined in Eq. (2.1); the polar turning point meansthat Θ( θ m ) = 0, where Θ( θ ) is defined in Eq. (2.2). Werequire that these conditions hold as the system evolvesdue to backreaction, which means that we require ddt R ( r a ) = 0 , ddt R ( r p ) = 0 , ddt Θ( θ m ) = 0 . (B1)Expanding these total time derivatives, we find that theresults can be written as a matrix equation, J Er a J L z r a J Qr a J Er p J L z r p J Qr p J Ex I J L z x I J Qx I · dE/dtdL z /dtdQ/dt = dr a /dtdr p /dtdx I /dt (B2)Computing the various matrix elements J ab , we find sim-ple expressions for dr a , p /dt and dx I /dt . Examine thechange to the inclination first: dx I dt = (1 − x I ) 2 (cid:112) − x I (cid:112) Q − a (1 − E )(1 − x I )( dL z /dt ) − x I ( dQ/dt ) − x I (1 − x I ) a E ( dE/dt )2[ Q − a (1 − E )(1 − x I ) ] (B3)= 2 x I (1 − x I ) L z ( dL z /dt ) − x I [( dQ/dt ) + 2(1 − x I ) a E ( dE/dt )]2[ L z + x I a (1 − E )] . (B4)Notice that Eq. (B3) is singular when Q → | x I | → L z → | x I | → | x I | ≤ .
5, and(B4) when x I > .
5. The two expressions yield identicalresults except right at their singular points.To present our results for dr a , p /dt , we first define D ( r ) ≡ M [ Q + ( L z − aE ) ] − r [ L z + Q − a (1 − E )]+6 M r − r (1 − E ) . (B5)Using this, we have J Er a , p ≡ aM ( L z − aE ) r a , p − Er , p ( a + r , p ) D ( r a , p ) ,J L z r a , p ≡ M ( aE − L z ) r a , p + 2 L z r , p D ( r a , p ) ,J Qr a , p ≡ r , p − M r a , p + a D ( r a , p ) . (B6) We then find dr a , p dt = J Er a , p dEdt + J L z r a , p dL z dt + J Qr a , p dQdt , (B7)Once dr a , p /dt are known, it is simple to compute dp/dt and de/dt : dpdt = (1 − e ) dr a dt + (1 + e ) dr p dt , (B8) dedt = (1 − e )2 p (cid:20) (1 − e ) dr a dt − (1 + e ) dr p dt (cid:21) . (B9)In the spherical limit e = 0, r a = r p ≡ r o , and theconditions R ( r a ) = 0, R ( r p ) = 0 are redundant. Theequations we have derived here do not work in that limit.Spherical orbits are instead governed by the conditions R ( r o ) = 0, R (cid:48) ( r o ) = 0, where R (cid:48) ≡ ∂R/∂r . Past work[54, 55] long ago proved that adiabatic dissipative self9interaction evolves a spherical orbit into a new sphericalorbit; the first derivation of dQ/dt [40] proved that thisresult respected the “spherical goes to spherical” con-straint. Given ( dE/dt, dL z /dt, dQ/dt ) from a sphericalorbit, one can infer dr o /dt by enforcing the condition ddt R (cid:48) ( r o ) = 0 . (B10)(Alternately, one can enforce the conditions dR ( r o ) /dt =0, dR (cid:48) ( r o ) /dt = 0. Doing so yields solutions for dr o /dt and dQ/dt given dE/dt and dL z /dt . This is how back-reaction on spherical orbits was computed [34] before dQ/dt was fully understood [40].) Implementing Eq. (B10), we find dr o dt = J c Er o dEdt + J c L z r o dL z dt + J c Qr o dQdt , (B11)where the Jacobian elements for spherical orbits are J c Er = 2 aM ( L z − aE ) − a Er − Er D c ( r ) , (B12) J c L z r = − M ( L z − aE ) + 2 L z r D c ( r ) , (B13) J c Qr = r − M D c ( r ) , (B14)with D c ( r ) = − [ a (1 − E ) + L z + Q ] + 6 M r − − E ) r . (B15) [1] A. Le Tiec, The overlap of numerical relativity, pertur-bation theory and post-newtonian theory in the binaryblack hole problem, Int. J. Mod. Phys. D , 1430022(2014), arXiv:1408.5505 [gr-qc].[2] S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sop-uerta, C. P. Berry, E. Berti, P. Amaro-Seoane, A. Pe-titeau, and A. Klein, Science with the space-based inter-ferometer LISA. V: Extreme mass-ratio inspirals, Phys.Rev. D , 103012 (2017), arXiv:1703.09722 [gr-qc].[3] L. Barack and A. Pound, Self-force and radiation reac-tion in general relativity, Rept. Prog. Phys. , 016904(2019), arXiv:1805.10385 [gr-qc].[4] T. Hinderer and E. E. Flanagan, Two timescale analysisof extreme mass ratio inspirals in Kerr. I. Orbital motion,Phys. Rev. D , 064028 (2008), arXiv:0805.3337 [gr-qc].[5] M. van de Meent, Gravitational self-force on genericbound geodesics in Kerr spacetime, Phys. Rev. D ,104033 (2018), arXiv:1711.09607 [gr-qc].[6] A. Pound, B. Wardell, N. Warburton, and J. Miller,Second-order self-force calculation of gravitational bind-ing energy in compact binaries, Phys. Rev. Lett. ,021101 (2020), arXiv:1908.07419 [gr-qc].[7] S. A. Hughes, S. Drasco, E. E. Flanagan, and J. Franklin,Gravitational radiation reaction and inspiral waveformsin the adiabatic limit, Phys. Rev. Lett. , 221101 (2005),arXiv:gr-qc/0504015.[8] L. Barack and C. Cutler, LISA capture sources: Approx-imate waveforms, signal-to-noise ratios, and parameterestimation accuracy, Phys. Rev. D , 082005 (2004),arXiv:gr-qc/0310125.[9] S. Babak, H. Fang, J. R. Gair, K. Glampedakis, andS. A. Hughes, ’Kludge’ gravitational waveforms for atest-body orbiting a Kerr black hole, Phys. Rev. D ,024005 (2007), [Erratum: Phys.Rev.D 77, 04990 (2008)],arXiv:gr-qc/0607007.[10] A. J. Chua and J. R. Gair, Improved analytic extreme-mass-ratio inspiral model for scoping out eLISA dataanalysis, Class. Quant. Grav. , 232002 (2015),arXiv:1510.06245 [gr-qc].[11] S. A. Hughes, Evolution of circular, nonequatorial orbits of Kerr black holes due to gravitational wave emission. II.Inspiral trajectories and gravitational wave forms, Phys.Rev. D , 064004 (2001), [Erratum: Phys.Rev.D 88,109902 (2013)], arXiv:gr-qc/0104041.[12] A. J. Chua, M. L. Katz, N. Warburton, and S. A.Hughes, Rapid generation of fully relativistic extreme-mass-ratio-inspiral waveform templates for LISA dataanalysis, (2020), arXiv:2008.06071 [gr-qc].[13] S. Drasco, E. E. Flanagan, and S. A. Hughes, Computinginspirals in Kerr in the adiabatic regime. I. The scalarcase, Class. Quant. Grav. , S801 (2005), arXiv:gr-qc/0505075.[14] P. A. Sundararajan, G. Khanna, S. A. Hughes, andS. Drasco, Towards adiabatic waveforms for inspiral intoKerr black holes: II. Dynamical sources and generic or-bits, Phys. Rev. D , 024022 (2008), arXiv:0803.0317[gr-qc].[15] A. Zenginoglu and G. Khanna, Null infinity waveformsfrom extreme-mass-ratio inspirals in Kerr spacetime,Phys. Rev. X , 021017 (2011), arXiv:1108.1816 [gr-qc].[16] J. Miller and A. Pound, Two-timescale evolutionof extreme-mass-ratio inspirals: waveform generationscheme for quasicircular orbits in Schwarzschild space-time, (2020), arXiv:2006.11263 [gr-qc].[17] A. Pound and B. Wardell, Black hole perturbation the-ory and gravitational self-force, (2021), arXiv:2101.04592[gr-qc].[18] N. Warburton, T. Osburn, and C. R. Evans, Evolutionof small-mass-ratio binaries with a spinning secondary,Phys. Rev. D , 084057 (2017), arXiv:1708.03720 [gr-qc].[19] T. Osburn, N. Warburton, and C. R. Evans, Highly ec-centric inspirals into a black hole, Phys. Rev. D ,064024 (2016), arXiv:1511.01498 [gr-qc].[20] G. A. Piovano, A. Maselli, and P. Pani, Extreme massratio inspirals with spinning secondary: a detailed studyof equatorial circular motion, Phys. Rev. D , 024041(2020), arXiv:2004.02654 [gr-qc].[21] V. Skoup´y and G. Lukes-Gerakopoulos, Gravitationalwave templates from Extreme Mass Ratio Inspirals, (2021), arXiv:2101.04533 [gr-qc].[22] M. Van De Meent and N. Warburton, Fast self-forced inspirals, Class. Quant. Grav. , 144003 (2018),arXiv:1802.05281 [gr-qc].[23] R. Fujita and M. Shibata, Extreme mass ratio inspiralson the equatorial plane in the adiabatic order, Phys. Rev.D , 064005 (2020), arXiv:2008.13554 [gr-qc].[24] C. W. Misner, K. Thorne, and J. Wheeler, Gravitation (W. H. Freeman, San Francisco, 1973).[25] W. Schmidt, Celestial mechanics in Kerr space-time, Class. Quant. Grav. , 2743 (2002), arXiv:gr-qc/0202090.[26] S. Drasco and S. A. Hughes, Rotating black hole orbitfunctionals in the frequency domain, Phys. Rev. D ,044015 (2004), arXiv:astro-ph/0308479.[27] R. Fujita and W. Hikida, Analytical solutions of boundtimelike geodesic orbits in Kerr spacetime, Class. Quant.Grav. , 135002 (2009), arXiv:0906.1420 [gr-qc].[28] M. van de Meent, Analytic solutions for parallel transportalong generic bound geodesics in Kerr spacetime, Class.Quant. Grav. , 145007 (2020), arXiv:1906.05090 [gr-qc].[29] S. Drasco and S. A. Hughes, Gravitational wave snap-shots of generic extreme mass ratio inspirals, Phys.Rev. D , 024027 (2006), [Erratum: Phys.Rev.D 88,109905 (2013), Erratum: Phys.Rev.D 90, 109905 (2014)],arXiv:gr-qc/0509101.[30] S. A. Teukolsky, Perturbations of a rotating black hole. I.Fundamental equations for gravitational electromagneticand neutrino field perturbations, Astrophys. J. , 635(1973).[31] E. Newman and R. Penrose, An approach to gravitationalradiation by a method of spin coefficients, J. Math. Phys. , 566 (1962).[32] T. C. Quinn and R. M. Wald, Energy conservation forpoint particles undergoing radiation reaction, Phys. Rev.D , 064009 (1999), arXiv:gr-qc/9903014.[33] ˜A. E. Flanagan, S. A. Hughes, and U. Ruangsri,Resonantly enhanced and diminished strong-fieldgravitational-wave fluxes, Phys. Rev. D , 084028(2014), arXiv:1208.3906 [gr-qc].[34] S. A. Hughes, Evolution of circular, nonequatorial or-bits of Kerr black holes due to gravitational waveemission, Phys. Rev. D , 084004 (2000), [Erratum:Phys.Rev.D 63, 049902 (2001), Erratum: Phys.Rev.D 65,069902 (2002), Erratum: Phys.Rev.D 67, 089901 (2003),Erratum: Phys.Rev.D 78, 109902 (2008), Erratum:Phys.Rev.D 90, 109904 (2014)], arXiv:gr-qc/9910091.[35] S. Mano, H. Suzuki, and E. Takasugi, Analytic solutionsof the Teukolsky equation and their low frequency ex-pansions, Prog. Theor. Phys. , 1079 (1996), arXiv:gr-qc/9603020.[36] R. Fujita and H. Tagoshi, New numerical methods toevaluate homogeneous solutions of the Teukolsky equa-tion, Prog. Theor. Phys. , 415 (2004), arXiv:gr-qc/0410018.[37] R. Fujita and H. Tagoshi, New numerical methods toevaluate homogeneous solutions of the Teukolsky equa-tion II. Solutions of the continued fraction equation,Prog. Theor. Phys. , 1165 (2005), arXiv:0904.3818[gr-qc].[38] S. O’Sullivan and S. A. Hughes, Strong-field tidal dis- tortions of rotating black holes: Formalism and resultsfor circular, equatorial orbits, Phys. Rev. D , 124039(2014), [Erratum: Phys.Rev.D 91, 109901 (2015)],arXiv:1407.6983 [gr-qc].[39] S. Teukolsky and W. Press, Perturbations of a rotatingblack hole. III. Interaction of the hole with gravitationaland electromagnetic radiation, Astrophys. J. , 443(1974).[40] N. Sago, T. Tanaka, W. Hikida, K. Ganz, and H. Nakano,The adiabatic evolution of orbital parameters in the Kerrspacetime, Prog. Theor. Phys. , 873 (2006), arXiv:gr-qc/0511151.[41] C. Cutler, D. Kennefick, and E. Poisson, Gravita-tional radiation reaction for bound motion around aSchwarzschild black hole, Phys. Rev. D , 3816 (1994).[42] C. W. Lincoln and C. M. Will, Coalescing binary sys-tems of compact objects to (post)5/2 Newtonian order:Late time evolution and gravitational radiation emission,Phys. Rev. D , 1123 (1990).[43] A. Pound and E. Poisson, Osculating orbits inSchwarzschild spacetime, with an application to extrememass-ratio inspirals, Phys. Rev. D , 044013 (2008),arXiv:0708.3033 [gr-qc].[44] J. R. Gair, E. E. Flanagan, S. Drasco, T. Hinderer, andS. Babak, Forced motion near black holes, Phys. Rev. D , 044037 (2011), arXiv:1012.5111 [gr-qc].[45] L. C. Stein and N. Warburton, Location of the last sta-ble orbit in Kerr spacetime, Phys. Rev. D , 064007(2020), arXiv:1912.07609 [gr-qc].[46] P. A. Sundararajan, Transition from adiabatic inspiral togeodesic plunge for a compact object around a massiveKerr black hole: Generic orbits, Phys. Rev. D , 124050(2008), arXiv:0803.4482 [gr-qc].[47] A. Apte and S. A. Hughes, Exciting black hole modesvia misaligned coalescences: I. Inspiral, transition, andplunge trajectories using a generalized Ori-Thorne proce-dure, Phys. Rev. D , 084031 (2019), arXiv:1901.05901[gr-qc].[48] S. Datta et al. ,064013 (2020), arXiv:1912.09461 [gr-qc].[53] G. Pratten, S. Husa, C. Garcia-Quiros, M. Colleoni,A. Ramos-Buades, H. Estelles, and R. Jaume, Settingthe cornerstone for a family of models for gravitationalwaves from compact binaries: The dominant harmonicfor nonprecessing quasicircular black holes, Phys. Rev. D , 064001 (2020), arXiv:2001.11412 [gr-qc].[54] D. Kennefick and A. Ori, Radiation reaction inducedevolution of circular orbits of particles around Kerrblack holes, Phys. Rev. D , 4319 (1996), arXiv:gr-qc/9512018.[55] F. D. Ryan, Effect of gravitational radiation reaction onnonequatorial orbits around a Kerr black hole, Phys. Rev.D53