Improved gravitational radiation time-scales: significance for LISA and LIGO-Virgo sources

We compute a revised version of Peters' (1964) time-scale for the gravitational-wave (GW) induced decay of two point masses, by taking into account post-Newtonian (PN) perturbations of the orbital motion. At first PN order, the corrected time-scale can be approximated by multiplying Peters' estimate by the simple factor $Q = 1 + 5 (r_{\rm S}/p)$, where $p$ is the periapsis and $r_{\rm S}$ the Schwarzschild radius of the system. We apply the revised time-scale to a set of typical LIGO-Virgo and LISA sources at the onset of their GW driven decay. We argue that our more accurate model for the orbital evolution will affect current event and detection rate estimates for mergers of compact object binaries, with stronger deviations for LISA sources (EMRIs, IMRIs, and SMBH binaries). We propose the correction factor $Q$ as a simple analytical prescription to quantify decay time-scales more accurately in future population synthesis models.


INTRODUCTION
The era of gravitational-wave (GW) astronomy started with the detection of merging stellar-mass black holes (BHs) and inspiralling neutron stars by the ground-based Laser Interferometer GW Observatory (LIGO) and Virgo detectors (Abbott et al. 2016(Abbott et al. , 2017(Abbott et al. , 2019. These two detectors, along with the upcoming Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017; Barack et al. 2019), not only allow physicists to observe the Universe, but offer an incredible tool to actively test the current theory of gravity: general relativity (GR). Making a prediction for a gravitational signal requires precise modelling of the physics of the source. Since the first successful numerical integration of the evolution of binary BH spacetime (Pretorius 2005), we have been in the position of simulating relativistic mergers in GR. The required integration is, however, extremely demanding for current computational resources, which makes covering a substantial part of the source's parameter space E-mail: zwicklo@ics.uzh.ch by means of fully relativistic simulations impractical. Moreover, different types of GW sources interact very differently with their astrophysical surroundings (see, e.g. Barausse et al. 2014). These interactions must be modelled in order to produce realistic event and detection rate estimates. Fortunately, there have been major developments in the field of approximations to GR, which allow to model both the sources and their surroundingsmore efficiently than a fully relativistic simulation. Schemes such as the post-Newtonian (PN) expansion (see, e.g. Blanchet 2014) or the effective one body approach (see, e.g. Damour & Nagar 2016) are used to find approximate solutions of Einstein's field equations and the associated equations of motion for the compact object binary. They successfully describe the relevant features of fully relativistic orbits, such as perihelion precession and radiation of GWs, and can be used to cover a larger amount of parameter space at the loss of some precision. The latter process is especially relevant, as it determines how relativistic orbits shrink as they lose energy due to gravitational radiation.
The first successful quantitative description of the evo-lution of a binary's orbit by means of GW induced decay is due to Peters (1964), in which (among other things) the quadrupole formula is used to calculate how the semimajor axis and the eccentricity of a Keplerian ellipse evolve in time. The appeal of the resulting equation, often called "Peters' time-scale", lies in its simplicity. It can be used in a wide range of applications to model the GW induced decay of any orbit. The crucial assumption in Peters (1964) often goes overlooked, even though it is clearly stated in the title of the paper itself. The assumption is that the orbital motion of binaries is described by Keplerian orbits even at scales where gravitational radiation is relevant.
At the present time, however, gravitational detectors such as LIGO/Virgo (and, in the future, LISA) are probing orbits that definitely require at least a PN analysis of orbital dynamics. In this paper, we argue that not considering deviations from Keplerian orbits is actually inconsistent (though often a good enough approximation) from a PN point of view when modelling GW induced decay. We then present an analytical modification to the commonly used result that takes into account the largest deviation of a PN orbit from a purely Keplerian ellipse. We claim that it will provide a more accurate time-scale of gravitational decay while retaining the appealing simplicity of Peters' result. In Section 2, we present the theoretical background of the problem, focussing on the PN series and the parametrization of relativistic orbits. In Section 3, we derive and discuss an analytical correction to Peters' time-scale. In Section 4, we discuss our findings and give a few examples for the most interesting sources of GWs, i.e. binaries of supermassive BHs (SMBHs), binaries of stellar-mass BHs, and extreme and intermediate mass-ratio inspirals (EMRIs and IMRIs). We speculate on the effects of the correction on event-rate estimates for GW detectors such as LIGO-Virgo and LISA.

THEORETICAL BACKGROUND
In order to analytically describe the motion of two bodies in GR, it is necessary to approximate Einstein's field equations. An important and powerful approximation scheme is the Arnowitt-Deser-Misner formalism (hereafter ADM; Arnowitt et al. 1959; see also Arnowitt et al. 2008). It offers a way to describe the motion of a binary system through the familiar concept of a Hamiltonian and Hamilton's equations for the separation and momentum vectors of the bodies. The ADM Hamiltonian, HADM, is organised into a PN series. Here the PN order is denoted with the index Hi and HN denotes the Newtonian contribution: where c is the speed of light in vacuum. Different gravitational effects (such as precession, spin-orbit coupling, etc.; see Schäfer & Jaranowski 2018 and references therein) enter the picture at their respective PN orders. Whereas all orders depend on the canonical variables, the 2.5 and all subsequent half-integer orders also explicitly depend on time. This is because half-integer orders describe energy dissipation due to GWs. The 2.5-term of the ADM Hamiltonian is determined by the quadrupole moment tensor of a Newtonian ellipse N: where G is the gravitational constant. The explicit time dependence in the half-integer terms of the ADM Hamiltonian precludes the existence of a general analytic solution to Hamilton's equations. The integer-order part of the ADM Hamiltonian is often called "conservative part" because it admits the definition of a conserved reduced (i.e. per unit mass) energy E and, in certain cases (e.g. when the spin is zero), even a conserved reduced angular momentum h (Schäfer & Jaranowski 2018). 1 These conserved quantities have been used in Damour & Deruelle (1985), Wex (1995), and finally in Memmesheimer et al. (2004) to produce parametric solutions to the PN Hamilton's equations that are exact up to third order for binaries with no spin (for another approach see, e.g. Boetzel et al. 2017). The solutions describe elliptical orbits undergoing periapsis precession; a simple way to describe them is to use the familiar concepts of semimajor axis, a, and eccentricity, e. However, some care is needed because the simple Newtonian relations between energy, angular momentum, and the orbital parameters are no longer valid in the PN model. Rather, the definitions of the orbital parameters are now given as PN series: where M = m1 + m2 (with m1 ≥ m2) is the combined mass of the two bodies. These orbital parameters (along with the periapsis precession) can be used to describe the thirdorder PN behaviour of the orbit in absence of the dissipative gravitational radiation term (2.5 PN). There is, however, a way to reintroduce dissipation effects in the description of the orbit. To understand it, we briefly return to the case of a Newtonian orbit.
The energy E and angular momentum h are conserved quantities of a Newtonian binary orbit. They completely determine its shape and size by means of the formulae On the other hand, Einstein's quadrupole formula (Einstein 1916) predicts that the (non-reduced) energy E and angular momentum L will be radiated from a gravitational system in the form of GWs: Improved gravitational radiation time-scales 3 Here the matrix M contains the components of the mass quadrupole moment tensor of the binary for an angular momentum pointing in the z direction. There is a clear inconsistency between the conserved quantities of a purely Newtonian description and the dissipation equations of the quadrupole formula. This inconsistency can be alleviated by noting that changes in E and L occur on time-scales much longer than the typical orbital (radial) period. Indeed, the time-scale over which the (non-reduced) energy changes scales as dt ∝ c 5 dE, whereas the typical orbital period is independent of c. This consideration is what allowed Peters (1964) to combine Keplerian celestial mechanics with the quadrupole formula to obtain the equations for the GW induced evolution of orbital parameters: where f (e) = 1 + 73 24 e 2 + 37 96 e 4 (1 − e 2 ) −7/2 (11) and q = m2/m1 ≤ 1 is the mass ratio. These equations describe how a Newtonian orbit slowly changes shape if its energy and angular momentum are radiated away according to the quadrupole formula. They predict that orbits decay over a particular time-scale, often referred to as Peters' time-scale (hereafter tP). By manipulating Equations (9) and (10), it is possible to find an analytic expression for tP. It scales with the fourth power of the initial Keplerian semimajor axis, a0, weighted by the eccentricity enhancement function f (e0), e0 being the initial Keplerian eccentricity: The assumption used in Peters (1964) is that, despite the radiation of energy, the binary is still moving along the orbit predicted by Newtonian mechanics. Here we argue that this assumption is not necessary. Note how the results of the quadrupole formula (Equation 7) in the context of Peters' work and the 2.5 PN Hamiltonian term (Equation 2) in the context of the ADM formalism are identical. Indeed, there are no additional terms in the ADM Hamiltonian that describe energy radiation before the 3.5 PN order. This suggests that the quadrupole formula should be able to describe the orbital evolution not only of a Newtonian orbit, but also of a PN orbit up to third order.
To prove this statement, one can attempt to "improve" the results of the quadrupole formula. We start with a simple Newtonian orbit and insert it in Equation (7). We then obtain an evolution equation for the energy E that scales with c −5 . As noted before, this is already equivalent to the 2.5 PN ADM Hamiltonian. From here on, there are only two possibilities to improve the accuracy of the description of energy radiation via GWs: the contribution of a 1 PN perturbed orbit to quadrupolar radiation and the contribution of a Newtonian orbit to octupolar radiation. To improve it even further, it is necessary to add more and more contributions (see, e.g. Thorne 1980;Galley & Leibovich 2012, for more details on this topic). We show the process in a schematic form: where we denote the operation of calculating the 2 n pole contribution of an orbit with the functions "2 n pole[ ]". The functions scale with c −1−2n . We also denote the various perturbations to the orbital dynamics as "Newton, 1 PN, 2 PN, ..., n PN". The perturbations scale as c −2n . Note how any improvement to the results of the quadrupole formula will necessarily have at least a c −7 factor. Indeed, any improvement to the quadrupole formula is of 3.5 PN order. The situation is very different if we consider the assumption of Keplerian orbits: perturbations already arise at the first PN order. Another way to state this is the following: the assumption of Keplerian orbits fails at a lower PN order than the assumption of quadrupolar energy radiation. The idea of this paper is thus to apply the quadrupole formula to the energy of a PN orbit rather than a Keplerian orbit. The parametric solution to the conservative ADM Hamiltonian equations explicitly depends on some kind of conserved energy and angular momentum, analogous but not necessarily equal to their Newtonian counterparts. By taking these quantities and using them as initial conditions for the energy dissipation equations, we can describe the full 3 PN evolution of the "conserved quantities" of an orbit, without needing any higher-order (e.g. 3.5 PN, 4.5 PN) energy dissipation contributions. This process should improve the results derived in Peters (1964). Before proceeding with the calculations, we sum up our arguments. In Peters (1964), the Newtonian definitions (Equations 5 and 6) are used to relate the evolution of the energy and angular momentum to the evolution of the shape of a binary's orbit. This is done without considering the fact that those definitions are only valid at the lowest (zeroth) PN order. In the PN case, they are given by PN expressions that depend on the PN conserved quantities (Equations 3 and 4). The use of the Newtonian definitions is actually inconsistent from a PN point of view, since it combines c 0 information on the orbit's geometry with c −5 information on the evolution of the orbit's parameters. We propose to adopt the more consistent PN definition of semimajor axis and eccentricity instead. This will necessarily change the predicted duration of the orbit's decay, leading to a revised self-consistent version of Peters' time-scale. In Section 3 and in Appendix A, we compute how much these considerations alter the original prediction of Peters (1964) for the duration of a binary decay. In the former, we calculate the effect of the inclusion of first-order PN dynamical information analytically. In the latter, we present a simple scheme on how to derive the fully consistent 3 PN results by means of a few numerical operations.

ANALYTIC DERIVATION OF THE LOWEST ORDER CORRECTION
The goal of this section is to derive an explicit analytic correction to Peters' time-scale in the PN framework. In order for the correction to be analytic, we have to restrict our calculations to the lowest significant order, effectively combining the c −2 conservative information with the c −5 nonconservative information. Readers interested in the complete derivation of self-consistent corrections may refer to Appendix A.
In their original form, Equations (9)-(10) simultaneously describe the change of the geometry and that of the energy state of a binary system. However, when comparing a Newtonian orbit with a perturbed one, some caution is advised: the shape (semimajor axis, eccentricity) of a Newtonian and a PN orbit with the same energy and angular momentum is not identical. Analogously, if the shape of the orbit is known at a moment in time, then its energy and angular momentum will differ from the ones computed in a purely Newtonian approximation. Here we treat these two cases separately: in Section 3.1, we start from a Newtonian and a PN orbit with identical initial energy and angular momentum but different orbital shape, whereas the case of identical initial orbital shape but different initial energy and angular momentum is presented in Section 3.2.

Orbits with identical initial energy and angular momentum
We start by recalling the explicit first-order expressions for the semimajor axis, a1, and eccentricity, e1, of the perturbed orbit (that describe a more realistic orbital shape), in terms of the PN conserved energy E1 and angular momentum h1 (Memmesheimer et al. 2004): e 2 1 = 1 + 2E1h 2 1 + E1 5E1h 2 1 3q 2 + 5q + 3 + 2 6q 2 + 11q + 6 Since the underlying assumption of this section is that the initial energy and angular momentum in the Newtonian and PN framework are the same, we can replace E1 and h1 with their classical definitions (EN and hN): A useful product of this manipulation is the formula for the perturbed semimajor axis in terms of its Newtonian equivalent, which shows how the PN perturbations slightly reduce the orbit's size: In order to quantify the time-scale for the perturbed orbit's decay, we must manipulate Equation (18) to obtain a new variable am, defined as the value of the Newtonian semimajor axis when the binary described via the perturbed orbit can be said to have reached its coalescence. In other words, we seek the value of aN for which the perturbed periapsis p1 = a1(1 − e1) of the orbit has reached the effective Schwarzschild radius of the two-body system (hereafter, simply Schwarzschild radius), rS = 2GM c −2 . If the periapsis is not smaller than the Schwarzschild radius at the beginning of the evolution, an orbit is expected to circularise via the emission of GWs before coalescing (in absence of other dynamical effects), allowing us to replace p1 = a1(1 − e1) with a1 in the following equation: In order to compare the time-scales for the decay, we can simply integrate Equation (9) in the two different cases: the Newtonian orbit has to decay until the Schwarzschild radius is reached, whereas the perturbed orbit has to decay until the value am is reached. Since am is always larger than the Schwarzschild radius, we come to a first qualitative result: PN orbits decay more quickly than what Peters' formula predicts, when comparing binaries with identical initial energy and angular momentum.
In order to solve the decay Equation (9) analytically, we have to make the simplifying assumption of small initial eccentricity (e 2 N ≈ 0). For a given initial semimajor axis a0, the time evolution reads The decay process takes a binary from its initial semimajor axis a0 (and associated energy) to coalescence. To know the duration of this process, we set Equation (20) equal to the value of the last allowed semimajor axis before the coalescence is achieved, and solve for the variable t. We obtain two different time-scales -tP for a Newtonian orbit (Peters' time-scale) and tc for a perturbed orbit (corrected time-scale) -by respectively choosing the final semimajor axis values as 2GM c −2 and am. To improve readability, we express the time-scales in terms of multiples of the Schwarzschild radius, by defining the number n ≥ 1 via n = a0/rS. We obtain To account for large initial eccentricities, we note that most of the decay time is spent in the neighbourhood of the initial conditions (this is the same argument used in Peters 1964). Therefore, we can adjust the time-scales by dividing them by the eccentricity enhancement function f (e) evaluated at the initial eccentricity e0: .
The equations listed so far display a dependence on the binary mass ratio. For convenience, we report here the ratio νE = tc/tP and the difference δE = tc −tP in the two limiting cases of q → 0, and q → 1, for orbits whose initial energy and angular momentum are the same.

Orbits with identical initial shape
In this section, we investigate the difference in the decay time-scale when we compare Newtonian and PN orbits with the same initial shape. This comparison is the crucial point of the present paper as it has more practical implications than the one discussed in the previous section. This is because, due to the typically adopted expression for Peters' time-scale, one usually computes the orbital decay time by plugging in an initial (Keplerian) semimajor axis and an eccentricity, rather than the associated (Keplerian) energy and angular momentum. We set the perturbed orbital parameters a1 and e1 equal to the Newtonian ones. In the case of a first-order PN analysis, the equations that one obtains are with a1 and e1 given by Equations (14) and (15). By rearranging the previous expressions we obtain a set of conditions for the PN energy, E1, and the PN angular momentum, h1. We can interpret these conditions as the values that E1 and h1 must have in order for the perturbed orbit to be identical to a Newtonian orbit with semimajor axis aN and eccentricity eN. The conditions read Equations (30) and (31) show that the perturbed orbit has a slightly larger energy and squared angular momentum than the non-perturbed one. In other words, it takes more energy and more angular momentum to sustain a certain shape in a PN gravity field rather than in a Newtonian gravity field. Since excess energy and momentum take time to be radiated away, we come to a second qualitative result: PN orbits decay more slowly than what Peters' formula predicts, when comparing binaries with identical initial shape (i.e. semimajor axis and eccentricity).
To quantify the time-scale for the decay, it is useful to define two new parameters, an effective semimajor axis, a eff , and an effective eccentricity, e eff : e 2 eff = 1 + 2h 2 1 E1.
Note that these effective quantities do not carry information on the actual shape of the orbit; we introduce them as convenient variables for the analysis. Equations (36) and (34) imply that: (i) the slight increase of the perturbed energy E1 with respect to the Newtonian energy causes the effective semimajor axis to be larger than the actual one; (ii) the slight increase of the perturbed squared angular momentum h 2 1 causes the squared effective eccentricity to decrease; and (iii) the squared effective eccentricity may in principle become negative (indeed, it is always negative if the physical orbit is circular).
Summing up, the effective ellipse is larger and more circular than the actual (PN) orbit because it is constructed using the larger PN energy and angular momentum. However, we incur into a problem if the physical orbit is already very circular. The crucial point is that PN orbits may have more angular momentum than what the classical definition of eccentricity and semimajor axis allow for a given energy. In other words, if the physical orbital eccentricity is smaller than a critical value ec, no Newtonian orbit exists with corresponding energy and angular momentum; as implied by the possible appearance of a negative "squared eccentricity". In Figure 1, we show the effective versus the Newtonian eccentricity for orbits with different initial sizes. There exists, however, a simple solution to this apparent inconsistency, which extends the validity of the effective orbit construction to all possible perturbed orbits. Note that, if negative, the absolute value of the squared effective eccentricity is always at least smaller than a 1 PN term: Therefore, if we allow for a small adjustment (smaller than a 1 PN term) to the squared angular momentum of the effective ellipse, we can neglect the cases where the squared effective eccentricity would be negative and instead treat the effective ellipse as having zero eccentricity.
In analogy with the previous subsection, the time-scale for the orbit's decay can now be computed by integrating Equation (9) from the value a eff to a value am. This time, however, the value of am is by construction equal to the Schwarzschild radius (this is also true for higher orders): If we want to proceed analytically, we have to assume low eccentricity (e 2 N ≈ 0). With this assumption, we can compute the corrected time-scale, τc, for the GW decay of a perturbed orbit by setting aN(t) = am and replacing the initial semimajor axis with the effective semimajor axis in Equation (20), and subsequently solving for t. We then reintroduce the effect of large eccentricities by means of the enhancement function f (e), which must be evaluated respectively at e eff for the perturbed case and at e0 for the Newtonian case. The results are most readable when expressed in terms of multiples a0 = nrS of the Schwarzschild radius: τc = 5(q + 1) 2 n + 7q 2 +13q+7 As done in the previous section, we report the ratio ) against the Newtonian eccentricity for orbits with different initial sizes (a N = nr S ). The effective eccentricity is always slightly smaller than its Newtonian counterpart. The significance of this difference is determined by the orbit's distance to the Schwarzschild radius of the system: it decreases for larger orbits, as expected for any PN perturbation. In very strong gravity regimes, there is no way to describe the energy and angular momentum of the orbit as an effective Newtonian ellipse. In the plot, we see this effect when the effective eccentricity becomes zero. νg = τc/tP and difference δg = τc − tP between time-scales in the two limiting cases of q → 0, and q → 1, computed for the same initial orbital shapes.

Discussion of the analytic results
In the previous sections, we have defined four quantities (νE, δE, νg, and δg) that describe the deviations from Peters' formula when taking PN perturbations to the orbit into account. Now, let us take a closer look at the ratios of the corrected to the standard Peters' time-scales (νE, νg).  Figure 2. We show the behaviour of the corrected to Peters' time-scale ratios ν E and νg plotted against the initial semimajor axis measured in Schwarzschild radii (a 0 = nr S ). The ratio ν E (assuming orbits starting with the same initial energy and angular momentum) does not depend on the orbital eccentricity, contrarily to νg (assuming orbits starting with the same initial orbital shape); for this reason, the latter is plotted for different values of the initial eccentricity e 0 . The ratio ν E very quickly approaches unity; on the other hand, νg can still be significantly different from one for very large n if the initial orbit is very eccentric, suggesting that the magnitude of this correction is sensitive to the periapsis of the system. Both ratios approach unity for n approaching infinity. This is expected, because any PN perturbation to the binary orbit has to vanish as we increase the relative distance between the two bodies. However, a quick glance at Figure 2 shows that, as n increases, the ratio νE approaches unity extremely quickly, whereas the ratio νg remains quite large for a wide range of n. More specifically, the behaviour of νg (which is also a function of the initial orbital eccentricity) strongly depends on the value of the orbital periapsis: the smaller the periapsis, the more enhanced the deviation from Peters' time-scale.
We now focus on the value of νg. Its explicit analytical formula is still rather complex. It is possible to find a much simpler formula that can approximate it very well by taking the limit of very large n (large orbits) and e0 approaching unity (very eccentric orbits): νg → 1 + 7 6q 2 + 11q + 6 8 (1 − e0) n(q + 1) 2 .
It turns out that the above approximation actually covers most of the phase space without producing a significant loss of precision. This aspect can be appreciated by look- ), for a wide range of initial conditions. For very eccentric orbits (1 − e 0 < 0.1) the error is of the magnitude of a PN term r S /p, where p is the periapsis. This implies, that the absolute error that arises by using this approximation instead of νg is as large as second-order PN term and can be neglected. The relative error increases for orbits with intermediate eccentricity (1 − e 0 ≈ 0.5) but dramatically falls for almost circular orbits (1 − e 0 ≈ 0.8). It also shows has a complicated behaviour for very low eccentricities but always remains small.
ing at Figure 3, which shows the discrepancy between the approximated νg (Equation 43) and the exact 1 PN value (Equation 39). In short, the absolute error of the approximation with respect to the full 1 PN result is of the order of a 2 PN term, and can therefore be neglected. We name the approximated value for the ratio νg the "PN correction factor" Q: Q = 1 + 7 6q 2 + 11q + 6 rS 8 (1 − e0) a0(q + 1) 2 1 + 5 rS a0(1 − e0) ≥ 1. (44) The PN correction factor can be easily interpreted: it shows that the time-scale of gravitational decay is affected in a significant way by PN perturbations if the orbital periapsis p comes close to the Schwarzschild radius of the system. The magnitude of the correction scales as K(rS/p), with K(q) being a monotonic function of the mass ratio q that decreases from 5.25 (or 21/4) to approximately 5.03 (or 161/32) as q ranges from 0 to 1. The variation is very weak because the first-order perturbations to the semimajor axis and eccentricity only weakly depend on the mass ratio. By using this simple factor, we can restate the qualitative result of Section 3.2 in a more precise form, valid when comparing orbits with identical initial semimajor axis and eccentricity: Peters' time-scale is too short by a factor of Q that is relevant whenever the periapsis p of the binary comes close to the Schwarzschild radius rS. The factor is given by the formula Q = 1 + 5rSp −1 .
Let us now briefly discuss the time differences δE (Equation 25) and δg (Equation 40). An important feature of these quantities is that they both scale with a system-specific dimensional term: rS/(cq). Its magnitude may vary considerably depending on the two masses m1 and m2 in the binary system: For a compact object of 10 M interacting with a Milky Way-like SMBH of 4.3 × 10 6 M (Schödel et al. 2002;Gillessen et al. 2017), its value amounts to ∼10 yr. The above factor has to be multiplied by a quantity of the order of 2f (e0) −1 in the case of δE and by a fourth-order polynomial in the initial semimajor axis n in the case of δg.
As one may already have noted, the magnitude of the corrections described so far is much larger if one compares a Newtonian and a perturbed orbit with identical initial orbital shape, i.e. semimajor axis and eccentricity (Section 3.2). The reason for this lies in the fact that the construction of the effective ellipse "adds" some energy and angular momentum at the beginning of the evolution, when the orbit is large and gravitational radiation is weak. If radiation is weak, then a small addition of energy will take a long time to be dissipated. In the case of orbits that start with the same energy and angular momentum (Section 3.1), the opposite is true: even though the system has to radiate slightly less energy than expected to reach coalescence, the difference only comes into play when the orbit is small and gravitational radiation is strong. When radiation is strong, a small difference in the energy will not affect the decay time in a significant manner. A more quantitative formulation comes from noting that Peters' time-scale tP scales with the fourth power of the semimajor axis: a small difference ∆a will change it proportionally to the cubed power of a itself: tP ∝ a 4 =⇒ ∆t ∝ 4a 3 ∆a.
We achieve the same conclusion by replacing the semimajor axis a with an inverse energy E −1 and remembering that the energy of a bound system is negative. A small difference in the amount of energy to be dissipated only affects the duration of the orbit's decay if applied when the orbit is large.

DISCUSSION AND CONCLUSIONS
A useful way to visualise the effect of the newly found correction factor is to compute isochrone curves. This amounts to setting the time-scale of gravitational decay equal to a fixed time τ (i.e. QtP = τ or tP = τ ) and, for any given eccentricity, solving for the semimajor axis. We do it for the cases of a SMBH binary and a stellar-mass BH binary and plot the results in the (a, 1 − e) phase space in Figure 4. The plots would be qualitatively similar for other types of GW sources such as EMRIs and IMRIs. The effect of the correction factor is to shift the isochrones towards the lower-left corner of the phase-space plots. Equivalently, any point in phase space that used to lay on an isochrone labelled by the time tP (given by Peters' time-scale) now lays on an isochrone labelled with the corrected time QtP. We also show some specific numerical examples in Table 1.
Regardless of total mass and mass ratio, there always exists a large portion of phase space where the correction The case for an EMRI with m 1 = 10 6 and m 2 = 10 M is indistinguishable from the SMBH binary in the top panel, but with isochrones labels ranging from hundreds to thousands of years (i.e. replace τ with 10 6 τ ). In both panels, the solid coloured lines represent the isochrone curves (for different τ ) obtained with the correction factor (Qt P = τ ), whereas the dashed coloured lines are obtained with Peters' time-scale (t P = τ ). The black solid curve represents orbits whose periapsis is one Schwarzschild radius. The black dashed curve represents the region in phase space where the correction factor Q is good at describing the 1 PN behaviour: leftwards of the curve, its absolute error is smaller than a 2 PN term. There would be no appreciable difference in the plot if one were to use the ratio νg instead of the correction factor Q.
factor Q is a significant fraction of the total decay time. It is worth noting that the discrepancy between the corrected and classical Peters' formula is enhanced for binaries with a large eccentricity. As a consequence, the sources that are going to be affected the most by the correction are those that are intrinsically more promising for GW observations (Equation 12). A wide variety of physical processes are known to enhance the eccentricity of a candidate GW source (and thus promote the inspiral) at all mass scales. For instance, the eccentricity of stellar compact binaries may increase as a result of the supernova kicks experienced by stellar BHs and neutron stars at birth (Giacobbo & Mapelli 2019b,a); in addition, evolution in triplets can induce large eccentricities due to Kozai-Lidov (Kozai 1962;Lidov 1962) cycles and chaotic evolution for both stellar (e.g. Antonini & Rasio 2016;Bonetti et al. 2018) and supermassive binaries (e.g. Bonetti et al. 2016). As a notable example, Bonetti et al. (2019) showed that evolution in triplets can produce supermassive binaries whose eccentricity is > 0.9 when the source enters the LISA band. Finally, most EMRIs are expected to enter the GW phase via the classical two-body relaxation channel (Hils & Bender 1995;Sigurdsson 1997) with eccentricities very close to unity (Amaro-Seoane et al. 2007b). For these reasons, we claim that the correction factor for Peters' time-scale will have an effect on the predicted detection rates of the LIGO-Virgo and LISA observatories. Peters' time-scale is commonly used to estimate the likelihood of a population of compact objects to produce gravitational signals detectable by GW observatories. It is used both in the modelling of the astrophysics (Barausse 2012;O'Leary et al. 2009) and in signal analysis (Klein et al. 2016). A fractional increase in the time-scale of gravitational decay corresponds to a fractional increase in the time that any signal remains in the optimal frequency band of the observatory. Moreover, it will have some astrophysical implication on the rates at which sources GW are produced. The interplay between these two effects will affect how many and what potential sources of GW are considered promising for both LISA and LIGO-Virgo observatories. We are currently working on a quantitative formulation of this effect, with a focus on the upcoming detector LISA. Preliminary calculations show that it might increase the predicted rates because of the better signal to noise ratios arising from longer decay durations. However, this effect must be weighted against other astrophysical implications of an increased GW timescale. As opposed to the correction factor Q, the absolute difference between Peters' and the corrected time-scale strongly depends on the mass ratio and total mass of the system. In Figure 4, this can be appreciated by noting that the labels of the isochrone curves range from mere seconds (bottom panel) to days (top panel) or even thousands of years (not shown) depending on the chosen masses, even though the curves themselves lie in an equivalent region of phase space (equivalent meaning similar periapses with respect to the Schwarzschild radius of the system). For the special case of EMRIs (top panel with τ multiplied by 10 6 ), mere fractional corrections to the GW time-scale can be of the order of astrophysical time-scales describing other dynamical processes around a SMBH. The correction factor Q might therefore have significant astrophysical implications. Most EMRIs are believed to be generated by the scattering of a stellar-mass compact object to an orbit with very low angular momentum (Hils & Bender 1995;Sigurdsson 1997) via the two-body relaxation process. However, two body scatterings can also deflect the orbit of an inspiral before it has had time to complete its gravitational decay. In order to understand whether a compact object reaches its GW induced coalescence, one has to compare the time-scale of gravitational decay to the two-body relaxation time-scale, given by . (47) where MSMBH = 4.3 × 10 6 M is the SMBH mass, ms = 1 M is the mass of a typical field star, mo = 10 M is the mass of the compact object undergoing the process, N = 10 6 is the number of stars in a given unit sphere of radius of 1 pc, γ is the exponent of the number density profile, and Cγ is a constant of order unity. The values of the parameters introduced above are chosen to represent the relaxation process in a Milky Way-analogue (see, e.g. Bortolas & Mapelli 2019). In Figure 5, we plot the curves that equate the twobody relaxation time-scale to the time-scale of gravitational decay for different power-law stellar distributions, solving tP = tr or QtP = tr for the seminajor axis, for a given ec-centricity (Amaro-Seoane et al. 2007a;Merritt et al. 2011;Bortolas & Mapelli 2019). In very rough terms, objects that fall above such curves in (a, 1 − e) phase space are unlikely to undergo an EMRI as their orbit will be likely scattered in a 'safer' zone by two-body relaxation. The effect of the correction factor Q is again to shift the curves towards the lower-left corner of the phase-space plots. For very steep power laws (i.e. γ = 2; 2.5), the shift is significant for all possible orbits. Steep density profiles could be produced by the effect of strong mass segregation (Alexander & Hopman 2009;Preto & Amaro-Seoane 2010) aided by the slow natal kicks received by stellar BHs at birth (Bortolas et al. 2017). For shallower power laws (i.e. γ = 1; 1.5) expected in the case of weak mass segregation (Bahcall & Wolf 1976, 1977, the shift is significant for eccentric orbits. This suggests that the correction factor Q will have a direct effect on previously computed event rates for EMRIs that use Peters' time-scale to model gravitational radiation (Amaro-Seoane 2018; Amaro-Seoane et al. 2007a). We are currently working on a quantitative estimate on the production rate that takes Q into account, along with the associated change in expected EMRI detection rates for LISA.

Conclusion
In this paper, we argue that the quadrupole formula can be used to describe energy radiation in perturbed PN orbits up to third order. We then derive a revised form of Peters' time-scale for the GW induced decay of compact objects by relaxing the assumption of Keplerian orbits in the historical derivation by Peters & Mathews (1963) and Peters (1964). Our key result is represented by the correction factor Q = 1 + 5(rS/p), where p is the orbit periapsis and rS is the effective Schwarzschild radius, valid for 1 PN orbits with identical initial orbital parameters. This factor can be multiplied by the well-known Peters' time-scale to obtain a more accurate estimate of the duration of a GW-induced inspiral, valid for any kind of binary source. The difference between the uncorrected and the corrected time-scales depends solely on the periapsis of the orbit in question and can be of fractional order (Q − 1 ≈ 0.1-1) for a large range of orbital parameters. We claim that this difference has an effect on current event-rate predictions for GW sources, as it influences both the astrophysical modelling of gravitational radiation and the quality of detectable signals. The advantage of the formula Q = 1 + 5(rS/p) is that it captures the largest deviation of fully relativistic orbits from Keplerian ones while remaining algebraically simple. Therefore, we propose that the correction factor Q should be implemented in future event-and detection-rate estimates for LISA and LIGO-Virgo sources, in order to model the effects of gravitational radiation more accurately. . We show two regions of phase space in the semimajor axis and eccentricity, for an EMRI system with m 1 = 4.3 × 10 6 and m 2 = 10 M . The black solid line represents the effective innermost stable circular orbit (ISCO; 3r S ), whereas the red and blue solid (dashed) lines are the curves obtained by equating the two-body relaxation time-scale to the corrected (uncorrected) Peters' time-scale, for four different density power laws (γ = 2; 2.5, top panel, and γ = 1; 1.5, bottom panel). Note how the corrected lines diverge from the uncorrected ones for very eccentric orbits. This is because the slope of the dashed curves is less steep than the lines of constant periapsis (grey lines parallel to the ISCO line): if one follows the curves from right to left, the periapsis will slowly decrease, in turn increasing the magnitude of the PN correction. The orbits that live under the ISCO line are expected to decay too quickly to produce easily detectable gravitational signals and are therefore not relevant. Similarly, orbits above the red and blue lines are scattered before gravitational radiation has time to significantly affect their orbit. This leaves us with the phase-space volume between the black and coloured lines as a proxy of the region of potential EMRI candidates. The various grey lines indicate different values for the correction factor Q and are a measure of how strong PN effects are in that region. The black dashed curve represents the region in phase space where the correction factor Q is good at describing the 1 PN behaviour: leftwards of the curve, its absolute error is smaller than a 2 PN term. For shallower density profiles, differences are only noticeable for extremely eccentric orbits.

a3(em)
To obtain the duration ∆tE of the orbital decay, we can integrate the evolution equation of the Newtonian eccentricity from the initial condition e0 to its value at coalescence, em: ∆tE = (A6)

A2 Orbits with identical initial shape
Our underlying assumption in this section is that the geometrical 3 PN orbital features are equal to the unperturbed Newtonian ones. We can therefore set them equal to each other: By rearranging Equations (A7) and (A8), we obtain a set of conditions for the PN energy E3 and the PN angular momentum h3. We can interpret these conditions as the values that E3 and h3 must have in order for the perturbed orbit to be identical to the Newtonian orbit with semimajor axis aN and eccentricity eN: h 2 3 = h 2 3 (aN, eN).
In analogy with Section 3.2, We then build the "effective ellipse" of the system, i.e. a fictitious Newtonian orbit constructed from the perturbed energy and angular momentum. It is characterised by a semimajor axis a eff and an eccentricity e eff containing third-order PN information: e 2 eff = 1 + 2E3h 2 3 .
Once the explicit numerical values of the effective parameters are obtained, one can use the standard Peters' formula to compute the duration of the decay. It turns out that the higher-order corrections are very small with respect to the first-order correction derived in the main text. For this reason, we have omitted the self-consistent calculations from it. Generally, such high PN orders are not relevant for the applications in which Peters' formula is commonly used. This paper has been typeset from a T E X/L A T E X file prepared by the author.