Improved gravitational radiation time-scales II: spin-orbit contributions and environmental perturbations

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 the convenience of a simple formula outweighs the need for precision. However, many promising sources of the Laser Interferometer Space Antenna (LISA), such as supermassive black hole binaries and extreme mass-ratio inspirals (EMRIs), are expected to enter the LISA band with highly eccentric ($e \gtrsim$ 0.9) and highly relativistic orbits. These are exactly the two limits in which Peters' estimate performs the worst. In this work, we expand upon previous results and give simple analytical fits to quantify how the inspiral time-scale is affected by the relative 1.5 post-Newtonian (PN) hereditary fluxes and spin-orbit couplings. We discuss several cases that demand a more accurate GW time-scale. We show how this can have a major influence on quantities that are relevant for LISA event-rate estimates, such as the EMRI critical semi-major axis. We further discuss two types of environmental perturbations that can play a role in the inspiral phase: the gravitational interaction with a third massive body and the energy loss due to dynamical friction and torques from a surrounding gas medium ubiquitous in galactic nuclei. With the aid of PN corrections to the time-scale in vacuum, we find simple analytical expressions for the regions of phase space in which environmental perturbations are of comparable strength to the effects of any particular PN order, being able to qualitatively reproduce the results of much more sophisticated analyses.


