$$
\begin{equation}
\begin{aligned}
I_{ijkl} = \mu & \left\{ x_{\langle ijk l \rangle} \left[ 1 - 3v + \gamma \left( \frac{3}{110} - \frac{25}{22}v + \frac{69}{22}v^2 \right) \right. \right. \\
& \left. \left. + \gamma^2 \left( -\frac{126901}{200200} - \frac{58101}{2600}v + \frac{204153}{2860}v^2 + \frac{1149}{1144}v^3 \right) \right] \right. \\
& \left. + \frac{r}{c} x_{\langle ijk v l \rangle} \gamma^2 \left( \frac{976}{55}v - \frac{3104}{55}v^2 \right) \right. \\
& \left. + \frac{r^2}{c^2} x_{\langle ij v k l \rangle} \left[ \frac{78}{55} (1 - 5v + 5v^2) \right. \right. \\
& \left. \left. + \gamma \left( \frac{30583}{3575} - \frac{107039}{3575}v + \frac{8792}{715}v^2 - \frac{639}{715}v^3 \right) \right] \right. \\
& \left. + \frac{71}{715} \frac{r^4}{c^4} v_{\langle ijk l \rangle} (1 - 7v + 14v^2 - 7v^3) \right\} + \mathcal{O}\left(\frac{1}{c^6}\right),
\end{aligned}
\tag{454a}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
J_{ijk} = \mu & \left\{ \epsilon_{ab \langle i x_j k \rangle a} v_b \left[ 1 - 3v + \gamma \left( \frac{181}{90} - \frac{109}{18}v + \frac{13}{18}v^2 \right) \right. \right. \\
& \left. \left. + \gamma^2 \left( \frac{1469}{3960} - \frac{5681}{264}v + \frac{48403}{660}v^2 - \frac{559}{3960}v^3 \right) \right] \right. \\
& \left. + \frac{r}{c} \epsilon_{ab \langle i X_j a v_k \rangle b} \gamma^2 \left( \frac{472}{45}v - \frac{496}{15}v^2 \right) \right. \\
& \left. + \frac{r^2}{c^2} \epsilon_{ab \langle i X_a v_j k \rangle b} \left[ \frac{7}{45} (1 - 5v + 5v^2) \right. \right. \\
& \left. \left. + \gamma \left( \frac{1621}{990} - \frac{4879}{990}v + \frac{1084}{495}v^2 - \frac{259}{990}v^3 \right) \right] \right\} + \mathcal{O}\left(\frac{1}{c^6}\right).
\end{aligned}
\tag{454b}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
I_{ijklm} = -\mu\Delta & \left\{ x_{\langle i j k l m \rangle} \left[ 1 - 2v + \gamma \left( \frac{2}{39} - \frac{47}{39}v + \frac{28}{13}v^2 \right) \right. \right. \\
& \left. \left. + \gamma^2 \left( -\frac{716}{819} - \frac{15595}{546}v + \frac{32245}{546}v^2 + \frac{631}{819}v^3 \right) \right] \right. \\
& \left. + \frac{r^2}{c^2} x_{\langle i j k v l m \rangle} \left[ \frac{70}{39} (1 - 4v + 3v^2) \right. \right. \\
& \left. \left. + \gamma \left( \frac{1238}{117} - \frac{986}{39}v + \frac{290}{39}v^2 - \frac{70}{117}v^3 \right) \right] \right. \\
& \left. + \frac{85}{273} \frac{r^4}{c^4} x_{\langle i v j k l m \rangle} (1 - 6v + 10v^2 - 4v^3) \right\} + \mathcal{O}\left(\frac{1}{c^5}\right),
\end{aligned}
\tag{454c}
\end{equation}
$$
$$
\begin{equation}
\begin{aligned}
J_{ijkl} = -\mu\Delta & \left\{ \epsilon_{ab\langle iX_{jkl}\rangle a} v_b \left[ 1 - 2v + \gamma \left( \frac{20}{11} - \frac{155}{44}v + \frac{5}{11}v^2 \right) \right. \right. \\
& \left. \left. + \gamma^2 \left( -\frac{1103}{5005} - \frac{141333}{5720}v + \frac{69783}{1144}v^2 - \frac{83}{572}v^3 \right) \right] \right. \\
& \left. + \frac{r^2}{c^2} \epsilon_{ab\langle iX_{ja}v_{kl}\rangle b} \left[ \frac{4}{11} (1 - 4v + 3v^2) \right. \right. \\
& \left. \left. + \gamma \left( \frac{2583}{715} - \frac{21887}{2860}v + \frac{1411}{572}v^2 - \frac{28}{143}v^3 \right) \right] \right\} + \mathcal{O}\left(\frac{1}{c^5}\right).
\end{aligned}
\tag{454d}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
I_{ijklmn} = \mu & \left\{ x_{\langle ijklmn\rangle} \left[ 1 - 5v + 5v^2 + \gamma \left( \frac{1}{14} - \frac{3}{2}v + 6v^2 - \frac{11}{2}v^3 \right) \right] \right. \\
& \left. + \frac{15}{7} \frac{r^2}{c^2} x_{\langle ijk l v_{mn}\rangle} (1 - 7v + 14v^2 - 7v^3) \right\} + \mathcal{O}\left(\frac{1}{c^4}\right),
\end{aligned}
\tag{454e}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
J_{ijklm} = \mu & \left\{ \epsilon_{ab\langle iX_{jklm}\rangle a} v_b \left[ 1 - 5v + 5v^2 \right. \right. \\
& \left. \left. + \gamma \left( \frac{1549}{910} - \frac{1081}{130}v + \frac{107}{13}v^2 - \frac{29}{26}v^3 \right) \right] \right. \\
& \left. + \frac{54}{91} \frac{r^2}{c^2} \epsilon_{ab\langle iX_{jka}v_{lm}\rangle b} (1 - 7v + 14v^2 - 7v^3) \right\} + \mathcal{O}\left(\frac{1}{c^4}\right).
\end{aligned}
\tag{454f}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
I_{ijklmno} = -\mu\Delta & \left\{ x_{\langle ijklmno\rangle} \left[ 1 - 4v + 3v^2 + \gamma \left( \frac{3}{34} - \frac{26}{17}v + \frac{83}{17}v^2 - \frac{57}{17}v^3 \right) \right] \right. \\
& \left. + \frac{42}{17} \frac{r^2}{c^2} x_{\langle ijklm v_{no}\rangle} (1 - 6v + 10v^2 - 4v^3) \right\} + \mathcal{O}\left(\frac{1}{c^3}\right),
\end{aligned}
\tag{454g}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
J_{ijklmn} = -\mu\Delta & \left\{ \epsilon_{ab\langle iX_{jklmn}\rangle a} v_b \left[ 1 - 4v + 3v^2 \right. \right. \\
& \left. \left. + \gamma \left( \frac{13}{8} - \frac{25}{4}v + \frac{55}{12}v^2 - \frac{2}{3}v^3 \right) \right] \right. \\
& \left. + \frac{5}{6} \frac{r^2}{c^2} \epsilon_{ab\langle iX_{jka}v_{mn}\rangle b} (1 - 6v + 10v^2 - 4v^3) \right\} + \mathcal{O}\left(\frac{1}{c^3}\right).
\end{aligned}
\tag{454h}
\end{equation}
$$

All the other higher-order moments are required at the Newtonian order (there are no $1/c$ corrections in $I_{ijklmno}$ and $J_{ijklmno}$), at which order they are trivial to compute, with result ($\forall l \in \mathbb{N}$)

$$
I_L = \mu \sigma_L(v) x_{\langle L \rangle} + \mathcal{O}\left(\frac{1}{c}\right), \qquad (455a)
$$
$$J_{L-1} = \mu \sigma_{\ell}(v) \epsilon_{ab \langle i_{l-1} X_{L-2} \rangle a} v_b + \mathcal{O}\left(\frac{1}{c}\right). \quad (455b)$$

Here we use the notation $\sigma_{\ell}(v) \equiv X_2^{\ell-1} + (-)^{\ell} X_1^{\ell-1}$, with $X_1 = \frac{m_1}{m}$ and $X_2 = \frac{m_2}{m}$ such that $X_1 + X_2 = 1$, $X_1 - X_2 = \Delta$ and $X_1 X_2 = \nu$. More explicit expressions are $(k \in \mathbb{N})$:

$$\sigma_{2k}(v) = \sum_{p=0}^{k-1} (-)^p \frac{2k-1}{2k-1-p} \binom{2k-1-p}{p} v^p, \quad (456a)$$

$$\sigma_{2k+1}(v) = -\Delta \sum_{p=0}^{k-1} (-)^p \binom{2k-1-p}{p} v^p, \quad (456b)$$

where $\binom{n}{p}$ is the usual binomial coefficient.

Evidently the results (448)–(455) permit the control only of that part of the total energy flux which is generated by the source multipole moments $I_L$ and $J_L$; from (179a) we have up to 4.5PN order

$$\mathcal{F}|_{I_L, J_L} = \frac{G}{c^5} \left\{ \frac{1}{5} (I_{ij}^{(3)})^2 + \frac{1}{c^2} \left[ \frac{1}{189} (I_{ijk}^{(4)})^2 + \frac{16}{45} (J_{ij}^{(3)})^2 \right] \right. \\ + \frac{1}{c^4} \left[ \frac{1}{9072} (I_{ijkl}^{(5)})^2 + \frac{1}{84} (J_{ijk}^{(4)})^2 \right] \\ + \frac{1}{c^6} \left[ \frac{1}{594000} (I_{ijklm}^{(6)})^2 + \frac{4}{14175} (J_{ijkl}^{(5)})^2 \right] \\ \left. + \frac{1}{c^8} \left[ \frac{1}{52123500} (I_{ijklmn}^{(7)})^2 + \frac{1}{213840} (J_{ijklm}^{(6)})^2 \right] + \mathcal{O}\left(\frac{1}{c^{10}}\right) \right\} \quad (457)$$

The time derivatives of these source moments are computed by means of the quasi circular-orbit equations of motion given by Eq. (371) together with (373). Next we need to add to the flux the contributions coming from the relations between radiative and canonical moments in Sect. 2.4.5, and between canonical and source ones in Sect. 2.4.6.

When computing the tails we shall also require, in addition to the above source moments $I_L$ and $J_L$, the binary's mass monopole M or ADM mass. In a realistic model where the binary system has been formed as a close compact binary at a finite instant in the past, this mass is equal to the sum of the rest masses $m = m_1 + m_2$, plus the total (constant) binary's mass-energy $E/c^2$ given by Eq. (374). To 3PN order the mass reads

$$M = m \left[ 1 - \frac{v}{2}\gamma + \frac{v}{8}(7-v)\gamma^2 + \frac{v}{16}(7-49v-v^2)\gamma^2 + \mathcal{O}\left(\frac{1}{c^8}\right) \right]. \quad (458)$$

This is required to compute the leading tail at 4.5PN order (i.e. 3PN relative order). Notice that 3PN order in M corresponds to 2PN order in E.
### 3.4.2 Post-adiabatic evaluation of hereditary integrals

The most important problem we face is the computation of the hereditary tail and related non-linear integrals occuring at the 4PN order, for instance the dominant 1.5PN quadrupole tail (239). To compute the tail integrals, one needs to specify the orbit's behaviour of the compact binary system in the remote past, as the effect is not localized in time, but integrates over the whole past history of the source. In the case of quasi-circular orbits, and up to 3.5PN precision in the moments, an adiabatic approximation (considering the orbital elements $r$ and $\Omega$ to be constant in time) is sufficient, and the integrals can be computed following the method of Blanchet and Schäfer (1993), Arun et al. (2004). However, Eqs. (372) show that this adiabatic approximation is no longer valid at the relative 2.5PN precision. As the tail enters at 1.5PN order in the moments, the first “post-adiabatic” correction will affect the moments at the 4PN order, thus we need to properly evaluate it in order to consistently derive the 4PN flux and (2, 2) GW mode.

As shown by Arun et al. (2004), the tail integrals on quasi-circular orbits reduce to elementary integrals of the type

$$ \mathcal{I}_{a,n}(U) = \int_{0}^{+\infty} d\tau [x(U-\tau)]^a e^{-in\phi(U-\tau)} \ln\left(\frac{\tau}{\tau_0}\right), \quad (459) $$

where $U$ is the current time, $U - \tau$ is any time in the past like in (239), and the orbital phase $\phi$, the orbital frequency $\Omega = \dot{\phi}$ and the PN parameter $x$ defined by (375) are integrated over the binary's past evolution. In (459), $n$ is a non-zero integer (the case $n = 0$ does not appear at 4PN order), $a$ can be any power ($a = 5$ for the leading tail), and the constant $\tau_0$ denotes either $2r_0/c$ or $2b_0/c$, see for instance Eq. (243). The strategy is to notice that the integral (459) involves a fast oscillating exponential, which is derivable in the form of an asymptotic series which is nothing but the post-adiabatic approximation. In the application to 4PN order we need only the first post-adiabatic correction, relevant to the dominant 1.5PN tail term (239); the other tail integrals like those appearing in the tails-of-memory (243) can be computed in the adiabatic limit.

Remark first that the difference between the current phase $\phi(U)$ and the phase $\phi(U - \tau)$ in the past is of the order of the inverse of the radiation reaction scale, i.e., $\phi(U) - \phi(U - \tau)$ scales as $O(c^5)$. To describe the radiation-reaction scale, we introduce a dimensionless adiabatic parameter at current time $U$,

$$ \xi(U) \equiv \frac{\dot{\Omega}(U)}{\Omega^2(U)} = O\left(\frac{1}{c^5}\right), \quad (460) $$

and compute the integral (459) in the post-adiabatic limit where $\xi(U) \to 0$. To this end, let us pose

$$ \phi(U) - \phi(U - \tau) = \frac{\lambda}{\xi(U)}, \quad (461) $$

and change the integration variable from $\tau$ to $\lambda$, so that
$$ \mathcal{I}_{a,n} = \frac{Gm}{c^3} \frac{e^{-in\phi(U)}}{\xi(U)} \int_{0}^{+\infty} d\lambda [x(U - \tau(\lambda))]^{a-3/2} \ln\left(\frac{\tau(\lambda)}{\tau_0}\right) e^{in\lambda/\xi(U)}. \quad (462) $$

Here $\tau(\lambda)$ is obtained by inverting Eq. (461) for $\tau$ as a function of $\lambda$. The main point is that, in the limit $\xi(U) \to 0$, since the imaginary exponential in the integrand of (462) oscillates very rapidly, the integral is a sum of alternatively positive and negative contributions which essentially sum up to zero. However there are no oscillations when $\lambda = 0$, therefore the integral is essentially given by the contribution in the neighbourhood of the bound at $\lambda = 0$. In other words we are entitled to compute the integral by replacing the integrand by its formal expansion series when $\lambda \to 0$, which gives a formal (asymptotic) series in powers of $\xi(U)$. In the limit $\tau \to 0$ we can expand Eq. (461) as

$$ \lambda = \xi \sum_{m=0}^{+\infty} \frac{(-)^m}{(m+1)!} \Omega^{(m)} \tau^{n+1}, \quad (463) $$

where we now pose $\xi = \xi(U)$ and $\Omega^{(m)} = (d^m\Omega/dU^m)(U)$ for the quantities at the current time $U$. Since $\Omega^{(1)} \equiv \dot{\Omega} = \xi\Omega^2$ the expansion (463) clearly represents the post-adiabatic approximation in powers of $\xi$ and arbitrary high time derivatives of $\xi$ denoted $\xi^{(m)}$. Note that each time derivative of $\xi$ adds a radiation reaction scale, hence we have $\xi^{(m)} = \mathcal{O}(\xi^{m+1}) = \mathcal{O}(c^{-5m-5})$. To obtain the expansion of (462) we need to invert the series (463) and obtain $\tau$ as a power series in $\lambda$. This is given by the Lagrange inversion theorem as

$$ \tau = \sum_{p=0}^{+\infty} \frac{1}{(p+1)!} f_p \left(\frac{\lambda}{\xi}\right)^{p+1}, \quad (464) $$

where the general coefficients $f_p$ read

$$ \begin{aligned} f_p &= \left(\frac{d}{d\tau}\right)^p \left[ \left(\frac{\tau}{\phi(U) - \phi(U-\tau)}\right)^{p+1} \right]_{\tau=0} \\ &= \left(\frac{d}{d\tau}\right)^p \left[ \left(\sum_{m=0}^{+\infty} \frac{(-)^m}{(m+1)!} \Omega^{(m)} \tau^m \right)^{-p-1} \right]_{\tau=0}, \end{aligned} \quad (465) $$

where $\tau = 0$ is applied after the $p$ differentiations with respect to $\tau$. For instance, up to second order we obtain

$$ \tau = \frac{\lambda}{\xi \Omega} \left[ 1 + \frac{\lambda}{2} + \left( 1 - \frac{\dot{\xi}}{\xi^2 \Omega} \right) \frac{\dot{\lambda}^2}{6} + \mathcal{O}(\lambda^3) \right]. \quad (466) $$

The previous formulas show that the method can be straightforwardly extended to any order. But to compute the tail term at 4PN order, we need only the relative correction of order 2.5PN, which corresponds to the first post-adiabatic approximation, hence just $\tau = \frac{\lambda}{\xi\Omega}[1 + \frac{\lambda}{2} + \mathcal{O}(\lambda^2)]$. To this order the integral (462) becomes (Blanchet et al. 2023b)
$$ \mathcal{I}_{a,n} = \frac{Gm}{c^3} \frac{x^{a-3/2} e^{-in\phi}}{\xi} \times \int_{0}^{+\infty} d\lambda \left[ \left( 1 + \left(1 - \frac{2a}{3}\right)\lambda \right) \ln\left(\frac{\lambda}{\xi\Omega\tau_0}\right) + \frac{\lambda}{2} \right] e^{in\lambda/\xi}, \quad (467) $$

modulo $O(\lambda^2)$ terms. The remaining integral is computed as follows. We transform the complex exponential into a real one by performing the change of variable $\lambda = i\sigma$ for $n > 0$, and $\lambda = -i\sigma$ for $n < 0$. The integration now takes place along the imaginary axis, which can be remedied by resorting to the Cauchy theorem on the closed contour made of the three following pieces, to be considered in the limit $R \to +\infty$: (i) the path from $\sigma = 0$ to $\sigma = R$ on the real axis, (ii) the oriented quarter of circle of radius $|\sigma| = R$ from arg $\sigma = 0$ to arg $\sigma = -\pi/2$ (or arg $\sigma = \pi/2$ if $n < 0$) and (iii) the segment of the imaginary axis going from $\sigma = -iR$ ($\sigma = iR$ if $n < 0$) to $\sigma = 0$. This leads to

$$ \mathcal{I}_{a,n} = \frac{Gm}{c^3} \frac{x^{a-3/2} e^{-in\phi}}{i\xi} \times \int_{0}^{+\infty} d\sigma \left[ \left( -s(n) + \frac{2a-3}{3} i\sigma \right) \left( \ln\left(\frac{\sigma}{\xi\Omega\tau_0}\right) + i\frac{\pi}{2}s(n) \right) + \frac{i\sigma}{2} \right] e^{-|n|\sigma/\xi}, \quad (468) $$

modulo $O(\sigma^2)$ terms, where $s(n)$ is the sign function. In this form, $\mathcal{I}_{a,n}$ can be integrated, as it boils down to elementary integrals. At the first post-adiabatic order the result is

$$ \mathcal{I}_{a,n} = \frac{Gm}{c^3} \frac{x^{a-3/2} e^{-in\phi}}{in} \left\{ \left( 1 + \frac{2a-3}{3} \frac{\xi}{in} \right) \left[ \ln(|n|\Omega\tau_0) + \gamma_E - i \frac{\pi}{2} s(n) \right] - \frac{4a-9}{6} \frac{\xi}{in} + O(\xi^2) \right\}. \quad (469) $$

Recall that all quantities are evaluated at current time $U$ and that the adiabatic parameter $\xi = \dot{\Omega}/\Omega^2$ is easily computed with Eqs. (372). With this result at hand, we are able to derive the tail integral (239) with 2.5PN relative precision, which impacts the computation and final results at the 4PN order.

For completeness we give the formulas needed to handle any approximation in the post-adiabatic expansion for tail like integrals:

$$ \int_{0}^{+\infty} d\sigma \sigma^{j} e^{-|n|\sigma/\xi} = j! \left(\frac{\xi}{|n|}\right)^{j+1}, \quad (470a) $$

$$ \int_{0}^{+\infty} d\sigma \sigma^j \ln\left(\frac{\sigma}{\xi\Omega\tau_0}\right) e^{-|n|\sigma/\xi} = j! \left(\frac{\xi}{|n|}\right)^{j+1} \left[H_j - \gamma_E - \ln(|n|\Omega\tau_0)\right], \quad (470b) $$

where $H_j = \sum_{k=1}^j \frac{1}{k}$ denotes the harmonic number. For instance the case $j = 1$ yields the standard integral for adiabatic tails (coming back to $\lambda = \pm i\sigma$)
$$ \int_{0}^{+\infty} d\lambda \ln\left(\frac{\lambda}{\xi\Omega\tau_0}\right) e^{-in\lambda/\xi} = \frac{i\xi}{n} \left(\ln(|n|\Omega\tau_0) + \gamma_E + i\frac{\pi}{2}s(n)\right). \quad (471) $$

For iterated tails, i.e. tails-of-tails at 3PN order and tails-of-tails-of-tails at 4.5PN order (see Sect. 2.4.5), the adiabatic approximation is sufficient and the required formulas are available from Gradshteyn and Ryzhik (1980):

$$ \int_{0}^{+\infty} d\lambda \ln^2\left(\frac{\lambda}{\xi\Omega\tau_0}\right) e^{-in\lambda/\xi} = -\frac{i\xi}{n} \left[ \left( \ln(|n|\Omega\tau_0) + \gamma_E + i\frac{\pi}{2}s(n) \right)^2 + \frac{\pi^2}{6} \right], \quad (472a) $$

$$ \int_{0}^{+\infty} d\lambda \ln^3\left(\frac{\lambda}{\xi\Omega\tau_0}\right) e^{-in\lambda/\xi} = \frac{i\xi}{n} \left[ \left( \ln(|n|\Omega\tau_0) + \gamma_E + i\frac{\pi}{2}s(n) \right)^3 + \frac{\pi^2}{2} \left( \ln(|n|\Omega\tau_0) + \gamma_E + i\frac{\pi}{2}s(n) \right) + 2\xi(3) \right]. \quad (472b) $$

Note the presence in the latter formula of the Apéry constant $\zeta(3)$, where $\zeta$ is the Riemann function, but this constant disappears from final results.

Next we deal with hereditary memory type terms, such as the first term of Eq. (240)—leading memory at 2.5PN order—or the first term of Eq. (243)—tail-of-memory at 4PN order. As was shown in Sect. 2.4.4 the memory terms enter only in the mass-type radiative moments $U_L$. In the case of quasi-circular orbits, they reduce to a sum of terms of the form

$$ \mathcal{J}_{a,n}(U) = \int_{0}^{+\infty} d\tau [x(U-\tau)]^a e^{-in\phi(U-\tau)}. \quad (473) $$

Besides the absence of the logarithm, the main difference with the tail integral (459) is the possibility of having $n = 0$, which corresponds to persistent or direct-current (“DC”) terms, in addition to oscillatory alternative-current (“AC”) terms having $n \neq 0$. As we remind below, the memory effect properly speaking is associated with the DC terms.

We focus first on AC terms, having $n \neq 0$. As they only involve a simple, logarithmic-free, integration, such terms will enter in the flux as instantaneous contributions. In particular, since they arise at the half-integer approximations 2.5PN and 3.5PN, they do not contribute to the flux for quasi-circular orbits. The evaluation of the integrals (473) is thus only required for the derivation of the modes. Concerning the quadrupole moment, since the memory effect enters at 2.5PN order [see Eq. (240)], it is thus required at relative 1.5PN order for the mode (2, 2) and thus, one can safely compute it in the adiabatic approximation. As proved by Arun et al. (2004), for circular orbits this is equivalent to taking the frequency $x$ constant together with a linear phase (appropriate for the exact circular orbit), and keeping only the contribution of the integral due to the bound $\tau = 0$. One finds
$$ \mathcal{J}_{a,n} = x^a e^{-in\phi} \int_0^1 d\tau e^{in\Omega\tau} + \mathcal{O}(\xi) = -\frac{x^a e^{-in\phi}}{in\Omega} + \mathcal{O}(\xi). \quad (474) $$

The DC terms obviously do not contribute to the flux, and they only contribute to the modes which have $m = 0$, for example $(\ell, m) = (2, 0)$ for the mass quadrupole. The absence of the fast oscillating exponential in $\mathcal{J}_{a,0}$—in contrast to (459)—generates an inverse power of the adiabatic parameter $\xi$, thus degrading the precision by the radiation-reaction scale 2.5PN. This is the memory effect: as it starts at 2.5PN order in the waveform, it finally enters the (2, 0) mode at Newtonian order. Several methods are possible to evaluate $\mathcal{J}_{a,0}$ (Wiseman and Will 1991; Arun et al. 2004; Favata 2009). Here we rely on the change of integration variables from $\tau$ to $x' = x(U - \tau)$:

$$ \mathcal{J}_{a,0} = \int_{0}^{+\infty} d\tau [x(U - \tau)]^a = \int_{0}^{x(U)} d\tau' \frac{x'^a}{\dot{x}(U - \tau)}, \quad (475) $$

reading $\dot{x}(U - \tau)$ as a function of $x'$ from Eqs. (372), and supposing that $a > 4$, as is the case in practical computations. To 1.5PN relative order, we have

$$ \mathcal{J}_{a,0} = \frac{5Gm}{64c^3v} \frac{x^{a-4}}{a-4} \left[ 1 + \frac{a-4}{a-3} \left( \frac{743}{336} + \frac{11}{4}v \right) x - 8\pi \frac{a-4}{2a-5} x^{3/2} + \mathcal{O}(x^2) \right]. \quad (476) $$

It is interesting to observe that the tail effect in the flux directly influences the DC memory in (476). Finally we look at the interesting case of the tails-of-memory; strictly speaking this is the term in the first line of (243), namely

$$ \mathcal{K}_{ij} = \frac{8G^2M}{7c^8} \int_0^{+\infty} d\rho M_{a\langle i}^{(4)}(u-\rho) \int_0^{+\infty} d\tau M_{j\rangle a}^{(4)}(u-\rho-\tau) \ln\left(\frac{\tau}{\tau_0}\right), \quad (477) $$

where $c\tau_0 = 2r_0 e^{1613/270}$. This expression (Trestini and Blanchet 2023) agrees with the tail-of-memory directly computed from the radiative quadrupole moment at infinity (Favata 2009, 2011c), see Sect. 2.4.4. We first perform the tail-like integral over $\tau$; as we need to evaluate it at relative Newtonian order only, we can safely use the adiabatic approximation. Next we perform the integral over $\tau$, which is found to be a simple DC memory integral of the type (475) with $a = 13/2$. Hence we find

$$ \mathcal{K}_{ij} = \frac{128\pi v^2 c^5}{7G} \ell_{\langle i \ell_j \rangle} \int_0^{+\infty} d\rho [x(U-\rho)]^{13/2} = \frac{4\pi}{7}mv^2 c^2 x^{5/2} \ell_{\langle i \ell_j \rangle}, \quad (478) $$

where $\ell^i = L^i/|L|$ is the constant unit direction of the Newtonian angular momentum, orthogonal to the orbital plane. Interestingly, the tail-of-memory effect (478) gives a contribution to the (2, 0) mode which exactly cancels the 1.5PN correction to the ordinary memory effect, i.e. the term $\propto \pi$ in (476). Hence there is no tail term at 1.5PN order in Eq. (522a) below.
### 3.4.3 The 4.5PN energy flux for circular orbits

The flux (457) due to source multipole moments at 4PN order can now be completed by all contributions coming from the non-linear interactions between the moments. The flux for quasi-circular orbits will finally be expressed as a function of the frequency of gravitational waves. However we have to take into account an important effect, which is that the GW frequency is modified with respect to the orbital frequency because of the propagation of tails in the wave zone. Indeed the phase $\psi$ of the gravitational wave [actually half the phase of the leading (2,2) mode], differs from the orbital phase $\phi$ by a logarithmic, tail-induced phase modulation, as

$$ \psi = \phi - \frac{2GM\Omega}{c^3} \ln\left(\frac{\Omega}{\Omega_0}\right), \qquad (479) $$

