Improved gravitational radiation time-scales II: spin-orbit contributions and environmental perturbations
Lorenz Zwick, Pedro R. Capelo, Elisa Bortolas, Veronica Vazquez-Aceves, Lucio Mayer, Pau Amaro-Seoane
MMNRAS , 1–14 (2020) Preprint 2 February 2021 Compiled using MNRAS L A TEX style file v3.0
Improved gravitational radiation time-scales II: spin-orbitcontributions and environmental perturbations
Lorenz Zwick, (cid:63) Pedro R. Capelo, Elisa Bortolas, , Verónica Vázquez-Aceves, Lucio Mayer and Pau Amaro-Seoane , , , , Center for Theoretical Astrophysics and Cosmology, Institute for Computational Science, University of Zurich,Winterthurerstrasse 190, CH-8057 Zürich, Switzerland Dipartimento di Fisica “G. Occhialini”, Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy Institute of Applied Mathematics, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, 100190 Beijing, China Universitat Politècnica de València, IGIC, Spain Deutsches Elektronen Synchrotron DESY, Platanenallee 6, D-15738 Zeuthen, Germany Kavli Institute for Astronomy and Astrophysics at Peking University, 100871 Beijing, China Zentrum für Astronomie und Astrophysik, TU Berlin, Hardenbergstraße 36, D-10623 Berlin, Germany
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Peters’ formula is an analytical estimate of the time-scale of gravitational wave (GW)-induced coalescence of binary systems. It is used in countless applications, where theconvenience of a simple formula outweighs the need for precision. However, manypromising sources of the Laser Interferometer Space Antenna (LISA), such as super-massive black hole binaries and extreme mass-ratio inspirals (EMRIs), are expectedto enter the LISA band with highly eccentric (e ∼ > Key words: black hole physics – gravitational waves – methods: analytical.
With the launch of space-borne gravitational-wave (GW)detectors in sight (e.g. the Laser Interferometer Space An-tenna, LISA; Amaro-Seoane et al. 2013; Barack et al. 2019),many questions still have to be answered with regards tothe possible sources of signal. In particular, large uncertain-ties remain on the expected event rates for different types ofsources and for different formation channels (e.g. Klein et al.2016; Babak et al. 2017). An essential part of these estimatesis a model of gravitational radiation and its effects on the (cid:63)
E-mail: [email protected] orbital dynamics of a system. Very often, it is convenient toexpress this effect with a simple analytical time-scale.Already in the early 1960s, Peters & Mathews (1963)calculated the radiation reaction equations and extracted asimple analytical formula for the time-scale of GW-induceddecay of a binary. This formula, commonly known as Pe-ters’ time-scale, has since seen an incredible amount of use,and more recently has been employed in many event-rateestimates and other source-modelling applications (see, e.g.Farris et al. 2015; VanLandingham et al. 2016; Bortolas& Mapelli 2019, amongst many others). Peters’ formula isthought to suffice as a first approximation, especially be-cause it is often used in conjunction with other dynamicaltime-scales that are themselves only order-of-magnitude es- © 2020 The Authors a r X i v : . [ a s t r o - ph . GA ] J a n L. Zwick et al. timates (see, e.g. Sesana & Khan 2015). However, many ofLISA most promising sources, such as supermassive blackhole (SMBH) binaries and extreme/intermediate/extremelylarge mass-ratio inspirals (EMRIs/IMRIs/XMRIs), are ex-pected to exhibit highly eccentric and highly relativistic or-bits (see, e.g. Amaro-Seoane et al. 2007; Antonini & Rasio2016; Bonetti et al. 2016, 2018; Khan et al. 2018; Amaro-Seoane 2018; Bonetti et al. 2019; Giacobbo & Mapelli 2019,2020; Amaro-Seoane 2019), which are the two regimes inwhich Peters’ estimate performs the worst. The commonlyused formula does not account for the Newtonian evolutionof the eccentricity. Moreover, it is based on the quadrupoleradiation of a Keplerian binary and therefore fails to cap-ture the more complex post-Newtonian (PN) behaviour ofhighly relativistic orbits.In a recent paper (Zwick et al. 2020, hereafter Paper I),we sought to fix these issues, in service of reducing the un-necessary errors that might arise from using Peters’ formulain its common form. We produced an updated formula thattakes into account the eccentricity evolution and the lowest-order PN correction, while still being algebraically very sim-ple. In this paper, we wish to extend our analysis to the nextorder of the PN series, for the two following reasons. Firstly,it will serve as a proof of good convergence for the previousresults. Secondly, two qualitatively new effects emerge atthe next order: hereditary (i.e. tail) and spin-orbit contribu-tions. Hereditary terms in the PN flux equations are knownto be disproportionately important in the amount of cyclesthat some GW sources complete in the LISA band (see, e.g.Blanchet 2014). Moreover, a significant fraction of SMBHsare expected to have large ( ∼ > . ) dimensionless spin val-ues (see, e.g. Reynolds 2020), whose effects are known tosignificantly alter the event rates of some GW sources, es-pecially EMRIs (Amaro-Seoane et al. 2013). This, combinedwith the weak convergence of the PN series in the case of ex-treme mass ratios (see, e.g. Yunes & Berti 2008; Zhang et al.2011), suggests that there is still precision to be gained be-yond the lowest-order PN correction.At some point, however, one inevitably clashes againstthe fact that most GW sources do not evolve in a pure vac-uum. Even at very small separations, the inspiral processcan be perturbed by the presence of other forces and thevacuum time-scale can be a precise estimate only if theseenvironmental perturbations are weak. Therefore, we takea look at two types of environmental perturbations to theinspiral time-scale. Firstly, we investigate the gravitationalinfluence of a third body. Secondly, we compare the strengthof gas friction and gas torques for an inspiral embedded ina gaseous disc. In both cases, we express the results as aseries of characteristic orbital separations at which the en-vironmental effects influence the inspiral time-scale as muchas a particular PN order. Even though these results are onlyorder-of-magnitude estimates, they clearly illustrate the re-gions of validity and usefulness for the corrections we pro-pose, and more generally for gravitational radiation time-scales.With the goals of this paper set, we wish to refer anyreader not at least partly familiar with Peters’ formula tothe original work from 1963 and to Peters’ Ph.D. thesis (Pe-ters & Mathews 1963; Peters 1964), as well as to a generalreview of PN theory in Schäfer & Jaranowski (2018). We also refer to Paper I for a more in-depth discussion of thelimitations of Peters’ derivation and on the considerationsrequired to go beyond Keplerian orbits and quadrupolar ra-diation. Section 2 presents a brief but complete overviewof the concepts that were used in this paper. In Section 3,we present our derivation of the 1.5 PN correction factorsto Peters’ formula, and discuss how their inclusion can af-fect current event-rate estimates for EMRIs. In Section 4,we discuss environmental effects and derive several usefulcharacteristic radii. Finally, in Section 5, we summarize ourresults. When going to higher PN orders, calculations become noto-riously algebraically intensive. It is therefore convenient towork in natural units, wherein M = m + m = 1 , G = 1 ,and c = 1 , where m and m are the masses of the two bod-ies, G is the gravitational constant, and c is the speed of lightin vacuum. In these units, we can express all expansions ofthe PN series as Taylor expansions in the so-called PN pa-rameter x = − E , where E is the Newtonian orbital energyper unit reduced mass µ = m m /M . The PN parameter x is also a gauge-invariant quantity (see, e.g. Blanchet 2014),meaning that the following calculations are valid regardlessof gauge choice. With the units taken care of, two differentingredients are needed to fully describe the orbital evolutionof a two-body system in PN theory.Firstly, we require a parametrization of the orbital sep-aration r ( φ ) , which solves the PN equations of motion for aconstant orbital energy E and angular momentum vector J : r ( φ ) = r N ( φ ) + r ( φ ) + r . ( φ ) + r ( φ ) + ..., (1)where φ is the polar angle and the first term of the right-hand side is the solution to the familiar Kepler problemin Newtonian gravity. There are many equivalent possibili-ties to express such parametrizations, and in this paper weuse the so-called quasi-Keplerian (QK) one (see Damour &Deruelle 1985; Wex 1995; Memmesheimer et al. 2004), inwhich the radial separation is expressed in terms of twogeneralised orbital parameters, a semi-major axis a and aneccentricity e . These parameters depend on the quantitiesof the binary system, such as the energy E , the angular mo-mentum vector J , and, in the case of BHs, two spin vectors S and S . They serve to greatly simplify the notation, justas introducing a semi-major axis and an eccentricity simpli-fies the solution of the Newtonian Kepler problem.Secondly, we require a set of differential equations thatdescribe how the aforementioned orbital parameters evolvedue to the flux of energy and angular momentum caused by In this work, the magnitude of the spin vector of a BH of angu-lar momentum’s magnitude L and mass M is the dimensionlessquantity defined by ≤ S ≡ cL/ ( GM ) ≤ . Note that, if thesecondary is not a BH, the definition of S needs to be revisited.For the scope of this paper, this clarification is not necessary,since we will neglect the spin of the secondary altogether.MNRAS000
E-mail: [email protected] orbital dynamics of a system. Very often, it is convenient toexpress this effect with a simple analytical time-scale.Already in the early 1960s, Peters & Mathews (1963)calculated the radiation reaction equations and extracted asimple analytical formula for the time-scale of GW-induceddecay of a binary. This formula, commonly known as Pe-ters’ time-scale, has since seen an incredible amount of use,and more recently has been employed in many event-rateestimates and other source-modelling applications (see, e.g.Farris et al. 2015; VanLandingham et al. 2016; Bortolas& Mapelli 2019, amongst many others). Peters’ formula isthought to suffice as a first approximation, especially be-cause it is often used in conjunction with other dynamicaltime-scales that are themselves only order-of-magnitude es- © 2020 The Authors a r X i v : . [ a s t r o - ph . GA ] J a n L. Zwick et al. timates (see, e.g. Sesana & Khan 2015). However, many ofLISA most promising sources, such as supermassive blackhole (SMBH) binaries and extreme/intermediate/extremelylarge mass-ratio inspirals (EMRIs/IMRIs/XMRIs), are ex-pected to exhibit highly eccentric and highly relativistic or-bits (see, e.g. Amaro-Seoane et al. 2007; Antonini & Rasio2016; Bonetti et al. 2016, 2018; Khan et al. 2018; Amaro-Seoane 2018; Bonetti et al. 2019; Giacobbo & Mapelli 2019,2020; Amaro-Seoane 2019), which are the two regimes inwhich Peters’ estimate performs the worst. The commonlyused formula does not account for the Newtonian evolutionof the eccentricity. Moreover, it is based on the quadrupoleradiation of a Keplerian binary and therefore fails to cap-ture the more complex post-Newtonian (PN) behaviour ofhighly relativistic orbits.In a recent paper (Zwick et al. 2020, hereafter Paper I),we sought to fix these issues, in service of reducing the un-necessary errors that might arise from using Peters’ formulain its common form. We produced an updated formula thattakes into account the eccentricity evolution and the lowest-order PN correction, while still being algebraically very sim-ple. In this paper, we wish to extend our analysis to the nextorder of the PN series, for the two following reasons. Firstly,it will serve as a proof of good convergence for the previousresults. Secondly, two qualitatively new effects emerge atthe next order: hereditary (i.e. tail) and spin-orbit contribu-tions. Hereditary terms in the PN flux equations are knownto be disproportionately important in the amount of cyclesthat some GW sources complete in the LISA band (see, e.g.Blanchet 2014). Moreover, a significant fraction of SMBHsare expected to have large ( ∼ > . ) dimensionless spin val-ues (see, e.g. Reynolds 2020), whose effects are known tosignificantly alter the event rates of some GW sources, es-pecially EMRIs (Amaro-Seoane et al. 2013). This, combinedwith the weak convergence of the PN series in the case of ex-treme mass ratios (see, e.g. Yunes & Berti 2008; Zhang et al.2011), suggests that there is still precision to be gained be-yond the lowest-order PN correction.At some point, however, one inevitably clashes againstthe fact that most GW sources do not evolve in a pure vac-uum. Even at very small separations, the inspiral processcan be perturbed by the presence of other forces and thevacuum time-scale can be a precise estimate only if theseenvironmental perturbations are weak. Therefore, we takea look at two types of environmental perturbations to theinspiral time-scale. Firstly, we investigate the gravitationalinfluence of a third body. Secondly, we compare the strengthof gas friction and gas torques for an inspiral embedded ina gaseous disc. In both cases, we express the results as aseries of characteristic orbital separations at which the en-vironmental effects influence the inspiral time-scale as muchas a particular PN order. Even though these results are onlyorder-of-magnitude estimates, they clearly illustrate the re-gions of validity and usefulness for the corrections we pro-pose, and more generally for gravitational radiation time-scales.With the goals of this paper set, we wish to refer anyreader not at least partly familiar with Peters’ formula tothe original work from 1963 and to Peters’ Ph.D. thesis (Pe-ters & Mathews 1963; Peters 1964), as well as to a generalreview of PN theory in Schäfer & Jaranowski (2018). We also refer to Paper I for a more in-depth discussion of thelimitations of Peters’ derivation and on the considerationsrequired to go beyond Keplerian orbits and quadrupolar ra-diation. Section 2 presents a brief but complete overviewof the concepts that were used in this paper. In Section 3,we present our derivation of the 1.5 PN correction factorsto Peters’ formula, and discuss how their inclusion can af-fect current event-rate estimates for EMRIs. In Section 4,we discuss environmental effects and derive several usefulcharacteristic radii. Finally, in Section 5, we summarize ourresults. When going to higher PN orders, calculations become noto-riously algebraically intensive. It is therefore convenient towork in natural units, wherein M = m + m = 1 , G = 1 ,and c = 1 , where m and m are the masses of the two bod-ies, G is the gravitational constant, and c is the speed of lightin vacuum. In these units, we can express all expansions ofthe PN series as Taylor expansions in the so-called PN pa-rameter x = − E , where E is the Newtonian orbital energyper unit reduced mass µ = m m /M . The PN parameter x is also a gauge-invariant quantity (see, e.g. Blanchet 2014),meaning that the following calculations are valid regardlessof gauge choice. With the units taken care of, two differentingredients are needed to fully describe the orbital evolutionof a two-body system in PN theory.Firstly, we require a parametrization of the orbital sep-aration r ( φ ) , which solves the PN equations of motion for aconstant orbital energy E and angular momentum vector J : r ( φ ) = r N ( φ ) + r ( φ ) + r . ( φ ) + r ( φ ) + ..., (1)where φ is the polar angle and the first term of the right-hand side is the solution to the familiar Kepler problemin Newtonian gravity. There are many equivalent possibili-ties to express such parametrizations, and in this paper weuse the so-called quasi-Keplerian (QK) one (see Damour &Deruelle 1985; Wex 1995; Memmesheimer et al. 2004), inwhich the radial separation is expressed in terms of twogeneralised orbital parameters, a semi-major axis a and aneccentricity e . These parameters depend on the quantitiesof the binary system, such as the energy E , the angular mo-mentum vector J , and, in the case of BHs, two spin vectors S and S . They serve to greatly simplify the notation, justas introducing a semi-major axis and an eccentricity simpli-fies the solution of the Newtonian Kepler problem.Secondly, we require a set of differential equations thatdescribe how the aforementioned orbital parameters evolvedue to the flux of energy and angular momentum caused by In this work, the magnitude of the spin vector of a BH of angu-lar momentum’s magnitude L and mass M is the dimensionlessquantity defined by ≤ S ≡ cL/ ( GM ) ≤ . Note that, if thesecondary is not a BH, the definition of S needs to be revisited.For the scope of this paper, this clarification is not necessary,since we will neglect the spin of the secondary altogether.MNRAS000 , 1–14 (2020) mproved gravitational radiation time-scales II the emission of GWs. A convenient form is found when theorbital parameters a and e are expressed in terms of the PNparameter x and a dummy eccentricity parameter e d , thatessentially is a measure of the angular momentum of thesystem: a ( x, e d ) = 1 x ∞ (cid:88) n =0 C an/ ( e d ) x n/ , (2) e ( x, e d ) = ∞ (cid:88) n =0 C en/ ( e d ) x n/ , (3)where the subscripts of the coefficients denote the relativePN order and some of the coefficients can vanish. The ex-plicit coefficients used in this paper are listed in the Ap-pendix, and can be found in Memmesheimer et al. (2004)for the spin-less part and in Klein & Jetzer (2010) for anextension that covers the spin-orbit coupling. The evolutionequations for x and e d take a relatively simple form: ˙ x ( x, e d ) = x ∞ (cid:88) n =0 A n/ ( e d ) x n/ , (4) ˙ e d ( x, e d ) = x ∞ (cid:88) n =0 B n/ ( e d ) x n/ , (5)where, similarly to before, the subscripts of the coefficientsdenote the relative PN order and some of the coefficientscan vanish. The explicit coefficients used in this paper canbe found in Boetzel et al. (2017) for the spin-less parts and inKlein & Jetzer (2010) for parts describing the spin-orbit cou-pling. Note that these evolution equations are always orbit-averaged, the assumption being that changes to the orbitalelements occur on time-scales that are much longer than asingle orbital period, which is the same assumption that isrequired to produce Peters’ formula in the first place. At the lowest order in the PN parameter x , Equations (2)and (3) reduce to a = x − and e = e d , respectively, whichdescribe a binary moving along a Keplerian orbit. Equa-tions (4) and (5) reduce instead to ˙ x = A x f ( e ) and ˙ e = B x g ( e ) , respectively, where A and B can be foundin the appendix, and f ( e ) and g ( e ) are known as the eccen-tricity enhancement functions: f ( e ) = (cid:18) e + 3796 e (cid:19) (1 − e ) − / , (6) g ( e ) = (cid:18) e (cid:19) (1 − e ) − / . (7)They describe a binary that is radiating GWs in accor-dance with Einstein (1916)’s quadrupole formula. Because ofthe coupled nature of these equations and the complexity ofthe eccentricity enhancement functions, it is impossible tosolve them analytically for an arbitrary initial semi-majoraxis a and eccentricity e . Nevertheless, Peters and Math-ews gave some approximated solutions for a ( t ) and e ( t ) (see Peters & Mathews 1963; Peters 1964). In particular, the so-lutions for a ( t ) can be set equal to zero and solved for thetime t in order to yield a decay time-scale of the orbitalevolution.The most commonly used formula for Peters’ time-scale, t P , is derived by solving Equation (2) at lowest order whilealso neglecting the evolution of the eccentricity, thus setting ˙ e = 0 . The solution reads t P = 5(1 + q ) q a f ( e ) = 5 c (1 + q ) G M q a f ( e ) , (8)where q = m /m ≤ is the mass ratio, and we reintroducedthe physical constants. As stated in the introduction, thisformula is currently used in a wide variety of applications,even though it is known to fail for highly eccentric and forhighly relativistic orbits.In Paper I, we found two simple factors, R and Q f , thatcan be multiplied by Peters’ time-scale to resolve its issues: t P ( a , e ) → t P ( a , e ) R ( e ) Q f ( a , e ) . (9)The fitting function R ( e ) = 8 −√ − e is purely Newto-nian, and is needed to resolve the lack of eccentricity evolu-tion in the standard formula. Other fits with the same pur-pose have been used in some select cases (see, e.g. Bonettiet al. 2018), but tend to be much more complex than thefactor R and are therefore uncommon.The factor Q f ( a , e ) is a simple fit that roughly ac-counts for both the more complex 1 PN parametrizationgiven by Equations (2) and (3), and the relative 1 PN fluxesin Equations (4) and (5). The corrected time-scale reads t P RQ f = 5 c (1 + q ) G M q a f ( e ) (cid:124) (cid:123)(cid:122) (cid:125) Peters’ formula −√ − e exp (cid:18) . r S a (1 − e ) (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) eccentricity and PN correction , (10)where r S is the effective Schwarzschild (1916) radius of thetwo-body system, r S = 2 GMc − , and the coefficient 2.5 inthe exponential function was chosen to fit a large rangeof initial eccentricities. In Paper I, we showed that Equa-tion (10) can resolve errors of order 1–10 than can arise whenapplying Peter’s formula to highly eccentric and relativisticorbits. BH spin couplings first appear at 1.5 PN order in the PNexpansion. From then on, both spins have to be taken intoaccount, drastically increasing the number of degrees of free-dom and dynamical variables of the system. The spins caninteract with the orbital angular momentum at 1.5 PN orderand also with each other at 2 PN order. As is shown in Ta-ble 1, this increases the total number of dynamical variablesfrom 2 ( a and e ) to 10. Sampling an appropriate range ofparameter space with a grand total of ten initial conditionswould simply be an exercise in frustration.At the cost of some loss of generality, however, we can MNRAS , 1–14 (2020)
L. Zwick et al.
Order Parametrization Fluxes Dyn. Variables0 PN Keplerian Quadrupole a , e a , e a , e , ˆJ , S , S a , e , ˆJ , S , S Table 1.
This table summarises the new elements that enter thepicture at a given PN order. The acronyms QK, S-O, and S-Sstand for “quasi-Keplerian”, “spin-orbit”, and “spin-spin”, respec-tively. After the first PN order, one must in principle considerthe evolution of both BH spins, S and S , as well as the specificangular momentum vector ˆJ . reduce the amount of new dynamical variables to only one.LISA BH sources will consist of three main categories: EM-RIs, IMRIs, and SMBH binaries. These are mainly distin-guished by ranges of their mass ratios: ∼ − , ∼ − to ∼ − and ∼ − to , respectively. For EMRIs and IM-RIs, the mass ratio is by definition very small, but even forSMBH mergers only a small minority will have mass ratiosof order one (see, e.g. O’Leary et al. 2021, where we assumethat galaxy mass is a tracer for SMBH mass). In light ofapplications for LISA sources, it is very natural to restrictour analysis to small mass ratios, i.e. q ∼ < . . This turnsout to be extremely convenient, because a small mass ratiosuppresses most of the complex effects of spin-orbit and es-pecially spin-spin couplings. We decided to also neglect thestandard, non-Kerr 2 PN contributions to the parametriza-tion and the fluxes of the two-body system. This last as-sumption is justified for two connected reasons: firstly, ournumerical investigations showed that integer PN orders tendto converge much faster than half-integer tail contributions(see also Blanchet 2014), making the 2 PN correction verysmall; secondly, our goal is to produce a correction to Pe-ters’ formula that is sufficiently simple to be used in thebroader context of source modelling for GW observatories.Fitting the minute 2 PN deviations would necessarily requirea much more complex fitting function. A divergence-free PN parametrization of the orbital motionthat includes spin-orbit and spin-spin couplings can be foundin Klein & Jetzer (2010). Up to 1.5 PN order, it has the formof a standard QK parametrization: r = a (1 − e cos u ) (11)The parameters a and e describe the shape of the orbit,whereas the parameter u is responsible for orbital precession.They are expressed as series in the PN parameter x and theyread a ( x, e d ) = 1 x −
74 + 2 √ x (cid:112) − e (cid:12)(cid:12)(cid:12) ˆJ × S (cid:12)(cid:12)(cid:12) + O [ x ] , (12) e ( x, e d ) = e d + 3 e d x − √ x (cid:112) − e (cid:12)(cid:12)(cid:12) ˆJ × S (cid:12)(cid:12)(cid:12) + O [ x ] , (13)where the 1 PN term is taken from equations (20a) and (20b)in Memmesheimer et al. (2004) and simplified according tothe small mass-ratio limit. We wish to compare a Newtonianorbit and a PN orbit that start out with the same orbitalelements at a given moment in time ( t = 0 ) and are sub-sequently evolved with the appropriate energy and angularmomentum fluxes. To enforce this, we set Equations (12)and (13) equal to their Newtonian counterparts. The onlyway for these two orbits to be identical at t = 0 is for them tohave different values of the PN parameter x and the dummyeccentricity e d . We denote them as ( x N , e N ) for the Newto-nian orbit and as ( x PN , e PN ) for the PN orbit. To find thecorrect values for x PN and e PN , we set the Newtonian andPN orbital elements equal to each other: x N ! = a ( x PN , e PN ) , (14) e N ! = e ( x PN , e PN ) . (15)We can solve these equations order by order and findthe values of the PN parameters so that the two orbits areidentical at the beginning of their evolution: x PN = x N − x + 2 (cid:112) x (cid:112) − e (cid:12)(cid:12)(cid:12) ˆJ × S (cid:12)(cid:12)(cid:12) + O [ x ] , (16) e PN = e N − e N x N + 2 (cid:112) x (cid:112) − e (cid:12)(cid:12)(cid:12) ˆJ × S (cid:12)(cid:12)(cid:12) + O [ x ] . (17)As noted before, the energy and angular momentumfluxes can be averaged into secular evolution equations forthe parameters x and e . The equations are organised intoPN series. Schematically, they read ˙ x = x (cid:16) A N + xA + x / ( A t + A SO ) (cid:17) + O [ x ] , (18) ˙ e d = x (cid:16) B N + xB + x / ( B t + B SO ) (cid:17) + O [ x ] , (19)where the different subscripts denote the Newtonian and1 PN contributions, as well as the lowest-order hereditarytail and the spin-orbit coupling contribution. The full ex-pressions for the coefficients can be found in Klein & Jetzer(2010) and Ebersold et al. (2019) or in the Appendix.In principle, we are missing various sets of evolutionequations for the direction of the orbital angular momentumvector and the two spin vectors. However, at lowest order inthe spin-orbit coupling and the mass ratio, the orbital an-gular momentum vector J only undergoes a Lense–Thirring(Lense & Thirring 1918) precession. Furthermore, all termsthat contain ˆJ are of the form MNRAS000
74 + 2 √ x (cid:112) − e (cid:12)(cid:12)(cid:12) ˆJ × S (cid:12)(cid:12)(cid:12) + O [ x ] , (12) e ( x, e d ) = e d + 3 e d x − √ x (cid:112) − e (cid:12)(cid:12)(cid:12) ˆJ × S (cid:12)(cid:12)(cid:12) + O [ x ] , (13)where the 1 PN term is taken from equations (20a) and (20b)in Memmesheimer et al. (2004) and simplified according tothe small mass-ratio limit. We wish to compare a Newtonianorbit and a PN orbit that start out with the same orbitalelements at a given moment in time ( t = 0 ) and are sub-sequently evolved with the appropriate energy and angularmomentum fluxes. To enforce this, we set Equations (12)and (13) equal to their Newtonian counterparts. The onlyway for these two orbits to be identical at t = 0 is for them tohave different values of the PN parameter x and the dummyeccentricity e d . We denote them as ( x N , e N ) for the Newto-nian orbit and as ( x PN , e PN ) for the PN orbit. To find thecorrect values for x PN and e PN , we set the Newtonian andPN orbital elements equal to each other: x N ! = a ( x PN , e PN ) , (14) e N ! = e ( x PN , e PN ) . (15)We can solve these equations order by order and findthe values of the PN parameters so that the two orbits areidentical at the beginning of their evolution: x PN = x N − x + 2 (cid:112) x (cid:112) − e (cid:12)(cid:12)(cid:12) ˆJ × S (cid:12)(cid:12)(cid:12) + O [ x ] , (16) e PN = e N − e N x N + 2 (cid:112) x (cid:112) − e (cid:12)(cid:12)(cid:12) ˆJ × S (cid:12)(cid:12)(cid:12) + O [ x ] . (17)As noted before, the energy and angular momentumfluxes can be averaged into secular evolution equations forthe parameters x and e . The equations are organised intoPN series. Schematically, they read ˙ x = x (cid:16) A N + xA + x / ( A t + A SO ) (cid:17) + O [ x ] , (18) ˙ e d = x (cid:16) B N + xB + x / ( B t + B SO ) (cid:17) + O [ x ] , (19)where the different subscripts denote the Newtonian and1 PN contributions, as well as the lowest-order hereditarytail and the spin-orbit coupling contribution. The full ex-pressions for the coefficients can be found in Klein & Jetzer(2010) and Ebersold et al. (2019) or in the Appendix.In principle, we are missing various sets of evolutionequations for the direction of the orbital angular momentumvector and the two spin vectors. However, at lowest order inthe spin-orbit coupling and the mass ratio, the orbital an-gular momentum vector J only undergoes a Lense–Thirring(Lense & Thirring 1918) precession. Furthermore, all termsthat contain ˆJ are of the form MNRAS000 , 1–14 (2020) mproved gravitational radiation time-scales II P ]10 a , p [ G M / c ] Newtonian1PN1.5PN, s = 01.5PN, s = -11.5PN, s = 1 0 1 2 3 4 5 6t [t P ]10 a , p [ G M / c ] P ]10 a , p [ G M / c ] Figure 1.
Three examples of the time evolution of the orbital semi-major axis a (full lines) and periapsis p (dotted lines), for differentinitial eccentricities (0.999, 0.9, 0.3). The different colours represent different PN orders and spin values for the primary BH. The spinvalues s = − (red) and s = 1 (purple) are chosen to show the two extreme cases of an equatorial co-rotating and counter-rotating orbitaround an extremal BH. The blue line ( s = 0 ) represents the Schwarzschild case with hereditary PN terms. Note how the higher-orderPN corrections converge around the first order one much more strongly for eccentric orbits. (cid:12)(cid:12)(cid:12) ˆJ × S (cid:12)(cid:12)(cid:12) = S cos θ = s , (20)where θ is the angle between the primary BH spin vectorand the orbital angular momentum vector. The magnitudeof the BH spin is assumed to be constant, because the evo-lution occurs in vacuum. Furthermore, the Lense–Thirringprecession does not affect the angle θ , hence the expres-sion in Equation (20) is constant in time. It follows thenthat Equations (18) and (19) are a closed system of differen-tial equations, which can be integrated numerically, startingfrom the correct initial conditions given by Equations (16)and (17).As we have noted before, there are two separate ef-fects that influence the inspiral time-scale when consideringhigher-order PN perturbations. On the one hand, the evo-lution equations for the orbital parameters gain additionalterms, which originate from the additional energy and an-gular momentum fluxes and, generally, tend to shorten theinspiral time-scale. On the other hand, some adjustments tothe initial conditions are required to ensure that the orbitalelements of a PN orbit and those of a Newtonian orbit areidentical at the beginning of the evolution, which tend to in-crease the overall time-scale. The reason for this increase isthat a binary in a PN gravity field requires slightly more en-ergy and angular momentum to support a given orbital sep-aration. This extra energy must be radiated away and there-fore tends to increase the overall inspiral time-scale. To get afeel for how different orders affect the overall behaviour, weshow the evolution of a highly eccentric ( e = 0 . ), mod-erately eccentric ( e = 0 . ) and low eccentricity ( e = 0 . )orbit in Figure 1. In all cases, Peters’ formula underestimatesthe inspiral time-scale, especially so for eccentric orbits. Theeffect of the 1 PN correction is also to further increase thetime-scale, while higher order corrections cluster around the1 PN result. In order to compute the inspiral time-scale, we performednumerous numerical integrations of the orbital evolution,sampling a wide range of initial orbital separations, eccen-tricities, and spin parameter values. We define the mergerevent as the moment at which the periapsis of the orbit, p = a (1 − e ) , reaches the effective Schwarzschild radius. Weexpress the corrections to Peters’ time-scale as ratios R and Q i , where the index i stands for the several different effectsthat we want to describe. The most complete and accurateformula for the time-scale is obtained by multiplying Petersformula t P by three separate correction ratios, t P ( a , e ) → t P ( a , e ) R ( e ) Q h ( a , e ) Q s ( a , e , s ) , (21)where the ratio R ( e ) accounts for the secular eccentricityevolution at quadrupole order, Q h ( a , e ) accounts for the1 PN and 1.5 PN spinless effects, and Q s ( a , e , s ) accountsfor the inclusion of the spin-orbit coupling. We calculatethese ratios numerically by integrating the appropriate evo-lution equations, and subsequently find various fits for theresults. The numerical results can be seen in the left panelof Figure 2, where we plot the ratio of the 1.5 PN time-scale to the Newtonian time-scale for different values of theperi-apsis, spin and initial eccentricity. As shown in PaperI, the correction generally grows for smaller periapsis, whilethe two edge cases of co- and counter-rotating orbits aroundan extremal BH envelop the spin-less results. In the circularcase, the spin parameter s is degenerate with the periapsisof the orbit by ∼ > r S , while it has much less of an influencein the eccentric case. Our proposal for the appropriate fitsto these ratios reads MNRAS , 1–14 (2020)
L. Zwick et al. e
10 %25 %50 %75 %100 % t s / t N p = 6 r S p = 9 r S p = 12 r S p = 18 r S e F i t R e l a t i v e E rr o r p = 6 r S p = 9 r S p = 12 r S p = 18 r S Figure 2.
In the left-hand panel, we plot the relative increase in the time-scale of GW-induced decay due to the 1 PN and 1.5 PNcontributions with respect to the Newtonian prediction. The solid lines represent the spinless case, whereas the dotted and dash-dottedlines represent the minimum ( − ) and maximum (1) values for the spin parameter s , respectively. These curves must be multiplied bythe correction factor R ( e ) to compute the complete correction to Peters’ formula. In the right-hand panel, we show the relative errorof the fits we propose for the same range of parameter space. For the sake of clarity, we avoid plotting the curves with spin, but notethat the fits perform similarly to the spinless case (with the single exception of a circular, equatorial and prograde orbit at r S arounda maximally spinning SMBH, where the fit fails by no more than 30 per cent). For a periapsis p ∼ > r S , the errors lie well below 1 percent. Even for the extreme case of p = 6 r S the fits perform rather well, considering the fact that the PN series itself is known to breakdown below those separations. R ( e ) = 8 −√ − e , (22) Q h ( p , e ) = exp (cid:18) . r S p (cid:19) f ( p , e ) , (23) f ( p , e ) = (cid:0) A ( p ) (1 − e ) + B ( p ) (1 − e ) (cid:1) ,A ( p ) = (cid:32) − (cid:18) . r S p (cid:19) + (cid:18) . r S p (cid:19) / (cid:33) ,B ( p ) = − (cid:18) . r S p (cid:19) / ,Q s ( p , e , s ) = exp (cid:16) A s ( p , e ) s + B s ( p , e ) | s | / (cid:17) , (24) A s ( p , e ) = e . r S p + (1 − e ) / (cid:18) . r S p (cid:19) / ,B s ( p , e ) = (cid:18) e . r S p (cid:19) / + (1 − e ) (cid:18) . r S p (cid:19) , valid for arbitrary values of the initial semi-major axis, ec-centricity, and (primary BH) spin parameter.These general correction factors turn out to be rathercomplicated and unwieldy. The reason for this is the compli-cated behaviour of the time-scale in a region of parameterspace between . ∼ < e ∼ < . , visible in Figure 2. The likelycause for this are some terms in the higher-PN fluxes switch-ing sign at eccentricities of ∼ e (cid:29) . by scattering compact ob-jects onto almost radial orbits (see e.g. Amaro-Seoane et al.2007, 2013); on the other hand, some EMRIs can be pro-duced when a binary containing a compact object is dis-rupted by the tidal influence of the SMBH (see e.g. Amaro-Seoane 2019). These would likely enter the GW-dominatedphase with very low eccentricities. For this reason, the mostgeneral form for the fits of the correction ratios is often notrequired. Rather, one can freely choose an appropriate ec-centricity range or limit that is expected for a particularformation channel. Here we will discuss the two most obvi-ous limits of e → and e → . In these limits, the fitsreduce to much more manageable formulas that satisfy ourgoal to be as simple as possible. If the binary enters the GW-dominated phase with a verylow eccentricity ( e ∼ < . ), the fits for the correction ratiossimplify to the following formulas: R → , (25) Q h Q s → exp (cid:32) . r S p + s (cid:18) . r S p (cid:19) / + | s | / (cid:18) . r S p (cid:19) (cid:33) . (26) MNRAS000
In the left-hand panel, we plot the relative increase in the time-scale of GW-induced decay due to the 1 PN and 1.5 PNcontributions with respect to the Newtonian prediction. The solid lines represent the spinless case, whereas the dotted and dash-dottedlines represent the minimum ( − ) and maximum (1) values for the spin parameter s , respectively. These curves must be multiplied bythe correction factor R ( e ) to compute the complete correction to Peters’ formula. In the right-hand panel, we show the relative errorof the fits we propose for the same range of parameter space. For the sake of clarity, we avoid plotting the curves with spin, but notethat the fits perform similarly to the spinless case (with the single exception of a circular, equatorial and prograde orbit at r S arounda maximally spinning SMBH, where the fit fails by no more than 30 per cent). For a periapsis p ∼ > r S , the errors lie well below 1 percent. Even for the extreme case of p = 6 r S the fits perform rather well, considering the fact that the PN series itself is known to breakdown below those separations. R ( e ) = 8 −√ − e , (22) Q h ( p , e ) = exp (cid:18) . r S p (cid:19) f ( p , e ) , (23) f ( p , e ) = (cid:0) A ( p ) (1 − e ) + B ( p ) (1 − e ) (cid:1) ,A ( p ) = (cid:32) − (cid:18) . r S p (cid:19) + (cid:18) . r S p (cid:19) / (cid:33) ,B ( p ) = − (cid:18) . r S p (cid:19) / ,Q s ( p , e , s ) = exp (cid:16) A s ( p , e ) s + B s ( p , e ) | s | / (cid:17) , (24) A s ( p , e ) = e . r S p + (1 − e ) / (cid:18) . r S p (cid:19) / ,B s ( p , e ) = (cid:18) e . r S p (cid:19) / + (1 − e ) (cid:18) . r S p (cid:19) , valid for arbitrary values of the initial semi-major axis, ec-centricity, and (primary BH) spin parameter.These general correction factors turn out to be rathercomplicated and unwieldy. The reason for this is the compli-cated behaviour of the time-scale in a region of parameterspace between . ∼ < e ∼ < . , visible in Figure 2. The likelycause for this are some terms in the higher-PN fluxes switch-ing sign at eccentricities of ∼ e (cid:29) . by scattering compact ob-jects onto almost radial orbits (see e.g. Amaro-Seoane et al.2007, 2013); on the other hand, some EMRIs can be pro-duced when a binary containing a compact object is dis-rupted by the tidal influence of the SMBH (see e.g. Amaro-Seoane 2019). These would likely enter the GW-dominatedphase with very low eccentricities. For this reason, the mostgeneral form for the fits of the correction ratios is often notrequired. Rather, one can freely choose an appropriate ec-centricity range or limit that is expected for a particularformation channel. Here we will discuss the two most obvi-ous limits of e → and e → . In these limits, the fitsreduce to much more manageable formulas that satisfy ourgoal to be as simple as possible. If the binary enters the GW-dominated phase with a verylow eccentricity ( e ∼ < . ), the fits for the correction ratiossimplify to the following formulas: R → , (25) Q h Q s → exp (cid:32) . r S p + s (cid:18) . r S p (cid:19) / + | s | / (cid:18) . r S p (cid:19) (cid:33) . (26) MNRAS000 , 1–14 (2020) mproved gravitational radiation time-scales II In this case, the eccentricity correction factor does nothave a significant effect, but the PN correction factors taketheir maximal values for the coefficients. This occurs be-cause binaries on circular orbits spend the whole durationof the orbit at small separations, maximising the influence ofhigher-order PN effects. In this limit, the fits are accurate to ∼
10 per cent of the relative value if p ∼ > r S (see Figure 2). If the binary enters the GW dominated phase with a veryhigh eccentricity ( e ∼ > . ), the fits for the correction ratiossimplify to the following formulas: R = 8 −√ − e (27) Q h Q s → exp (cid:32) . r S p + s (cid:18) . r S p (cid:19) + | s | / (cid:18) . r S p (cid:19) / (cid:33) . (28)In this case, the eccentricity correction factor domi-nates, whereas the PN correction factors are slightly sup-pressed with respect to the circular case. Nonetheless, bothfactors can compound on each other to reach values of ∼ for large regions of initial parameter space. In this limit,the fits are accurate to ∼
10 per cent of the relative value if p ∼ > r S , and to ∼ p ∼ > r S (see Figure 2). The EMRI critical semi-major axis, a crit , marks the separa-tion between the dynamical-evolution regime and the GW-emission regime, and is an essential ingredient for the cal-culation of EMRI event rates. The value of a crit is found byequating the inspiral time-scale to the dynamical relaxationtime-scale for a periapsis equal to the primary’s innermoststable circular orbit (ISCO). In Figure 3, we plot a crit as de-rived in Amaro-Seoane (2019), this time including the cor-rection factors to the GW-inspiral time-scale, for a systemof a 10 M (cid:12) BH inspiralling into a × M (cid:12) SMBH.For eccentric orbits, we find that the value of the criticalsemi-major axis decreases by an order of magnitude fromprevious results that use Peter’s formula as a simple modelof the inspiral time-scale. In the circular orbit case, the valueof a crit is generally decreased by a factor of a few, but theeffects of the spin correction factor Q s are clearly visible inthe steeper scaling of the counter-rotating case ( S < ), andin the slightly flattened scaling for the co-rotating case ( S > ). The difference between the circular and the eccentriccases is also increased by almost an order of magnitude,which suggests different relative efficiencies for formationchannels that preferentially create circular versus eccentricEMRIs. Note that the ISCO of a BH in a prograde orbit canbe as low as ∼ > r S , which is lower than the value where ourfits (and PN theory in general) are accurate. In these cases,we evaluate the correction factors at r S for highly eccentricorbits and at r S for circular orbits, to ensure that the resultsdo not diverge. This would represent a conservative estimateof the real change of the critical semi-major axis.The expected event rates for EMRIs are then computed a c r i t [ p c ] a crit ; e = 0not corrected a crit ; e = 0.999not corrected Figure 3.
We plot the corrected critical semi-major axis a crit (solid lines) as a function of the SMBH spin for an EMRI com-posed of a × M (cid:12) SMBH and a BH of 10 M (cid:12) , for two initialeccentricities. The results are shown for orbits with a nominal in-clination of θ = 0 . with respect to the SMBH equatorial plane,chosen to represent a typical case. The effect of the correctionfactors to Peters’ time-scale is to shift the curves downwards, andalso to change the scaling with the SMBH spin. by performing an integral over phase space, in which a crit isone upper limit (Amaro-Seoane et al. 2013; Amaro-Seoane2019). The total phase-space volume of integration scales as a , which means that even moderate changes to a crit canpotentially have a large effect on the expected event rates.For this reason, we are currently working on a precise re-derivation of EMRI and XMRI (Amaro-Seoane 2019) eventrates that include the corrections to the inspiral time-scale.Preliminary results show that the rates can decrease by afactor of ∼
20 for eccentric EMRIs and a factor of ∼ In this section, we consider the case in which a distant mas-sive body perturbs the orbit of an inspiralling LISA source.In other words, we consider a hierarchical three-body prob-lem in which the hard binary is made of compact objects atvery small orbital separations, and can exhibit relativisticbehaviour. Following Will (2014), we express the energy ofsuch a system in the following way: E = E N + E P + E + O (cid:104) m c (cid:105) , (29) MNRAS , 1–14 (2020)
L. Zwick et al. where E N and E are the Newtonian and 1 PN energycontributions, whereas E P is the change in energy of thehard binary of a given osculating semi-major axis and ec-centricity due to the influence of the distant third body.The explicit expression for E P , in the case of a perturber ina circular orbit, reads E P = − Gµm p R (cid:0) − ( I ) sin ( ω ) (cid:1) = − Gµm p R G ( I , ω ) (30)where the function G describes the angle dependence of theenergy shift, µ = m m /M , m is the mass of the perturber, R its distance to the centre of mass of the hard binary, I the inclination with respect to the orbital plane of the per-turber, and ω the angle of the periapsis with respect to theascending or descending nodes of the orbital plane of theperturber. The remarkable fact proven in Will (2014) is thatthis expression remains constant over time-scales at whichperihelion precession can significantly change the orientationof the orbit with respect to the plane of the perturber. Theterms that would normally vary in ω are cancelled by crossterms between the PN dynamics of the inner binary and thetidal forces induced by the perturber. Indeed, the second sin term is evaluated at an initial perihelion orientation,and therefore is not affected by perihelion precession. Thisresult is very convenient for our purposes, because it is es-sentially of the same form as the change in initial conditionscaused by a PN perturbation.We can estimate how the change in energy of the hardbinary caused by the perturber will affect the inspiral time-scale by readjusting the initial conditions of the GW-induceddecay. We define a fictitious semi-major axis a P and eccen-tricity e P , constructed from both the Newtonian energy ofthe hard binary, E N , and the constant shift induced by theperturber, E P : a P = − GM E N + E P ) ≈ a − m µa MR (1 − e ) G ( I , ω ) , (31) e ≈ e − E P E N (1 − e ) , (32)An estimate for the change in the inspiral time-scalecan by found by evaluating a particular vacuum time-scale t GW at the fictitious semi-major axis and dividing it by thecorresponding vacuum time-scale of the physical orbit: t GW ( a P , e P ) t GW ( a, e ) ≈ (cid:18) E N + E P E N (cid:19) − ≈ − E P E N F ( e ) ≈ − m p G ( I , ω ) MR (1 − e ) F ( e ) (33)where, in the first step, we assumed quadrupole radiationor, equivalently, the evolution equations of Peters & Math-ews (1963). The function F can be computed from the ec-centricity dependence of the inspiral time-scale, and ranges from 1/5 to 1/2 for most values of the eccentricity. Fromthe results of Paper I and previous sections, we know thatPN corrections to the Newtonian vacuum time-scale canbe described by factors of the form exp (cid:0) C PN r S p − (cid:1) ≈ C PN r S p − + O [( r S p − ) n ] , where C PN is a coefficient oforder unity and n is the PN order. We can therefore estimatethe radius at which the environmental perturbation becomesof Newtonian (or PN) order by solving the equation m p G ( I , ω ) MR (1 − e ) F ( e ) = C PN (cid:18) r S p (cid:19) n (34)and find the characteristic radii R P / PN at which the effectof the perturbation becomes as large as that of a given PNorder: R P / PN = (cid:18) G ( I , ω ) m p n (1 − e ) Mr n S (cid:19) / , (35)where we approximated F ( e ) / ≈ and C / ≈ . For theNewtonian order ( n = 0 and C PN = 1 ), the radius reads R P/ N ≈ G ( I , ω ) / (cid:16) m M (cid:17) / a (1 − e ) / , (36)which can only be of order of or smaller than the semi-majoraxis of the hard binary itself, breaking the assumption ofa hierarchical triplet. In general, every additional PN orderonly adds a factor of a few over the previous value and, sincewe are considering hierarchical triplets (where R (cid:29) a ), wecan conclude that the perturbation has very little influenceon the inspiral time-scale.On the other hand, if the accumulated signal-to-noiseratio (SNR) of the source is high enough, LISA might beable to detect deviations from very high PN orders. In Fig-ure 4, we show the minimum order of a PN-waveform tem-plate that is required to detect a shift in the time-scale fora perturber of a given mass and separation, in the case ofan inner binary composed of an SMBH of M (cid:12) and acompact object of 10 M (cid:12) . We compute it by solving Equa-tion (35) for the PN order n and adding one. If we assumethat the SNR is sufficient to reach current state-of-the-artwaveform templates with spin (depicted in Figure 4; Fujita2015), we find values that are similar to the results in Yuneset al. (2011), wherein it is shown that detecting the dephas-ing in the signal of such a system is in principle possible,but only realistic for very massive ( ∼ M (cid:12) ) and close(sub-parsec) perturbers.The shift in energy is not the only way in which a thirdbody can perturb the hard binary. For orbits that are suf-ficiently misaligned, Kozai–Lidov (KL; Kozai 1962; Lidov1962) oscillations can induce a secular variation in the ec-centricity of the orbit. If the inspiral time-scale is longerthan the time-scale of KL oscillations, t KL , the evolutionequations for the orbital parameters of the hard binary canbe significantly affected by the perturber. The detectabilityof KL oscillations in the GW signal of LISA sources has beenvery recently studied in Randall & Xianyu (2019) and Deme MNRAS000
L. Zwick et al. where E N and E are the Newtonian and 1 PN energycontributions, whereas E P is the change in energy of thehard binary of a given osculating semi-major axis and ec-centricity due to the influence of the distant third body.The explicit expression for E P , in the case of a perturber ina circular orbit, reads E P = − Gµm p R (cid:0) − ( I ) sin ( ω ) (cid:1) = − Gµm p R G ( I , ω ) (30)where the function G describes the angle dependence of theenergy shift, µ = m m /M , m is the mass of the perturber, R its distance to the centre of mass of the hard binary, I the inclination with respect to the orbital plane of the per-turber, and ω the angle of the periapsis with respect to theascending or descending nodes of the orbital plane of theperturber. The remarkable fact proven in Will (2014) is thatthis expression remains constant over time-scales at whichperihelion precession can significantly change the orientationof the orbit with respect to the plane of the perturber. Theterms that would normally vary in ω are cancelled by crossterms between the PN dynamics of the inner binary and thetidal forces induced by the perturber. Indeed, the second sin term is evaluated at an initial perihelion orientation,and therefore is not affected by perihelion precession. Thisresult is very convenient for our purposes, because it is es-sentially of the same form as the change in initial conditionscaused by a PN perturbation.We can estimate how the change in energy of the hardbinary caused by the perturber will affect the inspiral time-scale by readjusting the initial conditions of the GW-induceddecay. We define a fictitious semi-major axis a P and eccen-tricity e P , constructed from both the Newtonian energy ofthe hard binary, E N , and the constant shift induced by theperturber, E P : a P = − GM E N + E P ) ≈ a − m µa MR (1 − e ) G ( I , ω ) , (31) e ≈ e − E P E N (1 − e ) , (32)An estimate for the change in the inspiral time-scalecan by found by evaluating a particular vacuum time-scale t GW at the fictitious semi-major axis and dividing it by thecorresponding vacuum time-scale of the physical orbit: t GW ( a P , e P ) t GW ( a, e ) ≈ (cid:18) E N + E P E N (cid:19) − ≈ − E P E N F ( e ) ≈ − m p G ( I , ω ) MR (1 − e ) F ( e ) (33)where, in the first step, we assumed quadrupole radiationor, equivalently, the evolution equations of Peters & Math-ews (1963). The function F can be computed from the ec-centricity dependence of the inspiral time-scale, and ranges from 1/5 to 1/2 for most values of the eccentricity. Fromthe results of Paper I and previous sections, we know thatPN corrections to the Newtonian vacuum time-scale canbe described by factors of the form exp (cid:0) C PN r S p − (cid:1) ≈ C PN r S p − + O [( r S p − ) n ] , where C PN is a coefficient oforder unity and n is the PN order. We can therefore estimatethe radius at which the environmental perturbation becomesof Newtonian (or PN) order by solving the equation m p G ( I , ω ) MR (1 − e ) F ( e ) = C PN (cid:18) r S p (cid:19) n (34)and find the characteristic radii R P / PN at which the effectof the perturbation becomes as large as that of a given PNorder: R P / PN = (cid:18) G ( I , ω ) m p n (1 − e ) Mr n S (cid:19) / , (35)where we approximated F ( e ) / ≈ and C / ≈ . For theNewtonian order ( n = 0 and C PN = 1 ), the radius reads R P/ N ≈ G ( I , ω ) / (cid:16) m M (cid:17) / a (1 − e ) / , (36)which can only be of order of or smaller than the semi-majoraxis of the hard binary itself, breaking the assumption ofa hierarchical triplet. In general, every additional PN orderonly adds a factor of a few over the previous value and, sincewe are considering hierarchical triplets (where R (cid:29) a ), wecan conclude that the perturbation has very little influenceon the inspiral time-scale.On the other hand, if the accumulated signal-to-noiseratio (SNR) of the source is high enough, LISA might beable to detect deviations from very high PN orders. In Fig-ure 4, we show the minimum order of a PN-waveform tem-plate that is required to detect a shift in the time-scale fora perturber of a given mass and separation, in the case ofan inner binary composed of an SMBH of M (cid:12) and acompact object of 10 M (cid:12) . We compute it by solving Equa-tion (35) for the PN order n and adding one. If we assumethat the SNR is sufficient to reach current state-of-the-artwaveform templates with spin (depicted in Figure 4; Fujita2015), we find values that are similar to the results in Yuneset al. (2011), wherein it is shown that detecting the dephas-ing in the signal of such a system is in principle possible,but only realistic for very massive ( ∼ M (cid:12) ) and close(sub-parsec) perturbers.The shift in energy is not the only way in which a thirdbody can perturb the hard binary. For orbits that are suf-ficiently misaligned, Kozai–Lidov (KL; Kozai 1962; Lidov1962) oscillations can induce a secular variation in the ec-centricity of the orbit. If the inspiral time-scale is longerthan the time-scale of KL oscillations, t KL , the evolutionequations for the orbital parameters of the hard binary canbe significantly affected by the perturber. The detectabilityof KL oscillations in the GW signal of LISA sources has beenvery recently studied in Randall & Xianyu (2019) and Deme MNRAS000 , 1–14 (2020) mproved gravitational radiation time-scales II R [pc]0510152025 R e q u i r e d P N o r d e r Sago & Fujita (2015)Fujita (2015)Munna (2020) m = 10 M m = 10 M m = 10 M Figure 4.
We show the minimum required PN order of a wave-form template needed to detect a shift in the inspiral time-scaleof an EMRI with an orbital frequency in the LISA band ( f ≈ mHz) caused by the perturbation of a third body with mass m ,in the case of an inner binary composed of an SMBH of M (cid:12) and a compact object of 10 M (cid:12) . We calculate the value by solv-ing Equation (35) for the PN order n and adding one. The solidlines represent a quasi-circular inspiral, whereas the dotted linesrepresent orbits that have the same frequency but an eccentricityof 0.9. The horizontal lines stand for the state-of-the-art results inthe calculations of PN fluxes using different methods, for spinningBHs in eccentric (4 PN) and circular (11 PN) orbits (Fujita 2015;Sago & Fujita 2015), and also eccentric orbits for non-spinningBHs (Munna 2020). et al. (2020), where the original version of Peters’ time-scale is used as a simple model for the inspiral time-scale.An intermediate calculation in the aforementioned papersis the computation of the characteristic orbital separation, a KLO / GW , of the inner binary at which KL oscillations arequenched by GW emission. The radius is found by equatingthe GW and KL time-scales and solving for the semi-majoraxis: t GW ∼ R ( e ) Q h ( a, e ) Q s ( a, e, s ) t P = t KL = 2 π √ GMGm R a / (37) = ⇒ a KLO / GW ≈ (cid:18) πMqR f ( e )( GM ) / c m R ( e ) Q h ( a, e ) Q s ( a, e, s ) (cid:19) / , (38)where we model the corrections to Peters’ time-scale withthe correction factors R ( e ) Q h ( a, e ) Q s ( a, e, s ) . The generalsolution a KLO / GW must be computed numerically. However,we can see that (while suppressed by the small power) thecorrection factors can reduce the value of the quenching ra-dius a KLO / GW by a factor of up to ∼ a KLO / GW . An example would be the phase- space volume a / GW , which can easily change by a factorof 2–3. Gas drag and torques in accretion discs are expected to sig-nificantly influence the dynamics of an inspiral event andmight cause a detectable dephasing of the GW signal de-pending on the physical properties of the disc (see, e.g. Ba-rausse et al. 2014; Derdzinski et al. 2019). While severaldifferent disc models exist, many show a central region thatextends for ∼ × Schwarzschild radii, where the den-sity and temperature profiles are much more shallow thanin the external regions (see the profiles in, e.g. Shakura &Sunyaev 1973; Sirko & Goodman 2003; Thompson 2009).If we restrict our analysis to orbits within this region, it isa reasonable first-order approximation to describe the lo-cal density and temperature as constant. To describe theinfluence of the disc on the inspiral time-scale, we need tomodel the effects that have a direct influence on the orbitalelements of the secondary BH.One such effect is gas drag. While the subject of dynam-ical friction is complex, we will assume that, for BHs, gasdrag can be locally modelled by the Bondi–Hoyle–Lyttleton(BHL) drag (see, e.g. Hoyle & Lyttleton 1939; Bondi &Hoyle 1944; Bondi 1952; Ostriker 1999; Fabj et al. 2020).At every completed orbit, the energy of the secondary BHwill be reduced by an amount of work W BHL given by W BHL = (cid:90) S F BHL d s , (39)where S is the path of the orbit in space. The magnitude ofthe BHL drag reads F BHL = 4 πG m ρc + v , (40)where c s is the local speed of sound, ρ is the local disc den-sity, m is the mass of the inspiralling BH, and v rel the rel-ative velocity with respect to the Keplerian disc velocity. Ifthe orbit is very eccentric, the instantaneous velocity, v i , ofthe secondary will generally be misaligned with the Keple-rian disc velocity. In this case we can approximate v rel ≈ v i and v i (cid:29) c s , which reduces the BHL drag formula to thestandard high-velocity limit of dynamical friction (see, e.g.Mayer 2013, for a discussion of gas effects in highly eccentricorbits). Note that this might also be the more appropriatedescription, at least in a local sense, in discs that are veryturbulent, where gas velocity deviates from a Keplerian disc(but only if the sound speed is also smaller than the BHvelocity). We will consider both the BHL case and the high-velocity dynamical friction case, keeping in mind that theformer is more accurate for circular orbits and the latter ismore accurate for eccentric ones. The integral can only besolved numerically or as a series expansion in the eccentric-ity. We compute it for the two cases: MNRAS , 1–14 (2020) L. Zwick et al. W BHL = 8 aG M q π ρc W BHL ( e, n M ) , (41) W dyn = 8 a GMq π ρ W dyn ( e ) , (42)where W BHL is a function that depends on the eccentricity ofthe orbit and an effective average Mach number n M = v c − evaluated at a radius of one semi-major axis, where v d is theKeplerian disc velocity. It varies between one (for circularorbits) and several thousands, depending on the eccentric-ity, and converges very poorly for high Mach numbers. Thefunction W dyn only depends on eccentricity and can be com-puted as a polynomial expression in e : W dyn ( e ) = 1 + 34 e + 2164 e + 55256 e + O [10 − e ] . (43)The loss of energy due to drag results in a secular changeof the semi-major axis of the secondary BH, which we cancompare to similar effects caused by the PN fluxes. To obtainan expression for the orbit-averaged energy flux due to BHLdrag and dynamical friction, we simply divide Equation (41)by an orbital period T .Keeping this result in mind, we turn our attention toanother effect that directly influences the orbit of the sec-ondary BH, namely global torques. A realistic descriptionof torques in accretion discs is only possible through hydro-dynamical simulations. Nevertheless, we will assume thatthe simplest analytical models of Type I and Type II viscoustorques are an appropriate order-of-magnitude approxima-tion, when averaged over many orbits. We take the formulasby Tanaka et al. (2002) and Derdzinski et al. (2021), whichread ˙ L I = T − ρhr ( t ) q (cid:18) v d c s (cid:19) = ρhr ( t ) G m πa c , (44) ˙ L II = α πr ( t ) Ω c s ρh , (45)where Ω is the Keplerian frequency of the secondary and α is the viscous parameter, and we assume that the disc scaleheight, h , follows the standard scaling h ∼ c s Ω − , where Ω is the Keplerian disc orbital frequency. We average theexpressions over one orbit, yielding (cid:68) ˙ L i (cid:69) = 1 T (cid:90) T ˙ L i ˙ φ − dφ (46) = ⇒ (cid:68) ˙ L I (cid:69) = √ a G M q ρ π c s D I ( e ) (47) = ⇒ (cid:68) ˙ L II (cid:69) = α a / c ρ √ GM D II ( e ) (48)where D i ( e ) are functions of the eccentricity. The function D I ( e ) is approximated very well by the first four orders ofits polynomial expansion, whereas D II ( e ) has a closed formsolution: D I ( e ) = 1 + 9916 e + 34651024 e + 115516384 e + O [10 − e ] , (49) D II ( e ) = 1 + 152 e + 458 e + 516 e . (50)In the case of torques, the orbit-averaging procedure isnot completely legitimate: Type I torques are a consequenceof resonances between the orbital frequency of the secondaryand the disc, and as such it is unlikely that they would real-istically appear for eccentric orbits. Type II torques becomerelevant when the secondary is massive enough to open gapsin the disc structure (see, e.g. Goldreich & Tremaine 1980;Ward 1997). In a tight binary of SMBHs, the gap openingtorques result into the formation of a cavity, and of a cir-cumbinary disc around it (e.g. Cuadra et al. 2006; Roediget al. 2012; Tang et al. 2017). These torques depend on thedetailed thermodynamics and viscous properties of the discbecause both pressure forces and viscous forces oppose thegap-opening torque induced by the binary into the surround-ing disc gas. Therefore, one would in principle need to ac-count for the dependence of the torques on the many discproperties that determine the strength of the forces at play.However, here we take an explorative stand as we only wantto highlight order-of-magnitude effects, which could laterbe explored in a self-consistent model for the disc. In gen-eral, Type II torques will apply only to SMBH binaries asopposed to IMRIs/EMRIs, because there must be a massiveenough secondary to open a gap. The recent sub-pc scale hy-drodynamical simulations of Souza Lima et al. (2020) showthat even with a mass ratio of 20:1 the opening of a gapor cavity, is not always occurring, whereas Derdzinski et al.(2021) show that, for IMRIs, the orbital evolution proceedswith no gap opening. Additionally, we will also see that thefunctions D i do not play a significant role for the results ofthis section.Equations (47) and (48) are in the form of an orbit-averaged angular momentum flux, which we can comparedirectly to the PN fluxes. Note that, in the case of circularorbits, one can obtain results that are consistent with anyparticular disc model by simply replacing ρ and c s with thecorresponding profiles. For the eccentric case, the orbit aver-ages can be computed numerically for arbitrary disc models.In order to compare a large range of PN orders to theresults shown above, we have to find a simple way to char-acterise their scaling in the semi-major axis and eccentricityof the secondary’s orbit. By checking up to the first fourinteger-PN orders of both energy and angular momentumfluxes, we always find the following structure: In the following calculations, we will always assume that thespin of the primary BH does not evolve due to gas accretion, atleast over the time-scale of a single inspiral event. In this case, wecan use the vacuum solutions that were discussed in the previoussections. MNRAS000
We show the minimum required PN order of a wave-form template needed to detect a shift in the inspiral time-scaleof an EMRI with an orbital frequency in the LISA band ( f ≈ mHz) caused by the perturbation of a third body with mass m ,in the case of an inner binary composed of an SMBH of M (cid:12) and a compact object of 10 M (cid:12) . We calculate the value by solv-ing Equation (35) for the PN order n and adding one. The solidlines represent a quasi-circular inspiral, whereas the dotted linesrepresent orbits that have the same frequency but an eccentricityof 0.9. The horizontal lines stand for the state-of-the-art results inthe calculations of PN fluxes using different methods, for spinningBHs in eccentric (4 PN) and circular (11 PN) orbits (Fujita 2015;Sago & Fujita 2015), and also eccentric orbits for non-spinningBHs (Munna 2020). et al. (2020), where the original version of Peters’ time-scale is used as a simple model for the inspiral time-scale.An intermediate calculation in the aforementioned papersis the computation of the characteristic orbital separation, a KLO / GW , of the inner binary at which KL oscillations arequenched by GW emission. The radius is found by equatingthe GW and KL time-scales and solving for the semi-majoraxis: t GW ∼ R ( e ) Q h ( a, e ) Q s ( a, e, s ) t P = t KL = 2 π √ GMGm R a / (37) = ⇒ a KLO / GW ≈ (cid:18) πMqR f ( e )( GM ) / c m R ( e ) Q h ( a, e ) Q s ( a, e, s ) (cid:19) / , (38)where we model the corrections to Peters’ time-scale withthe correction factors R ( e ) Q h ( a, e ) Q s ( a, e, s ) . The generalsolution a KLO / GW must be computed numerically. However,we can see that (while suppressed by the small power) thecorrection factors can reduce the value of the quenching ra-dius a KLO / GW by a factor of up to ∼ a KLO / GW . An example would be the phase- space volume a / GW , which can easily change by a factorof 2–3. Gas drag and torques in accretion discs are expected to sig-nificantly influence the dynamics of an inspiral event andmight cause a detectable dephasing of the GW signal de-pending on the physical properties of the disc (see, e.g. Ba-rausse et al. 2014; Derdzinski et al. 2019). While severaldifferent disc models exist, many show a central region thatextends for ∼ × Schwarzschild radii, where the den-sity and temperature profiles are much more shallow thanin the external regions (see the profiles in, e.g. Shakura &Sunyaev 1973; Sirko & Goodman 2003; Thompson 2009).If we restrict our analysis to orbits within this region, it isa reasonable first-order approximation to describe the lo-cal density and temperature as constant. To describe theinfluence of the disc on the inspiral time-scale, we need tomodel the effects that have a direct influence on the orbitalelements of the secondary BH.One such effect is gas drag. While the subject of dynam-ical friction is complex, we will assume that, for BHs, gasdrag can be locally modelled by the Bondi–Hoyle–Lyttleton(BHL) drag (see, e.g. Hoyle & Lyttleton 1939; Bondi &Hoyle 1944; Bondi 1952; Ostriker 1999; Fabj et al. 2020).At every completed orbit, the energy of the secondary BHwill be reduced by an amount of work W BHL given by W BHL = (cid:90) S F BHL d s , (39)where S is the path of the orbit in space. The magnitude ofthe BHL drag reads F BHL = 4 πG m ρc + v , (40)where c s is the local speed of sound, ρ is the local disc den-sity, m is the mass of the inspiralling BH, and v rel the rel-ative velocity with respect to the Keplerian disc velocity. Ifthe orbit is very eccentric, the instantaneous velocity, v i , ofthe secondary will generally be misaligned with the Keple-rian disc velocity. In this case we can approximate v rel ≈ v i and v i (cid:29) c s , which reduces the BHL drag formula to thestandard high-velocity limit of dynamical friction (see, e.g.Mayer 2013, for a discussion of gas effects in highly eccentricorbits). Note that this might also be the more appropriatedescription, at least in a local sense, in discs that are veryturbulent, where gas velocity deviates from a Keplerian disc(but only if the sound speed is also smaller than the BHvelocity). We will consider both the BHL case and the high-velocity dynamical friction case, keeping in mind that theformer is more accurate for circular orbits and the latter ismore accurate for eccentric ones. The integral can only besolved numerically or as a series expansion in the eccentric-ity. We compute it for the two cases: MNRAS , 1–14 (2020) L. Zwick et al. W BHL = 8 aG M q π ρc W BHL ( e, n M ) , (41) W dyn = 8 a GMq π ρ W dyn ( e ) , (42)where W BHL is a function that depends on the eccentricity ofthe orbit and an effective average Mach number n M = v c − evaluated at a radius of one semi-major axis, where v d is theKeplerian disc velocity. It varies between one (for circularorbits) and several thousands, depending on the eccentric-ity, and converges very poorly for high Mach numbers. Thefunction W dyn only depends on eccentricity and can be com-puted as a polynomial expression in e : W dyn ( e ) = 1 + 34 e + 2164 e + 55256 e + O [10 − e ] . (43)The loss of energy due to drag results in a secular changeof the semi-major axis of the secondary BH, which we cancompare to similar effects caused by the PN fluxes. To obtainan expression for the orbit-averaged energy flux due to BHLdrag and dynamical friction, we simply divide Equation (41)by an orbital period T .Keeping this result in mind, we turn our attention toanother effect that directly influences the orbit of the sec-ondary BH, namely global torques. A realistic descriptionof torques in accretion discs is only possible through hydro-dynamical simulations. Nevertheless, we will assume thatthe simplest analytical models of Type I and Type II viscoustorques are an appropriate order-of-magnitude approxima-tion, when averaged over many orbits. We take the formulasby Tanaka et al. (2002) and Derdzinski et al. (2021), whichread ˙ L I = T − ρhr ( t ) q (cid:18) v d c s (cid:19) = ρhr ( t ) G m πa c , (44) ˙ L II = α πr ( t ) Ω c s ρh , (45)where Ω is the Keplerian frequency of the secondary and α is the viscous parameter, and we assume that the disc scaleheight, h , follows the standard scaling h ∼ c s Ω − , where Ω is the Keplerian disc orbital frequency. We average theexpressions over one orbit, yielding (cid:68) ˙ L i (cid:69) = 1 T (cid:90) T ˙ L i ˙ φ − dφ (46) = ⇒ (cid:68) ˙ L I (cid:69) = √ a G M q ρ π c s D I ( e ) (47) = ⇒ (cid:68) ˙ L II (cid:69) = α a / c ρ √ GM D II ( e ) (48)where D i ( e ) are functions of the eccentricity. The function D I ( e ) is approximated very well by the first four orders ofits polynomial expansion, whereas D II ( e ) has a closed formsolution: D I ( e ) = 1 + 9916 e + 34651024 e + 115516384 e + O [10 − e ] , (49) D II ( e ) = 1 + 152 e + 458 e + 516 e . (50)In the case of torques, the orbit-averaging procedure isnot completely legitimate: Type I torques are a consequenceof resonances between the orbital frequency of the secondaryand the disc, and as such it is unlikely that they would real-istically appear for eccentric orbits. Type II torques becomerelevant when the secondary is massive enough to open gapsin the disc structure (see, e.g. Goldreich & Tremaine 1980;Ward 1997). In a tight binary of SMBHs, the gap openingtorques result into the formation of a cavity, and of a cir-cumbinary disc around it (e.g. Cuadra et al. 2006; Roediget al. 2012; Tang et al. 2017). These torques depend on thedetailed thermodynamics and viscous properties of the discbecause both pressure forces and viscous forces oppose thegap-opening torque induced by the binary into the surround-ing disc gas. Therefore, one would in principle need to ac-count for the dependence of the torques on the many discproperties that determine the strength of the forces at play.However, here we take an explorative stand as we only wantto highlight order-of-magnitude effects, which could laterbe explored in a self-consistent model for the disc. In gen-eral, Type II torques will apply only to SMBH binaries asopposed to IMRIs/EMRIs, because there must be a massiveenough secondary to open a gap. The recent sub-pc scale hy-drodynamical simulations of Souza Lima et al. (2020) showthat even with a mass ratio of 20:1 the opening of a gapor cavity, is not always occurring, whereas Derdzinski et al.(2021) show that, for IMRIs, the orbital evolution proceedswith no gap opening. Additionally, we will also see that thefunctions D i do not play a significant role for the results ofthis section.Equations (47) and (48) are in the form of an orbit-averaged angular momentum flux, which we can comparedirectly to the PN fluxes. Note that, in the case of circularorbits, one can obtain results that are consistent with anyparticular disc model by simply replacing ρ and c s with thecorresponding profiles. For the eccentric case, the orbit aver-ages can be computed numerically for arbitrary disc models.In order to compare a large range of PN orders to theresults shown above, we have to find a simple way to char-acterise their scaling in the semi-major axis and eccentricityof the secondary’s orbit. By checking up to the first fourinteger-PN orders of both energy and angular momentumfluxes, we always find the following structure: In the following calculations, we will always assume that thespin of the primary BH does not evolve due to gas accretion, atleast over the time-scale of a single inspiral event. In this case, wecan use the vacuum solutions that were discussed in the previoussections. MNRAS000 , 1–14 (2020) mproved gravitational radiation time-scales II ˙ E n − PN ∝ ˙ E N (cid:16) r S a (cid:17) n (cid:18) C n ( e )(1 − e ) n (cid:19) , (51) ˙ L n − PN ∝ ˙ L N (cid:16) r S a (cid:17) n (cid:18) D n ( e )(1 − e ) n (cid:19) , (52)where the subscript n denotes the (integer) PN order and thesubscript N denotes the lowest-order (quadrupolar) fluxes,whereas the coefficients C i and D i are fractions of polynomi-als in e of order unity. The half-integer terms have slightlydifferent forms, but we will neglect them fur the purpose ofthe following calculations.To find a series of characteristic orbital separations, atwhich the effect of the gas on the evolution of the orbit isas strong as the n -th order PN fluxes, we equate the expres-sions and solve for the semi-major axis. Note that this isequivalent to equating the GW-induced inspiral time-scaleto a characteristic inspiral time-scale due to the effect ofthe gas. In the case of BHL drag, we can solve the equa-tion analytically only for the case of circular orbits, becausein general the function W BHL contains an arbitrary amountof powers of the semi-major axis a . For circular orbits, thecharacteristic semi-major axis a BHL / PN reads a e → / PN ≈ r S (cid:18) c Gr ρ (cid:19) n , (53)where we are able to neglect all the small numerical coeffi-cients because of the small exponent. For dynamical friction,the radii a dyn / PN read a dyn / PN ≈ r S (cid:18) c Gr ρ (cid:19) n (cid:18) f ( e )(1 − e ) − n W dyn ( e ) (cid:19) n , (54) a dyn / PN ≈ r S (cid:18) c Gr ρ (cid:19) n (cid:0) − e (cid:1) − n n , (55)where, in the last step, we approximated [ f ( e ) / W dyn ( e )] / (cid:39) . We will use the BHL resultsas the low-eccentricity limit and the dynamical frictionresults as the high-eccentricity limit.We repeat the same calculation for the case of Type Itorques and obtain similar characteristic separations a TI / PN as before. In this case, we are able to solve it analytically forarbitrary eccentricities, with the caveat that global torquesare unlikely to appear for non-circular orbits: a TI / PN ≈ r S (cid:18) c s cGr ρ (cid:19) n (cid:18) (1 − e ) − − n D I ( e ) (cid:19) n (56) ≈ r S (cid:18) c s cGr ρ (cid:19) n (cid:0) − e (cid:1) − n n . (57)Here we have a slightly different scaling in the physicalproperties of the disc, but otherwise the formulas look verysimilar. Equation (57) is a good approximation even for thewhole eccentricity range, because of the small range in values a [ r S ] t G W < y r s t G W < y r s t G W < y r s a TI/PN a TII/PN a BHL/PN , n M a dyn/PN ISCO
Figure 5.
Shown is a region of phase space in which the char-acteristic separations a TI / PN , a TII / PN , a BHL / PN , and a dyn / PN are plotted. From the highest to the lowest, curves of the samecolour separate regions of phase space in which the effect of gas isstronger than the 0th, first, second, and third PN order effects. Inother words, one can expect this particular GW source to decayin a time that is well described by n -th order PN theory only ifit starts from a location below the n -th set of coloured lines inphase space. The results shown are for an SMBH of M (cid:12) and asecondary of M (cid:12) embedded in an accretion disc with a centralmean density of ρ = 10 − g cm − and a viscous coefficient α of − , consistent with the models in Thompson (2009) . Changingthe density and viscosity parameter of the central part of the discby orders of magnitude only slightly shifts the coloured curvesup and down. The black line marks the location of the ISCO fora non-rotating central object, and the dotted lines denote theisochrone curves for different GW decay times. of the function D I ( e ) in the interval e ∈ [0 , . The sameis true for Type II torques, where the characteristic radii a TII / PN read a TII / PN ≈ r S (cid:18) c q αc Gr ρ (cid:19) n (cid:18) (1 − e ) − − n D II ( e ) (cid:19) n (58) ≈ r S (cid:18) c q αc Gr ρ (cid:19) n (cid:0) − e (cid:1) − n n . (59)With the aid of these characteristic radii, we can iden-tify different regions of phase space in which we can expectcertain PN orders to be accurate descriptions of the inspi-ral time-scale. The Newtonian ( n = 0 ) result is particularlytelling, since it separates the regions where gas effects dom-inate from regions where GW emission dominates.Figure 5 shows this procedure for the standard case ofa M (cid:12) SMBH and a M (cid:12) compact object. The resultsonly scale very weakly with disc properties such as centraldensity and temperature, suggesting that they are robust asan order-of-magnitude estimate. For example, one can usethe formulas for the characteristic radii to estimate devia- MNRAS , 1–14 (2020) L. Zwick et al. tions from the vacuum solutions as a function of the discproperties. As an example, let us take a M (cid:12) BH on acircular orbit at a radius of r S . If it is embedded in adisc with a local density of ρ = 10 − g cm − and viscosityparameter α = 10 − (Thompson 2009), we expect the effectof BHL drag (if present) to be almost as important as thelowest-order quadrupole radiation, while Type I or Type IItorques (if present) would only be as important as a 1 PNcorrection.Vice versa, one can also estimate the disc properties re-quired for a shift to be visible at a certain PN order, whichwe show in Figure 6 for varying values of the local disc den-sity. One can also use this type of analysis to speculate onthe likely environment of actual LISA sources. As an exam-ple, if we assume that the density of a typical circumnu-clear disc (CND) scales as a simple r − power law (Mestel1963) until it flattens out at a radius of ∼ × r S , wecan extrapolate some realistic densities from observationalresults for the total mass and extension of CNDs given inMedling et al. (2014), obtaining ρ ∼ − to − g cm − (a result more consistent with the disc models in Sirko &Goodman 2003). Note that these assumptions are also usedas initial conditions in state-of-the-art simulations of SMBHcoalescence (see, e.g. Souza Lima et al. 2017, 2020). Evenfor the wide range of observed CND masses ( ∼ M (cid:12) to ∼ × M (cid:12) ), the required PN order to detect a shift in theinspiral time-scale for an EMRI that is entering the LISAband (at ∼ ∼ > g cm − are present, a dephasing ofthe GWs is likely to be detectable by LISA. Note that thistype of calculation can be repeated with more sophisticateddensity profiles to yield more precise results. However, evenvery simple estimates such as assuming a r − density profileuntil × r S can yield correct order-of-magnitude esti-mates because of the weak scaling of the characteristic radiiwith disc properties. In this paper, we extended our analysis of the GW-induceddecay time-scale to include the 1.5 PN order. We produceda series of analytical fits that can be multiplied by Peter’sformula in order to model hereditary and spin-orbit effects(Equations 22 to 24). In the high-eccentricity (Equation 28)and low-eccentricity (Equation 26) limits, the fits reduce tosimple exponential functions that can be used to avoid er-rors of a factor ∼
10 and ∼
2, respectively. We show thatthese modifications to Peters’ formula are necessary in ap-plications that are highly sensitive to the inspiral time-scale,such as the calculation of the EMRI critical semi-major axis(Amaro-Seoane et al. 2013; Amaro-Seoane 2019; see alsoVázquez-Aceves et al., in prep.).We then used the results of Will (2014) in combinationwith our findings to compare the strength of arbitrarily highPN perturbations against gravitational perturbations in ahierarchical triple system. We found a series of characteris-tic radii at which the perturbations affect the inspiral time- [g/cm ]1234567 R e q u i r e d P N o r d e r BHL dragType I torqueType II torque M CND = 10 M M CND = 10 M M CND = 10 M Figure 6.
We show the required PN order needed to detect ashift from the vacuum time-scale in the cases of BHL drag, Type Iand Type II torques for different values of the local density. Thelines are computed for the standard case of a M (cid:12) SMBH anda 10 M (cid:12) secondary entering the LISA band at ∼ ∼ × cm s − . Thevertical lines represent a qualitative extrapolation of the localdensity from observed CND masses and extensions. The resultsonly scale very weakly with parameters such as the CND mass,local sound speed, or the mass ratio of the binary, but changesignificantly with the SMBH mass. This is, however, not due tothe properties of the disc, but rather to the fact that the efficiencyof GW emission at a fixed frequency strongly depends on the massof the SMBH. scale as much as any particular PN order (Equation 35),which we used to reproduce previously known results on thedetectability of the dephasing of GW signals (Yunes et al.2011). Furthermore, we briefly discussed the effects of KLoscillation and how corrections to Peters’ time-scale can af-fect recent calculations found in Randall & Xianyu (2019)and Deme et al. (2020) that depend on the typical radius atwhich KL oscillations are quenched by GW emission.Finally, we discussed the validity of vacuum time-scalesfor inspirals embedded in a gaseous medium, which is thelikely case for SMBHs sourrounded by accretion discs andCNDs. We analysed several types of drag forces (BHL dragand dynamical drag) and global torques (Type I torquesand Type II viscous torques), and compared them againstPN energy and angular momentum fluxes of various orders.We found several characteristic radii that separate the phasespace in different dynamical regions, where gas effects are asstrong as a particular PN order (Equations 53 to 59). Weused a simple prescription to extrapolate the gas densitiesthat would be likely encountered by a GW source enter-ing the LISA band in one of the several observed CNDs(Medling et al. 2014), and found that PN orders of 3–5 arerequired to detect a deviation from the vacuum case, whichis compatible with results from hydrodynamical simulations(Derdzinski et al. 2019, 2021). While many assumptions arerequired to derive the analytical formulas, the results scale MNRAS000
We show the required PN order needed to detect ashift from the vacuum time-scale in the cases of BHL drag, Type Iand Type II torques for different values of the local density. Thelines are computed for the standard case of a M (cid:12) SMBH anda 10 M (cid:12) secondary entering the LISA band at ∼ ∼ × cm s − . Thevertical lines represent a qualitative extrapolation of the localdensity from observed CND masses and extensions. The resultsonly scale very weakly with parameters such as the CND mass,local sound speed, or the mass ratio of the binary, but changesignificantly with the SMBH mass. This is, however, not due tothe properties of the disc, but rather to the fact that the efficiencyof GW emission at a fixed frequency strongly depends on the massof the SMBH. scale as much as any particular PN order (Equation 35),which we used to reproduce previously known results on thedetectability of the dephasing of GW signals (Yunes et al.2011). Furthermore, we briefly discussed the effects of KLoscillation and how corrections to Peters’ time-scale can af-fect recent calculations found in Randall & Xianyu (2019)and Deme et al. (2020) that depend on the typical radius atwhich KL oscillations are quenched by GW emission.Finally, we discussed the validity of vacuum time-scalesfor inspirals embedded in a gaseous medium, which is thelikely case for SMBHs sourrounded by accretion discs andCNDs. We analysed several types of drag forces (BHL dragand dynamical drag) and global torques (Type I torquesand Type II viscous torques), and compared them againstPN energy and angular momentum fluxes of various orders.We found several characteristic radii that separate the phasespace in different dynamical regions, where gas effects are asstrong as a particular PN order (Equations 53 to 59). Weused a simple prescription to extrapolate the gas densitiesthat would be likely encountered by a GW source enter-ing the LISA band in one of the several observed CNDs(Medling et al. 2014), and found that PN orders of 3–5 arerequired to detect a deviation from the vacuum case, whichis compatible with results from hydrodynamical simulations(Derdzinski et al. 2019, 2021). While many assumptions arerequired to derive the analytical formulas, the results scale MNRAS000 , 1–14 (2020) mproved gravitational radiation time-scales II only very weakly with disc properties, suggesting that theymight be a robust order-of-magnitude approximation of amore realistic treatment of gas-embedded inspirals. ACKNOWLEDGEMENTS
EB, PRC, LM, and LZ acknowledge support from theSwiss National Science Foundation under the Grant200020_178949. EB acknowledges support from the Euro-pean Research Council (ERC) under the European Union’sHorizon 2020 research and innovation program ERC-2018-COG under grant agreement N. 818691 (B Massive). PASacknowledges support from the Ramón y Cajal Programmeof the Ministry of Economy, Industry and Competitivenessof Spain, the COST Action GWverse CA16104, the NationalKey R&D Program of China (2016YFA0400702) and theNational Science Foundation of China (11721303).
DATA AVAILABILITY STATEMENT
The data underlying this article will be shared on reasonablerequest to the corresponding author.
APPENDIX A: PARAMETRISATION ANDEVOLUTION EQUATIONS
In this appendix, we list the PN parametrisation and evo-lution equations that were necessary for this work, adjustedto our notation.
First Order and Hereditary terms
From Ebersold et al. (2019), we take the 1 PN fluxes and1.5 PN hereditary fluxes radiation reaction equations: dxdt = 2 x (cid:16) X N + x X + x / X tail (cid:17) , (A1) de d dt = x (cid:16) E N + x E + x / E tail (cid:17) , (A2)where: X N = 1(1 − e ) / (cid:26)
965 + 292 e e (cid:27) , (A3) X = 1(1 − e ) / (cid:40) − − ν e (cid:18) − ν (cid:19) + e (cid:18) − ν (cid:19) + e (cid:18) − ν (cid:19) (cid:41) , (A4) X tail = x / φ ( e d ) , (A5) φ ( e d ) = 1 + 2335192 e + 42955768 e + 620464736864 e + O [ e ] , (A6) and E N = 1(1 − e ) / (cid:26) e (cid:27) , (A7) E = 1(1 − e ) / (cid:40) − − ν
45 + e (cid:18) − ν (cid:19) + e (cid:18) − ν (cid:19) (cid:41) , (A8) E tail = − (cid:18) − πx / ψ ( e d ) (cid:19) (A9) ψ ( e d ) = 192985 (cid:112) − e e (cid:20)(cid:113) − e φ ( e d ) − ˜ φ ( e d ) (cid:21) , (A10) ˜ φ ( e d ) = 1 + 20932 e + 2415128 e + 73075118432 e + O [ e ] . (A11) Spin-Orbit Coupling
From Klein & Jetzer (2010), we take the 1.5 PN spin-orbitcoupling parametrization equations, a = x − (cid:20) x / β (2 , √ − e (cid:21) , (A12) e = e d (cid:34) − x / β (2 , (cid:112) − e (cid:35) (A13)as well as the 1.5 PN spin-orbit coupling fluxes equations: dndt = νx / (1 − e ) / (cid:34) (cid:0)
96 + 292 e + 37 e (cid:1) − x /
10 (1 − e ) / β (3088 + 15 528 e + 7026 e + 195 e , e + 5982 e + 207 e ) (cid:35) , (A14) de dt = − νx (1 − e ) / (cid:34) e (cid:0)
304 + 121 e (cid:1) − e x /
15 (1 − e ) / β (13 048 + 12 000 e + 789 e , e + 835 e ) (cid:35) , (A15)where n = x / , (A16) β ( a, b ) = J J · ( aζ + bξ ) , (A17) ζ = S + S , (A18) ξ = m m S + m m S . (A19) MNRAS , 1–14 (2020) L. Zwick et al.
REFERENCES
Amaro-Seoane P., 2018, Living Reviews in Relativity, 21, 4Amaro-Seoane P., 2019, Physical Review D, 99Amaro-Seoane P., Gair J. R., Freitag M., Miller M. C., Mandel I.,Cutler C. J., Babak S., 2007, Classical and Quantum Gravity,24, R113Amaro-Seoane P., Sopuerta C. F., Freitag M. D., 2013, MNRAS,429, 3155Antonini F., Rasio F. A., 2016, ApJ, 831, 187Babak S., et al., 2017, Phys. Rev. D, 95, 103012Barack L., et al., 2019, Classical and Quantum Gravity, 36,143001Barausse E., Cardoso V., Pani P., 2014, Phys. Rev. D, 89, 104059Blanchet L., 2014, Living Reviews in Relativity, 17, 2Boetzel Y., Susobhanan A., Gopakumar A., Klein A., Jetzer P.,2017, Phys. Rev. D, 96, 044011Bondi H., 1952, MNRAS, 112, 195Bondi H., Hoyle F., 1944, MNRAS, 104, 273Bonetti M., Haardt F., Sesana A., Barausse E., 2016, MNRAS,461, 4419Bonetti M., Perego A., Capelo P. R., Dotti M., Miller M. C., 2018,Publ. Astron. Soc. Australia, 35, e017Bonetti M., Sesana A., Haardt F., Barausse E., Colpi M., 2019,MNRAS, 486, 4044Bortolas E., Mapelli M., 2019, MNRAS, 485, 2125Cuadra J., Nayakshin S., Springel V., Di Matteo T., 2006, MN-RAS, 366, 358Damour T., Deruelle N., 1985, Ann. Inst. Henri Poincaré Phys.Théor, 43, 107Deme B., Hoang B.-M., Naoz S., Kocsis B., 2020, ApJ, 901, 125Derdzinski A. M., D’Orazio D., Duffell P., Haiman Z., MacFadyenA., 2019, MNRAS, 486, 2754Derdzinski A., D’Orazio D., Duffell P., Haiman Z., MacFadyenA., 2021, MNRAS, 501, 3540Ebersold M., Boetzel Y., Faye G., Mishra C. K., Iyer B. R., JetzerP., 2019, Phys. Rev. D, 100, 084043Einstein A., 1916, Annalen der Physik, 354, 769Fabj G., Nasim S. S., Caban F., Ford K. E. S., McKernan B.,Bellovary J. M., 2020, MNRAS, 499, 2608Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2015, MN-RAS, 447, L80Fujita R., 2015, Progress of Theoretical and ExperimentalPhysics, 2015, 033E01Giacobbo N., Mapelli M., 2019, MNRAS, 482, 2234Giacobbo N., Mapelli M., 2020, ApJ, 891, 141Goldreich P., Tremaine S., 1980, ApJ, 241, 425Hoyle F., Lyttleton R. A., 1939, Proceedings of the CambridgePhilosophical Society, 35, 592Khan F. M., Capelo P. R., Mayer L., Berczik P., 2018, ApJ, 868,97Klein A., Jetzer P., 2010, Phys. Rev. D, 81, 124001Klein A., et al., 2016, Phys. Rev. D, 93, 024003Kozai Y., 1962, AJ, 67, 591Lense J., Thirring H., 1918, Physikalische Zeitschrift, 19, 156Lidov M. L., 1962, Planet. Space Sci., 9, 719Mayer L., 2013, Classical and Quantum Gravity, 30, 244008Medling A. M., et al., 2014, ApJ, 784, 70Memmesheimer R.-M., Gopakumar A., Schäfer G., 2004, Phys.Rev. D, 70, 104011Mestel L., 1963, MNRAS, 126, 553Munna C., 2020, Phys. Rev. D, 102, 124001O’Leary J. A., Moster B. P., Naab T., Somerville R. S., 2021,MNRAS, 501, 3215Ostriker E. C., 1999, ApJ, 513, 252Peters P. C., 1964, PhD thesis, California Institute of TechnologyPeters P. C., Mathews J., 1963, Physical Review, 131, 435Randall L., Xianyu Z.-Z., 2019, arXiv e-prints, p. arXiv:1902.08604Reynolds C. S., 2020, arXiv e-prints, p. arXiv:2011.08948Roedig C., Sesana A., Dotti M., Cuadra J., Amaro-Seoane P.,Haardt F., 2012, A&A, 545, A127Sago N., Fujita R., 2015, Progress of Theoretical and Experimen-tal Physics, 2015, 073E03Schäfer G., Jaranowski P., 2018, Living Reviews in Relativity, 21,7Schwarzschild K., 1916, Abh. Konigl. Preuss. Akad. Wis-senschaften Jahre 1906,92, Berlin,1907, 1916, 189Sesana A., Khan F. M., 2015, MNRAS, 454, L66Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33Sirko E., Goodman J., 2003, MNRAS, 341, 501Souza Lima R., Mayer L., Capelo P. R., Bellovary J. M., 2017,ApJ, 838, 13Souza Lima R., Mayer L., Capelo P. R., Bortolas E., Quinn T. R.,2020, ApJ, 899, 126Tanaka H., Takeuchi T., Ward W. R., 2002, ApJ, 565, 1257Tang Y., MacFadyen A., Haiman Z., 2017, MNRAS, 469, 4258Thompson T. A., 2009, in Wang W., Yang Z., Luo Z., ChenZ., eds, Astronomical Society of the Pacific ConferenceSeries Vol. 408, The Starburst-AGN Connection. p. 128( arXiv:0908.1756 )VanLandingham J. H., Miller M. C., Hamilton D. P., RichardsonD. C., 2016, ApJ, 828, 77Ward W. R., 1997, Icarus, 126, 261Wex N., 1995, Classical and Quantum Gravity, 12, 983Will C. M., 2014, Classical and Quantum Gravity, 31, 244001Yunes N., Berti E., 2008, Phys. Rev. D, 77, 124006Yunes N., Coleman Miller M., Thornburg J., 2011, Phys. Rev. D,83, 044030Zhang Z., Yunes N., Berti E., 2011, Phys. Rev. D, 84, 024029Zwick L., Capelo P. R., Bortolas E., Mayer L., Amaro-Seoane P.,2020, MNRAS, 495, 2321This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000