INTRODUCTION
With the launch of space-borne gravitational-wave (GW) detectors in sight (e.g. the Laser Interferometer Space Antenna, LISA; Amaro-Seoane et al. 2013; Barack et al. 2019), many questions still have to be answered with regards to the possible sources of signal. In particular, large uncertainties remain on the expected event rates for different types of sources and for different formation channels (e.g. Klein et al. 2016;Babak et al. 2017). An essential part of these estimates is a model of gravitational radiation and its effects on the E-mail: zwicklo@ics.uzh.ch orbital dynamics of a system. Very often, it is convenient to express this effect with a simple analytical time-scale.
Already in the early 1960s, Peters & Mathews (1963) calculated the radiation reaction equations and extracted a simple analytical formula for the time-scale of GW-induced decay of a binary. This formula, commonly known as Peters' time-scale, has since seen an incredible amount of use, and more recently has been employed in many event-rate estimates and other source-modelling applications (see, e.g. Farris et al. 2015;VanLandingham et al. 2016;Bortolas & Mapelli 2019, amongst many others). Peters' formula is thought to suffice as a first approximation, especially because it is often used in conjunction with other dynamical time-scales that are themselves only order-of-magnitude es-timates (see, e.g. Sesana & Khan 2015). However, many of LISA most promising sources, such as supermassive black hole (SMBH) binaries and extreme/intermediate/extremely large mass-ratio inspirals (EMRIs/IMRIs/XMRIs), are expected to exhibit highly eccentric and highly relativistic orbits (see, e.g. Amaro-Seoane et al. 2007a;Antonini & Rasio 2016;Bonetti et al. 2016Bonetti et al. , 2018Khan et al. 2018;Amaro-Seoane 2018;Bonetti et al. 2019;Giacobbo & Mapelli 2019; Amaro-Seoane 2019), which are the two regimes in which Peters' estimate performs the worst. The commonly used formula does not account for the Newtonian evolution of the eccentricity. Moreover, it is based on the quadrupole radiation of a Keplerian binary and therefore fails to capture the more complex post-Newtonian (PN) behaviour of highly relativistic orbits.
In a recent paper (Zwick et al. 2020, hereafter Paper I), we sought to fix these issues, in service of reducing the unnecessary errors that might arise from using Peters' formula in its common form. We produced an updated formula that takes into account the eccentricity evolution and the lowestorder PN correction, while still being algebraically very simple.
In this paper, we wish to extend our analysis to the next order of the PN series, for the two following reasons. Firstly, it will serve as a proof of good convergence for the previous results. Secondly, two qualitatively new effects emerge at the next order: hereditary (i.e. tail) and spin-orbit contributions. Hereditary terms in the PN flux equations are known to be disproportionately important in the amount of cycles that some GW sources complete in the LISA band (see, e.g. Blanchet 2014). Moreover, a significant fraction of SMBHs are expected to have large ( ∼ > 0.9) dimensionless spin values (see, e.g. Reynolds 2020), whose effects are known to significantly alter the event rates of some GW sources, especially EMRIs (Amaro-Seoane et al. 2013). This, combined with the weak convergence of the PN series in the case of extreme mass ratios (see, e.g. Yunes & Berti 2008;Zhang et al. 2011), suggests that there is still precision to be gained beyond the lowest-order PN correction.
At some point, however, one inevitably clashes against the fact that most GW sources do not evolve in a pure vacuum. Even at very small separations, the inspiral process can be perturbed by the presence of other forces and the vacuum time-scale can be a precise estimate only if these environmental perturbations are weak. Therefore, we take a look at two types of environmental perturbations to the inspiral time-scale. Firstly, we investigate the gravitational influence of a third body. Secondly, we compare the strength of gas friction and gas torques for an inspiral embedded in a gaseous disc. In both cases, we express the results as a series of characteristic orbital separations at which the environmental effects influence the inspiral time-scale as much as a particular PN order. Even though these results are only order-of-magnitude estimates, they clearly illustrate the regions of validity and usefulness for the corrections we propose, and more generally for gravitational radiation timescales.
With the goals of this paper set, we wish to refer any reader not at least partly familiar with Peters' formula to the original work from 1963 and to Peters' Ph.D. thesis (Peters & Mathews 1963;Peters 1964), as well as to a general review of PN theory in Schäfer & Jaranowski (2018). We also refer to Paper I for a more in-depth discussion of the limitations of Peters' derivation and on the considerations required to go beyond Keplerian orbits and quadrupolar radiation. Section 2 presents a brief but complete overview of the concepts that were used in this paper. In Section 3, we present our derivation of the 1.5 PN correction factors to Peters' formula, and discuss how their inclusion can affect current event-rate estimates for EMRIs. In Section 4, we discuss environmental effects and derive several useful characteristic radii. Finally, in Section 5, we summarize our results.
2 THEORETICAL BACKGROUND 2.1 Two-body problem in post-Newtonian theory When going to higher PN orders, calculations become notoriously algebraically intensive. It is therefore convenient to work in natural units, wherein M = m1 + m2 = 1, G = 1, and c = 1, where m1 and m2 are the masses of the two bodies, G is the gravitational constant, and c is the speed of light in vacuum. In these units, we can express all expansions of the PN series as Taylor expansions in the so-called PN parameter x = −2E, where E is the Newtonian orbital energy per unit reduced mass µ = m1m2/M . The PN parameter x is also a gauge-invariant quantity (see, e.g. Blanchet 2014), meaning that the following calculations are valid regardless of gauge choice. With the units taken care of, two different ingredients are needed to fully describe the orbital evolution of a two-body system in PN theory.
Firstly, we require a parametrization of the orbital separation r(φ), which solves the PN equations of motion for a constant orbital energy E and angular momentum vector J: where φ is the polar angle and the first term of the righthand side is the solution to the familiar Kepler problem in Newtonian gravity. There are many equivalent possibilities to express such parametrizations, and in this paper we use the so-called quasi-Keplerian (QK) one (see Damour & Deruelle 1985;Wex 1995;Memmesheimer et al. 2004), in which the radial separation is expressed in terms of two generalised orbital parameters, a semi-major axis a and an eccentricity e. These parameters depend on the quantities of the binary system, such as the energy E, the angular momentum vector J, and, in the case of BHs, two spin vectors S1 and S2. 1 They serve to greatly simplify the notation, just as introducing a semi-major axis and an eccentricity simplifies the solution of the Newtonian Kepler problem. Secondly, we require a set of differential equations that describe how the aforementioned orbital parameters evolve due to the flux of energy and angular momentum caused by the emission of GWs. A convenient form is found when the orbital parameters a and e are expressed in terms of the PN parameter x and a dummy eccentricity parameter e d , that essentially is a measure of the angular momentum of the system: where the subscripts of the coefficients denote the PN order and some of the coefficients can vanish. The explicit coefficients used in this paper are listed in the Appendix, and can be found in Memmesheimer et al. (2004) for the spinless part and in Klein & Jetzer (2010) for an extension that covers the spin-orbit coupling. The evolution equations for x and e d take a relatively simple form: where, similarly to before, the subscripts of the coefficients denote the relative PN order and some of the coefficients can vanish. 2 The explicit coefficients used in this paper can be found in Boetzel et al. (2017) for the spin-less parts and in Klein & Jetzer (2010) for parts describing the spin-orbit coupling. Note that these evolution equations are always orbitaveraged, the assumption being that changes to the orbital elements occur on time-scales that are much longer than a single orbital period, which is the same assumption that is required to produce Peters' formula in the first place. A different approach would be to use a Hamiltonian integrator that evolves the positions of the binary's components step by step (see, e.g. Berentzen et al. 2008;Lubich et al. 2010;Ireland et al. 2019). While a direct numerical integration of the equations of motion is in principle more precise, we chose to adopt the adiabatic approximation for two reasons. Firstly, we want to investigate a large range of parameter space in mass, mass ratio, orbital separation, eccentricity, and spin. The repeated integration of accurate equations of motion at high PN orders would simply be too costly. Secondly, IMRI/EMRI/XMRI sources complete thousands (if not millions) of orbits before coalescence, meaning that numerical errors, which accumulate orbit by orbit, would be very large. In those cases, the adiabatic approximation, whose PN terms are found through the Teukolsky formalism (Teukolsky 1973;Amaro-Seoane et al. 2007b), is more efficient and more accurate than orbital integrators, which are indeed often tested and calibrated on close to equal-mass ratio binaries.
They describe a binary that is radiating GWs in accordance with Einstein (1916)'s quadrupole formula. Because of the coupled nature of these equations and the complexity of the eccentricity enhancement functions, it is impossible to solve them analytically for an arbitrary initial semi-major axis a0 and eccentricity e0. Nevertheless, Peters and Mathews gave some approximated solutions for a(t) and e(t) (see Peters & Mathews 1963;Peters 1964). In particular, the solutions for a(t) can be set equal to zero and solved for the time t in order to yield a decay time-scale of the orbital evolution.
The most commonly used formula for Peters' time-scale, tP, is derived by solving Equation (2) at lowest order while also neglecting the evolution of the eccentricity, thus settinġ e = 0. The solution reads where q = m1/m2 ≤ 1 is the mass ratio, and we reintroduced the physical constants. As stated in the introduction, this formula is currently used in a wide variety of applications, even though it is known to fail for highly eccentric and for highly relativistic orbits. In Paper I, we found two simple factors, R and Q f , that can be multiplied by Peters' time-scale to resolve its issues: The fitting function R(e0) = 8 1− √ 1−e 0 is purely Newtonian, and is needed to resolve the lack of eccentricity evolution in the standard formula. Other fits with the same purpose have been used in some select cases (see, e.g. Bonetti et al. 2018), but tend to be much more complex than the factor R and are therefore uncommon.
The factor Q f (a0, e0) is a simple fit that roughly accounts for both the more complex 1 PN parametrization given by Equations (2) Table 1. This table summarises the new elements that enter the picture at a given PN order. The acronyms QK, S-O, and S-S stand for "quasi-Keplerian", "spin-orbit", and "spin-spin", respectively. After the first PN order, one must in principle consider the evolution of both BH spins, S 1 and S 2 , as well as the specific angular momentum vectorĴ.
where rS is the effective Schwarzschild (1916) radius of the two-body system, rS = 2GM c −2 , and the coefficient 2.5 in the exponential function was chosen to fit a large range of initial eccentricities. In Paper I, we showed that Equation (10) can resolve errors of order 1-10 than can arise when applying Peters' formula to highly eccentric and relativistic orbits.