where M is the constant ADM mass (i.e. $M = m + E/c^2$), and $\Omega_0$ is related to the length scale $b_0$ parametrizing the gauge transformation between radiative and harmonic coordinates, see Eq. (194), by $c\Omega_0^{-1} = 4b_0 e^{\gamma_E - 11/12}$. The phase modulation (479) was determined by Wiseman (1993), Blanchet and Schäfer (1993), and argued by Blanchet et al. (1996), Arun et al. (2004) to affect the waveform and flux at the 4PN order. Indeed, this phase modulation is responsible for a shift of the GW frequency, i.e. $\Omega_{GW} = d\psi/dt$, with respect to the orbital one, i.e. $\Omega = d\phi/dt$. Therefore the GW frequency, directly measurable from the waveform at infinity, reads

$$ \Omega_{GW} = \Omega - \frac{2GM\dot{\Omega}}{c^3} \left[ \ln\left(\frac{\Omega}{\Omega_0}\right) + 1 \right]. \qquad (480) $$

Using Eqs. (372) for the frequency chirp at the dominant order, and replacing M by $m$ at that order, yields

$$ \Omega_{GW} = \Omega \left\{ 1 - \frac{192}{5} v x^4 \left[ \ln\left(\frac{\Omega}{\Omega_0}\right) + 1 \right] + \mathcal{O}\left(\frac{1}{c^{10}}\right) \right\}, \qquad (481) $$

showing that this is a 4PN effect (i.e. $\propto x^4$), hence the fact that it did not affect previous computations of the phase (Blanchet et al. 1996; Arun et al. 2004). Therefore, once the flux has been obtained at 4PN order in terms of $x = (\frac{Gm\Omega}{c^3})^{2/3}$, we must reexpress it in terms of the measurable GW frequency using the inverse of (481), i.e.

$$ x = x_{GW} \left\{ 1 + \frac{192}{5} v x_{GW}^4 \left[ \ln\left(\frac{x_{GW}}{x_0}\right) + \frac{2}{3} \right] + \mathcal{O}\left(\frac{1}{c^{10}}\right) \right\}, \qquad (482) $$

where $x_{GW} = (\frac{Gm\Omega_{GW}}{c^3})^{2/3}$ and we pose $x_0 = (\frac{Gm\Omega_0}{c^3})^{2/3}$. With this procedure, once the flux is entirely expressed in terms of $x_{GW}$ (after required reexpansion at the 4PN order), we find that the constant $b_0$ cancels out as expected.

At this stage, it is appropriate to rename $x_{GW}$ as simply “$x$” in order to be consistent with the calculation of the binary’s equations of motion, and notably the
binary's energy E for circular orbits. Indeed the energy was expressed in Eq. (377) in terms of the orbital frequency x, which represents exactly the quantity which would be measured far away from the binary detector<fnref n="88" /> and actually agrees with the measurable GW frequency xGW obtained in (481). In particular we shall consistently use the balance equation (447), where both E(x) and F(x) are expressed in terms of the same frequency x, which is the orbital frequency measured far away in the left-hand side, and the GW frequency in the right-hand side.

At long last, adding all pieces together, i.e. Equation (457) plus contributions from tails, iterated tails and tails-of-memory in Sect. 2.4.5, we obtain the quasi-circular 4.5PN energy flux as (Blanchet et al. 2023c, b)—remind our rename of xGW as x:

$$ 
\begin{aligned}
\mathcal{F} = & \frac{32c^5}{5G} v^2 x^5 \left\{ 1 + \left( -\frac{1247}{336} - \frac{35}{12}v \right) x + 4\pi x^{3/2} \right. \\
& \qquad + \left( -\frac{44711}{9072} + \frac{9271}{504}v + \frac{65}{18}v^2 \right) x^2 + \left( -\frac{8191}{672} - \frac{583}{24}v \right) \pi x^{5/2} \\
& \qquad + \left[ \frac{6643739519}{69854400} + \frac{16}{3}\pi^2 - \frac{1712}{105}\gamma_E - \frac{856}{105}\ln(16x) + \left( -\frac{134543}{7776} + \frac{41}{48}\pi^2 \right)v \right. \\
& \qquad \left. - \frac{94403}{3024}v^2 - \frac{775}{324}v^3 \right] x^3 \\
& \qquad + \left( -\frac{16285}{504} + \frac{214745}{1728}v + \frac{193385}{3024}v^2 \right) \pi x^{7/2} \\
& \qquad + \left[ -\frac{323105549467}{3178375200} + \frac{232597}{4410}\gamma_E - \frac{1369}{126}\pi^2 + \frac{39931}{294}\ln 2 - \frac{47385}{1568}\ln 3 \right] \\
& \qquad + \frac{232597}{8820}\ln x \\
& \qquad + \left( -\frac{1452202403629}{1466942400} + \frac{41478}{245}\gamma_E - \frac{267127}{4608}\pi^2 + \frac{479062}{2205}\ln 2 + \frac{47385}{392}\ln 3 \right) v \\
& \qquad + \left( \frac{1607125}{6804} - \frac{3157}{384}\pi^2 \right) v^2 + \frac{6875}{504}v^3 + \frac{5}{6}v^4 \\
& \qquad + \left[ \frac{265978667519}{745113600} - \frac{6848}{105}\gamma_E - \frac{3424}{105}\ln(16x) + \left( \frac{2062241}{22176} + \frac{41}{12}\pi^2 \right)v \right. \\
& \qquad \left. - \frac{133112905}{290304}v^2 - \frac{3719141}{38016}v^3 \right] \pi x^{9/2} + \mathcal{O}\left( \frac{1}{c^{10}} \right)
\end{aligned}
\quad (483) 
$$

A technical check is that all arbitrary constants ($r_0$, $r'_0$ and $b_0$) which played a role in intermediate steps of the calculation finally cancel out. An important physical check

<footnote n="88">⁸⁸ For instance by analyzing EM radio pulses sent out to infinity by an observer equipped with a flashlight and sitting on one of the binary's particles (Detweiler 2005).</footnote>
is that the test-mass limit $v \to 0$ perfectly agrees with the result of linear black-hole perturbation theory by Tagoshi and Nakamura (1994), Sasaki (1994), Tagoshi and Sasaki (1994), Tanaka et al. (1996), Fujita (2012a, b). Since BH perturbations have recently been extended numerically to second order in the mass ratio $v$ (Wardell et al. 2023; Albertini et al. 2022a, b), it would be interesting to verify the consistency of this numerical result with the PN prediction (483). Note also that in the case of black holes, the contributions due to the absorption by the BH horizons are not included in the PN calculation, and should be added separately. The BH absorption is a 4PN effect for Schwarzschild black holes (Poisson and Sasaki 1995), and a 2.5PN effect for spinning ones (Tagoshi et al. 1997; Alvi 2001; Porto 2008; Chatzioioannou et al. 2013; Saketh et al. 2023). In the latter case, the flux has to be augmented by direct spin contributions, see Sect. 3.6.4. Note also that the flux has been confirmed by different groups up to 2PN order (Will and Wiseman 1996; Leibovich et al. 2020), and that the 4.5PN piece is in agreement with independent work by Messina and Nagar (2017).

### 3.4.4 The 4.5PN phase and frequency evolution

We shall now deduce the laws of variation with time of the orbital frequency and phase of an inspiralling compact binary from the energy balance equation (447). The center-of-mass energy $E$ is given by Eq. (377) and the total flux $\mathcal{F}$ by Eq. (483). For convenience we adopt the dimensionless time variable

$$ \tau \equiv \frac{v c^3}{5 G m} (t_0 - t), \qquad (484) $$

where $t$ is the coordinate time in the asymptotic radiative coordinate system [formerly denoted $U$, see Eq. (178)], and $t_0$ an integration constant. We transform the balance equation into an ordinary differential equation for the parameter $x$, which is immediately integrated. Then we use the freedom to redefine the integration constant $t_0$ as $t_0 + \zeta \frac{Gm}{c^3}$ where $\zeta$ is any constant, which amounts to the replacement of $\tau$ by $\tau[1 + \zeta \frac{v}{5\tau}]$. Although $t_0$ is not uniquely defined, it might be formally interpreted as the instant of coalescence, when $x \to +\infty$, and then it satisfies $t_0 - t = \mathcal{O}(c^5)$ in the PN regime. Thus we see that $\tau^{-1} = \mathcal{O}(c^{-8})$ is a small 4PN quantity. Using this fact, we conveniently adjust $\zeta$ so as to simplify as much as possible the result:
$$
\begin{equation}
\begin{aligned}
x = \frac{\tau^{-1/4}}{4} & \left\{ 1 + \left( \frac{743}{4032} + \frac{11}{48} v \right) \tau^{-1/4} - \frac{1}{5} \pi \tau^{-3/8} \right. \\
& + \left( \frac{19583}{254016} + \frac{24401}{193536} v + \frac{31}{288} v^2 \right) \tau^{-1/2} + \left( -\frac{11891}{53760} + \frac{109}{1920} v \right) \pi \tau^{-5/8} \\
& + \left[ -\frac{10052469856691}{6008596070400} + \frac{1}{6} \pi^2 + \frac{107}{420} \gamma_E - \frac{107}{3360} \ln\left(\frac{\tau}{256}\right) \right. \\
& \qquad \left. + \left( \frac{3147553127}{780337152} - \frac{451}{3072} v^2 \right) v - \frac{15211}{442368} v^2 + \frac{25565}{331776} v^3 \right] \tau^{-3/4} \\
& + \left( -\frac{113868647}{433520640} - \frac{31821}{143360} v + \frac{294941}{3870720} v^2 \right) \pi \tau^{-7/8} \\
& + \left[ -\frac{2518977598355703073}{3779358859513036800} + \frac{9203}{215040} \gamma_E + \frac{9049}{258048} \pi^2 \right. \\
& \qquad \left. + \frac{14873}{1128960} \ln 2 + \frac{47385}{1605632} \ln 3 - \frac{9203}{3440640} \ln \tau \right. \\
& \qquad \left. + \left( \frac{718143266031997}{576825222758400} + \frac{244493}{1128960} \gamma_E - \frac{65577}{1835008} \pi^2 \right. \right. \\
& \qquad \left. + \frac{15761}{47040} \ln 2 - \frac{47385}{401408} \ln 3 - \frac{244493}{18063360} \ln \tau \right) v \\
& + \left( -\frac{1502014727}{8323596288} + \frac{2255}{393216} \pi^2 \right) v^2 - \frac{258479}{33030144} v^3 + \frac{1195}{262144} v^4 \right] \tau^{-1} \ln \tau \\
& + \left[ -\frac{9965202491753717}{5768252227584000} + \frac{107}{600} \gamma_E + \frac{23}{600} \pi^2 - \frac{107}{4800} \ln\left(\frac{\tau}{256}\right) \right. \\
& \qquad \left. + \left( \frac{8248609881163}{2746786775040} - \frac{3157}{30720} \pi^2 \right) v - \frac{3590973803}{20808990720} v^2 \right. \\
& \qquad \left. - \frac{520159}{1634992128} v^3 \right] \pi \tau^{-9/8} + \mathcal{O}\left(\frac{1}{c^{10}}\right)
\end{aligned}
\tag{485}
\end{equation}
$$

We have $\dot{\psi} = \Omega$, where we recall that $\Omega \equiv \Omega_{GW}$ is the wave frequency, from which we determine the explicit expression of the time-domain GW half-phase $\psi$ in terms of $x = (Gm\Omega/c^3)^{2/3}$ as (Blanchet et al. 2023c, b)<fnref n="89" />

<footnote n="89">⁸⁹ This procedure for computing analytically the orbital phase corresponds to what is called the “Taylor T2 approximant”. We refer to Boyle et al. (2008), Buonanno et al. (2009a) for the definition and usefulness of several types of approximants for computing (in general numerically) the orbital phase.</footnote>
$$
\begin{equation}
\begin{aligned}
\psi = \psi_0 & - \frac{x^{-5/2}}{32v} \left\{ 1 + \left( \frac{3715}{1008} + \frac{55}{12}v \right) x - 10\pi x^{3/2} \right. \\
& + \left( \frac{15293365}{1016064} + \frac{27145}{1008}v + \frac{3085}{144}v^2 \right) x^2 + \left( \frac{38645}{1344} - \frac{65}{16}v \right) \pi x^{5/2} \ln x \\
& + \left[ \frac{12348611926451}{18776862720} - \frac{160}{3}\pi^2 - \frac{1712}{21}\gamma_E - \frac{856}{21}\ln(16x) \right. \\
& \qquad \left. + \left( -\frac{15737765635}{12192768} + \frac{2255}{48}\pi^2 \right) v + \frac{76055}{6912}v^2 - \frac{127825}{5184}v^3 \right] x^3 \\
& + \left( \frac{77096675}{2032128} + \frac{378515}{12096}v - \frac{74045}{6048}v^2 \right) \pi x^{7/2} \\
& + \left[ \frac{2550713843998885153}{2214468081745920} - \frac{9203}{126}\gamma_E - \frac{45245}{756}\pi^2 - \frac{252755}{2646}\ln 2 \right. \\
& \qquad \left. - \frac{78975}{1568}\ln 3 - \frac{9203}{252}\ln x \right. \\
& \qquad \left. + \left( -\frac{680712846248317}{337983528960} - \frac{488986}{1323}\gamma_E + \frac{109295}{1792}\pi^2 - \frac{1245514}{1323}\ln 2 \right. \right. \\
& \qquad \left. + \frac{78975}{392}\ln 3 - \frac{244493}{1323}\ln x \right) v \\
& + \left( \frac{7510073635}{24385536} - \frac{11275}{1152}\pi^2 \right) v^2 + \frac{1292395}{96768}v^3 - \frac{5975}{768}v^4 \right] x^4 \\
& + \left[ -\frac{93098188434443}{150214901760} + \frac{1712}{21}\gamma_E + \frac{80}{3}\pi^2 + \frac{856}{21}\ln(16x) \right. \\
& \qquad \left. + \left( \frac{1492917260735}{1072963584} - \frac{2255}{48}\pi^2 \right) v - \frac{45293335}{1016064}v^2 - \frac{10323755}{1596672}v^3 \right] \pi x^{9/2} \\
& + \left. \mathcal{O}\left(\frac{1}{c^{10}}\right) \right\},
\end{aligned}
\tag{486}
\end{equation}
$$

where the integration constant $\psi_0$ is determined by initial conditions, e.g. when the
wave frequency enters the detector’s band. The results (485)–(486) give the unique
prediction of Einstein’s general relativity for the GW frequency and phase chirp of
non-spinning circular compact binaries up to 4.5PN precision. Note that the con-
tributions of the quadrupole moments of compact objects which are induced by tidal
effects, are expected to come into play only at the 5PN order, see Sect. 1.1.3. On the
other hand, most of the frameworks for data analysis rely on the stationary phase
approximation (Tichy et al. 2000), for which the phase of the dominant quadrupole
mode reads
$$
\begin{aligned}
\psi_{\text{SPA}} = 2\pi F T_0 + \Psi_0 &+ \frac{3 v^{-5}}{128 v} \left\{ 1 + \left( \frac{3715}{756} + \frac{55}{9} v \right) v^2 - 16\pi v^3 \right. \\
&+ \left( \frac{15293365}{508032} + \frac{27145}{504} v + \frac{3085}{72} v^2 \right) v^4 + \left( \frac{38645}{252} - \frac{65}{3} v \right) \pi v^5 \ln v \\
&+ \left[ \frac{11583231236531}{4694215680} - \frac{640}{3} \pi^2 - \frac{6848}{21} \gamma_E - \frac{6848}{21} \ln(4v) \right. \\
&\left. + \left( -\frac{15737765635}{3048192} + \frac{2255}{12} \pi^2 \right) v + \frac{76055}{1728} v^2 - \frac{127825}{1296} v^3 \right] v^6 \\
&+ \left[ \frac{77096675}{254016} + \frac{378515}{1512} v - \frac{74045}{756} v^2 \right] \pi v^7 \\
&+ \left[ -\frac{2550713843998885153}{276808510218240} + \frac{90490}{189} \pi^2 + \frac{36812}{63} \gamma_E + \frac{1011020}{1323} \ln 2 \right. \\
&\left. + \frac{78975}{196} \ln 3 + \frac{18406}{63} \ln v \right. \\
&+ \left( \frac{680712846248317}{42247941120} - \frac{109295}{224} \pi^2 + \frac{3911888}{1323} \gamma_E + \frac{9964112}{1323} \ln 2 \right. \\
&\left. - \frac{78975}{49} \ln 3 + \frac{1955944}{1323} \ln v \right) v \\
&+ \left( -\frac{7510073635}{3048192} + \frac{11275}{144} \pi^2 \right) v^2 - \frac{1292395}{12096} v^3 + \frac{5975}{96} v^4 \right] v^8 \ln v \\
&+ \left[ \frac{105344279473163}{18776862720} - \frac{640}{3} \pi^2 - \frac{13696}{21} \gamma_E - \frac{13696}{21} \ln(4v) \right. \\
&\left. + \left( -\frac{1492917260735}{134120448} + \frac{2255}{6} \pi^2 \right) v + \frac{45293335}{127008} v^2 \right. \\
&\left. + \left[ \frac{10323755}{199584} v^3 \right] \pi v^9 + \mathcal{O}\left( \left( \frac{1}{c^{10}} \right) \right) \right\},
\end{aligned}
\quad (487)
$$

where $v \equiv (\pi GmF/c^3)^{1/3}$ with $F$ being the Fourier frequency, and where $T_0$ and $\Psi_0$ are two integration constants. Again we have adjusted $T_0$ in order to simplify the result. See e.g. Pratten et al. (2020) for application to gravitational-wave data analysis.

As a rough estimate of the relative importance of the various PN terms, we provide in Table 2 numerical estimates for the number of accumulated GW cycles in the frequency band of current and future detectors. Such an estimate is only indicative, because a full treatment would require the knowledge of the detector's noise, and a complete simulation of the parameter estimation using matched filtering techniques (Cutler and Flanagan 1994; Poisson and Will 1995; Królak et al. 1995;
<table><caption>**Table 2** Contribution of each PN order to the total number of accumulated cycles (488) inside the detector's frequency band, for typical (non-spinning) quasi-circular compact binaries observed by current and future detectors</caption><thead><tr><th rowspan="2">Detector<br>Masses (M<sub>☉</sub>)</th><th colspan="2">LIGO/Virgo</th><th colspan="2">ET</th><th colspan="2">LISA</th></tr><tr><th>1.4 × 1.4</th><th>10 × 10</th><th>1.4 × 1.4</th><th>500 × 500</th><th>10<sup>5</sup> × 10<sup>5</sup></th><th>10<sup>7</sup> × 10<sup>7</sup></th></tr></thead><tbody><tr><td>PN order</td><td colspan="6">Cumulative number of cycles N<sub>cycle</sub></td></tr><tr><td>Newtonian</td><td>2 562.599</td><td>95.502</td><td>744 401.36</td><td>37.90</td><td>28 095.39</td><td>9.534</td></tr><tr><td>1PN</td><td>143.453</td><td>17.879</td><td>4 433.85</td><td>9.60</td><td>618.31</td><td>3.386</td></tr><tr><td>1.5PN</td><td>-94.817</td><td>-20.797</td><td>-1 005.78</td><td>-12.63</td><td>-265.70</td><td>-5.181</td></tr><tr><td>2PN</td><td>5.811</td><td>2.124</td><td>23.94</td><td>1.44</td><td>11.35</td><td>0.677</td></tr><tr><td>2.5PN</td><td>-8.105</td><td>-4.604</td><td>-17.01</td><td>-3.42</td><td>-12.47</td><td>-1.821</td></tr><tr><td>3PN</td><td>1.858</td><td>1.731</td><td>2.69</td><td>1.43</td><td>2.59</td><td>0.876</td></tr><tr><td>3.5PN</td><td>-0.627</td><td>-0.689</td><td>-0.93</td><td>-0.59</td><td>-0.91</td><td>-0.383</td></tr><tr><td>4PN</td><td>-0.107</td><td>-0.064</td><td>-0.12</td><td>-0.04</td><td>-0.12</td><td>-0.013</td></tr><tr><td>4.5PN</td><td>0.098</td><td>0.118</td><td>0.14</td><td>0.10</td><td>0.14</td><td>0.065</td></tr></tbody></table>

We have approximated the frequency bands of LIGO/Virgo, Einstein Telescope (ET) and LISA with step functions, respectively between [$f_{entry}, f_{exit}$] = [30 Hz, $10^3$ Hz], [1 Hz, $10^4$ Hz] and [$10^{-4}$ Hz, $10^{-1}$ Hz]. When the merger occurs within the frequency band of the detector, the exit frequency is taken to be the Schwarzschild ISCO, $f_{exit} = c^3/(6^{3/2}\pi Gm)$

Arun et al. 2005). The number of cycles between entry and exit frequencies $f_{entry}$ and $f_{exit}$ (where $f \equiv \Omega/\pi$) is

$$N_{\text{cycle}} \equiv \frac{\psi_{\text{exit}} - \psi_{\text{entry}}}{\pi}. \qquad (488)$$

As we see from Table 2, the PN expansion seems to converge well. For all the typical compact binaries, we find that the 4PN and 4.5PN orders amount to about a tenth of a cycle (less than 1 radian). This suggests that systematic errors due to the PN modeling may be dominated by statistical errors and negligible for LISA. However, this should be confirmed by detailed investigations along the lines of Owen et al. (2023).

### 3.4.5 Spherical modes for data analysis and numerical relativity

Besides the chirp described by Eqs. (485)–(486), it is also important to compute the wave amplitude, in view of the data analysis of LISA (Hellings and Moore 2003; Moore and Hellings 2002). For instance it has been shown that using the full waveform information for LISA will yield substantial improvements (with respect to the “restricted post-Newtonian waveform”) of the angular resolution and the estimation of the luminosity distance of super-massive black hole binaries (Arun et al. 2007a, b; Trias and Sintes 2008). We also require the wave amplitude for high-accuracy comparisons with numerical relativity (see e.g. Boyle et al. 2007; Bernuzzi et al. 2012; Borhanian et al. 2020).
The two polarization waveforms defined by Eq. (183)<fnref n="90" /> are decomposed, at leading order in the distance $R$ to the binary (in radiative coordinates), onto the basis of spin-weighted spherical harmonics following Eq. (185). The spherical modes are the coefficients in this decomposition and can be evaluated either from the angular integration (186), or from the relations (187)–(188) giving them directly in terms of the radiative moments $U_L$ and $V_L$.

A useful fact to remember is that for non-spinning binaries, the mode $h^{lm}$ is entirely given by the mass multipole moment $U_L$ when $l+m$ is even, and by the current one $V_L$ when $l+m$ is odd. This is valid in general for non-spinning binaries, regardless of the orbit being quasi-circular or elliptical. The important point is only that the motion of the two particles must be planar, i.e., takes place in a fixed plane. This is the case if the particles are non-spinning, but this will also be the case if, more generally, the spins are aligned or anti-aligned with the orbital angular momentum, since there is no orbital precession in this case. Thus, for any “planar” binaries, Eq. (187) splits into (see Faye et al. 2012 for a proof)

$$h^{lm} = -\frac{G}{\sqrt{2}Rc^{l+2}} U^{lm} \quad (\text{when } l+m \text{ is even}), \quad (489a)$$

$$h^{lm} = \frac{G}{\sqrt{2}Rc^{l+3}} i V^{lm} \quad (\text{when } l+m \text{ is odd}). \quad (489b)$$

Let us factorize out in all the modes an overall coefficient including the appropriate phase factor $e^{-im\psi}$, where we recall that $\psi$ denotes the gravitational-wave half phase (486), and such that the dominant mode with $(l, m) = (2, 2)$ conventionally starts with one at the Newtonian order. We thus pose

$$h^{lm} = \frac{2G\mu x}{Rc^2} \sqrt{\frac{16\pi}{5}} H^{lm} e^{-im\psi}. \quad (490)$$

The dominant mode is $H^{22}$, due to the mass quadrupole radiation, and is primarily important for numerical relativity comparisons. This mode has been computed to 4PN order by Blanchet et al. (2023b, c), extending previous calculations to 3PN order by Kidder (2008), Blanchet et al. (2008) and to 3.5PN by Faye et al. (2012). We have

<footnote n="90">⁹⁰ For the choice of polarization vectors we adopt the following: $P$ and $Q$ lie along the major and minor axis, respectively, of the projection onto the plane of the sky of the binary’s orbit. Thus $P$ is oriented toward the orbit’s *ascending node*—namely the point $N$ at which the orbit intersects the plane of the sky and the bodies are moving *toward* the observer located in the direction $N$. The ascending node for the particle 1 is chosen for the origin of the phase $\psi$. The inclination angle $i$ is between the direction of the detector $N$ as seen from the binary’s center-of-mass, and the normal to the orbital plane (supposed to be right-handed with respect to the sense of motion, so that $0 \le i \le \pi$).</footnote>
$$
\begin{equation}
\begin{aligned}
H^{22} = & 1 + x \left(-\frac{107}{42} + \frac{55}{42}v\right) + 2\pi x^{3/2} \\
& + x^2 \left(-\frac{2173}{1512} - \frac{1069}{216}v + \frac{2047}{1512}v^2\right) + x^{5/2} \left[-\frac{107\pi}{21} + \left(\frac{34\pi}{21} - 24i\right)v\right] \\
& + x^3 \left[\frac{27027409}{646800} - \frac{856}{105}\gamma_E + \frac{428i\pi}{105} + \frac{2\pi^2}{3} + \left(-\frac{278185}{33264} + \frac{41\pi^2}{96}\right)v\right. \\
& \left. - \frac{20261}{2772}v^2 + \frac{114635}{99792}v^3 - \frac{428}{105}\ln(16x)\right] \\
& + x^{7/2} \left[-\frac{2173\pi}{756} + \left(-\frac{2495\pi}{378} + \frac{14333i}{162}\right)v + \left(\frac{40\pi}{27} - \frac{4066i}{945}\right)v^2\right] \\
& + x^4 \left[-\frac{846557506853}{12713500800} + \frac{45796}{2205}\gamma_E - \frac{22898}{2205}i\pi - \frac{107}{63}\pi^2 + \frac{22898}{2205}\ln(16x)\right. \\
& \left. + \left(-\frac{336005827477}{4237833600} + \frac{15284}{441}\gamma_E - \frac{219314}{2205}i\pi - \frac{9755}{32256}\pi^2 + \frac{7642}{441}\ln(16x)\right)v\right. \\
& \left. + \left(\frac{256450291}{7413120} - \frac{1025}{1008}\pi^2\right)v^2 - \frac{81579187}{15567552}v^3 + \frac{26251249}{31135104}v^4\right] \\
& + O\left(\frac{1}{c^9}\right).
\end{aligned}
\tag{491}
\end{equation}
$$

This result is in perfect agreement with linear black-hole perturbation theory in the limit when $v \to 0$ (Tagoshi and Sasaki 1994).

We notice that the mode (491) contains some imaginary parts. These can clearly be viewed as corrections to the phase rather than corrections in the amplitude. However, as corrections to the phase they are to be compared to the leading phase evolution $\psi(x)$ given by (486). Thus we see that their real order of magnitude is a factor $x^{5/2} = O(c^{-5})$ smaller than what is read from the amplitude. Since the imaginary parts in (491) start at 2.5PN order, they correspond to very small 5PN modulations in the phase, and therefore can be neglected (the same conclusion also applies to the other modes below).

The mass octupole modes $H^{33}$ and $H^{31}$ were computed to 3.5PN order by Faye et al. (2015), and the current quadrupole mode $H^{21}$ to 3.5PN order by Henry et al. (2021). Note that 3.5PN refers here to the precision of the waveform; for that one needs to control the mass octupole and current quadrupole radiative moments ($U_{ijk}$ and $V_{ij}$) to 3PN order. Therefore these modes are comparatively easier than the mass quadrupole mode (491) which required the full 4PN precision in the mass quadrupole moment $U_{ij}$. All the modes have been completed up to the 3.5PN level in the waveform by Henry (2023)—the overlap agrees with previous calculations (Blanchet et al. 2008). We have$^{91}$

<footnote>$^{91}$ We assume $m \ge 0$; the modes having $m<0$ are deduced using $H^{l,-m} = (-)^l \tilde{H}^{lm}$.</footnote>
$$
\begin{equation}
\begin{aligned}
H^{21} = \frac{1}{3} i \Delta & \left[ x^{1/2} + x^{3/2} \left( -\frac{17}{28} + \frac{5v}{7} \right) + x^2 \left( \pi + i \left( -\frac{1}{2} - 2 \ln 2 \right) \right) \right. \\
& + x^{5/2} \left( -\frac{43}{126} - \frac{509v}{126} + \frac{79v^2}{168} \right) + x^3 \left( -\frac{17\pi}{28} + \frac{3\pi v}{14} \right. \\
& \left. + i \left( \frac{17}{56} + v \left( -\frac{353}{28} - \frac{3\ln 2}{7} \right) + \frac{17\ln 2}{14} \right) \right) \\
& + x^{7/2} \left( \frac{15223771}{1455300} + \frac{\pi^2}{6} - \frac{214}{105} \gamma_E - \frac{107}{105} \ln(4x) - \ln 2 - 2(\ln 2)^2 \right. \\
& \left. + v \left( -\frac{102119}{2376} + \frac{205}{128} \pi^2 \right) - \frac{4211}{8316} v^2 + \frac{2263}{8316} v^3 + i \pi \left( \frac{109}{210} - 2 \ln 2 \right) \right] \\
& + \mathcal{O} \left( \frac{1}{c^8} \right),
\end{aligned}
\tag{492}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
H^{33} ={}& -\frac{3}{4} i \sqrt{\frac{15}{14}} \Delta \left[ x^{1/2} + x^{3/2}(-4+2v) + x^2 \left( 3\pi + i \left( -\frac{21}{5} + 6\ln(3/2) \right) \right) \right. \\
& \left. \qquad + x^{5/2} \left( \frac{123}{110} - \frac{1838v}{165} + \frac{887v^2}{330} \right) + x^3 \left( -12\pi + \frac{9\pi v}{2} \right. \right. \\
& \left. \left. \qquad + i \left( \frac{84}{5} - 24\ln(3/2) + v \left( -\frac{48103}{1215} + 9\ln(3/2) \right) \right) \right) \right] \\
& + x^{7/2} \left( \frac{19388147}{280280} + \frac{492}{35} \ln(3/2) - 18\ln^2(3/2) - \frac{78}{7}\gamma_E + \frac{3}{2}\pi^2 \right. \\
& \left. \qquad + 6i\pi \left( -\frac{41}{35} + 3\ln(3/2) \right) \right. \\
& \left. \qquad + \frac{v}{8} \left( -\frac{7055}{429} + \frac{41}{8}\pi^2 \right) - \frac{318841}{17160} v^2 + \frac{8237}{2860} v^3 - \frac{39}{7}\ln(16x) \right) \right] \\
& + \mathcal{O}\left(\frac{1}{c^8}\right),
\end{aligned}
\tag{493}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
H^{32} ={}& \frac{1}{3} \sqrt{\frac{5}{7}} \left[ x(1-3v) + x^2 \left( -\frac{193}{90} + \frac{145v}{18} - \frac{73v^2}{18} \right) \right. \\
& \left. \qquad + x^{5/2} \left( 2\pi - 6\pi v + i \left( -3 + \frac{66v}{5} \right) \right) \right. \\
& \left. \qquad + x^3 \left( -\frac{1451}{3960} - \frac{17387v}{3960} + \frac{5557v^2}{220} - \frac{5341v^3}{1320} \right) + x^{7/2} \left( \frac{193}{30}i - \frac{193}{45}\pi \right) \right. \\
& \left. \qquad + v \left( -\frac{258929}{5400}i + \frac{136}{9}\pi \right) + v^2 \left( \frac{33751}{450}i - \frac{46}{9}\pi \right) \right] \\
& + \mathcal{O}\left(\frac{1}{c^8}\right),
\end{aligned}
\tag{494}
\end{equation}
$$
$$
\begin{aligned}
H^{31} ={}& \frac{i \Delta}{12\sqrt{14}} \left[ x^{1/2} + x^{3/2} \left(-\frac{8}{3} - \frac{2v}{3}\right) + x^2 \left(\pi + i \left(-\frac{7}{5} - 2\ln 2\right)\right) \right. \\
& \left. \qquad + x^{5/2} \left(\frac{607}{198} - \frac{136v}{99} - \frac{247v^2}{198}\right) + x^3 \left(-\frac{8\pi}{3} - \frac{7\pi v}{6}\right) \right. \\
& \left. \qquad + i \left(\frac{56}{15} + \frac{16 \ln 2}{3} + v \left(-\frac{1}{15} + \frac{7 \ln 2}{3}\right)\right) \right] \\
& + x^{7/2} \left(\frac{10753397}{1513512} - 2 \ln 2 \left(\frac{212}{105} + \ln 2\right) \right. \\
& \left. \qquad - \frac{26}{21} \gamma_E + \frac{\pi^2}{6} - 2i\pi \left(\frac{41}{105} + \ln 2\right) \right. \\
& \left. \qquad + \frac{v}{8} \left(-\frac{1738843}{19305} + \frac{41}{8}\pi^2\right) + \frac{327059}{30888} v^2 - \frac{17525}{15444} v^3 - \frac{13}{21} \ln x \right] \\
& + \mathcal{O}\left(\frac{1}{c^8}\right),
\end{aligned}
\qquad (495)
$$

$$
H^{30} = -\frac{2}{5} i \sqrt{\frac{6}{7}} v \left[ x^{5/2} + x^{7/2} \left( -\frac{5017}{1296} - \frac{25}{108} v \right) \right] + \mathcal{O}\left( \frac{1}{c^8} \right), \quad (496)
$$

$$
\begin{aligned}
H^{44} ={}& -\frac{8}{9} \sqrt{\frac{5}{7}} \left[ x(1-3v) + x^2 \left( -\frac{593}{110} + \frac{1273v}{66} - \frac{175v^2}{22} \right) \right. \\
& \left. \qquad + x^{5/2} \left( 4\pi - 12\pi v + i \left( -\frac{42}{5} + v \left( \frac{1193}{40} - 24\ln 2 \right) + 8\ln 2 \right) \right) \right. \\
& \left. \qquad + x^3 \left( \frac{1068671}{200200} - \frac{1088119v}{28600} + \frac{146879v^2}{2340} - \frac{226097v^3}{17160} \right) \right. \\
& \left. \qquad + x^{7/2} \left( \frac{12453}{275}i - \frac{1186}{55}\pi - \frac{2372}{55}i\ln 2 \right. \right. \\
& \left. \qquad \left. + v \left( -\frac{31525499}{140800}i + \frac{2480}{33}\pi + \frac{4960}{33}i\ln 2 \right) \right. \right. \\
& \left. \qquad \left. + v^2 \left( \frac{4096237}{21120}i - \frac{284}{11}\pi - \frac{568}{11}i\ln 2 \right) \right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right),
\end{aligned}
\qquad (497)
$$

$$
\begin{aligned}
H^{43} ={}& -\frac{9i \Delta}{4\sqrt{70}} \left[ x^{3/2}(1-2v) + x^{5/2} \left( -\frac{39}{11} + \frac{1267v}{132} - \frac{131v^2}{33} \right) \right. \\
& \left. \qquad + x^3 \left( 3\pi - 6\pi v + i \left( -\frac{32}{5} + v \left( \frac{16301}{810} - 12\ln(3/2) \right) + 6\ln(3/2) \right) \right) \right. \\
& \left. \qquad \left. + x^{7/2} \left( \frac{7206}{5005} - \frac{82869}{5720}v + \frac{104839}{3432}v^2 - \frac{2987}{572}v^3 \right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right),
\end{aligned}
\qquad (498)
$$
$$
\begin{equation}
\begin{aligned}
H^{42} ={}& \frac{1}{63} \sqrt{5} \left[ x(1-3v) + x^2 \left( -\frac{437}{110} + \frac{805v}{66} - \frac{19v^2}{22} \right) \right. \\
& \left. \qquad + x^{5/2} \left( 2\pi - 6\pi v + i \left( -\frac{21}{5} + \frac{84v}{5} \right) \right) \right. \\
& \left. \qquad + x^3 \left( \frac{1038039}{200200} - \frac{606751v}{28600} + \frac{400453v^2}{25740} + \frac{25783v^3}{17160} \right) \right. \\
& \left. \qquad + x^{7/2} \left( \frac{9177}{550} i - \frac{437}{55} \pi + v \left( -\frac{83029}{880} i + \frac{772}{33} \pi \right) + v^2 \left( \frac{93081}{1100} i + \frac{14}{11} \pi \right) \right) \right] \\
& \left. \qquad + \mathcal{O}\left( \frac{1}{c^8} \right),
\end{aligned}
\tag{499}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
H^{41} ={}& \frac{i \Delta}{84\sqrt{10}} \left[ x^{3/2}(1-2v) + x^{5/2} \left( -\frac{101}{33} + \frac{337v}{44} - \frac{83v^2}{33} \right) \right. \\
& \left. \qquad + x^3 \left( \pi - 2\pi v + i \left( -\frac{32}{15} - 2\ln 2 + v \left( \frac{1661}{30} + 4\ln 2 \right) \right) \right) \right] \\
& \left. \qquad + x^{7/2} \left( \frac{42982}{15015} - \frac{513989}{51480} v + \frac{196957}{10296} v^2 - \frac{1195}{572} v^3 \right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right),
\end{aligned}
\tag{500}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
H^{55} ={}& \frac{625i \Delta}{96\sqrt{66}} \left[ x^{3/2}(1-2v) + x^{5/2} \left( -\frac{263}{39} + \frac{688v}{39} - \frac{256v^2}{39} \right) \right. \\
& \left. \qquad + x^3 \left( 5\pi - 10\pi v + i \left( -\frac{181}{14} + v \left( \frac{105834}{3125} - 20\ln(5/2) \right) + 10\ln(5/2) \right) \right) \right. \\
& \left. \qquad + x^{7/2} \left( \frac{9185}{819} - \frac{188765}{3276} v + \frac{54428}{819} v^2 - \frac{10567}{819} v^3 \right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right),
\end{aligned}
\tag{501}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
H^{54} ={}& -\frac{32}{9\sqrt{165}} \left[ x^2(1-5v+5v^2) + x^3 \left( -\frac{4451}{910} + \frac{3619v}{130} - \frac{521v^2}{13} + \frac{339v^3}{26} \right) \right. \\
& \left. \qquad + x^{7/2} \left( -\frac{52}{5}i + 4\pi + 8i \ln 2 + v \left( \frac{3351011}{53760}i - 20\pi - 40i \ln 2 \right) \right. \right. \\
& \left. \left. \qquad + v^2 \left( -\frac{10923}{128}i + 20\pi + 40i \ln 2 \right) \right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right),
\end{aligned}
\tag{502}
\end{equation}
$$
$$H^{53} = -\frac{9}{32}i\sqrt{\frac{3}{110}}\Delta \left[ x^{3/2}(1-2v) + x^{5/2}\left(-\frac{69}{13} + \frac{464v}{39} - \frac{88v^2}{39}\right) \right. \\ \left. \qquad + x^3\left(3\pi - 6\pi v + i\left(-\frac{543}{70} + v\left(\frac{83702}{3645} - 12\ln(3/2)\right) + 6\ln(3/2)\right)\right) \right. \\ \left. \qquad + x^{7/2}\left(\frac{12463}{1365} - \frac{56969}{1820}v + \frac{2172}{91}v^2 - \frac{365}{273}v^3\right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \qquad (503)$$

$$H^{52} = \frac{2}{27\sqrt{55}} \left[ x^2(1 - 5v + 5v^2) + x^3 \left( -\frac{3911}{910} + \frac{3079v}{130} - \frac{413v^2}{13} + \frac{231v^3}{26} \right) \right. \\ \left. \qquad + x^{7/2} \left( -\frac{26}{5}i + 2\pi + v \left( \frac{16237}{336}i - 10\pi \right) + v^2 \left( -\frac{1861}{20}i + 10\pi \right) \right) \right] \\ \left. \qquad + \mathcal{O}\left(\frac{1}{c^8}\right), \right] \qquad (504)$$

$$H^{51} = \frac{i \Delta}{288\sqrt{385}} \left[ x^{3/2}(1-2v) + x^{5/2}\left(-\frac{179}{39} + \frac{352v}{39} - \frac{4v^2}{39}\right) \right. \\ \left. \qquad + x^3\left(\pi - 2\pi v + i\left(-\frac{181}{70} - 2\ln 2 + v\left(\frac{626}{5} + 4\ln 2\right)\right)\right) \right. \\ \left. \qquad + x^{7/2}\left(\frac{5023}{585} - \frac{49447}{2340}v + \frac{68}{9}v^2 + \frac{287}{117}v^3\right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \qquad (505)$$

$$H^{50} = \frac{4117iv}{7560\sqrt{462}} x^{7/2} \left[ 1 - \frac{21588}{4117} v \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \qquad (506)$$

$$H^{66} = \frac{54}{5\sqrt{143}} \left[ x^2(1 - 5v + 5v^2) + x^3 \left( -\frac{113}{14} + \frac{91v}{2} - 64v^2 + \frac{39v^3}{2} \right) \right. \\ \left. \qquad + x^{7/2} \left( -\frac{249}{14}i + 6\pi + 12i \ln 3 + v \left( \frac{21787499}{217728}i - 30\pi - 60i \ln 3 \right) \right. \right. \\ \left. \left. \qquad + v^2 \left( -\frac{323903}{2592}i + 30\pi + 60i \ln 3 \right) \right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \qquad (507)$$

$$H^{65} = \frac{3125i \Delta}{504\sqrt{429}} \left[ x^{5/2}(1 - 4v + 3v^2) \right. \\ \left. \qquad + x^{7/2} \left( -\frac{149}{24} + \frac{349}{12}v - \frac{409}{12}v^2 + \frac{29}{3}v^3 \right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \qquad (508)$$
$$H^{64} = -\frac{128}{495} \sqrt{\frac{2}{39}} \left[ \begin{aligned} & x^2(1 - 5v + 5v^2) + x^3 \left( -\frac{93}{14} + \frac{71v}{2} - 44v^2 + \frac{19v^3}{2} \right) \\ & + x^{7/2} \left( -\frac{83}{7}i + 4\pi + 8i \ln 2 + v \left( \frac{3432215}{49152}i - 20\pi - 40i \ln 2 \right) \right) \\ & + v^2 \left( -\frac{382365}{4096}i + 20\pi + 40i \ln 2 \right) \end{aligned} \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (509)$$

$$H^{63} = -\frac{81i\Delta}{616\sqrt{65}} \left[ \begin{aligned} & x^{5/2}(1 - 4v + 3v^2) \\ & + x^{7/2} \left( -\frac{133}{24} + \frac{301}{12}v - \frac{329}{12}v^2 + 7v^3 \right) \end{aligned} \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (510)$$

$$H^{62} = \frac{2}{297\sqrt{65}} \left[ \begin{aligned} & x^2(1 - 5v + 5v^2) + x^3 \left( -\frac{81}{14} + \frac{59v}{2} - 32v^2 + \frac{7v^3}{2} \right) \\ & + x^{7/2} \left( -\frac{83}{14}i + 2\pi + v \left( \frac{799789}{13440}i - 10\pi \right) + v^2 \left( -\frac{19193}{160}i + 10\pi \right) \right) \end{aligned} \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (511)$$

$$H^{61} = \frac{i\Delta}{8316\sqrt{26}} \left[ \begin{aligned} & x^{5/2}(1 - 4v + 3v^2) \\ & + x^{7/2} \left( -\frac{125}{24} + \frac{277}{12}v - \frac{289}{12}v^2 + \frac{17}{3}v^3 \right) \end{aligned} \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (512)$$

$$H^{77} = -\frac{16807i\Delta}{1440\sqrt{858}} \left[ \begin{aligned} & x^{5/2}(1 - 4v + 3v^2) \\ & + x^{7/2} \left( -\frac{319}{34} + \frac{2225}{51}v - \frac{2558}{51}v^2 + \frac{230}{17}v^3 \right) \end{aligned} \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (513)$$

$$H^{76} = \frac{81}{35} \sqrt{\frac{3}{143}} x^3 \left[ 1 - 7v + 14v^2 - 7v^3 \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (514)$$

$$H^{75} = \frac{15625i\Delta}{26208\sqrt{66}} \left[ \begin{aligned} & x^{5/2}(1 - 4v + 3v^2) \\ & + x^{7/2} \left( -\frac{271}{34} + \frac{1793}{51}v - \frac{1838}{51}v^2 + \frac{134}{17}v^3 \right) \end{aligned} \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (515)$$

$$H^{74} = -\frac{128x^3}{1365} \sqrt{\frac{2}{33}} \left[ 1 - 7v + 14v^2 - 7v^3 \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (516)$$
$$H^{73} = -\frac{243i\Delta}{160160}\sqrt{\frac{3}{2}}\left[x^{5/2}(1-4v+3v^2) + x^{7/2}\left(-\frac{239}{34} + \frac{1505}{51}v - \frac{1358}{51}v^2 + \frac{70}{17}v^3\right)\right] + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (517)$$

$$H^{72} = \frac{x^3}{3003\sqrt{3}}(1 - 7v + 14v^2 - 7v^3) + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (518)$$

$$H^{71} = \frac{i \Delta}{864864\sqrt{2}} \left[ x^{5/2}(1 - 4v + 3v^2) + x^{7/2}\left(-\frac{223}{34} + \frac{1361}{51}v - \frac{1118}{51}v^2 + \frac{38}{17}v^3\right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (519)$$

$$H^{70} = \mathcal{O}\left(\frac{1}{c^8}\right). \quad (520)$$

With the 3.5PN approximation in the waveform all the modes with $ℓ ≥ 8$ can be considered as merely Newtonian. We give here the general Newtonian expressions of any mode with arbitrary $ℓ$ and non-zero $m$ (Kidder 2008):

$$H^{\ell m}|_{\text{leading}} = \frac{(-)^{(\ell-m+2)/2}}{2^{\ell+1}(\ell+m)!(\ell-m)!(2\ell-1)!!} \left( \frac{5(\ell+1)(\ell+2)(\ell+m)!(\ell-m)!}{\ell(\ell-1)(2\ell+1)} \right)^{1/2} \times \sigma_{\ell}(v) (im)^{\ell} x^{\ell/2-1} \quad (\text{for } \ell+m \text{ even}), \quad (521a)$$

$$H^{\ell m}|_{\text{leading}} = \frac{(-)^{(\ell-m-1)/2}}{2^{\ell-1}(\ell+m-1)!(\ell-m-1)!(2\ell+1)!!} \left( \frac{5(\ell+2)(2\ell+1)(\ell+m)!(\ell-m)!}{\ell(\ell-1)(\ell+1)} \right)^{1/2} \times \sigma_{\ell+1}(v) (im)^{\ell} x^{\ell-1/2} \quad (\text{for } \ell+m \text{ odd}), \quad (521b)$$

where we employ the function of the mass ratio $\sigma_{\ell}(v) = X_2^{\ell-1} + (-)^{\ell}X_1^{\ell-1}$ (with $X_a = m_a/m$), which is also given by Eqs. (456).

We have not yet written the modes $(\ell, 0)$ with $ℓ$ even. Indeed these modes with $m = 0$ arise from the DC memory integrals in the mass-type radiative moments $U_L$; hence they have $ℓ + m$ even, see (489). From Eq. (474) we know that these integrals exhibit the memory effect, i.e. their formal post-Newtonian order is enhanced by the inverse of the adiabatic parameter $\xi = \mathcal{O}(c^{-5})$. Therefore, starting with the explicit waveform at the 4PN order we can control the DC memory contributions only with the 1.5PN precision; the result is

$$H^{20} = -\frac{5}{14\sqrt{6}}\left[1 + x\left(-\frac{4075}{4032} + \frac{67v}{48}\right)\right] + \mathcal{O}\left(\frac{1}{c^4}\right), \quad (522a)$$
$$H^{40} = -\frac{1}{504\sqrt{2}} \left[ 1 + x \left( -\frac{180101}{29568} + \frac{27227}{1056} v \right) \right] + O\left(\frac{1}{c^4}\right), \quad (522b)$$

$$H^{60} = \frac{4195}{1419264\sqrt{273}} x \left[ 1 - \frac{3612}{839} v \right] + O\left(\frac{1}{c^4}\right), \quad (522c)$$

$$H^{80} = O\left(\frac{1}{c^4}\right). \quad (522d)$$

However, the DC modes can also be computed in a general way, using the general formalism of Sect. 2.4.4, in terms of the radiative moments, see Eqs. (225) and (237). Using this formalism the DC memory terms have been computed to 3PN order by Favata (2009); of course, the 1.5PN results (522) agree.

## 3.5 Eccentric compact binaries

Inspiring compact binaries are usually modelled as moving in quasi-circular orbits since gravitational radiation reaction circularizes the orbit towards the late stages of inspiral (Peters 1964), as discussed in Sect. 1.1.2. Nevertheless, there is an increased interest in inspiralling binaries moving in quasi-eccentric orbits. Astrophysical scenarios currently exist which lead to binaries with non-zero eccentricity in the gravitational-wave detector bandwidth, both terrestrial and space-based. For instance, inner binaries of hierarchical triplets undergoing ZLK (von Zeipel 1910; Lidov 1962; Kozai 1962) oscillations could not only merge due to gravitational radiation reaction but a fraction of them should have non negligible eccentricities when they enter the sensitivity band of advanced ground based interferometers (Wen 2003). On the other hand the population of stellar mass binaries in globular clusters is expected to have a thermal distribution of eccentricities (Benacquista and Downing 2013). In a study of the growth of intermediate black holes (Gultekin et al. 2004) in globular clusters it was found that the binaries have eccentricities between 0.1 and 0.2 in the LISA bandwidth. Though, supermassive black hole binaries are powerful gravitational-wave sources for LISA, it is not known if they would be in quasi-circular or quasi-eccentric orbits. If a ZLK mechanism is at work, these supermassive black hole binaries could be in highly eccentric orbits and merge within the Hubble time (Blaes et al. 2002). Sources of the kind discussed above provide the prime motivation for investigating higher post-Newtonian order modelling for quasi-eccentric binaries.

### 3.5.1 Doubly periodic structure of the motion of eccentric binaries

In Sect. 3.2.4 we have investigated the equations of motion of non-spinning compact binary systems in the frame of the center-of-mass for general orbits. We shall now (in this section and the next one) look for the explicit solution to those equations. In particular, let us discuss the general “doubly-periodic” structure of the post-Newtonian solution, closely following Damour (1983a, b), Damour and Deruelle (1985).
The 4PN equations of motion admit, when neglecting the radiation reaction terms at 2.5PN and 3.5PN orders, ten first integrals of the motion corresponding to the conservation of energy, angular momentum, linear momentum, and center of mass position. When restricted to the frame of the center of mass, the equations admit four first integrals associated with the energy E and the angular momentum vector **J**.

The motion takes place in the plane orthogonal to **J**. Denoting by $r = |\mathbf{x}|$ the binary's orbital separation in that plane, and by $\mathbf{v} = \mathbf{v}_1 - \mathbf{v}_2$ the relative velocity, we find that E and **J** are functions of $r$, $\dot{r}^2$, $\mathbf{v}^2$ and $\mathbf{x} \times \mathbf{v}$. We adopt polar coordinates $(r, \phi)$ in the orbital plane, and express E and the norm $\mathbf{J} = |\mathbf{J}|$, thanks to $\mathbf{v}^2 = \dot{r}^2 + r^2 \dot{\phi}^2$, as some explicit functions of $r$, $\dot{r}^2$ and $\dot{\phi}$. The latter functions can be inverted by means of a straightforward post-Newtonian iteration to give $\dot{r}^2$ and $\dot{\phi}$ in terms of $r$ and the constants of motion E and J. Hence,

$$ \dot{r}^2 = \mathcal{R}[r; E, \mathbf{J}], \qquad (523a) $$

$$ \dot{\phi} = \mathcal{S}[r; E, \mathbf{J}], \qquad (523b) $$

where $\mathcal{R}$ and $\mathcal{S}$ denote certain polynomials in $1/r$, the degree of which depends on the post-Newtonian approximation in question; for instance it is seventh degree for both $\mathcal{R}$ and $\mathcal{S}$ at 3PN order (Memmesheimer et al. 2004). At high PN orders one must use coordinate systems which avoid the presence of the logarithm of $r$. The various coefficients of the powers of $1/r$ are themselves polynomials in E and J, and also, of course, depend on the total mass $m$ and symmetric mass ratio $v$. In the case of bounded elliptic-like motion, one can prove (Damour 1983b) that the function $\mathcal{R}$ admits two real roots, say $r_p$ and $r_a$ such that $r_p \le r_a$, which admit some non-zero finite Newtonian limits when $c \to \infty$, and represent respectively the radii of the orbit's periastron ($p$) and apastron ($a$). The other roots are complex and tend to zero when $c \to \infty$.

Let us consider a given binary's orbital configuration, fully specified by some values of the integrals of motion E and J corresponding to quasi-elliptic motion.<fnref n="92" /> The binary's orbital period, or time of return to the periastron, is obtained by integrating the radial motion as

$$ P = 2 \int_{r_p}^{r_a} \frac{dr}{\sqrt{\mathcal{R}[r]}}. \qquad (524) $$

We introduce the fractional angle (i.e., the angle divided by $2\pi$) of the advance of the periastron *per* orbital revolution,

$$ K = \frac{1}{\pi} \int_{r_p}^{r_a} dr \frac{\mathcal{S}[r]}{\sqrt{\mathcal{R}[r]}}, \qquad (525) $$

which is such that the precession of the periastron *per* period is given by

<footnote n="92">⁹² The dependence on E and J will no longer be indicated but is always understood as implicit in what follows.</footnote>
$\Delta\Phi = 2\pi(K - 1)$. As $K$ tends to one in the limit $c \to \infty$ (as is easily checked from the usual Newtonian solution), it is often convenient to pose $k \equiv K - 1$, which will then entirely describe the *relativistic precession*.

Let us then define the mean anomaly $\ell$ and the mean motion $n$ by

$$ \ell \equiv n(t - t_p), \qquad (526a) $$

$$ n \equiv \frac{2\pi}{P}. \qquad (526b) $$

Here $t_p$ denotes the instant of passage to the periastron. For a given value of the mean anomaly $\ell$, the orbital separation $r$ is obtained by inversion of the integral equation

$$ \ell = n \int_{r_p}^{r} \frac{ds}{\sqrt{R[s]}}. \qquad (527) $$

This defines the function $r(\ell)$ which is a periodic function in $\ell$ with period $2\pi$. The orbital phase $\phi$ is then obtained in terms of the mean anomaly $\ell$ by integrating the angular motion as

$$ \phi = \phi_p + \frac{1}{n} \int_0^\ell d\ell S[r(\ell)], \qquad (528) $$

where $\phi_p$ denotes the value of the phase at the instant $t_p$. We may define the origin of the orbital phase at the ascending node $N$ with respect to some observer. In the particular case of a circular orbit, $r = \text{const}$, the phase evolves linearly with time, $\dot{\phi} = S[r] = \Omega$, where $\Omega$ is the orbital frequency of the circular orbit given by

$$ \Omega = Kn = (1 + k)n. \qquad (529) $$

In the general case of a non-circular orbit it is convenient to keep that definition $\Omega = Kn$ and to explicitly introduce the linearly growing part of the orbital phase (528) by writing it in the form

$$ \begin{aligned} \phi &= \phi_p + \Omega(t - t_p) + W(\ell) \\ &= \phi_p + K\ell + W(\ell). \end{aligned} \qquad (530) $$

Here $W(\ell)$ denotes a certain function of the mean anomaly which is periodic in $\ell$ with period $2\pi$, hence periodic in time with period $P$. According to Eq. (528) this function is given in terms of the mean anomaly $\ell$ by

$$ W(\ell) = \frac{1}{n} \int_{0}^{\ell} dl (S[r(l)] - \Omega). \qquad (531) $$

Finally, the decomposition (530) exhibits clearly the nature of the compact binary motion, which may be called doubly periodic in that the mean anomaly $\ell$ is periodic with period $2\pi$, and the periastron advance $K\ell$ is periodic with period $2\pi K$. Notice however that, though standard, the term “doubly periodic” is misleading since the
motion in physical space is not periodic in general. The radial motion $r(t)$ is periodic with period $P$ while the angular motion $\phi(t)$ is periodic [modulo $2\pi$] with period $P/k$ where $k = K - 1$. Only when the two periods are commensurable, i.e., when $k = 1/N$ where $N \in \mathbb{N}$, is the motion periodic in physical space (with period $NP$).

### 3.5.2 Quasi-Keplerian representation of the motion

The quasi-Keplerian (QK) representation of the motion of compact binaries is an elegant formulation of the solution of the 1PN equations of motion parametrized by the eccentric anomaly $u$ (entering a specific generalization of Kepler's equation) and depending on various orbital elements, such as several types of eccentricities. It was introduced by Damour and Deruelle (1985, 1986) to study the problem of binary pulsar timing data including relativistic corrections at the 1PN order, where notably the relativistic periastron precession complicates the simpler Keplerian solution.

In the QK representation the radial motion is given as

$$r = a_r(1 - e_r \cos u) + O\left(\frac{1}{c^4}\right), \qquad (532)$$

which is the standard form parametrized by the eccentric anomaly $u$, with $a_r$ and $e_r$ denoting two constants representing the semi-major axis of the orbit and its eccentricity. However, these constants are labelled after the radial coordinate $r$ to remember that they enter (by definition) into the radial equation; in particular $e_r$ will differ from other kinds of eccentricities $e_t$ and $e_\phi$. The “time” eccentricity $e_t$ enters the Kepler equation which at the 1PN order also takes the usual form

$$\ell = u - e_t \sin u + O\left(\frac{1}{c^4}\right), \qquad (533)$$

where the mean anomaly $\ell$, proportional to the time elapsed since the instant $t_p$ of passage at the periastron, and the mean motion $n$ are given by (526). The “angular” eccentricity $e_\phi$ enters the equation for the angular motion at 1PN order which is written in the form

$$\frac{\phi - \phi_p}{K} = v + O\left(\frac{1}{c^4}\right), \qquad (534)$$

where we define<fnref n="93" />

$$v \equiv 2 \arctan \left[ \left( \frac{1 + e_{\phi}}{1 - e_{\phi}} \right)^{1/2} \tan \frac{u}{2} \right]. \qquad (535)$$

The constant $K$ is the advance of periastron per orbital revolution defined by

$$v = \ell + \frac{W(\ell)}{K} + O\left(\frac{1}{c^4}\right).$$

<footnote n="93">⁹³ Comparing with Eqs. (530) we have also</footnote>
Eq. (525); it may be written as $K = \frac{\Phi}{2\pi}$ where $\Phi$ is the angle of return to the periastron.

Crucial to the QK formalism are the explicit expressions for the orbital elements $n$, $K$, $a_r$, $e_r$, $e_t$ and $e_\phi$ in terms of the conserved energy $E$ and angular momentum $J$ of the orbit. For convenience we introduce two dimensionless parameters directly linked to $E$ and $J$ by

$$ \epsilon \equiv -\frac{2E}{\mu c^2}, \qquad (536a) $$

$$ j \equiv -\frac{2E h^2}{\mu^3}, \qquad (536b) $$

where $\mu = mv$ is the reduced mass with $m$ the total mass (recall that $E<0$ for bound orbits) and we have used the intermediate standard notation $h \equiv \frac{1}{Gm}$. The equations to follow will then appear as expansions in powers of the small post-Newtonian parameter $\epsilon = \mathcal{O}(1/c^2)$,<fnref n="94" /> with coefficients depending on $j$ and the dimensionless reduced mass ratio $v$; notice that the parameter $j$ is of Newtonian order: $j = \mathcal{O}(1/c^0)$. We have (Damour and Deruelle 1985)

$$ n = \frac{\epsilon^{3/2} c^3}{Gm} \left\{ 1 + \frac{\epsilon}{8} (-15 + v) + \mathcal{O}\left(\frac{1}{c^4}\right) \right\}, \qquad (537a) $$

$$ K = 1 + \frac{3\epsilon}{j} + \mathcal{O}\left(\frac{1}{c^4}\right), \qquad (537b) $$

$$ a_r = \frac{Gm}{\epsilon c^2} \left\{ 1 + \frac{\epsilon}{4}(-7 + v) + \mathcal{O}\left(\frac{1}{c^4}\right) \right\}, \qquad (537c) $$

$$ e_r = \sqrt{1-j} + \frac{\epsilon}{8\sqrt{1-j}} \left[ 24 - 4v + 5j(-3+v) \right] + \mathcal{O}\left(\frac{1}{c^4}\right), \qquad (537d) $$

$$ e_t = \sqrt{1-j} + \frac{\epsilon}{8\sqrt{1-j}} \left[ -8 + 8v + j(17-7v) \right] + \mathcal{O}\left(\frac{1}{c^4}\right), \qquad (537e) $$

$$ e_{\phi} = \sqrt{1-j} + \frac{\epsilon}{8\sqrt{1-j}} \left[ 24 + j(-15+v) \right] + \mathcal{O}\left(\frac{1}{c^4}\right). \qquad (537f) $$

The dependence of such relations on the coordinate system will be discussed later. Notice the interesting point that there is no dependence of the mean motion $n$ and the radial semi-major axis $a_r$ on the angular momentum $J$ up to the 1PN order; such dependence starts at 2PN order $\propto \epsilon^2$, see e.g. (539a).

The QK representation of the compact binary motion has been generalized at the 2PN order by Damour and Schäfer (1988), Schäfer and Wex (1993), Wex (1995), at

<footnote n="94">⁹⁴ Note that the post-Newtonian parameter $\epsilon$ is precisely defined by Eq. (536a), and should not be confused with $\epsilon$ in Eq. (1) which is intended to represent a post-Newtonian *estimate*.</footnote>
3PN order by Memmesheimer et al. (2004), and at 4PN order by Cho et al. (2022b) with, however, an incomplete treatment of the 4PN tail term.<fnref n="95" /> The construction of a generalized QK representation at 3PN order exploits the fact that the radial equation given by Eq. (523a) is a polynomial in $1/r$. However, this is true only in coordinate systems avoiding the appearance of terms with the logarithm $\ln r$; the presence of logarithms in the standard harmonic (SH) coordinates at 3PN order will obstruct the construction of the QK parametrization. Therefore Memmesheimer et al. (2004) obtained it in the ADM coordinate system—free of logarithms—and also in the modified harmonic (MH) coordinate system, for which the coefficients in the equations of motion are given by Eqs. (362).

At the 3PN order the radial equation in ADM or MH coordinates is still given by Eq. (532). However, the Kepler equation (533) and angular equation (534) acquire extra contributions. The equations become

$$r = a_r(1 - e_r \cos u) + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (538a)$$

$$\ell = u - e_t \sin u + f_t \sin v + g_t (v - u) + i_t \sin 2v + h_t \sin 3v + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (538b)$$

$$\frac{\phi - \phi_p}{K} = v + f_{\phi} \sin 2v + g_{\phi} \sin 3v + i_{\phi} \sin 4v + h_{\phi} \sin 5v + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (538c)$$

in which $v$ is exactly given by Eq. (535). The new orbital elements $f_t, f_\phi, g_t, g_\phi$ are composed of 2PN and 3PN terms, while $i_t, i_\phi, h_t, h_\phi$ start only at 3PN order. All the orbital elements are now to be related, similarly to Eqs. (537), to the constants $\epsilon$ and $j$ with 3PN accuracy in a given coordinate system.

Let us make clear that in different coordinate systems such as MH and ADM coordinates, the QK representation takes exactly the same form as given by Eqs. (538). But, the relations linking the various orbital elements $a_r, e_r, e_t, e_\phi, f_t, f_\phi, \dots$ to E and J or $\epsilon$ and $j$, are different, with the notable exceptions of $n$ and $K$. Indeed, an important point related to the use of gauge invariant variables in the elliptical orbit case is that the functional forms of the mean motion $n$ and periastron advance $K$ in terms of the gauge invariant variables $\epsilon$ and $j$ are identical in different coordinate systems like the MH and ADM coordinates (Damour and Schäfer 1988). Their explicit expressions at 3PN order read

$$n = \frac{\epsilon^{3/2} c^3}{G m} \left\{ \begin{aligned} & 1 + \frac{\epsilon}{8} (-15 + v) + \frac{\epsilon^2}{128} \left[ 555 + 30v + 11v^2 + \frac{192}{j^{1/2}} (-5 + 2v) \right] \\ & + \frac{\epsilon^3}{3072} \left[ -29385 - 4995v - 315v^2 + 135v^3 + \frac{5760}{j^{1/2}} (17 - 9v + 2v^2) \right] \\ & + \frac{16}{j^{3/2}} \left( -10080 + (13952 - 123\pi^2)v - 1440v^2 \right) \end{aligned} \right\} + \mathcal{O}\left(\frac{1}{c^8}\right). \quad (539a)$$

<footnote n="95">⁹⁵ See also Trestini (2024) for the QK representation in scalar-tensor theory at 2PN order.</footnote>
$$ K = 1 + \frac{3\epsilon}{j} + \frac{\epsilon^2}{4} \left[ \frac{3}{j}(-5 + 2v) + \frac{15}{j^2}(7 - 2v) \right] \\ + \frac{\epsilon^3}{128} \left[ \frac{24}{j}(5 - 5v + 4v^2) + \frac{1}{j^2} \left( -10080 + (13952 - 123\pi^2)v - 1440v^2 \right) \right. \\ \left. + \frac{5}{j^3} \left( 7392 + (-8000 + 123\pi^2)v + 336v^2 \right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right). \tag{539b} $$

Because of their gauge invariant meaning, it is natural to use $n$ and $K$ as two independent gauge-invariant variables in the general orbit case. Actually, instead of working with the mean motion $n$ it is often preferable to use the orbital frequency $\Omega$ which has been defined for general quasi-elliptic orbits in Eq. (529). Moreover we can pose

$$ x = \left( \frac{Gm\Omega}{c^3} \right)^{2/3} \quad (\text{with } \Omega = Kn), \tag{540} $$

which constitutes the obvious generalization of the gauge invariant variable $x$ used in the circular orbit case. The use of $x$ as an independent parameter will thus facilitate the straightforward reading out and check of the circular orbit limit. The parameter $x$ is related to the energy and angular momentum variables $\epsilon$ and $j$ up to 3PN order by

$$ \begin{aligned} x = \epsilon \Biggl\{ & 1 + \epsilon \left[-\frac{5}{4} + \frac{1}{12}v + \frac{2}{j}\right] \\ & + \epsilon^2 \left[\frac{5}{2} + \frac{5}{24}v + \frac{1}{18}v^2 + \frac{1}{j^{1/2}}(-5+2v) + \frac{1}{j}\left(-5+\frac{7}{6}v\right) + \frac{1}{j^2}\left(\frac{33}{2}-5v\right)\right] \\ & + \epsilon^3 \left[-\frac{235}{48} - \frac{25}{24}v - \frac{25}{576}v^2 + \frac{35}{1296}v^3 + \frac{1}{j}\left(\frac{35}{4} - \frac{5}{3}v + \frac{25}{36}v^2\right)\right. \\ & \left. + \frac{1}{j^{1/2}}\left(\frac{145}{8} - \frac{235}{24}v + \frac{29}{12}v^2\right) + \frac{1}{j^{3/2}}\left(-45 + \left(\frac{472}{9} - \frac{41}{96}\pi^2\right)v - 5v^2\right) \right. \\ & \left. + \frac{1}{j^2}\left(-\frac{565}{8} + \left(\frac{1903}{24} - \frac{41}{64}\pi^2\right)v - \frac{95}{12}v^2\right) \right. \\ & \left. + \frac{1}{j^3}\left(\frac{529}{3} + \left(-\frac{610}{3} + \frac{205}{64}\pi^2\right)v + \frac{35}{4}v^2\right)\right] + \mathcal{O}\left(\frac{1}{c^8}\right) \Biggr\}. \end{aligned} \tag{541} $$

Besides the very useful gauge-invariant quantities $n$, $K$ and $x$, the other orbital elements $a_r$, $e_r$, $e_t$, $e_\phi$, $f_t$, $g_t$, $i_t$, $h_t$, $f_\phi$, $g_\phi$, $i_\phi$, $h_\phi$ parametrizing Eqs. (532) and (538) are not gauge invariant; their expressions in terms of $\epsilon$ and $j$ depend on the coordinate system in use. We refer to Memmesheimer et al. (2004), Arun et al. (2008a) for the full expressions of all the orbital elements at 3PN order in both MH and ADM coordinate systems. Here, for future use, we only give the expression of the time eccentricity $e_t$ (squared) in MH coordinates:
$$
\begin{equation}
\begin{aligned}
e_t^2 = {}& 1 - j + \frac{\epsilon}{4} \left[ -8 + 8v + j(17 - 7v) \right] \\
& + \frac{\epsilon^2}{8} \left[ 12 + 72v + 20v^2 + j(-112 + 47v - 16v^2) + 24j^{1/2}(5 - 2v) \right. \\
& \left. \qquad + \frac{16}{j}(4 - 7v) + \frac{24}{j^{1/2}}(-5 + 2v) \right] \\
& + \frac{\epsilon^3}{6720} \left[ 23520 - 464800v + 179760v^2 + 16800v^3 \right. \\
& \left. \qquad + 525j(528 - 200v + 77v^2 - 24v^3) \right] \\
& + 2520j^{1/2}(-265 + 193v - 46v^2) \\
& + \frac{6}{j} \left( -73920 + (260272 - 4305\pi^2)v - 61040v^2 \right) \\
& + \frac{70}{j^{1/2}} \left( 16380 + (-19964 + 123\pi^2)v + 3240v^2 \right) \\
& + \frac{70}{j^{3/2}} \left( -10080 + (13952 - 123\pi^2)v - 1440v^2 \right) \\
& \left. + \frac{8}{j^2} \left( 53760 + (-176024 + 4305\pi^2)v + 15120v^2 \right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right).
\end{aligned}
\tag{542}
\end{equation}
$$

Again, with our notation (536), this appears as a post-Newtonian expansion in the small parameter $\epsilon \to 0$ with fixed “Newtonian” parameter $j$.

### 3.5.3 Averaged energy and angular momentum fluxes

The gravitational-wave energy and angular momentum fluxes from a system of two point masses in elliptic motion was first computed by Peters and Mathews (1963), Peters (1964) at Newtonian level. The 1PN and 1.5PN corrections to the fluxes were provided by Wagoner and Will (1976), Blanchet and Schäfer (1989), Junker and Schäfer (1992), Blanchet and Schäfer (1993), Rieth and Schäfer (1997) and used to study the associated secular evolution of orbital elements under gravitational radiation reaction using the QK representation of the binary's orbit at 1PN order (Damour and Deruelle 1985). These results were extended to 2PN order by Gopakumar and Iyer (1997, 2002) for the instantaneous terms (leaving aside the tails) using the generalized QK representation (Damour and Schäfer 1988; Schäfer and Wex 1993; Wex 1995); the energy flux and waveform were in agreement with those of Will and Wiseman (1996) obtained using a different method. Arun et al. (2008a, 2009a, b) have generalized the results at 3PN order, including all tails and related hereditary contributions (Arun et al. 2008b; Loutrel and Yunes 2017), by computing the averaged energy and angular momentum fluxes for quasi-elliptical orbits using the QK representation at 3PN order (Memmesheimer et al. 2004), and deriving the secular evolution of the orbital elements under 3PN gravitational radiation reaction. The waveform of eccentric binaries has been derived up to the
2PN order in the Fourier domain by Tessmer and Schäfer (2010, 2011). The gravitational-wave modes, including the non-linear memory and tail contributions, have been derived at the 3PN order by Ebersold et al. (2019), Boetzel et al. (2019). The gravitational energy and angular momentum fluxes of eccentric orbits in the small mass ratio limit (case of EMRI systems) have been derived to high 9PN order by Forseth et al. (2016), Munna et al. (2020). Furthermore the expansion in powers of the eccentricity has been extended to high order using the method of singular eccentricity factors, allowing the flux to be accurately determined when $e \to 1$ (Forseth et al. 2016).

The secular evolution of orbital elements under gravitational radiation reaction is in principle only the starting point for constructing templates for eccentric binary orbits. To go beyond the secular evolution one needs to include in the evolution of the orbital elements, besides the averaged contributions in the fluxes, the post-adiabatic terms rapidly oscillating at the orbital period. An analytic approach, based on an improved method of variation of constants, has been introduced by Damour et al. (2004) at the leading post-adiabatic approximation, i.e. the level of 2.5PN radiation reaction.

The generalized QK representation of the motion discussed in Sect. 3.5.2 plays a crucial role in the procedure of averaging the energy and angular momentum fluxes $\mathcal{F}$ and $\mathcal{G}_i$ over one orbit.<fnref n="96" /> Actually the averaging procedure applies to the “instantaneous” parts of the fluxes, while the “hereditary” parts are treated separately for technical reasons (Arun et al. 2008b; Loutrel and Yunes 2017). We pose $\mathcal{F} = \mathcal{F}_{\text{inst}} + \mathcal{F}_{\text{hered}}$ where the hereditary part of the energy flux is composed of tails and tail-of-tails (limiting ourselves to 3PN). For the angular momentum flux one needs also to include a contribution from the memory (Arun et al. 2009a; Loutrel and Yunes 2017). We thus have to compute for the instantaneous part

$$ \langle \mathcal{F}_{\text{inst}} \rangle = \frac{1}{P} \int_{0}^{P} dt \, \mathcal{F}_{\text{inst}} = \frac{1}{2\pi} \int_{0}^{2\pi} du \, \frac{d\ell}{du} \mathcal{F}_{\text{inst}}, \qquad (543) $$

and similarly for the instantaneous part of the angular momentum flux $\mathcal{G}_i$.

Thanks to the QK representation, we can express $\mathcal{F}_{\text{inst}}$, which is initially a function of the natural variables $r, \dot{r}$ and $v^2$, as a function of the varying eccentric anomaly $u$, and depending on two constants: The frequency-related parameter $x$ defined by (540), and the “time” eccentricity $e_t$ given by (542). To do so one must select a particular coordinate system—the MH coordinates for instance. The choice of $e_t$ rather than $e_r$ (say) is a matter of convenience; since $e_t$ appears in the Kepler-like equation (538b) at leading order, it will directly be dealt with when averaging over one orbit. We note that in the expression of the energy flux at the 3PN order there are some logarithmic terms of the type $\ln(r/r_0)$ even in MH coordinates. Indeed, as we have seen in Sect. 3.2.4, the MH coordinates permit the removal of the logarithms $\ln(r/r'_0)$ in the equations of motion, where $r'_0$ is the UV scale in Eq. (357); however there are still some logarithms $\ln(r/r_0)$ which involve the IR constant $r_0$ entering the definition of the multipole moments for general sources, see Theorem 7

<footnote n="96">⁹⁶ Recall that the fluxes are defined in a general way, for any matter system, in terms of the radiative multipole moments by the expressions (179).</footnote>
and the regularization factor (72). As a result we find that the general structure of $\mathcal{F}_{\text{inst}}$ (and similarly for $\mathcal{G}_{\text{inst}}$, the norm of the angular momentum flux) consists of a finite sum of terms of the type

$$ \mathcal{F}_{\text{inst}} = \frac{du}{d\ell} \sum_{l} \frac{\alpha_l(x, e_l) + \beta_l(x, e_l) \sin u + \gamma_l(x, e_l) \ln(1 - e_l \cos u)}{(1 - e_l \cos u)^{l+1}}. \quad (544) $$

The factor $du/d\ell$ has been inserted to prepare for the orbital average (543). The coefficients $\alpha_l$, $\beta_l$ and $\gamma_l$ are straightforwardly computed using the QK parametrization as functions of $x$ and the time eccentricity $e_l$. The $\beta_l$'s correspond to 2.5PN radiation-reaction terms and will play no role, while the $\gamma_l$'s correspond to the logarithmic terms $\ln(r/r_0)$ arising at the 3PN order. For convenience the dependence on the constant $\ln r_0$ has been included into the coefficients $\alpha_l$'s. To compute the average we dispose of the following integration formulas ($l \in \mathbb{N}$)$^{97}$

$$ \int_{0}^{2\pi} \frac{du}{2\pi} \frac{\sin u}{(1 - e_l \cos u)^{l+1}} = 0, \quad (545a) $$

$$ \int_{0}^{2\pi} \frac{du}{2\pi} \frac{1}{(1 - e_l \cos u)^{l+1}} = \frac{(-1)^l}{l!} \left( \frac{d^l}{dz^l} \left[ \frac{1}{\sqrt{z^2 - e_l^2}} \right] \right) \Big|_{z=1}, \quad (545b) $$

$$ \int_{0}^{2\pi} \frac{du}{2\pi} \frac{\ln(1 - e_l \cos u)}{(1 - e_l \cos u)^{l+1}} = \frac{(-1)^l}{l!} \left( \frac{d^l}{dz^l} \left[ \frac{Z(z, e_l)}{\sqrt{z^2 - e_l^2}} \right] \right) \Big|_{z=1}. \quad (545c) $$

In the right-hand sides of Eqs. (545b) and (545c) we have to differentiate $l$ times with respect to the intermediate variable $z$ before applying $z = 1$. The equation (545c), necessary for dealing with the logarithmic terms, contains the non-trivial function

$$ Z(z, e_l) = \ln \left[ \frac{\sqrt{1 - e_l^2 + 1}}{2} \right] + 2 \ln \left[ 1 + \frac{\sqrt{1 - e_l^2 - 1}}{z + \sqrt{z^2 - e_l^2}} \right]. \quad (546) $$

From Eq. (545a) we see that there will be no radiation-reaction terms at 2.5PN order in the final result; the 2.5PN contribution is proportional to $\dot{r}$ and vanishes after averaging since it involves only odd functions of $u$.

Finally, after implementing all the above integrations, the averaged instantaneous energy flux in MH coordinates at the 3PN order is obtained as

$$ \langle \mathcal{F}_{\text{inst}} \rangle = \frac{32c^5}{5G} v^2 x^5 (\mathcal{I}_0 + x \mathcal{I}_1 + x^2 \mathcal{I}_2 + x^3 \mathcal{I}_3), \quad (547) $$

where we recall that the post-Newtonian parameter $x$ is defined by (540). The various

$$ \int_{0}^{2\pi} \frac{du}{2\pi} \frac{1}{(1 - e_l \cos u)^{l+1}} = \frac{1}{(1 - e_l^2)^{\frac{l+1}{2}}} P_l \left( \frac{1}{\sqrt{1 - e_l^2}} \right). $$

<footnote>$^{97}$ The second of these formulas can alternatively be written with the standard Legendre polynomial $P_l$ as</footnote>
instantaneous post-Newtonian pieces depend on the symmetric mass ratio $v$ and the time eccentricity $e_t$ in MH coordinates as (Arun et al. 2008a)

$$ \mathcal{I}_0 = \frac{1}{(1 - e_t^2)^{7/2}} \left( 1 + \frac{73}{24} e_t^2 + \frac{37}{96} e_t^4 \right), \qquad (548a) $$

$$ \mathcal{I}_1 = \frac{1}{(1 - e_t^2)^{9/2}} \left\{ -\frac{1247}{336} - \frac{35}{12} v + e_t^2 \left( \frac{10475}{672} - \frac{1081}{36} v \right) + e_t^4 \left( \frac{10043}{384} - \frac{311}{12} v \right) + e_t^6 \left( \frac{2179}{1792} - \frac{851}{576} v \right) \right\}, \qquad (548b) $$

$$ \begin{aligned} \mathcal{I}_2 = & \frac{1}{(1 - e_t^2)^{11/2}} \left\{ -\frac{203471}{9072} + \frac{12799}{504} v + \frac{65}{18} v^2 \right. \\ & + e_t^2 \left( -\frac{3807197}{18144} + \frac{116789}{2016} v + \frac{5935}{54} v^2 \right) \\ & + e_t^4 \left( -\frac{268447}{24192} - \frac{2465027}{8064} v + \frac{247805}{864} v^2 \right) \\ & + e_t^6 \left( \frac{1307105}{16128} - \frac{416945}{2688} v + \frac{185305}{1728} v^2 \right) \\ & + e_t^8 \left( \frac{86567}{64512} - \frac{9769}{4608} v + \frac{21275}{6912} v^2 \right) \\ & \left. + \sqrt{1 - e_t^2} \left[ \frac{35}{2} - 7v + e_t^2 \left( \frac{6425}{48} - \frac{1285}{24} v \right) \right. \right. \\ & \left. \left. + e_t^4 \left( \frac{5065}{64} - \frac{1013}{32} v \right) + e_t^6 \left( \frac{185}{96} - \frac{37}{48} v \right) \right] \right\}, \end{aligned} \qquad (548c) $$

$$ \begin{aligned} \mathcal{I}_3 = & \frac{1}{(1 - e_t^2)^{13/2}} \left\{ \frac{2193295679}{9979200} + \left[ \frac{8009293}{54432} - \frac{41\pi^2}{64} \right] v - \frac{209063}{3024} v^2 - \frac{775}{324} v^3 \right. \\ & + e_t^2 \left( \frac{20506331429}{19958400} + \left[ \frac{649801883}{272160} + \frac{4879\pi^2}{1536} \right] v - \frac{3008759}{3024} v^2 - \frac{53696}{243} v^3 \right) \\ & + e_t^4 \left( -\frac{3611354071}{13305600} + \left[ \frac{755536297}{136080} - \frac{29971\pi^2}{1024} \right] v - \frac{179375}{576} v^2 - \frac{10816087}{7776} v^3 \right) \\ & + e_t^6 \left( \frac{4786812253}{26611200} + \left[ \frac{1108811471}{1451520} - \frac{84501\pi^2}{4096} \right] v + \frac{87787969}{48384} v^2 - \frac{983251}{648} v^3 \right) \\ & + e_t^8 \left( \frac{21505140101}{141926400} + \left[ -\frac{32467919}{129024} - \frac{4059\pi^2}{4096} \right] v + \frac{79938097}{193536} v^2 - \frac{4586539}{15552} v^3 \right) \end{aligned} $$
$$ 
\begin{aligned}
&+e_t^{10} \left( -\frac{8977637}{11354112} + \frac{9287}{48384} v + \frac{8977}{55296} v^2 - \frac{567617}{124416} v^3 \right) \\
&+ \sqrt{1-e_t^2} \left[ -\frac{14047483}{151200} + \left[ -\frac{165761}{1008} + \frac{287\pi^2}{192} \right] v + \frac{455}{12} v^2 \right. \\
&\qquad \left. + e_t^2 \left( \frac{36863231}{100800} + \left[ -\frac{14935421}{6048} + \frac{52685\pi^2}{4608} \right] v + \frac{43559}{72} v^2 \right) \right. \\
&\qquad \left. + e_t^4 \left( \frac{759524951}{403200} + \left[ -\frac{31082483}{8064} + \frac{41533\pi^2}{6144} \right] v + \frac{303985}{288} v^2 \right) \right. \\
&\qquad \left. + e_t^6 \left( \frac{1399661203}{2419200} + \left[ -\frac{40922933}{48384} + \frac{1517\pi^2}{9216} \right] v + \frac{73357}{288} v^2 \right) \right. \\
&\qquad \left. + e_t^8 \left( \frac{185}{48} - \frac{1073}{288} v + \frac{407}{288} v^2 \right) \right] \\
&\qquad + \left( \frac{1712}{105} + \frac{14552}{63} e_t^2 + \frac{553297}{1260} e_t^4 \right. \\
&\qquad \left. + \frac{187357}{1260} e_t^6 + \frac{10593}{2240} e_t^8 \right) \ln \left[ \frac{x(1+\sqrt{1-e_t^2})}{z_0(2(1-e_t^2))} \right] \right\}.
\end{aligned} 
\quad (548d) $$

We recall that $e_t$ corresponds to MH coordinates and is given by Eq. (542). On the contrary $x$ is gauge invariant—since both $n$ and $K$ are gauge invariant. In the small mass ratio limit all the coefficients in (548) have been obtained independently by Forseth et al. (2016).

The Newtonian coefficient in (548) is nothing but the “enhancement” function of eccentricity which enters in the orbital radiation decay of the binary pulsar, see Eq. (12). It was computed by Peters and Mathews (1963) from the average of the usual quadrupole formula in the time domain, but also in the Fourier domain, using the Fourier decomposition of the Keplerian motion. We have $\mathcal{I}_0 \equiv f$ with

$$ f(e_t) = \sum_{p=1}^{+\infty} g(p, e_t) = \frac{1 + \frac{73}{24} e_t^2 + \frac{37}{96} e_t^4}{(1 - e_t^2)^{7/2}}, \quad (549) $$

where $g(p, e_t)$ is directly given by the discrete Fourier coefficients of the Newtonian STF quadrupole moment $\hat{Q}_{ij}$ (normalized by $\mu a^2$) as

$$ g(p, e_t) = \frac{p^6}{16} |\hat{Q}_{ij}|^2, \quad (550) $$

and admits the following expression in terms of Bessel functions:
$$g(p, e_t) = \frac{1}{2} p^2 \left\{ \left[ -\frac{4}{e_t^3} - 3e_t + \frac{7}{e_t} \right] p J_p(pe_t) J'_p(pe_t) \right. \\ \left. + \left[ \frac{1}{e_t^4} - \frac{1}{e_t^2} + \frac{1}{3} + \left( \frac{1}{e_t^4} - e_t^2 - \frac{3}{e_t^2} + 3 \right) p^2 \right] (J_p(pe_t))^2 \right. \\ \left. + \left[ \left( e_t^2 + \frac{1}{e_t^2} - 2 \right) p^2 + \frac{1}{e_t^2} - 1 \right] (J'_p(pe_t))^2 \right\}. \quad (551)$$

The last term in the 3PN coefficient of (548) is proportional to some logarithm which directly arises from the integration formula (545c). Inside the logarithm we have posed $z_0 \equiv Gm/(c^2r_0)$, exhibiting the dependence upon the arbitrary length scale $r_0$ introduced in the formalism through Eq. (72). Only after adding the hereditary contribution to the 3PN energy flux can we check the required cancellation of $r_0$. The hereditary part at 3PN order is of the form

$$\langle \mathcal{F}_{\text{hered}} \rangle = \frac{32c^5}{5G} v^2 x^5 \left( x^{3/2} \mathcal{K}_{3/2} + x^{5/2} \mathcal{K}_{5/2} + x^3 \mathcal{K}_3 \right), \quad (552)$$

where the 1.5PN and 2.5PN terms are due to tails, and the 3PN term is due to tails-of-tails. We have (Arun et al. 2008b)

$$\mathcal{K}_{3/2} = 4\pi \varphi(e_t), \quad (553a)$$

$$\mathcal{K}_{5/2} = -\frac{8191}{672} \pi \psi(e_t) - \frac{583}{24} v \pi \zeta(e_t), \quad (553b)$$

$$\mathcal{K}_3 = -\frac{116761}{3675} \kappa(e_t) \\ + \left[ \frac{16}{3} \pi^2 - \frac{1712}{105} \gamma_E - \frac{1712}{105} \ln\left(\frac{4x^{3/2}}{z_0}\right) \right] F(e_t), \quad (553c)$$

where $\varphi(e_t)$, $\psi(e_t)$, $\zeta(e_t)$, $\kappa(e_t)$ and $F(e_t)$ denote certain “enhancement” functions of the eccentricity, with normalization such that they equal 1 for circular orbits, when $e_t = 0$. The quadrupole moment (550) also determines the function $F(e_t)$ in factor of the logarithm in the 3PN piece, which turns out to admit a closed analytic form:

$$F(e_t) = \frac{1}{4} \sum_{p=1}^{+\infty} p^2 g(p, e_t) \\ = \frac{1}{(1 - e_t^2)^{13/2}} \left[ 1 + \frac{85}{6} e_t^2 + \frac{5171}{192} e_t^4 + \frac{1751}{192} e_t^6 + \frac{297}{1024} e_t^8 \right]. \quad (554)$$

This result is important for checking that the arbitrary constant $z_0$ disappears from the final result. Indeed this is immediately seen from comparing the last term in (548d) with (553c).

On the other hand, the four functions $\varphi(e_t)$, $\psi(e_t)$, $\zeta(e_t)$ and $\kappa(e_t)$ defined in Eqs. (553) do not admit analytic closed-form expressions. They have been derived by Arun et al. (2008b) (extending Blanchet and Schäfer 1993) in the form of infinite
series made out of quadratic products of Bessel functions, and computed numerically. For instance the numerical graph of the function $\varphi(e_t)$ is given in Fig. 8. The leading expansions when $e_t \to 0$ has been computed:

$$ \varphi(e_t) = 1 + \frac{2335}{192} e_t^2 + \mathcal{O}(e_t^4), \qquad (555a) $$

$$ \psi(e_t) = 1 - \frac{22988}{8191} e_t^2 + \mathcal{O}(e_t^4), \qquad (555b) $$

$$ \zeta(e_t) = 1 + \frac{1011565}{48972} e_t^2 + \mathcal{O}(e_t^4), \qquad (555c) $$

$$ \kappa(e_t) = 1 + \left( \frac{62}{3} - \frac{4613840}{350283} \ln 2 + \frac{24570945}{1868176} \ln 3 \right) e_t^2 + \mathcal{O}(e_t^4). \qquad (555d) $$

Important improvements on the analytic determination of these functions have been made by Forseth et al. (2016), using the method of “singular eccentricity factors”.<fnref n="98" /> Consider the function $\varphi(e_t)$ which enters the dominant 1.5PN tail term in Eq. (552). It is obtained as

$$ \varphi(e_t) = \sum_{p=1}^{+\infty} \frac{p}{2} g(p, e_t). \qquad (556) $$

The method of singular factors consists of determining a factor $(1 - e_t^2)^{-q}$ blowing up when $e_t \to 1$, such that $(1 - e_t^2)^q \varphi(e_t)$ admits a convergent expansion in powers of $e_t^2$ taking a finite value when $e_t = 1$. So, the divergence of the function when $e_t \to 1$ is factorized out into some $(1 - e_t^2)^{-q}$. Forseth et al. (2016) find $q = 5$ in the case of the function (556) and are able to determine

$$ \varphi(e_t) = \frac{1}{(1 - e_t^2)^5} \left( 1 + \frac{1375}{192} e_t^2 + \frac{3935}{768} e_t^4 + \frac{10007}{36864} e_t^6 + \dots \right), \qquad (557) $$

with many more coefficients in the series shown in Eq. (4.17) of Forseth et al. (2016). As checked numerically the power series which is left after the factorization (557) converges at $e_t = 1$. Finally all the functions can be pushed all the way to $e_t = 1$ which permits to extend analytically the validity of the flux in the regime of EMRIs, where $e_t$ is close to one. This is also useful for studying gravitational waves from binary black holes with moderately high eccentricities, such as those formed in globular clusters (Gultekin et al. 2004). On the other hand one can retrieve from (557) the expansion when $e_t \to 0$ to very high order,

$$ \varphi(e_t) = 1 + \frac{2335}{192} e_t^2 + \frac{42955}{756} e_t^4 + \frac{6204647}{36864} e_t^6 + \mathcal{O}(e_t^8). \qquad (558) $$

The averaged energy flux was presented above using $x$—the gauge invariant variable (540)—and the time eccentricity $e_t$ which however is gauge dependent. Of course it

<footnote n="98">⁹⁸ See also the resummation of asymptotic enhancement factors (Loutrel and Yunes 2017).</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 8 The enhancement factor $\varphi(e)$ as function of the eccentricity $e$. This function agrees with the numerical calculation of Blanchet and Schäfer (1993) modulo a trivial rescaling with the Peters–Mathews function (548a). The inset graph is a zoom of the function at a smaller scale. The dots represent the numerical computation and the solid line is a fit to the numerical points. In the circular orbit limit we have $\varphi(0) = 1$</figcaption></figure>

is possible to provide a fully gauge invariant formulation of the energy flux. The most natural choice is to express the result in terms of the conserved energy $E$ and angular momentum $J$, or, rather, in terms of the pair of rescaled variables $(\epsilon, j)$ defined by Eqs. (536). To this end it suffices to replace $e_t$ by its MH-coordinate expression (542) and to use Eq. (541) to re-express $x$ in terms of $\epsilon$ and $j$. However, there are other possible choices for a couple of gauge invariant quantities. As we have seen the mean motion $n$ and the periastron precession $K$ are separately gauge invariant so we may define the pair of variables to be $(x, t)$, where $x$ is given by (540) and we pose

$$t \equiv \frac{3x}{K-1}. \qquad (559)$$

Such choice would be motivated by the fact that $t$ reduces to the angular-momentum related variable $j$ in the limit $\epsilon \to 0$. Note however that with the latter choices $(\epsilon, j)$ or $(x, t)$ of gauge-invariant variables, the circular-orbit limit is not directly readable from the result; this is why we have preferred to present it in terms of the gauge dependent couple of variables $(x, e_t)$.

As we are interested in the phasing of binaries moving in quasi-eccentric orbits in the adiabatic approximation, we require the orbital averages not only of the energy flux $\mathcal{F}$ but also of the angular momentum flux $\mathcal{G}_j$. Since the quasi-Keplerian orbit is planar, we only need to average the magnitude $\mathcal{G}$ of the angular momentum flux. The complete computation thus becomes a generalisation of the previous computation of
the averaged energy flux requiring similar steps. The angular momentum flux is split into instantaneous $G_{inst}$ and hereditary $G_{hered}$ contributions; the instantaneous part is averaged using the QK representation in either MH or ADM coordinates; the hereditary part is evaluated separately and defined by means of several types of enhancement functions of the time eccentricity $e_i$; finally these are obtained numerically or analytically using the method leading to (557). At this stage we dispose of both the averaged energy and angular momentum fluxes $\langle \mathcal{F} \rangle$ and $\langle \mathcal{G} \rangle$.

The procedure to compute the secular evolution of the orbital elements under gravitational radiation-reaction is straightforward. Differentiating the orbital elements with respect to time, and using the heuristic balance equations, we equate the decreases of energy and angular momentum to the corresponding averaged fluxes $\langle \mathcal{F} \rangle$ and $\langle \mathcal{G} \rangle$ at 3PN order (Arun et al. 2009a). This extends earlier analyses at previous orders: Newtonian (Peters 1964) as we have reviewed in Sect. 1.1.2; 1PN (Blanchet and Schäfer 1989; Junker and Schäfer 1992); 1.5PN (Blanchet and Schäfer 1993; Rieth and Schäfer 1997) and 2PN (Gopakumar and Iyer 1997; Damour et al. 2004). Let us take the example of the mean motion $n$. From Eq. (539a) together with the definitions (536) we know the function $n(E, J)$ at 3PN order, where E and J are the orbit's constant energy and angular momentum. Thus,

$$ \frac{dn}{dt} = \frac{\partial n}{\partial E} \frac{dE}{dt} + \frac{\partial n}{\partial J} \frac{dJ}{dt}. \qquad (560) $$

The usual balance equations for energy and angular momentum

$$ \left\langle \frac{dE}{dt} \right\rangle = -\langle \mathcal{F} \rangle, \qquad (561a) $$

$$ \left\langle \frac{dJ}{dt} \right\rangle = -\langle \mathcal{G} \rangle, \qquad (561b) $$

have already been used at Newtonian order in Eqs. (10). Although heuristically assumed at 3PN order, they have been proved through 1.5PN order in Eq. (177). With the averaged fluxes known through 3PN order, we obtain the 3PN averaged evolution equation as

$$ \left\langle \frac{dn}{dt} \right\rangle = - \frac{\partial n}{\partial E} \langle \mathcal{F} \rangle - \frac{\partial n}{\partial J} \langle \mathcal{G} \rangle. \qquad (562) $$

We recall that this gives only the slow secular evolution under gravitational radiation reaction for eccentric orbits. The complete evolution includes also, superimposed on the averaged adiabatic evolution, some fast but smaller post-adiabatic oscillations at the orbital time scale (Damour et al. 2004; Königsdörffer and Gopakumar 2006).

## 3.6 Spinning compact binaries

The post-Newtonian templates have been developed so far for compact binary systems which can be described with high precision by point masses without spins. Here by spin, we mean the intrinsic (classical) angular momentum $S$ of the individual
compact body. Including the effects of spins is essential, as the astrophysical evidence indicates that stellar-mass black holes (Abramowicz and Kluźniak 2001; Strohmayer 2001; McClintock et al. 2006; Gou et al. 2011; Nowak et al. 2012) and supermassive black holes (Fabian and Miniutti 2009; Brenneman and Reynolds 2006; Brenneman et al. 2011) (see Reynolds 2014 for a review) can be generically close to maximally spinning. However the recent gravitational-wave observations have disfavored high spin magnitude (García-Bellido et al. 2021), although a fraction of binary black hole systems contain at least one black hole with a non-negligible spin (Adamcewicz et al. 2024).

The spins affect the orbital dynamics of the binary, in particular leading to orbital plane precession if they are not aligned with the orbital angular momentum (Cutler and Flanagan 1994; Apostolatos et al. 1994), and thereby to strong modulations in the observed signal frequency and phase. An important effort has been undertaken to compute spin effects to high post-Newtonian order in the dynamics and radiation of compact binaries:

1. *Dynamics*. The goal is to obtain the equations of motion and related conserved integrals of the motion, the equations of precession of the spins, and the post-Newtonian metric in the near zone. For this step we need a formulation of the conservative dynamics of particles with spins (either Lagrangian or Hamiltonian);

2. *Radiation*. The mass and current radiative multipole moments, including all hereditary effects, are to be computed for spinning binaries. One then deduces the gravitational waveform and the fluxes, from which we compute the secular evolution of the orbital phase. This step requires plugging the previous *dynamics* into the general wave generation formalism of Sect. 2.

We adopt a particular post-Newtonian counting for spin effects that actually refers to maximally spinning black holes. In this convention the two spin variables $S_a$ ($a = 1, 2$) have the dimension of an angular momentum multiplied by a factor $c$, and we pose

$$S_a = Gm_a^2 \chi_a, \qquad (563)$$

where $m_a$ is the mass of the compact body, and $\chi_a$ is the dimensionless spin parameter, which equals one for maximally spinning Kerr black holes. Thus the spins $S_a$ of the compact bodies can be considered as “Newtonian” quantities [there are no $c$'s in Eq. (563)], and all spin effects will carry (at least) an explicit $1/c$ factor with respect to non-spin effects. With this convention any post-Newtonian estimate is expected to be appropriate (i.e., numerically close to reality) in the case of maximal rotation. One should keep in mind that spin effects will be formally a factor $1/c$ smaller for non-maximally spinning objects such as neutron stars; thus in this case a given post-Newtonian computation will actually be a factor $1/c$ more accurate.

As usual we make a distinction between spin-orbit (SO) effects, which are linear in the spins, and spin-spin (SS) ones, which are quadratic. In the following sections we shall essentially focus on SO effects, as they play the most important role in
gravitational wave detection and parameter estimation; we shall give a short account of SS effects in Sect. 3.6.5.

The SO effects have been known at the leading level since the seminal works of Papapetrou (1951b), Corinaldesi and Papapetrou (1951), Tulczyjew (1957, 1959), Barker and O'Connell (1975, 1979), Kidder et al. (1993b), Kidder (1995). With our post-Newtonian counting such leading level corresponds to the 1.5PN order. The SO terms have been computed to the next-to-leading level which corresponds to 2.5PN order by Tagoshi et al. (2001), Faye et al. (2006), Damour et al. (2008), Levi (2010a), Porto (2010), Hartung and Steinhoff (2011a) for the equations of motion or dynamics, and by Blanchet et al. (2006), Blanchet et al. (2011) for the gravitational radiation field. Note that Tagoshi et al. (2001), Faye et al. (2006), Damour et al. (2008), Hartung and Steinhoff (2011a) employ traditional post-Newtonian methods (both harmonic-coordinates and Hamiltonian), but that Levi (2010a), Porto (2010) are based on the effective field theory (EFT) approach. The next-to-next-to-leading SO level corresponding to 3.5PN order has been obtained by Hartung and Steinhoff (2011b), Hartung et al. (2013) using the Hamiltonian method for the equations of motion, by Levi and Steinhoff (2015b) using the EFT, and by Marsat et al. (2013b), Bohé et al. (2013b) using the harmonic-coordinates method. Here we shall focus on the harmonic-coordinates approach (Marsat et al. 2013b; Bohé et al. 2013a, b; Marsat et al. 2013a) which is in fact well formulated using a Lagrangian, see Sect. 3.6.1. With this approach the next-to-next-to-leading SO level was derived not only for the equations of motion including precession, but also for the radiation field (energy flux and orbital phasing) by Bohé et al. (2013a), Marsat et al. (2013a). An analytic solution for the SO precession effects will be presented in Sect. 3.6.2. Concerning the radiation field the highest known SO level actually contains specific tail-induced contributions at 3PN and 4PN orders (Blanchet et al. 2011; Marsat et al. 2013a; Cho et al. 2022a). To next-to-next-to-next-to-leading order the SO terms in the equations of motion have been derived by Antonelli et al. (2020b), Kim et al. (2023b), Mandal et al. (2023b).

### 3.6.1 Lagrangian formalism for spinning point particles

Some necessary material for constructing a Lagrangian for a spinning point particle in curved spacetime is presented here. The formalism is issued from early works by Hanson and Regge (1974), Bailey and Israel (1975) and has also been developed in the context of the EFT approach by Porto (2006). Variants and alternatives (most importantly the associated Hamiltonian formalism) have been developed by Steinhoff et al. (2008c), Steinhoff (2011), Barausse et al. (2009). The formalism yields for the equations of motion of spinning particles and the equations of precession of the spins the classic results known in general relativity (Tulczyjew 1957, 1959; Mathisson 1937, 2010; Papapetrou 1951b; Corinaldesi and Papapetrou 1951; Trautman 2002; Dixon 1979).

Let us consider a single spinning point particle moving in a given curved background metric $g_{\alpha\beta}(x)$. The particle follows the worldline $y^\alpha(\tau)$, with tangent four-velocity $u^\alpha = dy^\alpha/d\tau$, where $\tau$ is a parameter along the representative worldline.
In a first stage we do not require that the four-velocity be normalized; thus $\tau$ needs not be the proper time elapsed along the worldline. To describe the internal degrees of freedom associated with the particle's spin, we introduce a moving orthonormal tetrad $e_A^\alpha(\tau)$ along the trajectory, which defines a “body-fixed” frame.<fnref n="99" /> The rotation tensor $\omega^{\alpha\beta}$ associated with the tetrad is defined by

$$ \frac{\mathrm{D}e_A^\alpha}{\mathrm{d}\tau} = -\omega^{\alpha\beta} e_{A\beta}, \qquad (564) $$

where $D/\mathrm{d}\tau \equiv u^\beta \nabla_\beta$ is the covariant derivative with respect to the parameter $\tau$ along the worldline; equivalently, we have

$$ \omega^{\alpha\beta} = e^{A\alpha} \frac{\mathrm{D}e_A^\beta}{\mathrm{d}\tau}. \qquad (565) $$

Because of the normalization of the tetrad the rotation tensor is antisymmetric: $\omega^{\alpha\beta} = -\omega^{\beta\alpha}$.

We look for an action principle for the spinning particle. Following Hanson and Regge (1974), Porto (2006) we require the following symmetries to hold:

1. The action is a covariant scalar, i.e., behaves as a scalar with respect to general space-time diffeomorphisms;

2. It is a global Lorentz scalar, i.e., stays invariant under an arbitrary change of the tetrad vectors: $e_A^\alpha(\tau) \to \Lambda_A^B e_B^\alpha(\tau)$ where $\Lambda_A^B$ is a constant Lorentz matrix;

3. It is parametrization-invariant, i.e., its form is independent of the parameter $\tau$ used to follow the particle's worldline.<fnref n="100" />

In addition to these symmetries we need to specify the dynamical degrees of freedom: These are chosen to be the particle's position $y^\alpha$ and the tetrad $e_A^\alpha$. Furthermore we restrict ourselves to a Lagrangian depending only on the four-velocity $u^\alpha$, the rotation tensor $\omega^{\alpha\beta}$, and the metric $g_{\alpha\beta}$. Thus, the postulated action is of the type

$$ S[y^\alpha, e_A^\alpha] = \int_{-\infty}^{+\infty} \mathrm{d}\tau L(u^\alpha, \omega^{\alpha\beta}, g_{\alpha\beta}). \qquad (566) $$

These assumptions confine the formalism to a “pole-dipole” model and to terms

$$ \mathrm{d}\tau = \mathrm{d}\bar{\tau}\sqrt{-\bar{u}^2}, \qquad u^\alpha = \frac{\bar{u}^\alpha}{\sqrt{-\bar{u}^2}}, \qquad L = \frac{\bar{L}}{\sqrt{-\bar{u}^2}}, $$

where $\bar{u}^\alpha = dy^\alpha/\mathrm{d}\bar{\tau}$ and $\bar{u}^2 = g_{\alpha\beta}\bar{u}^\alpha\bar{u}^\beta$. Then the action $S = \int d\bar{\tau}\bar{L}$ will be in parametrization-invariant form.

<footnote n="99">⁹⁹ The tetrad is orthonormal in the sense that $g_{\alpha\beta} e_A^\alpha e_B^\beta = \eta_{AB}$, where $\eta_{AB} = \mathrm{diag}(-1, 1, 1, 1)$ denotes a Minkowski metric. The indices $AB \cdots = 0, 1, 2, 3$ are the internal Lorentz indices, while as usual $\alpha\beta \cdots \mu\nu \cdots = 0, 1, 2, 3$ are the space-time covariant indices. The inverse dual tetrad $e_A^A$, defined by $e_A^\beta e_A^\alpha = \delta_\alpha^\beta$, satisfies $\eta_{AB} e_A^A e_B^B = g_{\alpha\beta}$. We have also the completeness relation $e_A^\beta e_\beta^\alpha = \delta_A^\beta$.</footnote>

<footnote n="100">¹⁰⁰ This condition can always be realized. Suppose that $\tau$ initially denotes the proper time so that $g_{\alpha\beta}u^\alpha u^\beta = -1$. For an arbitrary parametrization $\bar{\tau}$ we pose</footnote>
linear in the spins. An important point is that such a model is *universal* in the sense that it can be used for black holes as well as neutrons stars. Indeed, the internal structure of the spinning body appears only at the quadratic order in the spins, through the rotationally induced quadrupole moment.

As it is written in (566), i.e., depending only on Lorentz scalars, $L$ is automatically a Lorentz scalar. By performing an infinitesimal coordinate transformation, one easily sees that the requirement that the Lagrangian be a covariant scalar specifies its dependence on the metric to be such that (see e.g. Bailey and Israel 1975)

$$2 \frac{\partial L}{\partial g_{\alpha\beta}} = p^{\alpha} u^{\beta} + S_{\gamma}^{\alpha} \omega^{\beta\gamma}. \qquad (567)$$

We have defined the conjugate linear momentum $p^\alpha$ and the antisymmetric spin tensor $S^{\alpha\beta}$ by the partial derivatives

$$p_{\alpha} \equiv \left. \frac{\partial L}{\partial u^{\alpha}} \right|_{\omega,g}, \qquad S_{\alpha\beta} \equiv 2 \left. \frac{\partial L}{\partial \omega^{\alpha\beta}} \right|_{u,g}. \qquad (568)$$

Note that the right-hand side of Eq. (567) is necessarily symmetric by exchange of the indices $\alpha$ and $\beta$. Finally, imposing the invariance of the action (566) by reparametrization of the worldline, we find that the Lagrangian must be a homogeneous function of degree one in the velocity $u^\alpha$ and rotation tensor $\omega^{\alpha\beta}$. Applying Euler's theorem to the function $L(u^\alpha, \omega^{\alpha\beta})$ immediately gives

$$L = p_{\alpha} u^{\alpha} + \frac{1}{2} S_{\alpha\beta} \omega^{\alpha\beta}, \qquad (569)$$

where the functions $p_\alpha(u, \omega)$ and $S_{\alpha\beta}(u, \omega)$ must be parametrization invariant. Note that, at this stage, their explicit expressions are not known. They will be specified only when a spin supplementary condition is imposed, see (576).

We now investigate the unconstrained variations of the action (566) with respect to the dynamical variables $e_A^\alpha$, $y^\alpha$ and the metric. First, we vary it with respect to the tetrad $e_A^\alpha$ while keeping the position $y^\alpha$ fixed. A worry is that we must have a way to distinguish intrinsic variations of the tetrad from variations which are induced by a change of the metric $g_{\alpha\beta}$. This is conveniently solved by decomposing the variation $\delta e_A^\beta$ according to

$$\delta e_A^\beta = e_{A\alpha} \left( \delta\theta^{\alpha\beta} + \frac{1}{2} \delta g^{\alpha\beta} \right), \qquad (570)$$

in which we have introduced the antisymmetric tensor $\delta\theta^{\alpha\beta} \equiv e^{A[\alpha} \delta e_A^{\beta]}$, and where the corresponding symmetric part is simply given by the variation of the metric, i.e. $e^{A[\alpha} \delta e_A^{\beta]} \equiv \frac{1}{2} \delta g^{\alpha\beta}$. Then we can consider the independent variations $\delta\theta^{\alpha\beta}$ and $\delta g^{\alpha\beta}$. Varying with respect to $\delta\theta^{\alpha\beta}$, but holding the metric fixed, gives the equation of spin precession which is found to be
$$ \frac{DS_{\alpha\beta}}{d\tau} = \omega_{\alpha}^{\gamma} S_{\beta\gamma} - \omega_{\beta}^{\gamma} S_{\alpha\gamma}, \qquad (571) $$

or, alternatively, using the fact that the right-hand side of (567) is symmetric,

$$ \frac{DS_{\alpha\beta}}{d\tau} = p_{\alpha}u_{\beta} - p_{\beta}u_{\alpha}. \qquad (572) $$

We next vary with respect to the particle's position $y^{\alpha}$ while holding the tetrad $e_A^{\alpha}$ fixed. Operationally, this means that we have to parallel-transport the tetrad along the displacement vector, i.e., to impose

$$ \delta y^{\beta} \nabla_{\beta} e_{A}^{\alpha} = 0. \qquad (573) $$

A simple way to derive the result is to use locally inertial coordinates, such that the Christoffel symbols $\Gamma^{\alpha}_{\beta\gamma} = 0$ along the particle's worldline $y^{\alpha}(\tau)$; then, Eq. (573) gives $\delta e_{A}^{\alpha} = \delta y^{\beta} \partial_{\beta} e_{A}^{\alpha} = -\delta y^{\beta} \Gamma^{\alpha}_{\beta\gamma} e_{A}^{\gamma} = 0$. The variation leads then to the well-known (Mathisson 1937, 2010; Papapetrou 1951b; Corinaldesi and Papapetrou 1951) equation of motion

$$ \frac{Dp_{\alpha}}{d\tau} = -\frac{1}{2} u^{\beta} R_{\alpha\beta\mu\nu} S^{\mu\nu}, \qquad (574) $$

which involves the famous coupling of the spin tensor to the Riemann curvature.<fnref n="101" /> With a little more work, the equation of motion (574) can also be derived using an arbitrary coordinate system, making use of the parallel transport equation (573).

Finally, varying with respect to the metric while keeping $\delta\theta^{\alpha\beta} = 0$, gives the stress-energy tensor of the spinning particle. We must again take into account the scalarity of the action, as imposed by Eq. (567). We obtain the standard pole-dipole result (Tulczyjew 1957, 1959; Mathisson 1937, 2010; Papapetrou 1951b; Corinaldesi and Papapetrou 1951; Trautman 2002; Dixon 1979):

$$ T^{\alpha\beta} = \int_{-\infty}^{+\infty} d\tau p^{(\alpha} u^{\beta)} \frac{\delta_{(4}}(x-y)}{\sqrt{-g}} - \nabla_{\gamma} \int_{-\infty}^{+\infty} d\tau S^{\gamma(\alpha} u^{\beta)} \frac{\delta_{(4}}(x-y)}{\sqrt{-g}}, \qquad (575) $$

where $\delta_{(4}}(x-y)$ denotes the four-dimensional Dirac function. It can easily be checked that the covariant conservation law $\nabla_{\beta} T^{\alpha\beta} = 0$ holds as a consequence of the equation of motion (574) and the equation of spin precession (572).

Up to now we have considered unconstrained variations of the action (566), describing the particle's internal degrees of freedom by the six independent components of the tetrad $e_A^{\alpha}$ (namely a $4 \times 4$ matrix subject to the 10 constraints $g_{\alpha\beta} e_A^{\alpha} e_B^{\beta} = \eta_{AB}$). To correctly account for the number of degrees of freedom associated with the spin, we must impose three *supplementary spin conditions* (SSC).

Several choices are possible for a sensible SSC. Notice that in the case of extended bodies the choice of a SSC corresponds to the choice of a central worldline inside the

<footnote n="101">¹⁰¹ Our conventions for the Riemann tensor $R_{\alpha\beta\mu\nu}$ follow those of Misner et al. (1973).</footnote>
body with respect to which the spin angular momentum is defined (see Kidder 1995 for a discussion). Here we adopt the Tulczyjew (1957, 1959) covariant SSC

$$S^{\alpha\beta}p_{\beta} = 0. \qquad (576)$$

As shown by Hanson and Regge (1974) in the flat space-time case, it is possible to specify the Lagrangian in our original action (566) in such a way that the constraints (576) are directly the consequence of the equations derived from that Lagrangian. Here, for simplicity's sake, we shall simply impose the constraints (576) in the space of solutions of the Euler-Lagrange equations. From Eq. (576) we can introduce the covariant spin vector $S_μ$ associated with the spin tensor by<fnref n="102" />

$$S^{\alpha\beta} \equiv \frac{1}{m} \varepsilon^{\alpha\beta\mu\nu} p_{\mu} S_{\nu}, \qquad (577)$$

where we have defined the mass of the particle by $m^2 \equiv -g^{\mu\nu}p_\mu p_\nu$. By contracting Eq. (572) with $p^β$ and using the equation of motion (574), one obtains

$$p_{\alpha}(pu) + m^2 u_{\alpha} = \frac{1}{2} u^{\nu} R_{\gamma\mu\nu}^{\beta} S^{\mu\nu} S_{\alpha\beta}, \qquad (578)$$

where we denote $(pu) \equiv p_\mu u^\mu$. By further contracting Eq. (578) with $u^\alpha$ we obtain an explicit expression for $(pu)$, which can then be substituted back into (578) to provide the relation linking the four-momentum $p_\alpha$ to the four-velocity $u_\alpha$. It can be checked using (576) and (578) that the mass of the particle is constant along the particle's trajectory: $dm/d\tau = 0$. Furthermore the four-dimensional magnitude $s$ of the spin defined by $s^2 \equiv g^{\mu\nu}S_\mu S_\nu$ is also conserved: $ds/d\tau = 0$.

Henceforth we shall restrict our attention to SO interactions, linear in the spins (see Sect. 3.6.5 for SS interactions). We shall also adopt for the parameter $\tau$ along the particle's worldline the proper time $d\tau \equiv \sqrt{-g_{\mu\nu}}dy^\mu dy^\nu$, so that $g_{\mu\nu}u^\mu u^\nu = -1$. Neglecting quadratic spin-spin (SS) and higher-order interactions, the linear momentum is simply proportional to the normalized four-velocity: $p_\alpha = m u_\alpha + O(S^2)$. Hence, from Eq. (572) we deduce that $DS_{\alpha\beta}/d\tau = O(S^2)$. The equation for the spin covariant vector $S_\alpha$ then reduces at linear order to

$$\frac{DS_{\alpha}}{d\tau} = O(S^2). \qquad (579)$$

Thus the spin covector is parallel transported along the particle's trajectory at linear order in spin. For the spin covector $S_\alpha$ itself, we choose a four-vector which is purely spatial in the particle's instantaneous rest frame, where $u^\alpha = (1, \mathbf{0})$, hence the components of $S_\alpha$ are $(0, S)$ in that frame. Therefore, in any frame,

$$S_{\alpha} u^{\alpha} = 0. \qquad (580)$$

In applications (viz. the construction of gravitational-wave templates for the compact

<footnote n="102">¹⁰² The four-dimensional Levi-Civita tensor is defined by $\epsilon_{\alpha\beta\mu\nu} \equiv \sqrt{-g} \epsilon_{\alpha\beta\mu\nu}$ and $\epsilon^{\alpha\beta\mu\nu} \equiv -\epsilon^{\alpha\beta\mu\nu}/\sqrt{-g}$; here $\epsilon_{\alpha\beta\mu\nu} \equiv \epsilon^{\alpha\beta\mu\nu}$ denotes the completely anti-symmetric Levi-Civita symbol such that $\epsilon_{0123} = \epsilon^{0123} = 1$. For convenience we often pose $c = 1$.</footnote>
binary inspiral) it is very useful to introduce new spin variables that are designed to have a conserved three-dimensional Euclidean norm (numerically equal to s). Using conserved-norm spin vector variables is indeed the most natural choice when considering the dynamics of compact binaries reduced to the frame of the center of mass or to circular orbits (Bohé et al. 2013b). Indeed the evolution equations of such spin variables reduce, by construction, to ordinary precession equations, and these variables are secularly constant (see Will 2005).

A standard, general procedure to define a (Euclidean) conserved-norm spin spatial vector consists of projecting the spin covector $S_\alpha$ onto an orthonormal tetrad $e_A^\alpha$, which leads to the four scalar components ($A = 0, 1, 2, 3$)

$$S_A = e_A^\alpha S_\alpha. \qquad (581)$$

If we choose for the time-like tetrad vector the four-velocity itself, $e_0^\alpha = u^\alpha$,<fnref n="103" /> the time component tetrad projection $S_0$ vanishes because of the orthogonality condition (580). We have seen that $S_\alpha S^\alpha = s^2$ is conserved along the trajectory; because of (580) we can rewrite this as $\gamma^{\alpha\beta}S_\alpha S_\beta = s^2$, in which we have introduced the projector $\gamma^{\alpha\beta} = g^{\alpha\beta} + u^\alpha u^\beta$ onto the spatial hypersurface orthogonal to $u^\alpha$. From the orthonormality of the tetrad and our choice $e_0^\alpha = u^\alpha$, we have $\gamma^{\alpha\beta} = \delta^{ab}e_a^\alpha e_b^\beta$ in which $a, b = 1, 2, 3$ refer to the spatial values of the tetrad indices, i.e., $A = (0, a)$ and $B = (0, b)$. Therefore the conservation law $\gamma^{\alpha\beta}S_\alpha S_\beta = s^2$ becomes

$$\delta^{ab}S_a S_b = s^2, \qquad (582)$$

which is indeed the relation defining a Euclidean conserved-norm spin variable $S_a$.<fnref n="104" /> However, note that the choice of the spin variable $S_a$ is still somewhat arbitrary, since a rotation of the tetrad vectors can freely be performed. We refer to Damour et al. (2008), Bohé et al. (2013b) for the definition of some “canonical” choice for the tetrad in order to fix this residual freedom. Such choice presents the advantage of providing a unique determination of the conserved-norm spin variable in a given gauge. This canonical choice will be the one adopted in all explicit results presented in Sect. 3.6.4.

The evolution equation (579) for the original spin variable $S_\alpha$ now translates into an ordinary precession equation for the tetrad components $S_a$:

$$\frac{dS_a}{dt} = \Omega_a^b S_b, \qquad (583)$$

where the precession tensor $\Omega_{ab}$ is related to the tetrad components $\omega_{AB}$ of the rotation tensor defined in (565) by $\Omega_{ab} = z \omega_{ab}$ where we pose $z \equiv d\tau/dt$, remembering the redshift variable (427). The antisymmetric character of the matrix $\Omega_{ab}$ guarantees that $S_a$ satisfies the Euclidean precession equation

<footnote n="103">¹⁰³ Because of this choice, it is better to consider that the tetrad is not the same as the one we originally employed to construct the action (566).</footnote>

<footnote n="104">¹⁰⁴ Beware that here we employ the usual slight ambiguity in the notation when using the same carrier letter $S$ to denote the tetrad components (581) and the original spin covector. Thus, $S_a$ should not be confused with the spatial components $S_i$ (with $i = 1, 2, 3$) of the covariant vector $S_\alpha$.</footnote>
$$ \frac{d\mathbf{S}}{dt} = \mathbf{\Omega} \times \mathbf{S}, \qquad (584) $$

where we denote $\mathbf{S} = (S_a)$, and $\mathbf{\Omega} = (\Omega_a)$ with $\Omega_a = -\frac{1}{2}\epsilon_{abc}\mathbf{\Omega}^{bc}$. As a consequence of (584) the spin has a conserved Euclidean norm: $\mathbf{S}^2 = s^2$. From now on we shall no longer make any distinction between the spatial tetrad indices $ab \cdots$ and the ordinary spatial indices $ij \cdots$ which are raised and lowered with the Kronecker metric. Explicit results for the equations of motion and gravitational wave templates will be given in Sects. 3.6.2 and 3.6.4 using the canonical choice for the conserved-norm spin variable $\mathbf{S}$.

### 3.6.2 Equations of motion and precession for circular binaries

The previous formalism can be generalized to self-gravitating systems consisting of two spinning point particles. The metric generated by the system of particles, interacting only through gravitation, is solution of the Einstein field equations (47) with stress-energy tensor given by the sum of the individual stress-energy tensors (575) for each particles. The equations of motion of the particles are the Mathisson-Papapetrou equations (574) with the metric generated by the self-gravitating system evaluated at the location of the particles thanks to a regularization procedure (see Sect. 3.1). The precession equations of each of the spins are given by

$$ \frac{d\mathbf{S}_a}{dt} = \mathbf{\Omega}_a \times \mathbf{S}_a, \qquad (585) $$

where $a = 1, 2$ labels the particles. The spin variables $\mathbf{S}_a$ are the conserved-norm spins defined in Sect. 3.6.1. In the following it is convenient to introduce two combinations of the individuals spins defined by$^{105}$

$$ \mathbf{S} \equiv \mathbf{S}_1 + \mathbf{S}_2, \qquad (586a) $$

$$ \mathbf{\Sigma} \equiv \frac{\mathbf{S}_2}{X_2} - \frac{\mathbf{S}_1}{X_1}. \qquad (586b) $$

We shall investigate the case where the binary's orbit is quasi-circular, i.e., whose radius is constant apart from small perturbations induced by the spins (as usual we neglect the gravitational radiation damping effects). We denote by $\mathbf{x} = \mathbf{y}_1 - \mathbf{y}_2$ and $\mathbf{v} = d\mathbf{x}/dt$ the relative position and velocity.$^{106}$ We introduce an orthonormal moving triad $\{\mathbf{n}, \mathbf{\lambda}, \mathbf{\ell}\}$ defined by the unit separation vector $\mathbf{n} = \mathbf{x}/r$ (with $r = |\mathbf{x}|$) and the

$$ \mathbf{S}_1 = X_1 \mathbf{S} - \mathbf{v} \mathbf{\Sigma}, $$

$$ \mathbf{S}_2 = X_2 \mathbf{S} + \mathbf{v} \mathbf{\Sigma}. $$

<footnote>$^{105}$ This convenient notation is adopted by Kidder (1995); the inverse formulas read (recall $X_a \equiv m_a/m$ and $\mathbf{v} \equiv X_1\mathbf{X}_2$)</footnote>

<footnote>$^{106}$ Note that the individual particle's positions $\mathbf{y}_a$ in the frame of the center-of-mass (defined by the cancellation of the center-of-mass integral: $G = 0$) are related to $\mathbf{x}$ and $\mathbf{v}$ by expressions similar to (354), but containing spin effects starting at the order 1.5PN, see Sect. 3 of Bohé et al. (2013b).</footnote>
<figure><img src="image_10.png" /><figcaption>Fig. 9 Geometric definitions for the precessional motion of spinning compact binaries. We show (i) the source frame defined by the fixed orthonormal basis {$x, y, z$}; (ii) the instantaneous orbital plane which is described by the orthonormal basis {$x_\ell, y_\ell, \ell$}; (iii) the moving triad {$n, \lambda, \ell$} and the associated three Euler angles $\alpha, \iota$ and $\Phi$; (iv) the direction of the total angular momentum $J$ which coincides with the $z$-direction. Dashed lines show projections into the $x-y$ plane</figcaption></figure>

unit normal $\ell$ to the instantaneous orbital plane given by $\ell = n \times v / |n \times v|$; the orthonormal triad is then completed by $\lambda = \ell \times n$. Those vectors are represented on Fig. 9, which shows the geometry of the system. The orbital frequency $\Omega$ is defined for general orbits, not necessarily circular, by $v = \dot{r}n + r\Omega\lambda$ where $\dot{r} = n \cdot v$ is the derivative of $r$ with respect to the coordinate time $t$. The general expression for the relative acceleration $a \equiv dv/dt$ decomposed in the moving basis {$n, \lambda, \ell$} is

$$ a = (\ddot{r} - r\Omega^2)n + (r\dot{\Omega} + 2\dot{r}\Omega)\lambda + r\varpi\Omega\ell. \quad (587) $$

Here we have introduced the orbital plane precession $\varpi$ of the orbit defined by $\varpi \equiv -\lambda \cdot d\ell/dt$. Next we impose the restriction to quasi-circular precessing orbits which is defined by the conditions $\dot{r} = \dot{\Omega} = \mathcal{O}(1/c^5)$ and $\ddot{r} = \mathcal{O}(1/c^{10})$ so that $v^2 = r^2\Omega^2 + \mathcal{O}(1/c^{10})$; see Eqs. (372). Then $\lambda$ represents the direction of the velocity, and the precession frequency $\varpi$ is proportional to the variation of $\ell$ in the direction of the velocity. In this way we find that the equations of the relative motion in the frame of the center-of-mass are
$$ \mathbf{a} = -r \Omega^2 \mathbf{n} + r \bar{\omega} \Omega \ell + \mathcal{O}\left(\frac{1}{c^5}\right). \qquad (588) $$

Since we neglect the radiation reaction damping there is no component of the acceleration along $\lambda$. This equation represents the generalization of Eq. (371) for spinning quasi-circular binaries with no radiation reaction. The orbital frequency $\Omega$ will contain spin effects in addition to the non-spin terms given by (373), while the precessional frequency $\bar{\omega}$ will entirely be due to spins.

Here we report the results for the spin-orbit (SO) contributions into these quantities at the next-to-next-to-leading level corresponding to 3.5PN order (Marsat et al. 2013b; Bohé et al. 2013b). We project out the spins on the moving orthonormal basis, defining $\mathbf{S} = S_n\mathbf{n} + S_\lambda\mathbf{\lambda} + S_\ell\mathbf{\ell}$ and similarly for $\mathbf{\Sigma}$. We have

$$ \begin{aligned} \Omega_{\text{SO}}^2 = \frac{\gamma^{3/2}}{m r^3} \Biggl\{ & -5S_{\ell} - 3\Delta\Sigma_{\ell} + \gamma \left[ \left(\frac{45}{2} - \frac{27}{2}v\right)S_{\ell} + \Delta\left(\frac{27}{2} - \frac{13}{2}v\right)\Sigma_{\ell} \right] \\ & + \gamma^2 \left[ \left(-\frac{495}{8} - \frac{561}{8}v - \frac{51}{8}v^2\right)S_{\ell} + \Delta\left(-\frac{297}{8} - \frac{341}{8}v - \frac{21}{8}v^2\right)\Sigma_{\ell} \right] \Biggr\} \\ & + \mathcal{O}\left(\frac{1}{c^8}\right), \end{aligned} \qquad (589) $$

which has to be added to the non-spin terms (373) up to 3.5PN order. We recall that the ordering post-Newtonian parameter is $\gamma = \frac{Gm}{rc^2}$. The next-to-next-to-leading SO effects into the precessional frequency read

$$ \begin{aligned} \bar{\omega} = \frac{c^3 x^3}{G^2 m^3} \Biggl\{ & 7S_n + 3\Delta\Sigma_n + x \left[ (-3 - 12v)S_n + \Delta\left(-3 - \frac{11}{2}v\right)\Sigma_n \right] \\ & + x^2 \left[ \left(-\frac{3}{2} - \frac{59}{2}v + 9v^2\right)S_n + \Delta\left(-\frac{3}{2} - \frac{77}{8}v + \frac{13}{3}v^2\right)\Sigma_n \right] \Biggr\} \end{aligned} \qquad (590) \\ & + \mathcal{O}\left(\frac{1}{c^8}\right), \end{aligned} $$

where this time the ordering post-Newtonian parameter is $x \equiv (\frac{Gm\Omega}{c^3})^{2/3}$. In order to complete the evolution equations for quasi-circular orbits we need also the precession vectors $\Omega_a$ of the two spins as defined by Eq. (585). These are given by

$$ \begin{aligned} \Omega_1 = \frac{c^3 x^{5/2}}{Gm} \ell \Biggl\{ & \left[ \frac{3}{4} + \frac{1}{2}v - \frac{3}{4}\Delta + x \left[ \frac{9}{16} + \frac{5}{4}v - \frac{1}{24}v^2 + \Delta \left( -\frac{9}{16} + \frac{5}{8}v \right) \right] \right] \\ & + x^2 \left[ \frac{27}{32} + \frac{3}{16}v - \frac{105}{32}v^2 - \frac{1}{48}v^3 + \Delta \left( -\frac{27}{32} + \frac{39}{8}v - \frac{5}{32}v^2 \right) \right] \Biggr\} \end{aligned} \qquad (591) \\ & + \mathcal{O}\left(\frac{1}{c^8}\right). \end{aligned} $$
We obtain $\Omega_2$ by exchanging the masses, $\Delta \to -\Delta$. At the linear SO level the precession vectors $\Omega_a$ are independent of the spins.<fnref n="107" /> Finally we give the SO contributions to next-to-next-to-leading order in the binary's conservative energy for circular orbits as (Marsat et al. 2013b; Bohé et al. 2013b)

$$ E_{\text{SO}} = -\frac{\nu c^2 x^{5/2}}{2Gm} \left\{ \frac{14}{3} S_{\ell} + 2\Delta\Sigma_{\ell} \right. \\ \qquad + x \left[ \left( 11 - \frac{61}{9}v \right) S_{\ell} + \Delta \left( 3 - \frac{10}{3}v \right) \Sigma_{\ell} \right] \\ \qquad + x^2 \left[ \left( \frac{135}{4} - \frac{367}{4}v + \frac{29}{12}v^2 \right) S_{\ell} + \Delta \left( \frac{27}{4} - 39v + \frac{5}{4}v^2 \right) \Sigma_{\ell} \right] \\ \qquad \left. + O\left(\frac{1}{c^8}\right). \right\} \quad (592) $$

We recall that the non-spin contributions in E were given in Eq. (377).

### 3.6.3 Orbital precession at the spin-orbit level

We now investigate an analytical solution for the dynamics of compact spinning binaries on quasi-circular orbits, including the effects of spin precession (Blanchet et al. 2011; Marsat et al. 2013a). This solution will be valid whenever the radiation reaction effects can be neglected, and is restricted to the linear SO level.

In the following, we will extensively employ the total angular momentum of the system, that we denote by **J**, and which is conserved when radiation-reaction effects are neglected,

$$ \frac{d\mathbf{J}}{dt} = 0. \quad (593) $$

It is customary to decompose the conserved total angular momentum **J** as the sum of the orbital angular momentum **L** and of the two spins,<fnref n="108" />

$$ \mathbf{J} = \mathbf{L} + \frac{\mathbf{S}}{c}. \quad (594) $$

This split between **L** and **S** = **S**₁ + **S**₂ is specified by our choice of spin variables, here the conserved-norm spins defined in Sect. 3.6.1. Although **L** is called the “orbital” angular momentum, it actually includes both non-spin and spin contributions. We refer to Eq. (4.7) of Bohé et al. (2013b) for the expression of **L** at the next-to-next-to-leading SO level for quasi-circular orbits.

<footnote n="107">¹⁰⁷ Beware of our inevitably slightly confusing notation: $\Omega$ is the binary's orbital frequency and $\Omega_{\text{SO}}$ refers to the spin-orbit terms therein; $\Omega_a$ is the precession frequency of the a-th spin while $\tau$ is the precession frequency of the orbital plane; and $\omega_a$ defined earlier in Eqs. (395) and (436) is the rotation frequency of the a-th black hole. Such different notions nicely mix up in the first law of spinning binary black holes in Sect. 3.3.3; see Eq. (434) and the corotation condition (437).</footnote>

<footnote n="108">¹⁰⁸ Recall from Eq. (563) that in our convention the spins have the dimension of an angular momentum times c.</footnote>
Our solution will consist of some explicit expressions for the moving triad {$n, \lambda, \ell$} at the SO level in the conservative dynamics for quasi-circular orbits. With the previous definitions of the orbital frequency $\Omega$ and the precessional frequency $\omega$ we have the following system of equations for the time evolution of the triad vectors,

$$ \frac{dn}{dt} = \Omega \lambda, \qquad (595a) $$

$$ \frac{d\lambda}{dt} = -\Omega n + \omega \ell, \qquad (595b) $$

$$ \frac{d\ell}{dt} = -\omega \lambda. \qquad (595c) $$

Equivalently, introducing the orbital rotation vector $\Omega \equiv \Omega \ell$ and spin precession vector $\omega \equiv \omega n$, these equations can be elegantly written as

$$ \frac{dn}{dt} = \Omega \times n, \qquad (596a) $$

$$ \frac{d\lambda}{dt} = (\Omega + \omega) \times \lambda, \qquad (596b) $$

$$ \frac{d\ell}{dt} = \omega \times \ell. \qquad (596c) $$

Next we introduce a fixed (inertial) orthonormal basis {$x, y, z$} as follows: $z$ is defined to be the normalized value $J/J$ of the total angular momentum $J$; $y$ is orthogonal to the plane spanned by $z$ and the direction $N = X/R$ of the detector as seen from the source (notation of Sect. 2.4.1) and is defined by $y = z \times N/|z \times N|$; and $x$ completes the triad—see Fig. 9. Then, we introduce the standard spherical coordinates $(\alpha, \iota)$ of the vector $\ell$ measured in the inertial basis {$x, y, z$}. Since $\iota$ is the angle between the total and orbital angular momenta, we have

$$ \sin \iota = \frac{|J \times \ell|}{J}, \qquad (597) $$

where $J \equiv |J|$. The angles $(\alpha, \iota)$ are referred to as the *precession angles*.

We now derive from the time evolution of our triad vectors those of the precession angles $(\alpha, \iota)$, and of an appropriate phase $\Phi$ that specifies the position of $n$ with respect to some reference direction in the orbital plane denoted $\mathbf{x}_\ell$. Following Arun et al. (2009b), we pose

$$ \mathbf{x}_{\ell} = \frac{\mathbf{J} \times \ell}{|\mathbf{J} \times \ell|}, \qquad (598a) $$

$$ \mathbf{y}_{\ell} = \ell \times \mathbf{x}_{\ell}, \qquad (598b) $$

such that {$\mathbf{x}_\ell, \mathbf{y}_\ell, \ell$} forms an orthonormal basis. The motion takes place in the
instantaneous orbital plane spanned by **n** and **λ**, and the phase angle **Φ** is such that
(see Fig. 9):

$$
\mathbf{n} = \cos \Phi \mathbf{x}_{\ell} + \sin \Phi \mathbf{y}_{\ell}, \qquad (599a)
$$

$$
\boldsymbol{\lambda} = -\sin \Phi \boldsymbol{x}_{\ell} + \cos \Phi \boldsymbol{y}_{\ell}, \qquad (599\text{b})
$$

from which we deduce

$$
e^{-i\Phi} = \mathbf{x}_{\ell} \cdot (\mathbf{n} + i\mathbf{\lambda}) = \frac{J_{\lambda} - iJ_{n}}{\sqrt{J_{n}^{2} + J_{\lambda}^{2}}}. \quad (600)
$$

Combining Eqs. (600) with (597) we also get

$$
\sin t e^{-i\Phi} = \frac{J_{\lambda} - iJ_{n}}{J}. \qquad (601)
$$

By identifying the right-hand sides of (595) with the time-derivatives of the relations
(599) we obtain the following system of equations for the variations of $\alpha$, $t$ and $\Phi$,

$$
\frac{d\alpha}{dt} = \omega \frac{\sin \Phi}{\sin t}, \qquad (602a)
$$

$$
\frac{d t}{d t}=\omega \cos \Phi,
\qquad(602 \mathrm{~b})
$$

$$
\frac{d\Phi}{dt} = \Omega - \omega \frac{\sin \Phi}{\tan t}. \qquad (602c)
$$

On the other hand, using the decompositon of the total angular momentum (594)
together with the fact that the components of **L** projected along **n** and **λ** are of the
order O(S), Eq. (4.7) of Bohé et al. (2013b), we deduce that sin t is itself a small
quantity of order O(S). Since we also have ϖ = O(S), we conclude by direct
integration of the sum of Eqs. (602a) and (602c) that

$$
\Phi + \alpha = \phi + \mathcal{O}(S^2), \tag{603}
$$

in which we have defined the “carrier” phase as

$$
\phi \equiv \int \Omega \, dt = \Omega(t - t_0) + \phi_0, \qquad (604)
$$

with $\phi_0$ the value of the carrier phase at some arbitrary initial time $t_0$. An important
point we have used when integrating (604) is that the orbital frequency $\Omega$ is constant
at linear order in the spins. Indeed, from Eq. (589) we see that only the components
of the conserved-norm spin vectors along $\ell$ can contribute to $\Omega$ at linear order. As we
show in Eq. (607c) below, these components are in fact constant at linear order in
spins. Thus we can treat $\Omega$ as a constant for our purpose.
The combination $\Phi + \alpha$ being known by Eq. (603), we can further express the precession angles $\iota$ and $\alpha$ at linear order in spins in terms of the components $J_n$ and $J_\lambda$; from Eqs. (597) and (601),

$$ \sin \iota = \frac{\sqrt{J_n^2 + J_\lambda^2}}{L_{NS}} + \mathcal{O}(S^2), \qquad (605a) $$

$$ e^{i\alpha} = \frac{J_{\lambda} - iJ_{n}}{\sqrt{J_{n}^{2} + J_{\lambda}^{2}}} e^{i\phi} + \mathcal{O}(S^{2}), \qquad (605b) $$

where we denote by $L_{NS}$ the norm of the non-spin (NS) part of the orbital angular momentum $L$.

It remains to obtain the explicit time variation of the components of the two individual spins $S_n^a$, $S_\lambda^a$ and $S_\ell^a$ (with $a = 1, 2$). Using Eqs. (605) together with the decomposition (594) and the explicit expression of $L$ (see Bohé et al. 2013b), we shall then be able to obtain the explicit time variation of the precession angles $(\alpha, \iota)$ and phase $\Phi$. Combining (585) and (595) we obtain

$$ \frac{dS_n^a}{dt} = (\Omega - \Omega_a)S_\lambda^a, \qquad (606a) $$

$$ \frac{dS_\lambda^a}{dt} = -(\Omega - \Omega_a)S_n^a + \omega S_\ell^a, \qquad (606b) $$

$$ \frac{dS_\ell^a}{dt} = -\omega S_\lambda^a, \qquad (606c) $$

where $\Omega_a$ is the norm of the precession vector of the a-th spin as given by (591), and the precession frequency $\omega$ is explicitly given by (590). At linear order in spins these equations translate into

$$ \frac{dS_n^a}{dt} = (\Omega - \Omega_a)S_\lambda^a, \qquad (607a) $$

$$ \frac{dS_\lambda^a}{dt} = -(\Omega - \Omega_a)S_n^a + \mathcal{O}(S^2), \qquad (607b) $$

$$ \frac{dS_\ell^a}{dt} = \mathcal{O}(S^2). \qquad (607c) $$

We see that, as stated before, the spin components along $\ell$ are constant, and so is the orbital frequency $\Omega$ given by (589). At the linear SO level, the equations (607) can be decoupled and integrated as

$$ S_n^a = S_\perp^a \cos \psi_a, \qquad (608a) $$

$$ S_\lambda^a = -S_\perp^a \sin \psi_a, \qquad (608b) $$
$$S_{\ell}^{a} = S_{\parallel}^{a}. \qquad (608c)$$

Here $S_{\perp}^a$ and $S_{\parallel}^a$ denote two quantities for each spins, that are constant up to terms $O(S^2)$. The phase of the projection perpendicular to the direction $\ell$ of each of the spins is given by

$$\psi_a = (\Omega - \Omega_a)(t - t_0) + \psi_a^0, \qquad (609)$$

where $\psi_a^0$ is the constant initial phase at the reference time $t_0$.

Finally we can give in an explicit way, to linear SO order, the moving triad $\{n, \lambda, \ell\}$ in terms of the reference triad $\{n_0, \lambda_0, \ell_0\}$ at the reference time $t_0$ in Eqs. (609) and (604). The best way to express the result is to introduce the complex null vector $m = \frac{1}{\sqrt{2}}(n + i\lambda)$ and its complex conjugate $\bar{m}$; the normalization is chosen so that $m \cdot \bar{m} = 1$. We obtain

$$m = e^{-i(\phi - \phi_0)} m_0 + \frac{i}{\sqrt{2}} (\sin \ell e^{i\alpha} - \sin \ell_0 e^{i\alpha_0}) e^{-i\phi} \ell_0 + O(S^2), \qquad (610a)$$

$$\ell = \ell_0 + \left[ \frac{i}{\sqrt{2}} (\sin \ell e^{-i\alpha} - \sin \ell_0 e^{-i\alpha_0}) e^{i\phi_0} m_0 + \text{c.c.} \right] + O(S^2). \qquad (610b)$$

The precession effects in the dynamical solution for the evolution of the basis vectors $\{n, \lambda, \ell\}$ are given by the second terms in these equations. They depend only in the combination $\sin \ell e^{i\alpha}$ and its complex conjugate $\sin \ell e^{-i\alpha}$, which follows from Eqs. (605) and the known spin and non-spin contributions to the total angular momentum $J$. One can check that precession effects in the above dynamical solution (610) for the moving triad start at order $O(1/c^3)$.

### 3.6.4 Spin-orbit effects in the gravitational-wave flux and orbital phase

Like in Sect. 3.4 our main task is to control up to high post-Newtonian order the mass and current radiative multipole moments $U_L$ and $V_L$ which parametrize the asymptotic waveform and gravitational fluxes far away from the source, cf. Equations (178)–(179). The radiative multipole moments are in turn related to the source multipole moments $I_L$ and $J_L$ through complicated relationships involving tails and related effects, and described in Sects. 2.4.5 and 2.4.6.

The source moments have been expressed in Eqs. (7) in terms of some source densities $\Sigma$, $\Sigma_i$ and $\Sigma_{ij}$ defined from the components of the post-Newtonian expansion of the pseudo-tensor, denoted $\bar{\tau}^{\alpha\beta}$. To lowest order the (PN expansion of the) pseudo-tensor reduces to the matter tensor $T^{\alpha\beta}$ which has compact support, and the source densities $\Sigma$, $\Sigma_i$, $\Sigma_{ij}$ reduce to the compact support quantities $\sigma$, $\sigma_i$, $\sigma_{ij}$ given by Eqs. (160). Now, computing spin effects, the matter tensor $T^{\alpha\beta}$ has been found to be given by (575) in the framework of the pole-dipole approximation suitable for SO couplings (and sufficient also for SS interactions between different spins). Here, to give a flavor of the computation, we present the lowest order SO contributions to the general mass and current source multipole moments ($\forall \ell \in \mathbb{N}$):
$$
I_L^{SO} = \frac{2\ell v}{c^3(\ell + 1)} \left\{ \begin{aligned}
& \ell [\sigma_{\ell}(v)(\mathbf{v} \times \mathbf{S})^{\langle \ell \rangle} - \sigma_{\ell+1}(v)(\mathbf{v} \times \mathbf{\Sigma})^{\langle \ell \rangle}] x^{L-1} \\
& - (\ell - 1) [\sigma_{\ell}(v)(\mathbf{x} \times \mathbf{S})^{\langle \ell \rangle} - \sigma_{\ell+1}(v)(\mathbf{x} \times \mathbf{\Sigma})^{\langle \ell \rangle}] v^{\ell-1} x^{L-2} \\
& + \mathcal{O}\left(\frac{1}{c^5}\right),
\end{aligned} \right\} \quad (611a)
$$

$$
J_L^{SO} = \frac{(\ell + 1)v}{2c} \left[ \sigma_{\ell-1}(v)S^{\langle \ell \rangle} - \sigma_{\ell}(v)\Sigma^{\langle \ell \rangle} \right] x^{L-1} + \mathcal{O}\left(\frac{1}{c^3}\right). \quad (611b)
$$

Paralleling the similar expressions (455) for the Newtonian approximation to the
source moments in the non-spin case, we posed $\sigma_{\ell}(v) \equiv X_2^{\ell-1} + (-)^{\ell} X_1^{\ell-1}$, see also
Eqs. (456). In Eqs. (611) we employ the notation (586) for the two spins and the
ordinary cross product $\times$ of Euclidean vectors. Thus, the dominant level of spins is at
the 1.5PN order in the mass-type moments $I_L$, but only at the 0.5PN order in the
current-type moments $J_L$. It is then evident that the spin part of the current-type
moments will always dominate over that of the mass-type moments. We refer to
Blanchet et al. (2006), Bohé et al. (2013a) for higher order post-Newtonian
expressions of the source moments. If we insert the expressions (611) into tail
integrals like (192), we find that some spin contributions originate from tails starting
at the 3PN order (Blanchet et al. 2011; Marsat et al. 2013a; the equivalent calculation
has been done with EFT techniques: Cho et al. 2022a).

We are left with the following high-order result for the SO contributions to the
gravitational-wave energy flux, with 4PN level precision (Blanchet et al. 2006; Bohé
et al. 2013a; Marsat et al. 2013a; Cho et al. 2022a):

$$
\begin{equation}
\begin{aligned}
\mathcal{F}_{\text{SO}} = \frac{32v^2 c^5 x^{13/2}}{5G^2 m^2} & \left\{ -4S_{\ell} - \frac{5}{4} \Delta \Sigma_{\ell} \right. \\
& \left. + x \left[ \left( -\frac{9}{2} + \frac{272}{9} v \right) S_{\ell} + \left( -\frac{13}{16} + \frac{43}{4} v \right) \Delta \Sigma_{\ell} \right] \right. \\
& \left. + x^{3/2} \left[ -16\pi S_{\ell} - \frac{31\pi}{6} \Delta \Sigma_{\ell} \right] \right. \\
& \left. + x^2 \left[ \left( \frac{476645}{6804} + \frac{6172}{189} v - \frac{2810}{27} v^2 \right) S_{\ell} + \left( \frac{9535}{336} + \frac{1849}{126} v - \frac{1501}{36} v^2 \right) \Delta \Sigma_{\ell} \right] \right. \\
& \left. + x^{5/2} \left[ \left( -\frac{3485\pi}{96} + \frac{13879\pi}{72} v \right) S_{\ell} + \left( -\frac{7163\pi}{672} + \frac{130583\pi}{2016} v \right) \Delta \Sigma_{\ell} \right] \right. \\
& \left. + \mathcal{O}\left(\frac{1}{c^6}\right) \right\}.
\end{aligned}
\tag{612}
\end{equation}
$$

We recall that $S_\ell \equiv \ell \cdot S$ and $\Sigma_\ell \equiv \ell \cdot \Sigma$, with $S$ and $\Sigma$ denoting the combinations (586), and the individual spins are the specific conserved-norm spins that have been constructed in Sect. 3.6.1. The result (612) superposes to the non-spin contributions given by Eq. (483). Satisfyingly it is in complete agreement in the test-mass limit where $v \to 0$ with the result of black-hole perturbation theory on a Kerr background obtained by Tagoshi et al. (1996).

Finally we can compute the spin effects in the time evolution of the binary’s
orbital frequency Ω. As in Sect. 3.4 we rely on Eq. (447) balancing the total emitted
energy flux $\mathcal{F}$ with the variation of the binary's center-of-mass energy E. The SO part of the energy E have been provided in Eq. (592). Using E and $\mathcal{F}$ expressed as functions of the orbital frequency $\Omega$ (through $x$) and of the spin variables (through $S_\ell$ and $\Sigma_\ell$), we transform the balance equation into

$$ \dot{\Omega}_{\text{SO}} = - \left( \frac{\mathcal{F}}{\text{d}E/\text{d}\Omega} \right)_{\text{SO}} \qquad (613) $$

However, in writing the latter equation it is important to justify that the spin quantities $S_\ell$ and $\Sigma_\ell$ are secularly constant, i.e., do not evolve on a gravitational radiation reaction time scale so we can neglect their variations when taking the time derivative of Eq. (592). Fortunately, this is the case of the conserved-norm spin variables, as proved by Will (2005) up to relative 1PN order, i.e., considering radiation reaction effects up to 3.5PN order. Furthermore this can also be shown from the following structural general argument valid at linear order in spins (Blanchet et al. 2011; Bohé et al. 2013a). In the center-of-mass frame, the only vectors at our disposal, except for the spins, are $n$ and $v$. Recalling that the spin vectors are pseudovectors regarding parity transformation, we see that the only way SO contributions can enter scalars such as the energy E or the flux $\mathcal{F}$ is through the mixed products $(n, v, S_a)$, i.e., through the components $S_\ell^a$. Now, the same reasoning applies to the precession vectors $\Omega_a$ in Eqs. (585): They must be pseudovectors, and, at linear order in spin, they must only depend on $n$ and $v$; so that they must be proportional to $\ell$, as can be explicitly seen for instance in Eq. (591). Now, the time derivative of the components along $\ell$ of the spins are given by $dS_\ell^a/dt = S_a \cdot (d\ell/dt + \ell \times \Omega_a)$. The second term vanishes because $\Omega_a \propto \ell$, and since $d\ell/dt = \mathcal{O}(S)$, we obtain that $S_\ell^a$ is constant at linear order in the spins. We have already met an instance of this important fact in Eq. (607c). This argument is valid at any post-Newtonian order and for general orbits, but is limited to spin-orbit terms; furthermore it does not specify any time scale for the variation, so it applies to short time scales such as the orbital and precessional periods, as well as to the long gravitational radiation reaction time scale (see also Gergely 2000 and references therein for related discussions).

We are then allowed to apply Eq. (613) with conserved-norm spin variables at the SO level. We thus obtain the secular evolution of $\Omega$ and from that we deduce by a further integration (following the Taylor approximant T2) the secular evolution of the carrier phase $\phi \equiv \int \Omega dt$:

$$ \phi_{\text{SO}} = -\frac{x^{-1}}{32 G m^2 v} \left\{ \begin{aligned} & \frac{235}{6} S_{\ell} + \frac{125}{8} \Delta \Sigma_{\ell} \\ & + x \ln x \left[ \left( -\frac{554345}{2016} - \frac{55}{8} v \right) S_{\ell} + \left( -\frac{41745}{448} + \frac{15}{8} v \right) \Delta \Sigma_{\ell} \right] \\ & + x^{3/2} \left[ \frac{940\pi}{3} S_{\ell} + \frac{745\pi}{6} \Delta \Sigma_{\ell} \right] \\ & + x^2 \left[ \left( -\frac{8980424995}{6096384} + \frac{6586595}{6048} - v - \frac{305}{288} v^2 \right) S_{\ell} \right] \end{aligned} \right\} $$
$$ + \left( -\frac{170978035}{387072} + \frac{2876425}{5376} v + \frac{4735}{1152} v^2 \right) \Delta\Sigma_{\ell} + x^{5/2} \left[ \left( \frac{2388425\pi}{3024} - \frac{9925\pi}{36} v \right) S_{\ell} + \left( \frac{3237995\pi}{12096} - \frac{258245\pi}{2016} v \right) \Delta\Sigma_{\ell} \right] + O\left(\frac{1}{c^6}\right) \tag{614} $$

This expression, when added to the expression for the non-spin terms given by Eq. (486), and considering also the SS terms, constitutes the main theoretical input needed for the construction of templates for spinning compact binaries. However, recall that in the case of precessional binaries, for which the spins are not aligned or anti-aligned with the orbital angular momentum, we must subtract to the carrier phase $\phi$ the precessional correction $\alpha$ arising from the precession of the orbital plane. Indeed the physical phase variable $\Phi$ which is defined in Fig. 9, has been proved to be given by $\Phi = \phi - \alpha$ at linear order in spins, cf. Equation (603). The precessional correction $\alpha$ can be computed at linear order in spins from the results of Sect. 3.6.2.

### 3.6.5 Spin-spin effects

The SS effects are known at the leading level corresponding to 2PN order from Barker and O'Connell (1975, 1979) in the equations of motion (see Kidder 1995; Porto 2006; Buonanno et al. 2013 for subsequent derivations), and from Kidder et al. (1993b); Kidder (1995) in the radiation field. Next-to-leading SS contributions are at 3PN order and have been obtained with the Hamiltonian (Steinhoff et al. 2008a, b, c; Hergt et al. 2010; Hartung and Steinhoff 2011a), the EFT (Porto and Rothstein 2006, 2008a, b; Levi 2010b; Levi and Steinhoff 2015a) and the harmonic coordinates (Bohé et al. 2015) techniques. Bohé et al. (2015), Cho et al. (2021) obtained also the next-to-leading SS terms in the gravitational-wave flux. The next-to-next-to-next-to-leading order SS effects in the equations of motion have been obtained by Antonelli et al. (2020a), Kim et al. (2023a), Mandal et al. (2023a).

With SS effects in a compact binary system one must make a distinction between the spin squared terms, involving the coupling between the two same spins $S_1$ or $S_2$, and the interaction terms, involving the coupling between the two different spins $S_1$ and $S_2$. The spin-squared terms $S_1^2$ and $S_2^2$ arise due to the effects on the dynamics of the quadrupole moments of the compact bodies that are induced by their spins, as shown by Poisson (1997). They have been computed through 2PN order in the fluxes and orbital phase by Gergely (1999, 2000), Mikóczi et al. (2005). The interaction terms $S_1 \times S_2$ are simpler, as they can be computed with the pole-dipole formalism reviewed in Sect. 3.6.1. The interaction terms $S_1 \times S_2$ between different spins have been derived to next-to-next-to-leading 4PN order for the equations of motion by Levi (2012), Levi and Steinhoff (2016) (EFT) and by Hartung and Steinhoff (2011c) (Hamiltonian).
The SS interactions (spin-squared terms $S_1^2$ and $S_2^2$) can be computed by representing each compact object as a spinning point particle whose internal structure and finite-size effects are parametrized by an effective quadrupole moment. The dynamics of such particles with higher order multipoles was investigated by Dixon (1964, 1973, 1979), Bailey and Israel (1975) who constructed the effective stress-energy tensor and Lagrangian to encode information about the internal structure of the body. In the case of a spinning particle endowed with a quadrupolar structure the equations generalizing the Mathisson-Papapetrou equations (574) and (572) read (with $c = 1$)

$$ \frac{Dp_{\alpha}}{d\tau} = -\frac{1}{2}u^{\beta}R_{\alpha\beta\mu\nu}S^{\mu\nu} - \frac{1}{3}\nabla_{\lambda}R_{\alpha\beta\mu\nu}J^{\lambda\beta\mu\nu}, \qquad (615a) $$

$$ \frac{DS^{\alpha\beta}}{d\tau} = 2p^{[\alpha}u^{\beta]} + \frac{4}{3}R^{[\alpha}_{\gamma\mu\nu}J^{\beta]}\gamma^{\mu\nu}, \qquad (615b) $$

where $J^{\alpha\beta\mu\nu}$ is the Dixon quadrupolar tensor, which is at this stage only constrained to have the same symmetry properties as $R_{\alpha\beta\mu\nu}$. We account for the correct number of degrees of freedom of the spin by imposing the same covariant SSC (576) as before. Similarly for the quadrupole tensor $J^{\alpha\beta\mu\nu}$, not all of its degrees of freedom are physical, and we restrict these by relating the quadrupole tensor to the four-velocity $u^{\mu}$ and the spin tensor $S^{\mu\nu}$. In order to describe a physical body with a spin and a non-vanishing spin-induced quadrupole moment one has (Porto 2006; Porto and Rothstein 2008a; Steinhoff and Puetzfeld 2010; Steinhoff 2011; Buonanno et al. 2013; Bohé et al. 2015)

$$ J^{\alpha\beta\mu\nu} = \frac{3\kappa}{m} S^{\lambda[\alpha} u^{\beta]} S^{\mu]u}{}^{\nu}, \qquad (616) $$

where the mass parameter is defined by $m^2 = -g^{\mu\nu}p_{\mu}p_{\nu}$ as in Sect. 3.6.1, but is here constant only modulo higher-order $O(S^2)$ terms. The constant $\kappa$ characterizes the quadrupolar “polarisability” of the body; for a neutron star the numerical value of $\kappa$ depends on the equation of state of the nuclear matter.

As in Sect. 3.6.1 we apply the formalism to self-gravitating binary systems (strictly speaking the formalism is applied to spin-squared contributions $S_1^2$ and $S_2^2$; the pole-dipole formalism is sufficient for $S_1 \times S_2$). We present the complete SS contributions, including all $S_1^2$, $S_2^2$ and $S_1 \times S_2$ terms, to the conservative energy in the case of spins aligned with the orbital angular momentum and for circular orbits (Bohé et al. 2015):
$$
\begin{equation}
\begin{aligned}
E_{SS} = & -\frac{\nu c^2 x^3}{2G^2 m^3} \left\{ S_{\ell}^2 (-\kappa_{+} - 2) + S_{\ell} \Sigma_{\ell} (-\Delta\kappa_{+} - 2\Delta + \kappa_{-}) \right. \\
& \left. \qquad + \Sigma_{\ell}^2 \left( \left( \frac{\Delta\kappa_{-}}{2} - \frac{\kappa_{+}}{2} \right) + \nu(\kappa_{+} + 2) \right) \right. \\
& \left. \qquad + x \left[ S_{\ell}^2 \left( \left( -\frac{5\Delta\kappa_{-}}{3} - \frac{25\kappa_{+}}{6} + \frac{50}{9} \right) + \nu \left( \frac{5\kappa_{+}}{6} + \frac{5}{3} \right) \right) \right. \right. \\
& \left. \qquad + S_{\ell} \Sigma_{\ell} \left( \left( -\frac{5\Delta\kappa_{+}}{2} + \frac{25\Delta}{3} + \frac{5\kappa_{-}}{2} \right) + \nu \left( \frac{5\Delta\kappa_{+}}{6} + \frac{5\Delta}{3} + \frac{35\kappa_{-}}{6} \right) \right) \right. \\
& \left. \qquad + \Sigma_{\ell}^2 \left( \left( \frac{5\Delta\kappa_{-}}{4} - \frac{5\kappa_{+}}{4} + 5 \right) + \nu \left( \frac{5\Delta\kappa_{-}}{4} + \frac{5\kappa_{+}}{4} - 10 \right) \right. \right. \\
& \left. \qquad \left. + \nu^2 \left( -\frac{5\kappa_{+}}{6} - \frac{5}{3} \right) \right] \right\} + \mathcal{O}\left(\frac{1}{c^4}\right).
\end{aligned}
\tag{617}
\end{equation}
$$

This expression is to be added to the non-spin contributions (377) and to the SO ones given in Eq. (592). Here $\kappa_+ = \kappa_1 + \kappa_2$ and $\kappa_- = \kappa_1 - \kappa_2$ where $\kappa_a$ are the polar-isability parameters of the two bodies. The expression (617) can be shown to be in agreement, in the test-mass limit for one body, with the energy of a test particle in circular equatorial orbits around a Kerr black hole (Bardeen et al. 1972).

The energy flux for SS interactions—thus complementing the non-spin (483) and SO (612) terms—in the case of aligned spins (no orbital precession) and circular orbits, reads (Bohé et al. 2015; Cho et al. 2021)

$$
\begin{equation}
\begin{aligned}
\mathcal{F}_{SS} = & \frac{32\nu^2 c^5 x^7}{5G^3 m^4} \left\{ S_{\ell}^2 (2\kappa_+ + 4) + S_{\ell} \Sigma_{\ell} (2\Delta\kappa_+ + 4\Delta - 2\kappa_-) \right. \\
& \left. \qquad + \Sigma_{\ell}^2 \left( \left(-\Delta\kappa_- + \kappa_+ + \frac{1}{16}\right) + \nu(-2\kappa_+ - 4) \right) \right. \\
& \left. \qquad + x \left[ S_{\ell}^2 \left( \left(\frac{41\Delta\kappa_-}{16} - \frac{271\kappa_+}{112} - \frac{5239}{504}\right) + \nu\left(-\frac{43\kappa_+}{4} - \frac{43}{2}\right) \right) \right. \right. \\
& \left. \qquad + S_{\ell} \Sigma_{\ell} \left( \left(-\frac{279\Delta\kappa_+}{56} - \frac{817\Delta}{56} + \frac{279\kappa_-}{56}\right) \right. \right. \\
& \left. \qquad \left. + \nu\left(-\frac{43\Delta\kappa_+}{4} - \frac{43\Delta}{2} + \frac{\kappa_-}{2}\right) \right) \right. \\
& \left. \qquad + \Sigma_{\ell}^2 \left( \left(\frac{279\Delta\kappa_-}{112} - \frac{279\kappa_+}{112} - \frac{25}{8}\right) + \nu\left(\frac{45\Delta\kappa_-}{16} + \frac{243\kappa_+}{112} + \frac{344}{21}\right) \right) \right. \\
& \left. \qquad \left. + \nu^2 \left(\frac{43\kappa_+}{4} + \frac{43}{2}\right) \right] \right\} + \mathcal{O}\left(\frac{1}{c^3}\right).
\end{aligned}
\tag{618}
\end{equation}
$$

This result is in agreement, in the limit of a test particle orbiting a Kerr black hole, with the black hole perturbation calculation of Tagoshi et al. (1996). Using (618) together with the expression of the orbital energy (617), we can solve the energy balance equation to obtain the phase evolution of the binary for circular orbits and aligned spins. This yields (Bohé et al. 2015)
$$ \phi_{SS} = -\frac{x^{-1/2}}{32G^2m^4v} \left\{ S_{\ell}^2(-25\kappa_+ - 50) + S_{\ell}\Sigma_{\ell}(-25\Delta\kappa_+ - 50\Delta + 25\kappa_-) \right. \\ + \Sigma_{\ell}^2 \left( \left( \frac{25\Delta\kappa_-}{2} - \frac{25\kappa_+}{2} - \frac{5}{16} \right) + v(25\kappa_+ + 50) \right) \\ + x \left[ S_{\ell}^2 \left( \left( \frac{2215\Delta\kappa_-}{48} + \frac{15635\kappa_+}{84} - \frac{31075}{126} \right) + v(30\kappa_+ + 60) \right) \right. \\ + S_{\ell}\Sigma_{\ell} \left( \left( \frac{47035\Delta\kappa_+}{336} - \frac{9775\Delta}{42} - \frac{47035\kappa_-}{336} \right) \right. \\ + v \left( 30\Delta\kappa_+ + 60\Delta - \frac{2575\kappa_-}{12} \right) \left. \right) \\ + \Sigma_{\ell}^2 \left( \left( -\frac{47035\Delta\kappa_-}{672} + \frac{47035\kappa_+}{672} - \frac{410825}{2688} \right) \right. \\ + v \left( -\frac{2935\Delta\kappa_-}{48} - \frac{4415\kappa_+}{56} + \frac{23535}{112} \right) \\ + v^2(-30\kappa_+ - 60) \left. \right] + \mathcal{O}\left(\frac{1}{c^3}\right) \right\}. \tag{619} $$

Concerning the waveform the spin contributions for spin-aligned circular binaries have been included in the polarization modes at the 3.5PN order (Henry et al. 2022); for non-circular binaries (with aligned or anti-aligned spins) they are known to the third post-Newtonian order (Henry and Khalil 2023). We mention also the calculations of the subdominant cubic spin-spin-spin (SSS) and even quartic SSSS effects (Hergt and Schäfer 2008, 2008; Levi and Steinhoff 2015a; Marsat 2015; Vaidya 2015).

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Abbott BP et al (2016a) Binary black hole mergers in the first Advanced LIGO observing run. Phys Rev X 6:041015. https://doi.org/10.1103/PhysRevX.6.041015. arXiv:1606.04856 [gr-qc]

Abbott BP et al (2016b) Tests of general relativity with GW150914. Phys Rev Lett 116:221101. https://doi.org/10.1103/PhysRevLett.116.221101. arXiv:1602.03841 [gr-qc]

Abbott BP et al (2016c) Observation of gravitational waves from a binary black hole merger. Phys Rev Lett 116:061102. https://doi.org/10.1103/PhysRevLett.116.061102. arXiv:1602.03837 [gr-qc]

Abbott B et al (2017a) Multi-messenger observations of a binary neutron star merger. Astrophys J Lett 848: L12. https://doi.org/10.3847/2041-8213/aa91c9. arXiv:1710.05833 [gr-qc]

Abbott BP et al (2017b) GW170817: observation of gravitational waves from a binary neutron star inspiral. Phys Rev Lett 119:161101. https://doi.org/10.1103/PhysRevLett.119.161101. arXiv:1710.05832 [gr-qc]
Abbott BP et al (2019) Tests of general relativity with GW170817. Phys Rev Lett 123(1):011102. https://doi.org/10.1103/PhysRevLett.123.011102

Abdelsalhin T, Gualtieri L, Pani P (2018) Post-Newtonian spin-tidal couplings for compact binaries. Phys Rev D 98. https://doi.org/10.1103/PhysRevD.98.104046

Abramowicz MA, Kluźniak W (2001) A precise determination of black hole spin in GRO J1655-40. Astron Astrophys 374:L19-L20. https://doi.org/10.1051/0004-6361:20010791. arXiv:astro-ph/0105077

Adamcewicz C, Galaudage S, Lasky PD, Thrane E (2024) Which black hole is spinning? Probing the origin of black hole spin with gravitational waves. Astrophys J Lett 964(1):L6. https://doi.org/10.3847/2041-8213/ad2df2. arXiv:2311.05182 [astro-ph.HE]

Ajith P, Iyer BR, Robinson CAK, Sathyaprakash BS (2005) New class of post-Newtonian approximants to the waveform templates of inspiralling compact binaries: test mass in the Schwarzschild spacetime. Phys Rev D 71:044029 arXiv:gr-qc/0412033

Ajith P, Babak S, Chen Y, Hewitson M, Krishnan B, Sintes AM, Whelan JT, Brügmann B, Diener P, Dorband N, Gonzalez J, Hannam M, Husa S, Pollney D, Rezzolla L, Santamaría L, Sperhake U, Thornburg J (2008) Template bank for gravitational waveforms from coalescing binary black holes: nonspinning binaries. Phys Rev D 77:104017. https://doi.org/10.1103/PhysRevD.77.104017, Erratum: Phys Rev D, 79, 129901(E) (2009). arXiv:0710.2335 [gr-qc]

Albertini A, Nagar A, Pound A, Warburton N, Wardell B, Durkan L, Miller J (2022a) Comparing second-order gravitational self-force and effective one body waveforms from inspiralling, quasicircular and nonspinning black hole binaries. II. The large-mass-ratio case. Phys Rev D 106(8):084062. https://doi.org/10.1103/PhysRevD.106.084062. arXiv:2208.02055 [gr-qc]

Albertini A, Nagar A, Pound A, Warburton N, Wardell B, Durkan L, Miller J (2022b) Comparing second-order gravitational self-force, numerical relativity, and effective one body waveforms from inspiralling, quasicircular, and nonspinning black hole binaries. Phys Rev D 106(8):084061. https://doi.org/10.1103/PhysRevD.106.084061. arXiv:2208.01049 [gr-qc]

Almeida GL, Foffa S, Sturani R (2021) Gravitational multipole renormalization. Phys Rev D 104 (8):084095. https://doi.org/10.1103/PhysRevD.104.084095. arXiv:2107.02634 [gr-qc]

Almeida GL, Foffa S, Sturani R (2023) Gravitational radiation contributions to the two-body scattering angle. Phys Rev D 107(2):024020. https://doi.org/10.1103/PhysRevD.107.024020. arXiv:2209.11594 [gr-qc]

Alvi K (2001) Energy and angular momentum flow into a black hole in a binary. Phys Rev D 64:104020. https://doi.org/10.1103/PhysRevD.64.104020. arXiv:0107080 [gr-qc]

Anderson JL, DeCanio TC (1975) Equations of hydrodynamics in general relativity in the slow motion approximation. Gen Relativ Gravit 6:197-237. https://doi.org/10.1007/BF00769986

Anderson JL, Kates RE, Kegeles LS, Madonna RG (1982) Divergent integrals of post-Newtonian gravity: nonanalytic terms in the near-zone expansion of a gravitationally radiating system found by matching. Phys Rev D 25:2038-2048. https://doi.org/10.1103/PhysRevD.25.2038

Antonelli A, Kavanagh C, Khalil M, Steinhoff J, Vines J (2020a) Gravitational spin-orbit and aligned spin1-spin2 couplings through third-subleading post-Newtonian orders. Phys Rev D 102(12). https://doi.org/10.1103/physrevd.102.124024

Antonelli A, Kavanagh C, Khalil M, Steinhoff J, Vines J (2020b) Gravitational spin-orbit coupling through third-subleading post-Newtonian order: from first-order self-force to arbitrary mass ratios. Phys Rev Lett 125(1). https://doi.org/10.1103/physrevlett.125.011103

Apostolatos TA, Cutler C, Sussman GJ, Thorne KS (1994) Spin induced orbital precession and its modulation of the gravitational wave forms from merging binaries. Phys Rev D 49:6274-6297. https://doi.org/10.1103/PhysRevD.49.6274

Arun KG, Blanchet L, Iyer BR, Qusailah MSS (2004) The 2.5PN gravitational wave polarizations from inspiralling compact binaries in circular orbits. Class Quantum Grav 21:3771-3801. https://doi.org/10.1088/0264-9381/21/15/010, Erratum: Class Quantum Grav, 22:3115 (2005). arXiv:gr-qc/0404185

Arun KG, Iyer BR, Sathyaprakash BS, Sundararajan PA (2005) Parameter estimation of inspiralling compact binaries using 35 post-Newtonian gravitational wave phasing: the nonspinning case. Phys Rev D 71:084008. https://doi.org/10.1103/PhysRevD.71.084008. arXiv:gr-qc/0411146

Arun KG, Iyer BR, Qusailah MSS, Sathyaprakash BS (2006a) Probing the non-linear structure of general relativity with black hole binaries. Phys Rev D 74:024006. https://doi.org/10.1103/PhysRevD.74.024006. arXiv:gr-qc/0604067
Arun KG, Iyer BR, Qusailah MSS, Sathyaprakash BS (2006b) Testing post-Newtonian theory with gravitational wave observations. Class Quantum Grav 23:L37–L43. https://doi.org/10.1088/0264-9381/23/9/L01. arXiv:gr-qc/0604018

Arun KG, Iyer BR, Sathyaprakash BS, Sinha S (2007a) Higher harmonics increase LISA’s mass reach for supermassive black holes. Phys Rev D 75:124002. https://doi.org/10.1103/PhysRevD.75.124002. arXiv:0704.1086

Arun KG, Iyer BR, Sathyaprakash BS, Sinha S, Van Den Broeck C (2007b) Higher signal harmonics, LISA’s angular resolution, and dark energy. Phys Rev D 76:104016. https://doi.org/10.1103/PhysRevD.76.104016. arXiv:0707.3920

Arun KG, Blanchet L, Iyer BR, Qusailah MS (2008a) Inspiralling compact binaries in quasi-elliptical orbits: the complete 3PN energy flux. Phys Rev D 77:064035. https://doi.org/10.1103/PhysRevD.77.064035. arXiv:0711.0302

Arun KG, Blanchet L, Iyer BR, Qusailah MS (2008b) Tail effects in the 3PN gravitational wave energy flux of compact binaries in quasi-elliptical orbits. Phys Rev D 77:064034. https://doi.org/10.1103/PhysRevD.77.064034. arXiv:0711.0250

Arun KG, Blanchet L, Iyer BR, Sinha S (2009a) Third post-Newtonian angular momentum flux and the secular evolution of orbital elements for inspiralling compact binaries in quasi-elliptical orbits. Phys Rev D 80:124018. https://doi.org/10.1103/PhysRevD.80.124018. arXiv:0908.3854

Arun KG, Buonanno A, Faye G, Ochsner E (2009b) Higher-order spin effects in the amplitude and phase of gravitational waveforms emitted by inspiralling compact binaries: ready-to-use gravitational waveforms. Phys Rev D 79:104023. https://doi.org/10.1103/PhysRevD.79.104023. arXiv:0810.5336

Ashtekar A, Campiglia M, Laddha A (2018) Null infinity, the BMS group and infrared issues. Gen Relativ Gravit 50:1–23. https://doi.org/10.1007/s10714-018-2464-3

Bailey I, Israel W (1975) Lagrangian dynamics of spinning particles and polarized media in general relativity. Commun Math Phys 42:65. https://doi.org/10.1007/BF01609434

Baker JG, Centrella J, Choi DI, Koppitz M, van Meter J, Miller MC (2006a) Getting a kick out of numerical relativity. Astrophys J 653:L93–L96. https://doi.org/10.1086/510448. arXiv:astro-ph/0603204

Baker JG, Centrella J, Choi DI, Koppitz M, van Meter JR (2006b) Gravitational-wave extraction from an inspiralling configuration of merging black holes. Phys Rev Lett 96:111102. https://doi.org/10.1103/PhysRevLett.96.111102. arXiv:gr-qc/0511103

Banihashemi B, Vines J (2020) Gravitomagnetic tidal effects in gravitational waves from neutron star binaries. Phys Rev D 101(6):064003. https://doi.org/10.1103/PhysRevD.101.064003. arXiv:1805.07266 [gr-qc]

Barack L (2009) Gravitational self-force in extreme mass-ratio inspirals. Class Quantum Grav 26:213001. https://doi.org/10.1088/0264-9381/26/21/213001. arXiv:0908.1664 [gr-qc]

Barack L (2011) Computational methods for the self-force in black hole spacetimes. In: Blanchet L, Spallicci A, Whiting B (eds) Mass and motion in general relativity. Fundamental theories of physics, vol 162. Springer, New York, pp 327–366. https://doi.org/10.1007/978-90-481-3015-3_12

Barack L, Pound A (2018) Self-force and radiation reaction in general relativity. Rep Prog Phys 82 (1):016904. https://doi.org/10.1088/1361-6633/aae552. arXiv:1805.10385 [gr-qc]

Barack L, Sago N (2009) Gravitational self-force correction to the innermost stable circular orbit of a Schwarzschild black hole. Phys Rev Lett 102:191101. https://doi.org/10.1103/PhysRevLett.102.191101. arXiv:0902.0573

Barack L, Sago N (2011) Beyond the geodesic approximation: conservative effects of the gravitational self-force in eccentric orbits around a Schwarzschild black hole. Phys Rev D 83:084023. https://doi.org/10.1103/PhysRevD.83.084023. arXiv:1101.3331 [gr-qc]

Barack L, Damour T, Sago N (2010) Precession effect of the gravitational self-force in a Schwarzschild spacetime and the effective one-body formalism. Phys Rev D 82(8):084036. https://doi.org/10.1103/PhysRevD.82.084036

Barausse E, Racine E, Buonanno A (2009) Hamiltonian of a spinning test particle in curved spacetime. Phys Rev D 80:104025. https://doi.org/10.1103/PhysRevD.80.104025. arXiv:0907.4745 [gr-qc]

Bardeen JM, Press WH, Teukolsky SA (1972) Rotating black holes: locally nonrotating frames, energy extraction, and scalar synchrotron radiation. Astrophys J 178:347. https://doi.org/10.1086/151796

Bardeen JM, Carter B, Hawking SW (1973) The four laws of black hole mechanics. Commun Math Phys 31:161–170. https://doi.org/10.1007/BF01645742

Barker BM, O’Connell RF (1975) Gravitational two-body problem with arbitrary masses, spins, and quadrupole moments. Phys Rev D 12:329–335. https://doi.org/10.1103/PhysRevD.12.329
Barker BM, O'Connell RF (1979) The gravitational interaction: spin, rotation, and quantum effects—a review. Gen Relativ Gravit 11:149–175. https://doi.org/10.1007/BF00756587

Baumgarte TW (2000) Innermost stable circular orbit of binary black holes. Phys Rev D 62:024018. https://doi.org/10.1103/PhysRevD.62.024018

Bekenstein JD (1973) Gravitational radiation recoil and runaway black holes. Astrophys J 183:657–664. https://doi.org/10.1086/152255

Bel L, Damour T, Deruelle N, Ibáñez J, Martin J (1981) Poincaré-invariant gravitational field and equations of motion of two pointlike objects: the postlinear approximation of general relativity. Gen Relativ Gravit 13:963–1004. https://doi.org/10.1007/BF00756073

Benacquista MJ, Downing JMB (2013) Relativistic binaries in globular clusters. Living Rev Relativ 16:4. https://doi.org/10.12942/lrr-2013-4. arXiv:1110.4423

Bern Z, Cheung C, Roiban R, Shen CH, Solon MP, Zeng M (2019a) Black hole binary dynamics from the double copy and effective theory. J High Energy Phys 10:206. https://doi.org/10.1007/JHEP10(2019) 206. arXiv:1908.01493 [hep-th]

Bern Z, Cheung C, Roiban R, Shen CH, Solon MP, Zeng M (2019b) Scattering amplitudes and the conservative Hamiltonian for binary systems at third post-Minkowskian order. Phys Rev Lett 122 (20):201603. https://doi.org/10.1103/PhysRevLett.122.201603. arXiv:1901.04424 [hep-th]

Bern Z, Parra-Martinez J, Roiban R, Ruf MS, Shen CH, Solon MP, Zeng M (2021) Scattering amplitudes and conservative binary dynamics at $O(G^4)$. Phys Rev Lett 126(17):171601. https://doi.org/10.1103/PhysRevLett.126.171601. arXiv:2112.10750

Bernard L (2018) Dynamics of compact binary systems in scalar-tensor theories: equations of motion to the third post-Newtonian order. Phys Rev D 98:044004. https://doi.org/10.1103/PhysRevD.98.044004. arXiv:1802.10201 [gr-qc]

Bernard L (2019) Dynamics of compact binary systems in scalar-tensor theories. II. Center-of-mass and conserved quantities to 3PN order. Phys Rev D 99(4):044047. https://doi.org/10.1103/PhysRevD.99. 044047. arXiv:1812.04169 [gr-qc]

Bernard L, Blanchet L, Bohé A, Faye G, Marsat S (2016) Fokker action of non-spinning compact binaries at the fourth post-Newtonian approximation. Phys Rev D 93:084037. https://doi.org/10.1103/ PhysRevD.93.084037. arXiv:1512.02876 [gr-qc]

Bernard L, Blanchet L, Bohé A, Faye G, Marsat S (2017a) Dimensional regularization of the IR divergences in the Fokker action of point-particle binaries at the fourth post-Newtonian order. Phys Rev D 96:104043. https://doi.org/10.1103/PhysRevD.96.104043. arXiv:1706.08480 [gr-qc]

Bernard L, Blanchet L, Bohé A, Faye G, Marsat S (2017b) Energy and periastron advance of compact binaries on circular orbits at the fourth post-Newtonian order. Phys Rev D 95:044026. https://doi.org/ 10.1103/PhysRevD.95.044026. arXiv:1610.07934 [gr-qc]

Bernard L, Blanchet L, Faye G, Marchand T (2018) Center-of-mass equations of motion and conserved integrals of compact binary systems at the fourth post-Newtonian order. Phys Rev D 97:044037. https://doi.org/10.1103/PhysRevD.97.044037. arXiv:1711.00283 [gr-qc]

Bernard L, Blanchet L, Trestini D (2022) Gravitational waves in scalar-tensor theory to one-and-a-half post-Newtonian order. J Cosmol Astropart Phys 08:008. https://doi.org/10.1088/1475-7516/2022/08/ 008. arXiv:2201.10924 [gr-qc]

Bernuzzi S, Thierfelder M, Brügmann B (2012) Accuracy of numerical relativity waveforms from binary neutron star mergers and their comparison with post-Newtonian waveforms. Phys Rev D 85:104030. https://doi.org/10.1103/PhysRevD.85.104030. arXiv:1109.3611 [gr-qc]

Bertotti B, Plebański JF (1960) Theory of gravitational perturbations in the fast motion approximation. Ann Phys (NY) 11:169–200. https://doi.org/10.1016/0003-4916(60)90132-9

Bini D, Damour T (2013) Analytical determination of the two-body gravitational interaction potential at the fourth post-Newtonian approximation. Phys Rev D 87:121501. https://doi.org/10.1103/ PhysRevD.87.121501. arXiv:1305.4884 [gr-qc]

Bini D, Damour T (2014a) Analytic determination of the eight-and-a-half post-Newtonian self-force contributions to the two-body gravitational interaction potential. Phys Rev D 89:104047. https://doi. org/10.1103/PhysRevD.89.104047. arXiv:1403.2366 [gr-qc]

Bini D, Damour T (2014b) High-order post-Newtonian contributions to the two-body gravitational interaction potential from analytical gravitational self-force calculations. Phys Rev D 89:064063. https://doi.org/10.1103/PhysRevD.89.064063. arXiv:1312.2503 [gr-qc]
Bini D, Damour T, Faye G (2012) Effective action approach to higher-order relativistic tidal interactions in binary systems and their effective one body description. Phys Rev D 85:124034. https://doi.org/10.1103/PhysRevD.85.124034. arXiv:1202.3565 [gr-qc]

Bini D, Damour T, Geralico A (2016) Confirming and improving post-Newtonian and effective-one-body results from self-force computations along eccentric orbits around a Schwarzschild black hole. Phys Rev D 93:064023. https://doi.org/10.1103/PhysRevD.93.064023. arXiv:1511.04533 [gr-qc]

Bini D, Damour T, Geralico A (2019) Novel approach to binary dynamics: application to the fifth post-Newtonian level. Phys Rev Lett 123(23):231104. https://doi.org/10.1103/PhysRevLett.123.231104. arXiv:1909.02375 [gr-qc]

Bini D, Damour T, Geralico A (2020a) Binary dynamics at the fifth and fifth-and-a-half post-Newtonian orders. Phys Rev D 102(2):024062. https://doi.org/10.1103/PhysRevD.102.024062. arXiv:2003.11891 [gr-qc]

Bini D, Damour T, Geralico A (2020b) Sixth post-Newtonian local-in-time dynamics of binary systems. Phys Rev D 102(2):024061. https://doi.org/10.1103/PhysRevD.102.024061. arXiv:2004.05407 [gr-qc]

Bini D, Damour T, Geralico A (2020c) Sixth post-Newtonian nonlocal-in-time dynamics of binary systems. Phys Rev D 102(8):084047. https://doi.org/10.1103/PhysRevD.102.084047. arXiv:2007.11239 [gr-qc]

Bini D, Damour T, Geralico A, Laporta S, Mastrolia P (2021) Gravitational scattering at the seventh order in $G$: nonlocal contribution at the sixth post-Newtonian accuracy. Phys Rev D 103(4):044038. https://doi.org/10.1103/PhysRevD.103.044038. arXiv:2012.12918 [gr-qc]

Bini D, Damour T, Geralico A (2023) Comparing one-loop gravitational bremsstrahlung amplitudes to the multipolar-post-Minkowskian waveform. Phys Rev D 108(12):124052. https://doi.org/10.1103/PhysRevD.108.124052. arXiv:2309.14925 [gr-qc]

Bini D, Damour T, De Angelis S, Geralico A, Herderschee A, Roiban R, Teng F (2024) Gravitational waveform: a tale of two formalisms. arXiv e-prints arXiv:2402.06604 [hep-th]

Binnington T, Poisson E (2009) Relativistic theory of tidal Love numbers. Phys Rev D 80:084018. https://doi.org/10.1103/PhysRevD.80.084018. arXiv:0906.1366 [gr-qc]

Blaes O, Lee MH, Socrates A (2002) The Kozai mechanism and the evolution of binary supermassive black holes. Astrophys J 578:775-786. https://doi.org/10.1086/342655. arXiv:astro-ph/0203370

Blanchet L (1987) Radiative gravitational fields in general-relativity. II. Asymptotic-behaviour at future null infinity. Proc R Soc London Ser A 409:383-399. https://doi.org/10.1098/rspa.1987.0022

Blanchet L (1990) Contribution à l'étude du rayonnement gravitationnel émis par un système isolé. PhD thesis, Université Pierre et Marie Curie, Paris VI. http://www2.iap.fr/users/blanchet/TheseHabilitation1990.pdf, see Chapter VI, pp 205-214

Blanchet L (1993) Time asymmetric structure of gravitational radiation. Phys Rev D 47:4392-4420. https://doi.org/10.1103/PhysRevD.47.4392

Blanchet L (1995) Second-post-Newtonian generation of gravitational radiation. Phys Rev D 51:2559-2583. https://doi.org/10.1103/PhysRevD.51.2559. arXiv:gr-qc/9501030

Blanchet L (1996) Energy losses by gravitational radiation in inspiralling compact binaries to five halves post-Newtonian order. Phys Rev D 54:1417-1438. https://doi.org/10.1103/PhysRevD.54.1417, Erratum: Phys Rev D 71:129904(E) (2005). arXiv:gr-qc/9603048

Blanchet L (1997a) Gravitational radiation from relativistic sources. In: Marck JA, Lasota JP (eds) Relativistic gravitation and gravitational radiation. Cambridge University Press, Cambridge, Cambridge Contemporary Astrophysics, pp 33-66. arXiv:gr-qc/9609049

Blanchet L (1997b) Gravitational radiation reaction and balance equations to post-Newtonian order. Phys Rev D 55:714-732. https://doi.org/10.1103/PhysRevD.55.714. arXiv:gr-qc/9609049

Blanchet L (1998a) Gravitational-wave tails of tails. Class Quantum Grav 15:113-141. https://doi.org/10.1088/0264-9381/15/1/009, Erratum: Class Quantum Grav 22:3381 (2005). arXiv:gr-qc/9710038

Blanchet L (1998b) On the multipole expansion of the gravitational field. Class Quantum Grav 15:1971-1999. https://doi.org/10.1088/0264-9381/15/7/013. arXiv:gr-qc/9801101

Blanchet L (1998c) Quadrupole-quadrupole gravitational waves. Class Quantum Grav 15:89-111. https://doi.org/10.1088/0264-9381/15/1/008. arXiv:gr-qc/9710037

Blanchet L (2002) Innermost circular orbit of binary black holes at the third post-Newtonian approximation. Phys Rev D 65:124009. https://doi.org/10.1103/PhysRevD.65.124009. arXiv:gr-qc/0112056
Blanchet L (2011a) On the accuracy of the post-Newtonian approximation. In: Ciufolini I, Dominici D, Lusanna L (eds) 2001: A relativistic spacetime odyssey. World Scientific, pp 411–430. https://doi.org/10.1142/9789812791368_0022. arXiv:gr-qc/0207037

Blanchet L (2011b) Post-Newtonian theory and the two-body problem. In: Blanchet L, Spallicci A, Whiting B (eds) Mass and motion in general relativity. Fundamental theories of physics. Springer, Dordrecht, pp 125–166. https://doi.org/10.1007/978-90-481-3015-3_5. arXiv:0907.3596 [gr-qc]

Blanchet L (2014) Gravitational radiation from post-Newtonian sources and inspiralling compact binaries. Living Rev Relativ 17:2. https://doi.org/10.12942/lrr-2014-2. arXiv:1310.1528 [gr-qc]

Blanchet L, Damour T (1984) Multipolar radiation reaction in general relativity. Phys Lett A 104:82–86. https://doi.org/10.1016/0375-9601(84)90967-8

Blanchet L, Damour T (1986) Radiative gravitational fields in general relativity I. General structure of the field outside the source. Philos Trans R Soc London Ser A 320:379–430. https://doi.org/10.1098/rsta.1986.0125

Blanchet L, Damour T (1988) Tail-transported temporal correlations in the dynamics of a gravitating system. Phys Rev D 37:1410–1435. https://doi.org/10.1103/PhysRevD.37.1410

Blanchet L, Damour T (1989) Post-Newtonian generation of gravitational waves. Ann Inst Henri Poincare A 50:377–408

Blanchet L, Damour T (1992) Hereditary effects in gravitational radiation. Phys Rev D 46:4304–4319. https://doi.org/10.1103/PhysRevD.46.4304

Blanchet L, Faye G (2000a) Equations of motion of point-particle binaries at the third post-Newtonian order. Phys Lett A 271:58–64. https://doi.org/10.1016/S0375-9601(00)00360-1. arXiv:gr-qc/0004009

Blanchet L, Faye G (2000b) Hadamard regularization. J Math Phys 41:7675–7714. https://doi.org/10.1063/1.1308506. arXiv:gr-qc/0004008

Blanchet L, Faye G (2001a) General relativistic dynamics of compact binaries at the third post-Newtonian order. Phys Rev D 63:062005. https://doi.org/10.1103/PhysRevD.63.062005. arXiv:gr-qc/0007051

Blanchet L, Faye G (2001b) Lorentzian regularization and the problem of point-like particles in general relativity. J Math Phys 42:4391–4418. https://doi.org/10.1063/1.1384864. arXiv:gr-qc/0006100 [gr-qc]

Blanchet L, Faye G (2019) Flux-balance equations for linear momentum and center-of-mass position of self-gravitating post-Newtonian systems. Class Quantum Grav 36(8):085003. https://doi.org/10.1088/1361-6382/ab0d4f

Blanchet L, Fokas A (2018) Equations of motion of self-gravitating N-body systems in the first post-Minkowskian approximation. Phys Rev D 98:084005. https://doi.org/10.1103/PhysRevD.98.084005. arXiv:0812.4413 [gr-qc]

Blanchet L, Iyer BR (2003) Third post-Newtonian dynamics of compact binaries: Equations of motion in the center-of-mass frame. Class Quantum Grav 20:755. https://doi.org/10.1088/0264-9381/20/4/309. arXiv:gr-qc/0209089

Blanchet L, Iyer BR (2005) Hadamard regularization of the third post-Newtonian gravitational wave generation of two point masses. Phys Rev D 71:024004. https://doi.org/10.1103/PhysRevD.71.024004. arXiv:gr-qc/0409094

Blanchet L, Le Tiec A (2017) First law of compact binary mechanics with gravitational-wave tails. Class Quantum Grav 34:164001. https://doi.org/10.1088/1361-6382/aa79d7. arXiv:gr-qc/1702.06839

Blanchet L, Sathyaprakash BS (1994) Signal analysis of gravitational wave tails. Class Quantum Grav 11:2807–2831. https://doi.org/10.1088/0264-9381/11/11/020

Blanchet L, Sathyaprakash BS (1995) Detecting a tail effect in gravitational-wave experiments. Phys Rev Lett 74:1067–1070. https://doi.org/10.1103/PhysRevLett.74.1067

Blanchet L, Schäfer G (1989) Higher order gravitational radiation losses in binary systems. Mon Not R Astron Soc 239:845–867. https://doi.org/10.1093/mnras/239.3.845

Blanchet L, Schäfer G (1993) Gravitational wave tails and binary star systems. Class Quantum Grav 10:2699–2721. https://doi.org/10.1088/0264-9381/10/12/026

Blanchet L, Damour T, Iyer BR (1995a) Gravitational waves from inspiralling compact binaries: energy loss and wave form to second post-Newtonian order. Phys Rev D 51:5360–5386. https://doi.org/10.1103/PhysRevD.51.5360. arXiv:gr-qc/9501029

Blanchet L, Damour T, Iyer BR, Will CM, Wiseman AG (1995b) Gravitational-radiation damping of compact binary systems to second post-Newtonian order. Phys Rev Lett 74:3515–3518. https://doi.org/10.1103/PhysRevLett.74.3515. arXiv:gr-qc/9501027
Blanchet L, Iyer BR, Will CM, Wiseman AG (1996) Gravitational wave forms from inspiralling compact binaries to second-post-Newtonian order. Class Quantum Grav 13:575-584. https://doi.org/10.1088/0264-9381/13/4/002. arXiv:gr-qc/9602024

Blanchet L, Faye G, Ponsot B (1998) Gravitational field and equations of motion of compact binaries to 5/2 post-Newtonian order. Phys Rev D 58:124002. https://doi.org/10.1103/PhysRevD.58.124002. arXiv:gr-qc/9804079

Blanchet L, Faye G, Iyer BR, Joguet B (2002a) Gravitational-wave inspiral of compact binary systems to 7/2 post-Newtonian order. Phys Rev D 65:061501(R). https://doi.org/10.1103/PhysRevD.65.061501, Erratum: Phys Rev D, 71:129902(E) (2005). arXiv:gr-qc/0105099

Blanchet L, Iyer BR, Joguet B (2002b) Gravitational waves from inspiralling compact binaries: energy flux to third post-Newtonian order. Phys Rev D 65:064005. Erratum: Phys Rev D 71:129903(E) (2005). arXiv:gr-qc/0105098

Blanchet L, Damour T, Esposito-Farèse G (2004a) Dimensional regularization of the third post-Newtonian dynamics of point particles in harmonic coordinates. Phys Rev D 69:124007. https://doi.org/10.1103/PhysRevD.69.124007. arXiv:gr-qc/0311052

Blanchet L, Damour T, Esposito-Farèse G, Iyer BR (2004b) Gravitational radiation from inspiralling compact binaries completed at the third post-Newtonian order. Phys Rev Lett 93:091101. https://doi.org/10.1103/PhysRevLett.93.091101. arXiv:gr-qc/0406012

Blanchet L, Damour T, Esposito-Farèse G, Iyer BR (2005a) Dimensional regularization of the third post-Newtonian gravitational wave generation of two point masses. Phys Rev D 71:124004. https://doi.org/10.1103/PhysRevD.71.124004. arXiv:gr-qc/0503044

Blanchet L, Damour T, Iyer BR (2005b) Surface-integral expressions for the multipole moments of post-Newtonian sources and the boosted Schwarzschild solution. Class Quantum Grav 22:155. https://doi.org/10.1088/0264-9381/22/1/011. arXiv:gr-qc/0410021

Blanchet L, Faye G, Nissanke S (2005c) Structure of the post-Newtonian expansion in general relativity. Phys Rev D 72:044024. https://doi.org/10.1103/PhysRevD.72.044024

Blanchet L, Qusailah MS, Will CM (2005d) Gravitational recoil of inspiraling black-hole binaries to second post-Newtonian order. Astrophys J 635:508. https://doi.org/10.1086/497332. arXiv:astro-ph/0507692

Blanchet L, Buonanno A, Faye G (2006) Higher-order spin effects in the dynamics of compact binaries II. Radiation field. Phys Rev D 74:104034. https://doi.org/10.1103/PhysRevD.74.104034, Erratum: Phys Rev D 75:049903 (2007). arXiv:gr-qc/0605140

Blanchet L, Faye G, Iyer BR, Sinha S (2008) The third post-Newtonian gravitational wave polarisations and associated spherical harmonic modes for inspiralling compact binaries in quasi-circular orbits. Class Quantum Grav 25:165003. https://doi.org/10.1088/0264-9381/25/16/165003. arXiv:0802.1249

Blanchet L, Detweiler S, Le Tiec A, Whiting BF (2010a) Higher-order Post-Newtonian fit of the gravitational self-force for circular orbits in the Schwarzschild geometry. Phys Rev D 81:084033. https://doi.org/10.1103/PhysRevD.81.084033. arXiv:1002.0726 [gr-qc]

Blanchet L, Detweiler S, Le Tiec A, Whiting BF (2010b) Post-Newtonian and numerical calculations of the gravitational self-force for circular orbits in the Schwarzschild geometry. Phys Rev D 81:064004. https://doi.org/10.1103/PhysRevD.81.064004. arXiv:0910.0207 [gr-qc]

Blanchet L, Buonanno A, Faye G (2011) Tail-induced spin-orbit effect in the gravitational radiation of compact binaries. Phys Rev D 84:064041. https://doi.org/10.1103/PhysRevD.84.064041. arXiv:1104.5659 [gr-qc]

Blanchet L, Buonanno A, Le Tiec A (2013) First law of mechanics for black hole binaries with spins. Phys Rev D 87:024030. https://doi.org/10.1103/PhysRevD.87.024030. arXiv:1211.1060 [gr-qc]

Blanchet L, Faye G, Whiting BF (2014a) Half-integral conservative post-Newtonian approximations in the redshift factor of black hole binaries. Phys Rev D 89:064026. https://doi.org/10.1103/PhysRevD.89.064026. arXiv:1312.2975 [gr-qc]

Blanchet L, Faye G, Whiting BF (2014b) High-order half-integral conservative post-Newtonian coefficients in the redshift factor of black hole binaries. Phys Rev D 90:044017. https://doi.org/10.1103/PhysRevD.90.044017. arXiv:1405.5151 [gr-qc]

Blanchet L, Foffa S, Larrouturou F, Sturani R (2020) Logarithmic tail contributions to the energy function of circular compact binaries. Phys Rev D 101(8):084045. https://doi.org/10.1103/PhysRevD.101.084045. arXiv:1912.12359 [gr-qc]

Blanchet L, Compère G, Faye G, Oliveri R, Seraj A (2021) Multipole expansion of gravitational waves: from harmonic to Bondi coordinates. J High Energy Phys 2:29. https://doi.org/10.1007/JHEP02(2021)029. arXiv:2011.10000 [gr-qc]