Assumptions needed for higher orders
BH spin couplings first appear at 1.5 PN order in the PN expansion. From then on, both spins have to be taken into account, drastically increasing the number of degrees of freedom and dynamical variables of the system. The spins can interact with the orbital angular momentum at 1.5 PN order and also with each other at 2 PN order. As is shown in Table 1, this increases the total number of dynamical variables from 2 (a and e) to 10. Sampling an appropriate range of parameter space with a grand total of ten initial conditions would simply be an exercise in frustration. At the cost of some loss of generality, however, we can 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 distinguished by ranges of their mass ratios: ∼10 −6 , ∼10 −4 to ∼10 −2 , and ∼10 −2 to 10 0 , respectively. For EMRIs and IM-RIs, the mass ratio is by definition very small, but even for SMBH mergers only a small minority will have mass ratios of order one (see, e.g. O'Leary et al. 2021, where we assume that galaxy mass is a tracer for SMBH mass). In light of applications for LISA sources, it is very natural to restrict our analysis to small mass ratios, i.e. q ∼ < 0.1. This turns out to be extremely convenient, because a small mass ratio suppresses most of the complex effects of spin-orbit and especially spin-spin couplings. We decided to also neglect the standard, non-Kerr 2 PN contributions to the parametrization and the fluxes of the two-body system. This last assumption is justified for two connected reasons: firstly, our numerical investigations showed that integer-PN orders tend to converge much faster than half-integer tail contributions (see also Blanchet 2014), making the 2 PN correction very small; secondly, our goal is to produce a correction to Peters' formula that is sufficiently simple to be used in the broader context of source modelling for GW observatories. Fitting the minute 2 PN deviations would necessarily require a much more complex fitting function.

Initial conditions and orbital evolution
A divergence-free PN parametrization of the orbital motion that includes spin-orbit and spin-spin couplings can be found in Klein & Jetzer (2010). Up to 1.5 PN order, it has the form of a standard QK parametrization: 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 they read where the 1 PN term is taken from equations (20a) and (20b) in Memmesheimer et al. (2004) and simplified according to the small mass-ratio limit. We wish to compare a Newtonian orbit and a PN orbit that start out with the same orbital elements at a given moment in time (t = 0) and are subsequently evolved with the appropriate energy and angular momentum fluxes. To enforce this, we set Equations (12) and (13) equal to their Newtonian counterparts. The only way for these two orbits to be identical at t = 0 is for them to have different values of the PN parameter x and the dummy eccentricity e d . We denote them as (xN, eN) for the Newtonian orbit and as (xPN, ePN) for the PN orbit. To find the correct values for xPN and ePN, we set the Newtonian and PN orbital elements equal to each other: We can solve these equations order by order and find the values of the PN parameters so that the two orbits are identical at the beginning of their evolution: As noted before, the energy and angular momentum fluxes can be averaged into secular evolution equations for the parameters x and e. The equations are organised into PN series. Schematically, they read x where the different subscripts denote the Newtonian and 1 PN contributions, as well as the lowest-order hereditary tail and the spin-orbit coupling contribution. The full expressions 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 evolution equations for the direction of the orbital angular momentum vector and the two spin vectors. However, at lowest order in the spin-orbit coupling and the mass ratio, the orbital angular momentum vector J only undergoes a Lense-Thirring (Lense & Thirring 1918) precession. Furthermore, all terms that containĴ are of the form where θ is the angle between the primary BH spin vector and the orbital angular momentum vector. The magnitude of the BH spin is assumed to be constant, because the evolution occurs in vacuum. Furthermore, the Lense-Thirring precession does not affect the angle θ, hence the expression in Equation (20) is constant in time. It follows then that Equations (18) and (19) are a closed system of differential equations, which can be integrated numerically, starting from the correct initial conditions given by Equations (16) and (17).
As we have noted before, there are two separate effects that influence the inspiral time-scale when considering higher-order PN perturbations. On the one hand, the evolution equations for the orbital parameters gain additional terms, which originate from the additional energy and angular momentum fluxes and, generally, tend to shorten the inspiral time-scale. On the other hand, some adjustments to the initial conditions are required to ensure that the orbital elements of a PN orbit and those of a Newtonian orbit are identical at the beginning of the evolution, which tend to increase the overall time-scale. The reason for this increase is that a binary in a PN gravity field requires slightly more energy and angular momentum to support a given orbital separation. This extra energy must be radiated away and therefore tends to increase the overall inspiral time-scale. To get a feel for how different orders affect the overall behaviour, we show the evolution of a highly eccentric (e0 = 0.999), moderately eccentric (e0 = 0.9), and low-eccentricity (e0 = 0.3) orbit in Figure 1. In all cases, Peters' formula underestimates the inspiral time-scale, especially so for eccentric orbits. The effect of the 1 PN correction is also to further increase the time-scale, whereas higher-order corrections cluster around the 1 PN result.

Correction ratios
In order to compute the inspiral time-scale, we performed numerous numerical integrations of the orbital evolution, sampling a wide range of initial orbital separations, eccentricities, and spin parameter values. We define the merger event as the moment at which the periapsis of the orbit, p = a(1 − e), reaches the effective Schwarzschild radius. We express the corrections to Peters' time-scale as ratios R and Qi, where the index i stands for the several different effects that we want to describe. The most complete and accurate formula for the time-scale is obtained by multiplying Peters' formula tP by three separate correction ratios, where the ratio R(e0) accounts for the secular eccentricity evolution at quadrupole order, Q h (a0, e0) accounts for the 1 PN and 1.5 PN spinless effects, and Qs(a0, e0, s1) accounts for the inclusion of the spin-orbit coupling. We calculate these ratios numerically by integrating the appropriate evolution equations, and subsequently find various fits for the  (1) values for the spin parameter s 1 , respectively. In the right-hand panel, we show the relative error of the fits we propose for the same range of parameter space. For the sake of clarity, we avoid plotting the curves with spin, but note that the fits perform similarly to the spinless case (with the single exception of a circular, equatorial and prograde orbit at 6r S around a maximally spinning primary BH, where the fit fails by no more than 30 per cent). For a periapsis p 0 ∼ > 10r S , the errors lie well below 1 per cent. Even for the extreme case of p 0 = 6r S the fits perform rather well, considering the fact that the PN series itself is known to break down below those separations.
results. The numerical results can be seen in the left-hand panel of Figure 2, where we plot the ratio of the 1.5 PN timescale to the Newtonian time-scale for different values of the initial periapsis, spin, and initial eccentricity. As shown in Paper I, the correction generally grows for smaller periapsis, while the two edge cases of co-and counter-rotating orbits around an extremal BH envelop the spin-less results. In the circular case, the spin parameter s1 is degenerate with the periapsis of the orbit by ∼ > 5rS, while it has much less of an influence in the eccentric case. Our proposal for the appropriate fits to these ratios reads , Qs(p0, e0, s1) = exp As(p0, e0)s1 + Bs(p0, e0)|s1| 3/2 , (24) As(p0, e0) = e0 0.3rS p0 + (1 − e0) 3/2 3.7rS p0 3/2 , Bs(p0, e0) = e0 1.1rS p0 valid for arbitrary values of the initial semi-major axis, eccentricity, and (primary BH) spin parameter. As shown in Figure 2, we find these fits to be an accurate estimate of the numerical results for orbits whose initial periapsis is no lower than p0 ∼ 3rs(2 − e0). Below that point, the assumptions of the PN expansion break down, and we cannot expect our results to be accurate.
These general correction factors turn out to be rather complicated and unwieldy. The reason for this is the complicated behaviour of the time-scale in a region of parameter space between 0.1 ∼ < e0 ∼ < 0.8, visible in Figure 2. The likely cause for this are some terms in the higher-PN fluxes switching sign at eccentricities of ∼0.5, complicating the dynamics for orbits that start evolving from that region of phase space. With regards to LISA sources, however, most formation channels do not predict an isotropic distribution of eccentricities, but rather a strong preference for either highly eccentric or almost circular orbits. Let us take EMRIs as an example: on the one hand, dynamical processes tend to generate EMRIs with e0 0.9 by scattering compact objects onto almost radial orbits (see, e.g. Amaro-Seoane et al. 2007a, 2013; on the other hand, some EMRIs can be produced when a binary containing a compact object is disrupted by the tidal influence of the SMBH (see, e.g. Amaro-Seoane 2019). These would likely enter the GW-dominated phase with very low eccentricities. For this reason, the most general form for the fits of the correction ratios is often not required. Rather, one can freely choose an appropriate eccentricity range or limit that is expected for a particular formation channel. Here we will discuss the two most obvious limits of e0 → 0 and e0 → 1. In these limits, the fits reduce to much more manageable formulas that satisfy our goal to be as simple as possible.

Quasi-circular orbits
If the binary enters the GW-dominated phase with a very low eccentricity (e0 ∼ < 0.1), the fits for the correction ratios simplify to the following formulas: In this case, the eccentricity correction factor does not have a significant effect, but the PN correction factors take their maximal values for the coefficients. This occurs because binaries on circular orbits spend the whole duration of the orbit at small separations, maximising the influence of higher-order PN effects. In this limit, the fits are accurate to ∼10 per cent of the relative value if p0 ∼ > 6rS (see Figure 2).

Highly eccentric orbits
If the binary enters the GW dominated phase with a very high eccentricity (e0 ∼ > 0.8), the fits for the correction ratios simplify to the following formulas: In this case, the eccentricity correction factor dominates, whereas the PN correction factors are slightly suppressed with respect to the circular case. Nonetheless, both factors can compound on each other to reach values of ∼10 for large regions of initial parameter space. In this limit, the fits are accurate to ∼10 per cent of the relative value if p0 ∼ > 3rS, and to ∼1 per cent if p0 ∼ > 6rS (see Figure 2).

Application: EMRI critical semi-major axis
The EMRI critical semi-major axis, acrit, marks the separation between the dynamical-evolution regime and the GWemission regime, and is an essential ingredient for the calculation of EMRI event rates. The value of acrit is found by equating the inspiral time-scale to the dynamical relaxation time-scale for a periapsis equal to the primary's last stable parabolic orbit. In Figure 3, we plot acrit as derived in Amaro-Seoane (2019), this time including the correction factors to the GW-inspiral time-scale, for a system of a 10 M BH inspiralling into a 4 × 10 6 M SMBH. For eccentric orbits, we find that the value of the critical semi-major axis decreases by an order of magnitude from previous results that use Peters' formula as a simple model of the inspiral time-scale. In the circular orbit case, the value of acrit is generally decreased by a factor of a few, but the effects of the spin correction factor Qs are clearly visible in the steeper scaling of the counter-rotating case (S1 < 0), and in the slightly flattened scaling for the co-rotating case (S1 > 0). The difference between the circular and the eccentric cases is also increased by almost an order of magnitude, which suggests different relative efficiencies for formation channels that preferentially create circular versus eccentric EMRIs. Note that the ISCO of a BH in a prograde orbit can be as low as ∼ > 0.5rS, which is lower than the value where our fits (and PN theory in general) are accurate. In these cases, we evaluate the correction factors at 3rS for highly eccentric orbits and at 6rS for circular orbits, to ensure that the results do not diverge. This would represent a conservative estimate of the real change of the critical semi-major axis.
The expected event rates for EMRIs are then computed by performing an integral over phase space, in which acrit is one upper limit (Amaro-Seoane et al. 2013; Amaro-Seoane 2019). The total phase-space volume of integration scales as a 3 crit , which means that even moderate changes to acrit can potentially 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) event rates that include the corrections to the inspiral timescale. Preliminary results show that the rates decrease by more than one order of magnitude when switching from Peters' time-scale to the corrected formula. However, note that some event-rate estimates compute the critical semi-major axis using a version of the merger time-scale that focuses specifically on highly eccentric orbits (Amaro-Seoane 2019). In that case, our correction to the critical semi-major axis is only due to PN effects and the rates generally vary by a factor of a few.

Perturbations by a third body
In this section, we consider the case in which a distant massive body perturbs the orbit of an inspiralling LISA source.
In other words, we consider a hierarchical three-body problem in which the hard binary is made of compact objects at very small orbital separations, and can exhibit relativistic behaviour. Following Will (2014), we express the energy of such a system in the following way: where EN and E1 PN are the Newtonian and 1 PN energy contributions, whereas EP is the change in energy of the hard binary of a given osculating semi-major axis and eccentricity due to the influence of the distant third body. The explicit expression for EP, in the case of a perturber in a circular orbit, reads where the function G describes the angle dependence of the energy shift, µ = m1m2/M , m3 is the mass of the perturber, R its distance to the centre of mass of the hard binary, I0 the inclination with respect to the orbital plane of the perturber, and ω0 the angle of the periapsis with respect to the ascending or descending nodes of the orbital plane of the perturber. The remarkable fact proven in  is that this expression remains constant over time-scales at which perihelion precession can significantly change the orientation of the orbit with respect to the plane of the perturber. The terms that would normally vary in ω are cancelled by cross terms between the PN dynamics of the inner binary and the tidal forces induced by the perturber. Indeed, the second sin 2 term is evaluated at an initial perihelion orientation, and therefore is not affected by perihelion precession. This result is very convenient for our purposes, because it is essentially of the same form as the change in initial conditions caused by a PN perturbation. We can estimate how the change in energy of the hard binary caused by the perturber will affect the inspiral timescale by readjusting the initial conditions of the GW-induced decay. We define a fictitious semi-major axis aP and eccentricity eP, constructed from both the Newtonian energy of the hard binary, EN, and the constant shift induced by the perturber, EP: An estimate for the change in the inspiral time-scale can by found by evaluating a particular vacuum time-scale tGW at the fictitious semi-major axis and dividing it by the corresponding vacuum time-scale of the physical orbit: where, in the first step, we assumed quadrupole radiation or, equivalently, the evolution equations of Peters & Mathews (1963). The function F can be computed from the eccentricity dependence of the inspiral time-scale, and ranges from 1/5 to 1/2 for most values of the eccentricity. From the results of Paper I and previous sections, we know that PN corrections to the Newtonian vacuum time-scale can be described by factors of the form exp CPNrSp −1 where CPN is a coefficient of order unity and n is the PN order. We can therefore estimate the radius at which the environmental perturbation becomes of Newtonian (or PN) order by solving the equation and find the characteristic radii R P/PN at which the effect of the perturbation becomes as large as that of a given PN order: where we approximated F(e) 1/3 ≈ 1 and C 1/3 PN ≈ 1. For the Newtonian order (n = 0 and CPN = 1), the radius reads which can only be of order of or smaller than the semi-major axis of the hard binary itself, breaking the assumption of a hierarchical triplet. In general, every additional PN order only adds a factor of a few over the previous value and, since we are considering hierarchical triplets (where R a0), we can conclude that the perturbation has very little influence on the inspiral time-scale.  (2015) Fujita (2015) Munna (2020) m 3 = 10 1 M m 3 = 10 3 M m 3 = 10 6 M Figure 4. We show the minimum required PN order of a waveform template needed to detect a shift in the inspiral time-scale of an EMRI with an orbital frequency in the LISA band (f ≈ 1 mHz) caused by the perturbation of a third body with mass m 3 , in the case of an inner binary composed of an SMBH of 10 6 M and a compact object of 10 M . We calculate the value by solving Equation (35) for the PN order n and adding one. The solid lines represent a quasi-circular inspiral, whereas the dotted lines represent orbits that have the same frequency but an eccentricity of 0.9. The horizontal lines stand for the state-of-the-art results in the calculations of PN fluxes using different methods, for spinning BHs in eccentric (4 PN) and circular (11 PN) orbits (Fujita 2015;Sago & Fujita 2015), and also eccentric orbits for non-spinning BHs (Munna 2020).
On the other hand, if the accumulated signal-to-noise ratio (SNR) of the source is high enough, LISA might be able to detect deviations from very high PN orders. In Figure 4, we show the minimum order of a PN-waveform template that is required to detect a shift in the time-scale for a perturber of a given mass and separation, in the case of an inner binary composed of an SMBH of 10 6 M and a compact object of 10 M . We compute it by solving Equation (35) for the PN order n and adding one. If we assume that the SNR is sufficient to reach current state-of-the-art waveform templates with spin (depicted in Figure 4; Fujita 2015), we find values that are similar to the results in Yunes et al. (2011), wherein it is shown that detecting the dephasing in the signal of such a system is in principle possible, but only realistic for very massive (∼10 6 M ) and close (sub-parsec) perturbers.
The shift in energy is not the only way in which a third body can perturb the hard binary. For orbits that are sufficiently misaligned, Kozai-Lidov (KL;Kozai 1962;Lidov 1962) oscillations can induce a secular variation in the eccentricity of the orbit. If the inspiral time-scale is longer than the time-scale of KL oscillations, tKL, the evolution equations for the orbital parameters of the hard binary can be significantly affected by the perturber. The detectability of KL oscillations in the GW signal of LISA sources has been very recently studied in Randall & Xianyu (2019) and Deme et al. (2020), where the original version of Peters' timescale is used as a simple model for the inspiral time-scale. An intermediate calculation in the aforementioned papers is the computation of the characteristic orbital separation, a KLO/GW , of the inner binary at which KL oscillations are quenched by GW emission. The radius is found by equating the GW and KL time-scales and solving for the semi-major axis: tGW ∼ R(e)Q h (a, e)Qs(a, e, s1)tP = tKL = 2π =⇒ a KLO/GW ≈ 2 2 7 πM qR 3 f (e)(GM ) 5/2 5c 5 m3R(e)Q h (a, e)Qs(a, e, s1) where we model the corrections to Peters' time-scale with the correction factors R(e)Q h (a, e)Qs(a, e, s1). The general solution a KLO/GW must be computed numerically. However, we can see that (while suppressed by the small power) the correction factors can reduce the value of the quenching radius a KLO/GW by a factor of up to ∼1.5, which might be enough to significantly affect any further results that depend on powers of a KLO/GW . An example would be the phasespace volume a 3 KLO/GW , which can easily change by a factor of 2-3.

Gas drag and torques
Gas drag and torques in accretion discs are expected to significantly influence the dynamics of an inspiral event and might cause a detectable dephasing of the GW signal depending on the physical properties of the disc (see, e.g. Barausse et al. 2014;Derdzinski et al. 2019). While several different disc models exist, many show a central region that extends for ∼5 × 10 3 Schwarzschild radii, where the density and temperature profiles are much more shallow than in 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 is a reasonable first-order approximation to describe the local density and temperature as constant. To describe the influence of the disc on the inspiral time-scale, we need to model the effects that have a direct influence on the orbital elements of the secondary BH.
One such effect is gas drag. While the subject of dynamical friction is complex, we will assume that, for BHs, gas drag 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 BH will be reduced by an amount of work WBHL given by where S is the path of the orbit in space. The magnitude of the BHL drag reads where cs is the local speed of sound, ρ is the local disc density, m is the mass of the inspiralling BH, and v rel the relative velocity with respect to the Keplerian disc velocity. If the orbit is very eccentric, the instantaneous velocity, vi, of the secondary will generally be misaligned with the Keplerian disc velocity. In this case we can approximate v rel ≈ vi and vi cs, which reduces the BHL drag formula to the standard high-velocity limit of dynamical friction (see, e.g. Mayer 2013, for a discussion of gas effects in highly eccentric orbits). Note that this might also be the more appropriate description, at least in a local sense, in discs that are very turbulent, where gas velocity deviates from a Keplerian disc (but only if the sound speed is also smaller than the BH velocity). We will consider both the BHL case and the highvelocity dynamical friction case, keeping in mind that the former is more accurate for circular orbits and the latter is more accurate for eccentric ones. The integral can only be solved numerically or as a series expansion in the eccentricity. We compute it for the two cases: where WBHL is a function that depends on the eccentricity of the orbit and an effective average Mach number nM = v 2 d c −2 s evaluated at a radius of one semi-major axis, where v d is the Keplerian disc velocity. It varies between one (for circular orbits) and several thousands, depending on the eccentricity, and converges very poorly for high Mach numbers. The function W dyn only depends on eccentricity and can be computed as a polynomial expression in e 2 : W dyn (e) = 1 + 3 4 e 2 + 21 64 e 4 + 55 256 The loss of energy due to drag results in a secular change of the semi-major axis of the secondary BH, which we can compare to similar effects caused by the PN fluxes. To obtain an expression for the orbit-averaged energy flux due to BHL drag and dynamical friction, we simply divide Equation (41) by an orbital period T .
Keeping this result in mind, we turn our attention to another effect that directly influences the orbit of the secondary BH, namely global torques. A realistic description of torques in accretion discs is only possible through hydrodynamical simulations. Nevertheless, we will assume that the simplest analytical models of Type I and Type II viscous torques are an appropriate order-of-magnitude approximation, when averaged over many orbits. We take the formulas by Tanaka et al. (2002) and Derdzinski et al. (2021), which readL where Ω2 is the Keplerian frequency of the secondary and α is the viscous parameter, and we assume that the disc scale height, h, follows the standard scaling h ∼ csΩ −1 , where Ω is the Keplerian disc orbital frequency. We average the expressions over one orbit, yielding where Di(e) are functions of the eccentricity. The function DI(e) is approximated very well by the first four orders of its polynomial expansion, whereas DII(e) has a closed form solution: DI ( In the case of torques, the orbit-averaging procedure is not completely legitimate: Type I torques are a consequence of resonances between the orbital frequency of the secondary and the disc, and as such it is unlikely that they would realistically appear for eccentric orbits. Type II torques become relevant when the secondary is massive enough to open gaps in the disc structure (see, e.g. Goldreich & Tremaine 1980;Ward 1997). In a tight binary of SMBHs, the gap opening torques result into the formation of a cavity, and of a circumbinary disc around it (e.g. Cuadra et al. 2006;Roedig et al. 2012;Tang et al. 2017). These torques depend on the detailed thermodynamics and viscous properties of the disc because both pressure forces and viscous forces oppose the gap-opening torque induced by the binary into the surrounding disc gas. Therefore, one would in principle need to account for the dependence of the torques on the many disc properties that determine the strength of the forces at play. However, here we take an explorative stand as we only want to highlight order-of-magnitude effects, which could later be explored in a self-consistent model for the disc. In general, Type II torques will apply only to SMBH binaries as opposed to IMRIs/EMRIs, because there must be a massive enough secondary to open a gap. The recent sub-pc scale hydrodynamical simulations of Souza Lima et al. (2020) show that even with a mass ratio of 20:1 the opening of a gap or cavity, is not always occurring, whereas Derdzinski et al. (2021) show that, for IMRIs, the orbital evolution proceeds with no gap opening. Additionally, we will also see that the functions Di do not play a significant role for the results of this section.
Equations (47) and (48) are in the form of an orbitaveraged angular momentum flux, which we can compare directly to the PN fluxes. 3 Note that, in the case of circular orbits, one can obtain results that are consistent with any particular disc model by simply replacing ρ and cs with the corresponding profiles. For the eccentric case, the orbit averages can be computed numerically for arbitrary disc models.
In order to compare a large range of PN orders to the results shown above, we have to find a simple way to characterise their scaling in the semi-major axis and eccentricity of the secondary's orbit. By checking up to the first four integer-PN orders of both energy and angular momentum fluxes, we always find the following structure: where the subscript n denotes the (integer) PN order and the subscript N denotes the lowest-order (quadrupolar) fluxes, whereas the coefficients Ci and Di are fractions of polynomials in e 2 of order unity. The half-integer terms have slightly different forms, but we will neglect them fur the purpose of the following calculations.
To find a series of characteristic orbital separations, at which the effect of the gas on the evolution of the orbit is as strong as the n-th order PN fluxes, we equate the expressions and solve for the semi-major axis. Note that this is equivalent to equating the GW-induced inspiral time-scale to a characteristic inspiral time-scale due to the effect of the gas. In the case of BHL drag, we can solve the equation analytically only for the case of circular orbits, because in general the function WBHL contains an arbitrary amount of powers of the semi-major axis a. For circular orbits, the characteristic semi-major axis a BHL/PN reads where we are able to neglect all the small numerical coefficients because of the small exponent. For dynamical friction, the radii a dyn/PN read where, in the last step, we approximated [f (e)/W dyn (e)] 2/11 1. We will use the BHL results as the low-eccentricity limit and the dynamical friction results as the high-eccentricity limit.
We repeat the same calculation for the case of Type I torques and obtain similar characteristic separations a TI/PN as before. In this case, we are able to solve it analytically for arbitrary eccentricities, with the caveat that global torques are unlikely to appear for non-circular orbits: Here we have a slightly different scaling in the physical properties of the disc, but otherwise the formulas look very similar. Equation (57) With the aid of these characteristic radii, we can identify different regions of phase space in which we can expect certain PN orders to be accurate descriptions of the inspiral time-scale. The Newtonian (n = 0) result is particularly telling, since it separates the regions where gas effects dominate from regions where GW emission dominates. Figure 5 shows this procedure for the standard case of a 10 6 M SMBH and a 10 M compact object. The results only scale very weakly with disc properties such as central density and temperature, suggesting that they are robust as an order-of-magnitude estimate. For example, one can use the formulas for the characteristic radii to estimate deviations from the vacuum solutions as a function of the disc properties. As an example, let us take a 10 M BH on a circular orbit at a radius of 10 3 rS. If it is embedded in a disc with a local density of ρ = 10 −13 g cm −3 and viscosity parameter α = 10 −3 (Thompson 2009), we expect the effect of BHL drag (if present) to be almost as important as the lowest-order quadrupole radiation, while Type I or Type II torques (if present) would only be as important as a 1 PN correction.
Vice versa, one can also estimate the disc properties required for a shift to be visible at a certain PN order, which we show in Figure 6 for varying values of the local disc density. One can also use this type of analysis to speculate on the likely environment of actual LISA sources. As an example, if we assume that the density of a typical circumnuclear disc (CND) scales as a simple r −2 power law (Mestel 1963) until it flattens out at a radius of ∼5 × 10 3 rS, we can extrapolate some realistic densities from observational results for the total mass and extension of CNDs given in t G W < 1 0 2 y rs t G W < 1 0 4 y rs t G W < 1 0 8 y rs a TI/PN a TII/PN a BHL/PN , n M 1 a dyn/PN ISCO Figure 5. Shown is a region of phase space in which the characteristic 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 same colour separate regions of phase space in which the effect of gas is stronger than the 0th, first, second, and third PN order effects. In other words, one can expect this particular GW source to decay in a time that is well described by n-th order PN theory only if it starts from a location below the n-th set of coloured lines in phase space. The results shown are for an SMBH of 10 6 M and a secondary of 10 M embedded in an accretion disc with a central mean density of ρ = 10 −13 g cm −3 and a viscous coefficient α of 10 −3 , consistent with the models in Thompson (2009) . Changing the density and viscosity parameter of the central part of the disc by orders of magnitude only slightly shifts the coloured curves up and down. The black line marks the location of the ISCO for a non-rotating central object, and the dotted lines denote the isochrone curves for different GW decay times. Medling et al. (2014), obtaining ρ ∼ 10 −11 to 10 −9 g cm −3 (a result more consistent with the disc models in Sirko & Goodman 2003). Note that these assumptions are also used as initial conditions in state-of-the-art simulations of SMBH coalescence (see, e.g. Souza Lima et al. 2017Lima et al. , 2020. Even for the wide range of observed CND masses (∼10 8 M to ∼5×10 9 M ), the required PN order to detect a shift in the inspiral time-scale for an EMRI that is entering the LISA band (at ∼1 mHz) only ranges between 3 and 5, which is realistic given a source with sufficient SNR. These results also qualitatively agree with more sophisticated hydrodynamic simulations such as those in Derdzinski et al. (2019) and Derdzinski et al. (2021), where it is shown that, if surface densities of ∼ > 10 3 g cm −2 are present, a dephasing of the GWs is likely to be detectable by LISA. Note that this type of calculation can be repeated with more sophisticated density profiles to yield more precise results. However, even very simple estimates such as assuming a r −2 density profile until 5 × 10 3 rS can yield correct order-of-magnitude estimates because of the weak scaling of the characteristic radii with disc properties. . We show the required PN order needed to detect a shift from the vacuum time-scale in the cases of BHL drag, Type I and Type II torques for different values of the local density. The lines are computed for the standard case of a 10 6 M SMBH and a 10 M secondary entering the LISA band at ∼1 mHz (solid coloured lines) with a local sound speed of ∼ 3 × 10 7 cm s −1 . The vertical lines represent a qualitative extrapolation of the local density from observed CND masses and extensions. The results only scale very weakly with parameters such as the CND mass, local sound speed, or the mass ratio of the binary, but change significantly with the SMBH mass. This is, however, not due to the properties of the disc, but rather to the fact that the efficiency of GW emission at a fixed frequency strongly depends on the mass of the SMBH.

SUMMARY AND CONCLUSION
In this paper, we extended our analysis of the GW-induced decay time-scale to include the 1.5 PN order. We produced a series of analytical fits that can be multiplied by Peters' formula 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 to simple exponential functions that can be used to avoid errors of a factor ∼10 and ∼2, respectively. We show that these modifications to Peters' formula are necessary in applications 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 also Vázquez-Aceves et al., in prep.). We then used the results of  in combination with our findings to compare the strength of arbitrarily high PN perturbations against gravitational perturbations in a hierarchical triple system. We found a series of characteristic radii at which the perturbations affect the inspiral timescale as much as any particular PN order (Equation 35), which we used to reproduce previously known results on the detectability of the dephasing of GW signals . Furthermore, we briefly discussed the effects of KL oscillation and how corrections to Peters' time-scale can affect recent calculations found in Randall & Xianyu (2019) and Deme et al. (2020) that depend on the typical radius at which KL oscillations are quenched by GW emission.
Finally, we discussed the validity of vacuum time-scales for inspirals embedded in a gaseous medium, which is the likely case for SMBHs sourrounded by accretion discs and CNDs. We analysed several types of drag forces (BHL drag and dynamical drag) and global torques (Type I torques and Type II viscous torques), and compared them against PN energy and angular momentum fluxes of various orders. We found several characteristic radii that separate the phase space in different dynamical regions, where gas effects are as strong as a particular PN order (Equations 53 to 59). We used a simple prescription to extrapolate the gas densities that would be likely encountered by a GW source entering the LISA band in one of the several observed CNDs (Medling et al. 2014), and found that PN orders of 3-5 are required to detect a deviation from the vacuum case, which is compatible with results from hydrodynamical simulations (Derdzinski et al. 2019(Derdzinski et al. , 2021. While many assumptions are required to derive the analytical formulas, the results scale only very weakly with disc properties, suggesting that they might be a robust order-of-magnitude approximation of a more realistic treatment of gas-embedded inspirals.