$$ \delta \dot{P}_{\Psi} = 0, \qquad (264b) $$

$$ \delta \dot{R} = \sigma_0 \delta P_R, \qquad (264c) $$

$$ \delta \Omega = \rho_0 \delta R + \tau_0 \delta P_{\Psi}, \qquad (264d) $$

where the coefficients, which depend on the unperturbed orbit, are given by

$$ \pi_0 = \frac{\partial^2 \mathcal{H}}{\partial R^2} [R_0, 0, P_{\Psi}^0], \qquad (265a) $$

$$ \rho_0 = \frac{\partial^2 \mathcal{H}}{\partial R \partial P_{\Psi}} [R_0, 0, P_{\Psi}^0], \qquad (265b) $$

$$ \sigma_0 = \frac{\partial^2 \mathcal{H}}{\partial P_R^2} [R_0, 0, P_{\Psi}^0], \qquad (265c) $$

$$ \tau_0 = \frac{\partial^2 \mathcal{H}}{\partial P_{\Psi}^2} [R_0, 0, P_{\Psi}^0]. \qquad (265d) $$

By looking to solutions proportional to some $e^{i\sigma t}$ one obtains some real frequencies, and therefore one finds stable circular orbits, if and only if

$$ \hat{C}_0 \equiv \pi_0 \sigma_0 > 0. \qquad (266) $$

Using explicitly the Hamiltonian (223) we readily obtain

$$ \hat{C}_0 = \frac{Gm}{R_0^3} \left\{ 1 + \frac{Gm}{R_0 c^2} (-9 + \nu) + \frac{G^2 m^2}{R_0^2 c^4} \left( \frac{117}{4} + \frac{43}{8}\nu + \nu^2 \right) \right. \\ \left. + \frac{G^3 m^3}{R_0^3 c^6} \left( -61 + \left[ \frac{4777}{48} - \frac{325}{64}\pi^2 \right] \nu - \frac{31}{8}\nu^2 + \nu^3 \right) + \mathcal{O}\left(\frac{1}{c^8}\right) \right\}. \qquad (267) $$

This result does not look the same as our previous result (257), but this is simply due to the fact that it depends on the ADM radial separation $R_0$ instead of the harmonic one $r_0$. Fortunately we have derived in Section 7.2 the material needed to connect $R_0$ to $r_0$ with the 3PN accuracy. Indeed, with Eqs. (210) we have the relation valid for general orbits in an arbitrary frame between the separation vectors in both coordinate systems. Specializing that relation to circular orbits we find

$$ R_0 = r_0 \left\{ 1 + \frac{G^2 m^2}{r_0^2 c^4} \left( -\frac{1}{4} - \frac{29}{8}\nu \right) + \frac{G^3 m^3}{r_0^3 c^6} \left( \left[ \frac{3163}{1680} + \frac{21}{32}\pi^2 - \frac{22}{3}\ln\left(\frac{r_0}{r'_0}\right) \right] \nu + \frac{3}{8}\nu^2 \right) + \mathcal{O}\left(\frac{1}{c^8}\right) \right\}. \qquad (268) $$

Note that the difference between $R_0$ and $r_0$ starts only at 2PN order. That relation easily permits to perfectly reconcile both expressions (257) and (267).

Finally let us give to $\hat{C}_0$ an invariant meaning by expressing it with the help of the orbital frequency $\Omega_0$ of the circular orbit, or, more conveniently, of the frequency-related parameter $x_0 \equiv (Gm\Omega_0/c^3)^{2/3} - \text{cf. Eq. (230)}$. This allows us to write the criterion for stability as $C_0 > 0$, where $C_0 = \frac{G^2 m^2}{x_0^3} \hat{C}_0$ admits the gauge-invariant form

$$ C_0 = 1 - 6x_0 + 14\nu x_0^2 + \left( \left[ \frac{397}{2} - \frac{123}{16}\pi^2 \right] \nu - 14\nu^2 \right) x_0^3 + \mathcal{O}(x_0^4). \qquad (269) $$

This form is more interesting than the coordinate-dependent expressions (257) or (267), not only because of its invariant form, but also because as we see the 1PN term yields exactly the Schwarzschild
result that the innermost stable circular orbit or ISCO of a test particle (i.e., in the limit $\nu \to 0$) is located at $x_{\text{ISCO}} = \frac{1}{6}$. Thus we find that, at the 1PN order, but for any mass ratio $\nu$,

$$x_{\text{ISCO}}^{\text{1PN}} = \frac{1}{6}. \qquad (270)$$

One could have expected that some deviations of the order of $\nu$ already occur at the 1PN order, but it turns out that only from the 2PN order does one find the occurrence of some non-Schwarzschildian corrections proportional to $\nu$. At the 2PN order we obtain

$$x_{\text{ISCO}}^{\text{2PN}} = \frac{3}{14\nu} \left( 1 - \sqrt{1 - \frac{14\nu}{9}} \right). \qquad (271)$$

For equal masses this gives $x_{\text{ISCO}}^{\text{2PN}} \simeq 0.187$. Notice also that the effect of the finite mass corrections is to increase the frequency of the ISCO with respect to the Schwarzschild result, i.e., to make it more inward:<fnref n="58" />

$$x_{\text{ISCO}}^{\text{2PN}} = \frac{1}{6} \left[ 1 + \frac{7}{18}\nu + \mathcal{O}(\nu^2) \right]. \qquad (272)$$

Finally, at the 3PN order, for equal masses $\nu = \frac{1}{4}$, we find that according to our criterion all the circular orbits are stable. More generally, we find that at the 3PN order all orbits are stable when the mass ratio $\nu$ is larger than some critical value $\nu_c \simeq 0.183$.

The stability criterion (269) has been compared in great details to various other stability criteria by Favata [191] and shown to perform very well, and has also been generalized to spinning black hole binaries in Ref. [190]. Note that this criterion is based on the physical requirement that a stable perturbation should have a real frequency. It gives an innermost stable circular orbit, when it exists, which differs from the innermost circular orbit or ICO defined in Section 8.1; see Ref. [378] for a discussion on the difference between an ISCO and the ICO in the PN context. Note also that the criterion (269) is based on systematic post-Newtonian expansions, without resorting for instance to Padé approximants. Nevertheless, it performs better than other criteria based on various resummation techniques, as discussed in Ref. [191].

## 8.3 The first law of binary point-particle mechanics

In this section we shall review a very interesting relation for binary systems of point particles modelling black hole binaries and moving on circular orbits, known as the “first law of point-particle mechanics”. This law was obtained using post-Newtonian methods in Ref. [289], but is actually a particular case of a more general law, valid for systems of black holes and extended fluid balls, derived by Friedman, Uryū & Shibata [208].

Before tackling the problem it is necessary to make more precise the notion of circular orbits. These are obtained from the conservative part of the dynamics, neglecting the dissipative radiation-reaction force responsible for the gravitational-wave inspiral. In post-Newtonian theory this means neglecting the radiation-reaction force at 2.5PN and 3.5PN orders, i.e., considering only the conservative dynamics at the even-parity 1PN, 2PN and 3PN orders. We have seen in Sections 5.2 and 5.4 that this clean separation between conservative even-parity and dissipative odd-parity terms breaks at 4PN order, because of a contribution originating from gravitational-wave tails in the radiation-reaction force. We expect that at any higher order 4PN, 4.5PN, 5PN, etc. there will be a mixture of conservative and dissipative effects; here we assume that at any higher order we can neglect the radiation-reaction dissipation effects.

<footnote n="58">⁵⁸ This tendency is in agreement with numerical and analytical self-force calculations [24, 287].</footnote>
Consider a system of two compact objects moving on circular orbits. We examine first the case of non-spinning objects. With exactly circular orbits the geometry admits a helical Killing vector (HKV) field $K^\alpha$, satisfying the Killing equation $\nabla^\alpha K^\beta + \nabla^\beta K^\alpha = 0$. Imposing the existence of the HKV is the rigorous way to implement the notion of circular orbits. A Killing vector is only defined up to an overall constant factor. The helical Killing vector $K^\alpha$ extends out to a large distance where the geometry is essentially flat. There,

$$K^\alpha \partial_\alpha = \partial_t + \Omega \partial_\varphi, \qquad (273)$$

in any natural coordinate system which respects the helical symmetry [370]. We let this equality define the overall constant factor, thereby specifying the Killing vector field uniquely. In Eq. (273) $\Omega$ denotes the angular frequency of the binary's circular motion.

An observer moving with one of the particles (say the particle 1), while orbiting around the other particle, would detect no change in the local geometry. Thus the four-velocity $u_1^\alpha$ of that particle is tangent to the Killing vector $K^\alpha$ evaluated at the location of the particle, which we denote by $K_1^\alpha$. A physical quantity is then defined as the constant of proportionality $u_1^T$ between these two vectors, namely

$$u_1^\alpha = u_1^T K_1^\alpha. \qquad (274)$$

The four-velocity of the particle is normalized by $(g_{\mu\nu})_1 u_1^\mu u_1^\nu = -1$, where $(g_{\mu\nu})_1$ denotes the metric at the location of the particle. For a self-gravitating compact binary system, the metric at point 1 is generated by the two particles and has to be regularized according to one of the self-field regularizations discussed in Section 6. It will in fact be sometimes more convenient to work with the inverse of $u_1^T$, denoted $z_1 \equiv 1/u_1^T$. From Eq. (274) we get

$$z_1 = -(u_1 K_1), \qquad (275)$$

where $(u_1 K_1) = (g_{\mu\nu})_1 u_1^\mu K_1^\nu$ denotes the usual space-time dot product. Thus we can regard $z_1$ as the Killing energy of the particle that is associated with the HKV field $K^\alpha$. The quantity $z_1$ represents also the redshift of light rays emitted from the particle and received on the helical symmetry axis perpendicular to the orbital plane at large distances from it [176]. In the following we shall refer to $z_1$ as the redshift observable.

If we choose a coordinate system such that Eq. (273) is satisfied everywhere, then in particular $K_1^t = 1$, thus $u_1^T$ simply agrees with $u_1^t$, the t-component of the four-velocity of the particle. The Killing vector on the particle is then $K_1^\alpha = u_1^\alpha/u_1^t$, and simply reduces to the particle's ordinary coordinate velocity: $K_1^\alpha = v_1^\alpha/c$ where $v_1^\alpha = dy_1^\alpha/dt$ and $y_1^\alpha(t) = [ct, y_1(t)]$ denotes the particle's trajectory in that coordinate system. The redshift observable we are thus considering is

$$z_1 = \frac{1}{u_1^T} = \sqrt{-(g_{\mu\nu})_1 v_1^\mu v_1^\nu / c^2}. \qquad (276)$$

It is important to note that for circular orbits this quantity does not depend upon the choice of coordinates; in a perturbative approach in which the perturbative parameter is the particles' mass ratio $\nu \ll 1$, it does not depend upon the choice of perturbative gauge with respect to the background metric. We shall be interested in the invariant scalar function $z_1(\Omega)$, where $\Omega$ is the angular frequency of the circular orbit introduced when imposing Eq. (273).

We have obtained in Section 7.4 the expressions of the post-Newtonian binding energy $E$ and angular momentum $J$ for point-particle binaries on circular orbits. We shall now show that there are some differential and algebraic relations linking $E$ and $J$ to the redshift observables $z_1$ and $z_2$ associated with the two individual particles. Here we prefer to introduce instead of $E$ the total relativistic (ADM) mass of the binary system

$$M = m + \frac{E}{c^2}, \qquad (277)$$
where $m$ is the sum of the two post-Newtonian individual masses $m_1$ and $m_2$ – those which have been used up to now, for instance in Eq. (203). Note that in the spinning case such post-Newtonian masses acquire some spin contributions given, e.g., by Eqs. (243) – (246).

For point particles without spins, the ADM mass $M$, angular momentum $J$, and redshifts $z_a$, are functions of three independent variables, namely the orbital frequency $\Omega$ that is imposed by the existence of the HKV, and the individual masses $m_a$. For spinning point particles, we have also the two spins $S_a$ which are necessarily aligned with the orbital angular momentum. We first recover that the ADM quantities obey the “thermodynamical” relation already met in Eq. (235),

$$ \frac{\partial M}{\partial \Omega} = \Omega \frac{\partial J}{\partial \Omega}. \qquad (278) $$

Such relation is commonly used in post-Newtonian theory (see e.g., [160, 51]). It states that the gravitational-wave energy and angular momentum fluxes are strictly proportional for circular orbits, with $\Omega$ being the coefficient of proportionality. This relation is used in computations of the binary evolution based on a sequence of quasi-equilibrium configurations [228, 232, 133, 121], as discussed in Section 8.1.

The first law will be a thermodynamical generalization of Eq. (278), describing the changes in the ADM quantities not only when the orbital frequency $\Omega$ varies with fixed masses, but also when the individual masses $m_a$ of the particles vary with fixed orbital frequency. That is, one compares together different conservative dynamics with different masses but the same frequency. This situation is answered by the differential equations

$$ \frac{\partial M}{\partial m_a} - \Omega \frac{\partial J}{\partial m_a} = z_a \quad (a = 1, 2). \qquad (279) $$

Finally the three relations (278) – (279) can be summarized in the following result.

**Theorem 11.** The changes in the ADM mass and angular momentum of a binary system made of point particles in response to infinitesimal variations of the individual masses of the point particles, are related together by the first law of binary point-particle mechanics as [208, 289]

$$ \delta M - \Omega \delta J = \sum_a z_a \delta m_a. \qquad (280) $$

This law was proved in a very general way in Ref. [208] for systems of black holes and extended bodies under some arbitrary Killing symmetry. The particular form given in Eq. (280) is a specialization to the case of point particle binaries with helical Killing symmetry. It has been proved directly in this form in Ref. [289] up to high post-Newtonian order, namely 3PN order plus the logarithmic contributions occurring at 4PN and 5PN orders.

The first law of binary point-particle mechanics (280) is of course reminiscent of the celebrated first law of black hole mechanics $\delta M - \omega_H \delta J = \frac{\kappa}{5\pi} \delta A$, which holds for any non-singular, asymptotically flat perturbation of a stationary and axisymmetric black hole of mass $M$, intrinsic angular momentum (or spin) $J \equiv Ma$, surface area $A$, uniform surface gravity $\kappa$, and angular frequency $\omega_H$ on the horizon [26, 417]; see Ref. [289] for a discussion.

An interesting by-product of the first law (280) is the remarkably simple algebraic relation

$$ M - 2\Omega J = \sum_a z_a m_a, \qquad (281) $$

which can be seen as a first integral of the differential relation (280). Note that the existence of such a simple algebraic relation between the local quantities $z_1$ and $z_2$ on one hand, and the globally defined quantities $M$ and $J$ on the other hand, is not trivial.
Next, we report the result of a generalization of the first law applicable to systems of point particles with spins (moving on circular orbits).<fnref n="59" /> This result is valid through linear order in the spin of each particle, but holds also for the quadratic coupling between different spins (interaction spin terms $S_1 \times S_2$ in the language of Section 11). To be consistent with the HKV symmetry, we must assume that the two spins $S_a$ are aligned or anti-aligned with the orbital angular momentum. We introduce the total (ADM-like) angular momentum $J$ which is related to the orbital angular momentum $L$ by $J = L + \sum_a S_a$ for aligned or anti-aligned spins. The first law now becomes [55]

$$ \delta M - \Omega \delta J = \sum_a \left[ z_a \delta m_a + (\Omega_a - \Omega) \delta S_a \right], \qquad (282) $$

where $\Omega_a = |\Omega_a|$ denotes the precession frequency of the spins. This law has been derived in Ref. [55] from the canonical Hamiltonian formalism. The spin variables used here are the canonical spins $S_a$, that are easily seen to obey, from the algebra satisfied by the canonical variables, the usual Newtonian-looking precession equations $dS_a/dt = \Omega_a \times S_a$. These variables are identical to the “constant-in-magnitude” spins which will be defined and extensively used in Section 11. Similarly, to Eq. (281) we have also a first integral associated with the variational law (282):

$$ M - 2\Omega J = \sum_a \left[ z_a m_a - 2(\Omega - \Omega_a)S_a \right]. \qquad (283) $$

Notice that the relation (282) has been derived for point particles and arbitrary aligned spins. We would like now to derive the analogous relation for binary black holes. The key difference is that black holes are extended finite-size objects while point particles have by definition no spatial extension. For point particle binaries the spins can have arbitrary magnitude and still be compatible with the HKV. In this case the law (282) would describe also (super-extremal) naked singularities. For black hole binaries the HKV constraints the rotational state of each black hole and the binary system must be corotating.

Let us derive, in a heuristic way, the analogue of the first law (282) for black holes by introducing some “constitutive relations” $m_a(\mu_a, S_a, \dots)$ specifying the energy content of the bodies, i.e., the relations linking their masses $m_a$ to the spins $S_a$ and to some irreducible masses $\mu_a$. More precisely, we define for each spinning particle the analogue of an irreducible mass $\mu_a \equiv m_a^{\text{irr}}$ via the variational relation $\delta m_a = c_a \delta \mu_a + \omega_a \delta S_a$,<fnref n="60" /> in which the “response coefficient” $c_a$ of the body and its proper rotation frequency $\omega_a$ are associated with the internal structure:

$$ c_a = \left. \frac{\partial m_a}{\partial \mu_a} \right|_{S_a}, \qquad (284a) $$

$$ \omega_a = \left. \frac{\partial m_a}{\partial S_a} \right|_{\mu_a}. \qquad (284b) $$

For instance, using the Christodoulou mass formula (243) for Kerr black holes, we obtain the rotation frequency $\omega_a$ given by Eq. (244). On the other hand, the response coefficient $c_a$ differs from 1 only because of spin effects, and we can check that $c_a = 1 + \mathcal{O}(S_a^2)$.

Within the latter heuristic model a condition for the corotation of black hole binaries has been proposed in Ref. [55] as

$$ z_a \omega_a + \Omega_a = \Omega. \qquad (285) $$

This condition determines the value of the proper frequency $\omega_a$ of each black hole appropriate to the corotation state. When expanded to 2PN order the condition (285) leads to Eq. (247)

<footnote n="59">⁵⁹ The first law (280) has also been generalized for binary systems of point masses moving along generic stable bound (eccentric) orbits in Ref. [286].</footnote>

<footnote n="60">⁶⁰ In the case of extended material bodies, $\mu_a$ would represent the baryonic mass of the bodies.</footnote>
that we have already used in Section 8.1. With Eq. (285) imposed, the first law (282) simplifies considerably:

$$ \delta M - \Omega \delta J = \sum_a c_a z_a \delta \mu_a. \qquad (286) $$

This is almost identical to the first law for non-spinning binaries given by Eq. (280); indeed it simply differs from it by the substitutions $c_a \to 1$ and $\mu_a \to m_a$. Since the irreducible mass $\mu_a$ of a rotating black hole is the spin-independent part of its total mass $m_a$, this observation suggests that corotating binaries are very similar to non-spinning binaries, at least from the perspective of the first law. Finally we can easily reconcile the first law (286) for corotating systems with the known first law of binary black hole mechanics [208], namely

$$ \delta M - \Omega \delta J = \sum_a \frac{\kappa_a}{8\pi} \delta A_a. \qquad (287) $$

Indeed it suffices to make the formal identification in Eq. (286) of $c_a z_a$ with $4\mu_a \kappa_a$, where $\kappa_a$ denotes the constant surface gravity, and using the surface areas $A_a = 16\pi\mu_a^2$ instead of the irreducible masses of the black holes. This shows that the heuristic model based on the constitutive relations (284) is able to capture the physics of corotating black hole binary systems.

## 8.4 Post-Newtonian approximation versus gravitational self-force

The high-accuracy predictions from general relativity we have drawn up to now are well suited to describe the inspiralling phase of compact binaries, when the post-Newtonian parameter (1) is small independently of the mass ratio $q \equiv m_1/m_2$ between the compact bodies. In this section we investigate how well does the post-Newtonian expansion compare with another very important approximation scheme in general relativity: The gravitational self-force approach, based on black-hole perturbation theory, which gives an accurate description of extreme mass ratio binaries having $q \ll 1$ or equivalently $\nu \ll 1$, even in the strong field regime. The gravitational self-force analysis [317, 360, 178, 231] (see [348, 177, 23] for reviews) is thus expected to provide templates for extreme mass ratio inspirals (EMRI) anticipated to be present in the bandwidth of space-based detectors.

Consider a system of two (non-spinning) compact objects with $q = m_1/m_2 \ll 1$; we shall call the smaller mass $m_1$ the “particle”, and the larger mass $m_2$ the “black hole”. The orbit of the particle around the black hole is supposed to be exactly circular as we neglect the radiation-reaction effects. With circular orbits and no dissipation, we are considering the conservative part of the dynamics, and the geometry admits the HKV field (273). Note that in self-force theory there is a clean split between the dissipative and conservative parts of the dynamics (see e.g., [22]). This split is particularly transparent for an exact circular orbit, since the radial component (along $r$) is the only non-vanishing component of the conservative self-force, while the dissipative part of the self-force are the components along $t$ and $\varphi$.

The problem of the comparison between the post-Newtonian and perturbative self-force analyses in their common domain of validity, that of the slow-motion and weak-field regime of an extreme mass ratio binary, is illustrated in Figure 2. This problem has been tackled by Detweiler [176], who computed numerically within the self-force (SF) approach the redshift observable $u_1^T$ associated with the particle, and compared it with the 2PN prediction extracted from existing post-Newtonian results [76]. This comparison proved to be successful, and was later systematically implemented and extended to higher post-Newtonian orders in Refs. [68, 67]. In this section we review the works [68, 67] which have demonstrated an excellent agreement between the analytical post-Newtonian result derived through 3PN order, with inclusion of specific logarithmic terms at 4PN and 5PN orders, and the exact numerical SF result.
<figure><img src="image_6.png" /><figcaption>**Figure 2:** Different analytical approximation schemes and numerical techniques to study black hole binaries, depending on the mass ratio $q = m_1/m_2$ and the post-Newtonian parameter $\epsilon^2 \sim v^2/c^2 \sim Gm/(c^2r_{12})$. Post-Newtonian theory and perturbative self-force analysis can be compared in the post-Newtonian regime ($\epsilon \ll 1$ thus $r_{12} \gg Gm/c^2$) of an extreme mass ratio ($m_1 \ll m_2$) binary.</figcaption></figure>

For the PN-SF comparison, we require two physical quantities which are precisely defined in the context of each of the approximation schemes. The orbital frequency $\Omega$ of the circular orbit as measured by a distant observer is one such quantity and has been introduced in Eq. (273); the second quantity is the redshift observable $u_1^T$ (or equivalently $z_1 = 1/u_1^T$) associated with the smaller mass $m_1 \ll m_2$ and defined by Eqs. (274) or (275). The truly coordinate and perturbative-gauge independent properties of $\Omega$ and the redshift observable $u_1^T$ play a crucial role in this comparison. In the perturbative self-force approach we use Schwarzschild coordinates for the background, and we refer to “gauge invariance” as a property which holds within the restricted class of gauges for which (273) is a helical Killing vector. In all other respects, the gauge choice is arbitrary. In the post-Newtonian approach we work with harmonic coordinates and compute the explicit expression (276) of the redshift observable.

The main difficulty in the post-Newtonian calculation is the control to high PN order of the near-zone metric $(g_{\mu\nu})_1$ entering the definition of the redshift observable (276), and which has to be regularized at the location of the particle by means of dimensional regularization (see Sections 6.3–6.4). Up to 2.5PN order the Hadamard regularization is sufficient and the regularized metric has been provided in Eqs. (242). Here we report the end result of the post-Newtonian computation of the redshift observable including all terms up to the 3PN order, and augmented by the logarithmic contributions up to the 5PN order (and also the known Schwarzschild limit) [68, 67, 289]:

$$u_1^T = 1 + \left(\frac{3}{4} - \frac{3}{4}\Delta - \frac{\nu}{2}\right)x + \left(\frac{27}{16} - \frac{27}{16}\Delta - \frac{5}{2}\nu + \frac{5}{8}\Delta\nu + \frac{1}{24}\nu^2\right)x^2 + \left(\frac{135}{32} - \frac{135}{32}\Delta - \frac{37}{4}\nu - \frac{67}{16}\Delta\nu + \frac{115}{32}\nu^2 - \frac{5}{32}\Delta\nu^2 + \frac{1}{48}\nu^3\right)x^3$$
$$ 
\begin{aligned}
& + \left( \frac{2835}{256} - \frac{2835}{256} \Delta + \left[ -\frac{2183}{48} + \frac{41}{64} \pi^2 \right] \nu + \left[ \frac{12199}{384} - \frac{41}{64} \pi^2 \right] \Delta \nu \right. \\
& \left. + \left[ \frac{17201}{576} - \frac{41}{192} \pi^2 \right] \nu^2 - \frac{795}{128} \Delta \nu^2 - \frac{2827}{864} \nu^3 - \frac{25}{1728} \Delta \nu^3 + \frac{35}{10368} \nu^4 \right) x^4 \\
& + \left( \frac{15309}{512} - \frac{15309}{512} \Delta + \nu [u_4(\nu) + \Delta v_4(\nu)] + \left[ -\frac{32}{5} + \frac{32}{5} \Delta - \frac{32}{15} \nu \right] \nu \ln x \right) x^5 \\
& + \left( \frac{168399}{2048} - \frac{168399}{2048} \Delta + \nu [u_5(\nu) + \Delta v_5(\nu)] \right. \\
& \left. + \left[ \frac{478}{105} - \frac{478}{105} \Delta + \frac{10306}{105} \nu - 36 \Delta \nu - \frac{296}{15} \nu^2 \right] \nu \ln x \right) x^6 + \mathcal{O}\left( \frac{1}{c^{14}} \right).
\end{aligned} 
\quad (288) $$

We recall that $x$ denotes the post-Newtonian parameter (230), $\nu$ is the mass ratio (215), and $\Delta = (m_1 - m_2)/m$. The redshift observable of the other particle is deduced by setting $\Delta \to -\Delta$.

In Eq. (288) we denote by $u_4(\nu)$, $v_4(\nu)$ and $u_5(\nu)$, $v_5(\nu)$ some unknown 4PN and 5PN coefficients, which are however polynomials of the symmetric mass ratio $\nu$. They can be entirely determined from the related coefficients $e_4(\nu)$, $j_4(\nu)$ and $e_5(\nu)$, $j_5(\nu)$ in the expressions of the energy and angular momentum in Eqs. (233) and (234). To this aim it suffices to apply the differential first law (280) up to 5PN order; see Ref. [289] for more details.

The post-Newtonian result (288) is valid for any mass ratio, and for comparison purpose with the SF calculation we now investigate the small mass ratio regime $q \ll 1$. We introduce a post-Newtonian parameter appropriate to the small mass limit of the “particle”,

$$ y = \left( \frac{G m_2 \Omega}{c^3} \right)^{2/3} = x (1+q)^{-2/3}. \quad (289) $$

We express the symmetric mass ratio in terms of the asymmetric one: $\nu = q/(1+q)^2$, together with $\Delta = (q-1)/(q+1)$. Then Eq. (288), expanded through first order in $q$, which means including only the linear self-force level, reads

$$ u_1^T = u_{\text{Schw}}^T + q u_{\text{SF}}^T + \mathcal{O}(q^2). \quad (290) $$

The Schwarzschildian result is known in closed form as

$$ u_{\text{Schw}}^T = \frac{1}{\sqrt{1-3y}}, \quad (291) $$

and for the self-force contribution one obtains$^{61}$

$$ 
\begin{aligned}
u_{\text{SF}}^{T}(y) = & -y - 2y^{2} - 5y^{3} + \left( -\frac{121}{3} + \frac{41}{32} \pi^{2} \right) y^{4} \\
& + \left( \alpha_{4} - \frac{64}{5} \ln y \right) y^{5} + \left( \alpha_{5} + \frac{956}{105} \ln y \right) y^{6} + o(y^{6}).
\end{aligned} 
\quad (292) $$

The analytic coefficients were determined up to 2PN order in Ref. [176]; the 3PN term was computed in Ref. [68] making full use of dimensional regularization; the logarithmic contributions at the 4PN and 5PN orders were added in Refs. [67, 146].

The coefficients $\alpha_4$ and $\alpha_5$ represent some pure numbers at the 4PN and 5PN orders. By an analytic self-force calculation [36] the coefficient $\alpha_4$ has been obtained as

$$ \alpha_4 = -\frac{1157}{15} + \frac{677}{512} \pi^2 - \frac{256}{5} \ln 2 - \frac{128}{5} \gamma_E. \quad (293) $$

<footnote>$^{61}$ Since there are logarithms in this expansion we use the Landau $o$-symbol for remainders; see the footnote 20.</footnote>
Using the first law (280), we know how to deduce from the PN coefficients in the redshift variable the corresponding PN coefficients in the energy function (233). Thus, the result reported in Eq. (236) for the 4PN term in the energy function for circular orbits has been deduced from Eq. (293) by application of the first law.

On the self-force front the main problem is to control the numerical resolution of the computation of the redshift observable in order to distinguish more accurately the contributions of very high order PN terms. The comparison of the post-Newtonian expansion (292) with the numerical SF data has confirmed with high precision the determination of the 3PN coefficient [68, 67]: Witness Table 1 where the agreement with the analytical value involves 7 significant digits. Notice that such agreement provides an independent check of the dimensional regularization procedure invoked in the PN expansion scheme (see Sections 6.3–6.4). It is remarkable that such procedure is equivalent to the procedure of subtraction of the singular field in the SF approach [178].

<table><caption>**Table 1:** Numerically determined value of the 3PN coefficient for the SF part of the redshift observable defined by Eq. (292). The analytic post-Newtonian computation [68] is confirmed with many significant digits.</caption><thead><tr><th>3PN coefficient</th><th>SF value</th></tr></thead><tbody><tr><td>α<sub>3</sub> ≡ −<sup>121</sup>&frasl;<sub>3</sub> + <sup>41</sup>&frasl;<sub>32</sub>π<sup>2</sup> = −27.6879026 ... −27.6879034 ± 0.0000004</td><td></td></tr></tbody></table>

<table><caption>**Table 2:** Numerically determined values of higher-order PN coefficients for the SF part of the redshift observable defined by Eq. (292). The uncertainty in the last digit (or two last digits) is indicated in parentheses. The 4PN numerical value agrees with the analytical expression (293).</caption><thead><tr><th>PN coefficient</th><th>SF value</th></tr></thead><tbody><tr><td>α<sub>4</sub></td><td>−114.34747(5)</td></tr><tr><td>α<sub>5</sub></td><td>−245.53(1)</td></tr><tr><td>α<sub>6</sub></td><td>−695(2)</td></tr><tr><td>β<sub>6</sub></td><td>+339.3(5)</td></tr><tr><td>α<sub>7</sub></td><td>−5837(16)</td></tr></tbody></table>

Furthermore the PN-SF comparison has permitted to measure the coefficients $α_4$ and $α_5$ with at least 8 significant digits for the 4PN coefficient, and 5 significant digits for the 5PN one. In Table 2 we report the result of the analysis performed in Refs. [68, 67] by making maximum use of the analytical coefficients available at the time, i.e., all the coefficients up to 3PN order plus the logarithmic contributions at 4PN and 5PN orders. One uses a set of five basis functions corresponding to the unknown non-logarithmic 4PN and 5PN coefficients $α_4$, $α_5$ in Eq. (292), and augmented by the 6PN and 7PN non-logarithmic coefficients $α_6$, $α_7$ plus a coefficient $β_6$ for the logarithm at 6PN. A contribution $β_7$ from a logarithm at 7PN order is likely to confound with the $α_7$ coefficient. There is also the possibility of the contribution of a logarithmic squared at 7PN order, but such a small effect is not permitted in this fit.

Gladly we discover that the more recent analytical value of the 4PN coefficient, Eq. (293), matches the numerical value which was earlier measured in Ref. [67] (see Table 2). This highlights the predictive power of perturbative self-force calculations in determining numerically new post-Newtonian coefficients [176, 68, 67]. This ability is obviously due to the fact (illustrated in Figure 2) that perturbation theory is legitimate in the strong field regime of the coalescence of black hole binary systems, which is inaccessible to the post-Newtonian method. Of course, the limitation
of the self-force approach is the small mass-ratio limit; in this respect it is taken over by the post-Newtonian approximation.

More recently, the accuracy of the numerical computation of the self-force, and the comparison with the post-Newtonian expansion, have been drastically improved by Shah, Friedman & Whit- ing [383]. The PN coefficients of the redshift observable were obtained to very high 10.5PN order both numerically and also analytically, for a subset of coefficients that are either rational or made of the product of π with a rational. The analytical values of the coefficients up to 6PN order have also been obtained from an alternative self-force calculation [38, 37]. An interesting feature of the post-Newtonian expansion at high order is the appearance of half-integral PN coefficients (i.e., of the type n/2 PN where n is an odd integer) in the conservative dynamics of binary point particles, moving on exactly circular orbits. This is interesting because any instantaneous (non-tail) term at any half-integral PN order will be zero for circular orbits, as can be shown by a simple dimensional argument [77]. Therefore half-integral coefficients can appear only due to truly hereditary (tail) integrals. Using standard post-Newtonian methods it has been proved in Refs. [77, 78] that the dominant half-integral PN term in the redshift observable (292) occurs at the 5.5PN order (con- firming the earlier finding of Ref. [383]) and originates from the non-linear “tail-of-tail” integrals investigated in Section 3.2. The results for the 5.5PN coefficient in Eq. (292), and also for the next-to-leading 6.5PN and 7.5PN ones, are

$$
\alpha_{5.5} = -\frac{13696}{525}\pi, \quad \alpha_{6.5} = \frac{81077}{3675}\pi, \quad \alpha_{7.5} = \frac{82561159}{467775}\pi, \qquad (294)
$$

and fully agree between the PN and SF computations. We emphasize that the results (294) are achieved by the traditional PN approach, which is completely general (contrary to various analytical and numerical SF calculations [383, 38, 37, 268]), i.e., is not tuned to a particular type of source but is applicable to any extended PN source (see Part A). Note that Eqs. (294) represent the complete coefficients as there are no logarithms at these orders.

To conclude, the consistency of this “cross-cultural” comparison between the analytical post-Newtonian and the perturbative self-force approaches confirms the soundness of both approximations in describing the dynamics of compact binaries. Furthermore this interplay between PN and SF efforts (which is now rapidly growing [383]) is important for the synthesis of template waveforms of EMRs to be analysed by space-based gravitational wave detectors, and has also an impact on efforts of numerical relativity in the case of comparable masses.

9 Gravitational Waves from Compact Binaries

We pointed out that the 3.5PN equations of motion, Eqs. (203) or (219) – (220), are merely 1PN as regards the radiative aspects of the problem, because the radiation reaction force starts at the 2.5PN order. A solution would be to extend the precision of the equations of motion so as to include the full relative 3PN or 3.5PN precision into the radiation reaction force, but the equations of motion up to the 5.5PN or 6PN order are beyond the present state-of-the-art. The much better alternative solution is to apply the wave-generation formalism described in Part A, and to determine by its means the work done by the radiation reaction force directly as a total energy flux at future null infinity.<fnref n="62" /> In this approach, we replace the knowledge of the higher-order radiation reaction force by the computation of the total flux $\mathcal{F}$, and we apply the energy balance equation

$$
\frac{dE}{dt} = -\mathcal{F}. \qquad (295)
$$

Therefore, the result (232) that we found for the 3.5PN binary’s center-of-mass energy *E* constitutes only “half” of the solution of the problem. The second “half” consists of finding the rate of decrease

<footnote n="62">⁶² In addition, the wave generation formalism will provide the waveform itself, see Sections 9.4 and 9.5.</footnote>
$dE/dt$, which by the balance equation is nothing but the total gravitational-wave flux $\mathcal{F}$ at the relative 3.5PN order beyond the Einstein quadrupole formula (4).

Because the orbit of inspiralling binaries is circular, the energy balance equation is sufficient, and there is no need to invoke the angular momentum balance equation for computing the evolution of the orbital period $\dot{P}$ and eccentricity $\dot{e}$, see Eqs. (9)–(13) in the case of the binary pulsar. Furthermore the time average over one orbital period as in Eqs. (9) is here irrelevant, and the energy and angular momentum fluxes are related by $\mathcal{F} = \Omega G$. This all sounds good, but it is important to remind that we shall use the balance equation (295) at the very high 3.5PN order, and that at such order one is missing a complete proof of it (following from first principles in general relativity). Nevertheless, in addition to its physically obvious character, Eq. (295) has been verified by radiation-reaction calculations, in the cases of point-particle binaries [258, 259] and extended post-Newtonian fluids [43, 47], at the 1PN order and even at 1.5PN, the latter order being especially important because of the first appearance of wave tails; see Section 5.4. One should also quote here Refs. [260, 336, 278, 322, 254] for the 3.5PN terms in the binary’s equations of motion, fully consistent with the balance equations.

Obtaining the energy flux $\mathcal{F}$ can be divided into two equally important steps: Computing the source multipole moments $I_L$ and $J_L$ of the compact binary system with due account of a self-field regularization; and controlling the tails and related non-linear effects occurring in the relation between the binary’s source moments and the radiative ones $U_L$ and $V_L$ observed at future null infinity (cf. the general formalism of Part A).

## 9.1 The binary’s multipole moments

The general expressions of the source multipole moments given by Theorem 6, Eqs. (123), are worked out explicitly for general fluid systems at the 3.5PN order. For this computation one uses the formula (126), and we insert the 3.5PN metric coefficients (in harmonic coordinates) expressed in Eqs. (144) by means of the retarded-type elementary potentials (146)–(148). Then we specialize each of the (quite numerous) terms to the case of point-particle binaries by inserting, for the matter stress-energy tensor $T^{\alpha\beta}$, the standard expression made out of Dirac delta-functions. In Section 11 we shall consider spinning point particle binaries, and in that case the stress-energy tensor is given by Eq. (378). The infinite self-field of point-particles is removed by means of the Hadamard regularization; and, as we discussed in Section 6.4, dimensional regularization is used to fix the values of a few ambiguity parameters. This computation has been performed in Ref. [86] at the 1PN order, and in [64] at the 2PN order; we report below the most accurate 3PN results obtained in Refs. [81, 80, 62, 63] for the flux and [11, 74, 197] for the waveform.

A difficult part of the analysis is to find the closed-form expressions, fully explicit in terms of the particle’s positions and velocities, of many non-linear integrals. Let us give a few examples of the type of technical formulas that are employed in this calculation. Typically we have to compute some integrals like

$$Y_L^{(p,q)}(\mathbf{y}_1, \mathbf{y}_2) = -\frac{1}{2\pi} \mathcal{F}\mathcal{P} \int d^3\mathbf{x} \hat{x}_L r_1^p r_2^q, \qquad (296)$$

where $r_1 = |\mathbf{x} - \mathbf{y}_1|$ and $r_2 = |\mathbf{x} - \mathbf{y}_2|$. When $p > -3$ and $q > -3$, this integral is perfectly well-defined, since the finite part $\mathcal{F}\mathcal{P}$ deals with the IR regularization of the bound at infinity. When $p \le -3$ or $q \le -3$, we cure the UV divergencies by means of the Hadamard *partie finie* defined by Eq. (162); so a *partie finie* prescription $Pf$ is implicit in Eq. (296). An example of closed-form formula we get is

$$\dot{Y}_L^{(-1,-1)} = \frac{r_{12}}{\ell+1} \sum_{k=0}^{\ell} y_1^{(L-K)} y_2^K, \qquad (297)$$
where we pose $y_1^{L-K} \equiv y_1^{i_1} \cdots y_1^{i_{\ell-k}}$ and $y_2^K \equiv y_2^{i_{\ell-k+1}} \cdots y_2^{i_\ell}$, the brackets surrounding indices denoting the STF projection. Another example, in which the $\mathcal{FP}$ regularization is crucial, is (in the quadrupole case $\ell = 2$)

$$ (\begin{matrix} -2, & -1 \\ Y_{ij} & \end{matrix}) = y_1^{\langle ij \rangle} \left[ \frac{16}{15} \ln \left( \frac{r_{12}}{r_0} \right) - \frac{188}{225} \right] + y_1^{\langle i \rangle} y_2^{\langle j \rangle} \left[ \frac{8}{15} \ln \left( \frac{r_{12}}{r_0} \right) - \frac{4}{225} \right] + y_2^{\langle ij \rangle} \left[ \frac{2}{5} \ln \left( \frac{r_{12}}{r_0} \right) - \frac{2}{25} \right], \quad (298) $$

where the IR scale $r_0$ is defined in Eq. (42). Still another example, which necessitates both the $\mathcal{FP}$ and a UV partie finie regularization at the point 1, is

$$ Y_{ij}^{(-3,0)} = \left[ 2 \ln \left( \frac{s_1}{r_0} \right) + \frac{16}{15} \right] y_1^{\langle ij \rangle}, \quad (299) $$

where $s_1$ is the Hadamard-regularization constant introduced in Eq. (162).

The most important input for the computation of the waveform and flux is the mass quadrupole moment $I_{ij}$, since this moment necessitates the full post-Newtonian precision. Here we give the mass quadrupole moment complete to order 3.5PN, for non-spinning compact binaries on circular orbits, as

$$ I_{ij} = \mu \left( A x_{\langle ij \rangle} + B \frac{r^2}{c^2} v^{\langle ij \rangle} + \frac{48}{7} \frac{G^2 m^2 \nu}{c^5 r} C x_{\langle i \rangle} v_{\langle j \rangle} \right) + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (300) $$

where $x = y_1 - y_2 = (x_i)$ and $v = v_1 - v_2 = (v_i)$ are the orbital separation and relative velocity. The third term with coefficient $C$ is a radiation-reaction term, which will affect the waveform at orders 2.5PN and 3.5PN; however it does not contribute to the energy flux for circular orbits. The two conservative coefficients are $A$ and $B$. All those coefficients are [81, 74, 197]

$$ A = 1 + \gamma \left( -\frac{1}{42} - \frac{13}{14}\nu \right) + \gamma^2 \left( -\frac{461}{1512} - \frac{18395}{1512}\nu - \frac{241}{1512}\nu^2 \right) + \gamma^3 \left( \frac{395899}{13200} - \frac{428}{105}\ln\left(\frac{r_{12}}{r_0}\right) + \left[ \frac{3304319}{166320} - \frac{44}{3}\ln\left(\frac{r_{12}}{r_0}\right) \right]\nu + \frac{162539}{16632}\nu^2 + \frac{2351}{33264}\nu^3 \right), \quad (301a) $$

$$ B = \frac{11}{21} - \frac{11}{7}\nu + \gamma \left( \frac{1607}{378} - \frac{1681}{378}\nu + \frac{229}{378}\nu^2 \right) + \gamma^2 \left( -\frac{357761}{19800} + \frac{428}{105}\ln\left(\frac{r_{12}}{r_0}\right) - \frac{92339}{5544}\nu + \frac{35759}{924}\nu^2 + \frac{457}{5544}\nu^3 \right), \quad (301b) $$

$$ C = 1 + \gamma \left( -\frac{256}{135} - \frac{1532}{405}\nu \right). \quad (301c) $$

These expressions are valid in harmonic coordinates via the post-Newtonian parameter $\gamma = \frac{Gm}{rc^2}$ defined in Eq. (225). As we see, there are two types of logarithms at 3PN order in the quadrupole moment: One type involves the UV length scale $r'_0$ related by Eq. (221) to the two gauge constants $r'_1$ and $r'_2$ present in the 3PN equations of motion; the other type contains the IR length scale $r_0$ coming from the general formalism of Part A – indeed, recall that there is a $\mathcal{FP}$ operator in front of the source multipole moments in Theorem 6. As we know, the UV scale $r'_0$ is specific to the standard harmonic (SH) coordinate system and is pure gauge (see Section 7.3): It will disappear from our physical results at the end. On the other hand, we have proved that the multipole expansion outside a general post-Newtonian source is actually free of $r_0$, since the $r_0$'s present in the two terms of Eq. (105) cancel out. Indeed we have already found in Eqs. (93)–(94) that the constant $r_0$ present in $I_{ij}$ is compensated by the same constant coming from the non-linear wave “tails of tails” in the radiative moment $U_{ij}$. For a while, the expressions (301) contained the ambiguity parameters $\xi$, $\kappa$ and $\zeta$, which have now been replaced by their correct values (173).
Besides the 3.5PN mass quadrupole (300)-(301), we need also (for the 3PN waveform) the mass octupole moment $I_{ijk}$ and current quadrupole moment $J_{ij}$, both of them at the 2.5PN order; these are given for circular orbits by [81, 74]

$$ I_{ijk} = -\nu m \Delta \left\{ x_{\langle ijk \rangle} \left[ 1 - \gamma\nu - \gamma^2 \left( \frac{139}{330} + \frac{11923}{660}\nu + \frac{29}{110}\nu^2 \right) \right] \right. \\ \left. + \frac{r^2}{c^2} x_{\langle i v_{jk} \rangle} \left[ 1 - 2\nu - \gamma \left( -\frac{1066}{165} + \frac{1433}{330}\nu - \frac{21}{55}\nu^2 \right) \right] \right. \\ \left. + \frac{196}{15} \frac{r}{c} \gamma^2 \nu x_{\langle ij v_k \rangle} \right\} + \mathcal{O} \left( \frac{1}{c^6} \right), \quad (302a) $$

$$ J_{ij} = -\nu m \Delta \left\{ \epsilon_{ab\langle i x_j \rangle a} v_b \left[ 1 + \gamma \left( \frac{67}{28} - \frac{2}{7}\nu \right) + \gamma^2 \left( \frac{13}{9} - \frac{4651}{252}\nu - \frac{1}{168}\nu^2 \right) \right] \right. \\ \left. - \frac{188}{35} \frac{r}{c} \gamma^2 \nu \epsilon_{ab\langle i v_j \rangle a} x_b \right\} + \mathcal{O} \left( \frac{1}{c^6} \right). \quad (302b) $$

The list of required source moments for the 3PN waveform continues with the 2PN mass 2$^4$-pole and current 2$^3$-pole (octupole) moments, and so on. Here we give the most updated moments:$^{63}$

$$ I_{ijkl} = \nu m \left\{ x_{\langle i j k l \rangle} \left[ 1 - 3\nu + \gamma \left( \frac{3}{110} - \frac{25}{22}\nu + \frac{69}{22}\nu^2 \right) \right. \right. \\ \left. + \gamma^2 \left( -\frac{126901}{200200} - \frac{58101}{2600}\nu + \frac{204153}{2860}\nu^2 + \frac{1149}{1144}\nu^3 \right) \right] \\ \left. + \frac{r^2}{c^2} x_{\langle i j v_{kl} \rangle} \left[ \frac{78}{55}(1 - 5\nu + 5\nu^2) \right. \right. \\ \left. + \gamma \left( \frac{30583}{3575} - \frac{107039}{3575}\nu + \frac{8792}{715}\nu^2 - \frac{639}{715}\nu^3 \right) \right] \\ \left. + \frac{71}{715} \frac{r^4}{c^4} v_{\langle i j k l \rangle} (1 - 7\nu + 14\nu^2 - 7\nu^3) \right\} + \mathcal{O} \left( \frac{1}{c^5} \right), \quad (303a) $$

$$ J_{ijk} = \nu m \left\{ \epsilon_{ab\langle i x_{jk} \rangle a} v_b \left[ 1 - 3\nu + \gamma \left( \frac{181}{90} - \frac{109}{18}\nu + \frac{13}{18}\nu^2 \right) \right. \right. \\ \left. + \gamma^2 \left( \frac{1469}{3960} - \frac{5681}{264}\nu + \frac{48403}{660}\nu^2 - \frac{559}{3960}\nu^3 \right) \right] \\ \left. + \frac{r^2}{c^2} \epsilon_{ab\langle i x_a v_{jk} \rangle b} \left[ \frac{7}{45}(1 - 5\nu + 5\nu^2) + \gamma \left( \frac{1621}{990} - \frac{4879}{990}\nu + \frac{1084}{495}\nu^2 - \frac{259}{990}\nu^3 \right) \right] \right\} \\ \left. + \mathcal{O} \left( \frac{1}{c^5} \right) . \quad (303b) $$

$$ I_{ijklm} = -\nu m \Delta \left\{ x_{\langle i j k l m \rangle} \left[ 1 - 2\nu + \gamma \left( \frac{2}{39} - \frac{47}{39}\nu + \frac{28}{13}\nu^2 \right) \right] \right. \\ \left. + \frac{70}{39} \frac{r^2}{c^2} x_{\langle i j k v_{lm} \rangle} (1 - 4\nu + 3\nu^2) \right\} + \mathcal{O} \left( \frac{1}{c^4} \right), \quad (303c) $$

$$ J_{ijkl} = -\nu m \Delta \left\{ \epsilon_{ab\langle i x_{jk} \rangle a} v_b \left[ 1 - 2\nu + \gamma \left( \frac{20}{11} - \frac{155}{44}\nu + \frac{5}{11}\nu^2 \right) \right] \right. \\ \left. + \frac{4}{11} \frac{r^2}{c^2} \epsilon_{ab\langle i x_{ja} v_{kl} \rangle b} (1 - 4\nu + 3\nu^2) \right\} + \mathcal{O} \left( \frac{1}{c^4} \right). \quad (303d) $$

<footnote>$^{63}$ The STF projection $\langle \rangle$ applies only on “living” indices $ijkl...$ but not on the summed up indices $a$ and $b$.</footnote>
$$ I_{ijklmn} = \nu m \left\{ x_{\langleijklmn\rangle} \left[ 1 - 5\nu + 5\nu^2 + \gamma \left( \frac{1}{14} - \frac{3}{2}\nu + 6\nu^2 - \frac{11}{2}\nu^3 \right) \right] \right. \\ \left. + \frac{15}{7} \frac{r^2}{c^2} x_{\langleijklv_{mn}\rangle} (1 - 7\nu + 14\nu^2 - 7\nu^3) \right\} + \mathcal{O} \left( \frac{1}{c^4} \right), \qquad (303e) $$

$$ J_{ijklm} = \nu m \left\{ \epsilon_{ab\langle i x_{jklm} \rangle a v_b} \left[ 1 - 5\nu + 5\nu^2 + \gamma \left( \frac{1549}{910} - \frac{1081}{130}\nu + \frac{107}{13}\nu^2 - \frac{29}{26}\nu^3 \right) \right] \right. \\ \left. + \frac{54}{91} \frac{r^2}{c^2} \epsilon_{ab\langle i x_{jka} v_{lm} \rangle b} (1 - 7\nu + 14\nu^2 - 7\nu^3) \right\} + \mathcal{O} \left( \frac{1}{c^4} \right). \qquad (303f) $$

$$ I_{ijklmno} = -\nu m \Delta(1 - 4\nu + 3\nu^2) x_{\langleijklmno\rangle} + \mathcal{O} \left( \frac{1}{c^2} \right), \qquad (303g) $$

$$ J_{ijklmn} = -\nu m \Delta(1 - 4\nu + 3\nu^2) \epsilon_{ab\langle i x_{jklmn} \rangle a v_b} + \mathcal{O} \left( \frac{1}{c^2} \right). \qquad (303h) $$

All the other higher-order moments are required at the Newtonian order, at which they are trivial to compute with result ($\forall \ell \in \mathbb{N}$)

$$ I_L = \nu m \sigma_\ell(\nu) x_{\langle L \rangle} + \mathcal{O} \left( \frac{1}{c} \right), \qquad (304a) $$

$$ J_{L-1} = \nu m \sigma_{\ell}(\nu) \epsilon_{ab\langle i_{l-1} x_{L-2} \rangle a v_b} + \mathcal{O} \left( \frac{1}{c} \right). \qquad (304b) $$

Here we introduce the useful notation $\sigma_\ell(\nu) \equiv X_2^{\ell-1} + (-)^\ell X_1^{\ell-1}$, where $X_1 = \frac{m_1}{m}$ and $X_2 = \frac{m_2}{m}$ are 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}(\nu) = \sum_{p=0}^{k-1} (-)^p \frac{2k-1}{2k-1-p} \binom{2k-1-p}{p} \nu^p, \qquad (305a) $$

$$ \sigma_{2k+1}(\nu) = -\Delta \sum_{p=0}^{k-1} (-)^p \binom{2k-1-p}{p} \nu^p, \qquad (305b) $$

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

## 9.2 Gravitational wave energy flux

The results (300)–(304) permit the control of the instantaneous part of the total energy flux, by which we mean that part of the flux which is generated solely by the source multipole moments, i.e., not counting the hereditary tail and related integrals. The instantaneous flux $\mathcal{F}_{\text{inst}}$ is defined by the replacement into the general expression of $\mathcal{F}$ given by Eq. (68a) of all the radiative moments $U_L$ and $V_L$ by the corresponding $\ell$-th time derivatives of the source moments $I_L$ and $J_L$. Up to the 3.5PN order we have

$$ \mathcal{F}_{\text{inst}} = \frac{G}{c^5} \left\{ \frac{1}{5} I_{ij}^{(3)} I_{ij}^{(3)} + \frac{1}{c^2} \left[ \frac{1}{189} I_{ijk}^{(4)} I_{ijk}^{(4)} + \frac{16}{45} J_{ij}^{(3)} J_{ij}^{(3)} \right] + \frac{1}{c^4} \left[ \frac{1}{9072} I_{ijkm}^{(5)} I_{ijkm}^{(5)} + \frac{1}{84} J_{ijk}^{(4)} J_{ijk}^{(4)} \right] \right. \\ \left. + \frac{1}{c^6} \left[ \frac{1}{594000} I_{ijklmn}^{(6)} I_{ijklmn}^{(6)} + \frac{4}{14175} J_{ijklm}^{(5)} J_{ijklm}^{(5)} \right] + \mathcal{O} \left( \frac{1}{c^8} \right) \right\}. \qquad (306) $$

in which we insert the explicit expressions (300)–(304) for the moments. The time derivatives of these source moments are computed by means of the circular-orbit equations of motion given by
Eq. (226) together with (228). The net result is

$$ \mathcal{F}_{\text{inst}} = \frac{32c^5}{5G} \nu^2 \gamma^5 \left\{ 1 + \left( -\frac{2927}{336} - \frac{5}{4}\nu \right) \gamma + \left( \frac{293383}{9072} + \frac{380}{9}\nu \right) \gamma^2 \right. \\ \left. + \left[ \frac{53712289}{1108800} - \frac{1712}{105} \ln\left(\frac{r_{12}}{r_0}\right) \right. \right. \\ \left. \left. + \left( -\frac{50625}{112} + \frac{123}{64}\pi^2 + \frac{110}{3}\ln\left(\frac{r_{12}}{r'_0}\right) \right) \nu - \frac{383}{9}\nu^2 \right] \gamma^3 + \mathcal{O}\left(\frac{1}{c^8}\right) \right\}. \quad (307) $$

The Newtonian approximation agrees with the prediction of the Einstein quadrupole formula (4), as reduced for quasi-circular binary orbits by Landau & Lifshitz [285]. At the 3PN order in Eq. (307), there was some Hadamard regularization ambiguity parameters which have been replaced by their values computed with dimensional regularization.

To the instantaneous part of the flux, we must add the contribution of non-linear multipole interactions contained in the relationship between the source and radiative moments. The needed material has already been provided in Sections 3.3. Similar to the decomposition of the radiative quadrupole moment in Eq. (88), we can split the energy flux into the different terms

$$ \mathcal{F} = \mathcal{F}_{\text{inst}} + \mathcal{F}_{\text{tail}} + \mathcal{F}_{\text{tail-tail}}, \quad (308) $$

where $\mathcal{F}_{\text{inst}}$ has just been obtained in Eq. (307); $\mathcal{F}_{\text{tail}}$ is made of the quadratic multipolar tail integrals in Eqs. (90) and (95); $\mathcal{F}_{\text{tail-tail}}$ involves the square of the quadrupole tail in Eq. (90) and the quadrupole tail of tail given in Eq. (91).

We shall see that the tails play a crucial role in the predicted signal of compact binaries. It is quite remarkable that so small an effect as a “tail of tail” should be relevant to the data analysis of the current generation of gravitational wave detectors. By contrast, the non-linear memory effects, given by the integrals inside the 2.5PN and 3.5PN terms in Eq. (92), do not contribute to the gravitational-wave energy flux before the 4PN order in the case of circular-orbit binaries. Indeed the memory integrals are actually “instantaneous” in the flux, and a simple general argument based on dimensional analysis shows that instantaneous terms cannot contribute to the energy flux for circular orbits.<fnref n="64" /> Therefore the memory effect has rather poor observational consequences for future detections of inspiralling compact binaries.

Let us also recall that following the general formalism of Part A, the mass M which appears in front of the tail integrals of Sections 3.2 and 3.3 represents 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 binary’s mass-energy $E/c^2$ given for instance by Eq. (229). At 3.5PN order we need 2PN corrections in the tails and therefore 2PN also in the mass M, thus

$$ M = m \left[ 1 - \frac{\nu}{2}\gamma + \frac{\nu}{8}(7-\nu)\gamma^2 + \mathcal{O}\left(\frac{1}{c^6}\right) \right]. \quad (309) $$

Notice that 2PN order in M corresponds to 1PN order in E.

We give the two basic technical formulas needed when carrying out the reduction of the tail and tail-of-tail integrals to circular orbits (see e.g., [230]):

$$ \int_{0}^{+\infty} d\tau \ln\left(\frac{\tau}{2\tau_0}\right) e^{-i\omega\tau} = \frac{i}{\omega} \left( \ln(2\omega\tau_0) + \gamma_E + i\frac{\pi}{2} \right), \quad (310a) $$

<footnote n="64">⁶⁴ The same argument shows that the non-linear multipole interactions in Eq. (89) as well as those in Eqs. (97) and (98) do not contribute to the flux for circular orbits.</footnote>
$$ \int_{0}^{+\infty} d\tau \ln^2\left(\frac{\tau}{2\tau_0}\right) e^{-i\omega\tau} = -\frac{i}{\omega} \left[ \left(\ln(2\omega\tau_0) + \gamma_E + i\frac{\pi}{2}\right)^2 + \frac{\pi^2}{6} \right], \quad (310b) $$

where $\omega > 0$ is a strictly positive frequency (a multiple of the orbital frequency $\Omega$), where $\tau_0 = r_0/c$ and $\gamma_E$ is the Euler constant.

Notice the important point that the tail (and tail-of-tail) integrals can be evaluated, thanks to these formulas, for a fixed (i.e., non-decaying) circular orbit. Indeed it can be shown [60, 87] that the “remote-past” contribution to the tail integrals is negligible; the errors due to the fact that the orbit has actually evolved in the past, and spiraled in by emission of gravitational radiation, are of the order of the radiation-reaction scale $O(c^{-5})$,<fnref n="65" /> and do not affect the signal before the 4PN order. We then find, for the quadratic tails stricto sensu, the 1.5PN, 2.5PN and 3.5PN contributions

$$ \mathcal{F}_{\text{tail}} = \frac{32c^5}{5G} \gamma^5 \nu^2 \left\{ 4\pi\gamma^{3/2} + \left(-\frac{25663}{672} - \frac{125}{8}\nu\right) \pi\gamma^{5/2} + \left(\frac{90205}{576} + \frac{505747}{1512}\nu + \frac{12809}{756}\nu^2\right) \pi\gamma^{7/2} + \mathcal{O}\left(\frac{1}{c^8}\right) \right\}. \quad (311) $$

For the sum of squared tails and cubic tails of tails at 3PN, we get

$$ \mathcal{F}_{\text{tail-tail}} = \frac{32c^5}{5G} \gamma^5 \nu^2 \left\{ \left( -\frac{116761}{3675} + \frac{16}{3}\pi^2 - \frac{1712}{105}\gamma_E + \frac{1712}{105}\ln\left(\frac{r_{12}}{r_0}\right) - \frac{856}{105}\ln(16\gamma) \right) \gamma^3 + \mathcal{O}\left(\frac{1}{c^8}\right) \right\}. \quad (312) $$

By comparing Eqs. (307) and (312) we observe that the constants $r_0$ cleanly cancel out. Adding together these contributions we obtain

$$ \begin{aligned} \mathcal{F} = \frac{32c^5}{5G} \gamma^5 \nu^2 \left\{ \right. & 1 + \left(-\frac{2927}{336} - \frac{5}{4}\nu\right) \gamma + 4\pi\gamma^{3/2} \\ & + \left(\frac{293383}{9072} + \frac{380}{9}\nu\right) \gamma^2 + \left(-\frac{25663}{672} - \frac{125}{8}\nu\right) \pi\gamma^{5/2} \\ & + \left[ \frac{129386791}{7761600} + \frac{16\pi^2}{3} - \frac{1712}{105}\gamma_E - \frac{856}{105}\ln(16\gamma) \right. \\ & \qquad \left. + \left(-\frac{50625}{112} + \frac{110}{3}\ln\left(\frac{r_{12}}{r_0}\right) + \frac{123\pi^2}{64}\right) \nu - \frac{383}{9}\nu^2 \right] \gamma^3 \\ & + \left( \frac{90205}{576} + \frac{505747}{1512}\nu + \frac{12809}{756}\nu^2 \right) \pi\gamma^{7/2} + \mathcal{O}\left(\frac{1}{c^8}\right) \left. \right\}. \end{aligned} \quad (313) $$

The gauge constant $r'_0$ has not yet disappeared because the post-Newtonian expansion is still parametrized by $\gamma$ instead of the frequency-related parameter $x$ defined by Eq. (230) – just as for $E$ when it was given by Eq. (229). After substituting the expression $\gamma(x)$ given by Eq. (231), we find that $r'_0$ does cancel as well. Because the relation $\gamma(x)$ is issued from the equations of motion, the latter cancellation represents an interesting test of the consistency of the two computations, in harmonic coordinates, of the 3PN multipole moments and the 3PN equations of motion. At long last we obtain our end result:<fnref n="66" />

$$ \mathcal{F} = \frac{32c^5}{5G} \nu^2 x^5 \left\{ 1 + \left(-\frac{1247}{336} - \frac{35}{12}\nu\right) x + 4\pi x^{3/2} \right\}. $$

<footnote n="65">⁶⁵ Or, rather, $O(c^{-5} \ln c)$ as shown in the Appendix of Ref. [87].</footnote>

<footnote n="66">⁶⁶ See Section 10 for the generalization of the flux of energy to eccentric binary orbits.</footnote>
$$ + \left( -\frac{44711}{9072} + \frac{9271}{504}\nu + \frac{65}{18}\nu^2 \right) x^2 + \left( -\frac{8191}{672} - \frac{583}{24}\nu \right) \pi x^{5/2} \\ + \left[ \frac{6643739519}{69854400} + \frac{16}{3}\pi^2 - \frac{1712}{105}\gamma_E - \frac{856}{105}\ln(16x) \right. \\ + \left( -\frac{134543}{7776} + \frac{41}{48}\pi^2 \right) \nu - \frac{94403}{3024}\nu^2 - \frac{775}{324}\nu^3 \right] x^3 \\ + \left( -\frac{16285}{504} + \frac{214745}{1728}\nu + \frac{193385}{3024}\nu^2 \right) \pi x^{7/2} + \mathcal{O}\left( \frac{1}{c^8} \right) \right\} . \quad (314) $$

In the test-mass limit $\nu \to 0$ for one of the bodies, we recover exactly the result following from linear black-hole perturbations obtained by Tagoshi & Sasaki [395] (see also [393, 397]). In particular, the rational fraction $6643739519/69854400$ comes out exactly the same as in black-hole perturbations. On the other hand, the ambiguity parameters discussed in Section 6.2 were part of the rational fraction $-134543/7776$, belonging to the coefficient of the term at 3PN order proportional to $\nu$ (hence this coefficient cannot be computed by linear black-hole perturbations).

The effects due to the spins of the two black holes arise at the 1.5PN order for the spin-orbit (SO) coupling, and at the 2PN order for the spin-spin (SS) coupling, for maximally rotating black holes. Spin effects will be discussed in Section 11. On the other hand, the terms due to the radiating energy flowing into the black-hole horizons and absorbed rather than escaping to infinity, have to be added to the standard post-Newtonian calculation based on point particles as presented here; such terms arise at the 4PN order for Schwarzschild black holes [349] and at 2.5PN order for Kerr black holes [392].

## 9.3 Orbital phase 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 (295). The center-of-mass energy $E$ is given by Eq. (232) and the total flux $\mathcal{F}$ by Eq. (314). For convenience we adopt the dimensionless time variable$^{67}$

$$ \Theta \equiv \frac{\nu c^3}{5Gm} (t_c - t), \quad (315) $$

where $t_c$ denotes the instant of coalescence, at which the frequency formally tends to infinity, although evidently, the post-Newtonian method breaks down well before this point. We transform the balance equation into an ordinary differential equation for the parameter $x$, which is immediately integrated with the result

$$ x = \frac{1}{4}\Theta^{-1/4}\left\{1 + \left(\frac{743}{4032} + \frac{11}{48}\nu\right)\Theta^{-1/4} - \frac{1}{5}\pi\Theta^{-3/8}\right. \\ + \left(\frac{19583}{254016} + \frac{24401}{193536}\nu + \frac{31}{288}\nu^2\right)\Theta^{-1/2} + \left(-\frac{11891}{53760} + \frac{109}{1920}\nu\right)\pi\Theta^{-5/8} \\ + \left[ -\frac{10052469856691}{6008596070400} + \frac{1}{6}\pi^2 + \frac{107}{420}\gamma_E - \frac{107}{3360}\ln\left(\frac{\Theta}{256}\right) \right. \\ + \left( \frac{3147553127}{780337152} - \frac{451}{3072}\pi^2 \right)\nu - \frac{15211}{442368}\nu^2 + \frac{25565}{331776}\nu^3 \right] \Theta^{-3/4} \\ + \left( -\frac{113868647}{433520640} - \frac{31821}{143360}\nu + \frac{294941}{3870720}\nu^2 \right)\pi\Theta^{-7/8} + \mathcal{O}\left(\frac{1}{c^8}\right) \right\}. \quad (316) $$

<footnote>$^{67}$ Notice the “strange” post-Newtonian order of this time variable: $\Theta = \mathcal{O}(c^{+8})$.</footnote>
The orbital phase is defined as the angle $\phi$, oriented in the sense of the motion, between the separation of the two bodies and the direction of the ascending node (called $\mathcal{N}$ in Section 9.4) within the plane of the sky. We have $d\phi/dt = \Omega$, which translates, with our notation, into $d\phi/d\Theta = -5x^{3/2}/\nu$, from which we determine$^{68}$

$$ \phi = -\frac{1}{\nu}\Theta^{5/8} \left\{ 1 + \left(\frac{3715}{8064} + \frac{55}{96}\nu\right)\Theta^{-1/4} - \frac{3}{4}\pi\Theta^{-3/8} \right. \\ \qquad + \left(\frac{9275495}{14450688} + \frac{284875}{258048}\nu + \frac{1855}{2048}\nu^2\right)\Theta^{-1/2} + \left(-\frac{38645}{172032} + \frac{65}{2048}\nu\right)\pi\Theta^{-5/8}\ln\left(\frac{\Theta}{\Theta_0}\right) \\ \qquad + \left[ \frac{831032450749357}{57682522275840} - \frac{53}{40}\pi^2 - \frac{107}{56}\gamma_E + \frac{107}{448}\ln\left(\frac{\Theta}{256}\right) \right. \\ \qquad \qquad \left. + \left(-\frac{126510089885}{4161798144} + \frac{2255}{2048}\pi^2\right)\nu + \frac{154565}{1835008}\nu^2 - \frac{1179625}{1769472}\nu^3 \right] \Theta^{-3/4} \\ \qquad + \left(\frac{188516689}{173408256x} + \frac{488825}{516096}\nu - \frac{141769}{516096}\nu^2\right)\pi\Theta^{-7/8} + \mathcal{O}\left(\frac{1}{c^8}\right) \left. \right\}, \quad (317) $$

where $\Theta_0$ is a constant of integration that can be fixed by the initial conditions when the wave frequency enters the detector. Finally we want also to dispose of the important expression of the phase in terms of the frequency $x$. For this we get

$$ \phi = -\frac{x^{-5/2}}{32\nu} \left\{ 1 + \left(\frac{3715}{1008} + \frac{55}{12}\nu\right)x - 10\pi x^{3/2} \right. \\ \qquad + \left(\frac{15293365}{1016064} + \frac{27145}{1008}\nu + \frac{3085}{144}\nu^2\right)x^2 + \left(\frac{38645}{1344} - \frac{65}{16}\nu\right)\pi x^{5/2} \ln\left(\frac{x}{x_0}\right) \\ \qquad + \left[ \frac{12348611926451}{18776862720} - \frac{160}{3}\pi^2 - \frac{1712}{21}\gamma_E - \frac{856}{21}\ln(16x) \right. \\ \qquad \qquad \left. + \left(-\frac{15737765635}{12192768} + \frac{2255}{48}\pi^2\right)\nu + \frac{76055}{6912}\nu^2 - \frac{127825}{5184}\nu^3 \right] x^3 \\ \qquad + \left(\frac{77096675}{2032128} + \frac{378515}{12096}\nu - \frac{74045}{6048}\nu^2\right)\pi x^{7/2} + \mathcal{O}\left(\frac{1}{c^8}\right) \left. \right\}, \quad (318) $$

where $x_0$ is another constant of integration. With the formula (318) the orbital phase is complete up to the 3.5PN order for non-spinning compact binaries. Note that the contributions of the quadrupole moments of compact objects which are induced by tidal effects, are expected from Eq. (16) to come into play only at the 5PN order.

As a rough estimate of the relative importance of the various post-Newtonian terms, we give in Table 3 their contributions to the accumulated number of gravitational-wave cycles $\mathcal{N}_{cycle}$ in the bandwidth of ground-based detectors. Note that such an estimate is only indicative, because a full treatment would require the knowledge of the detector’s power spectral density of noise, and a complete simulation of the parameter estimation using matched filtering techniques [138, 350, 284]. We define $\mathcal{N}_{cycle}$ as

$$ \mathcal{N}_{cycle} \equiv \frac{\phi_{ISCO} - \phi_{seismic}}{\pi}. \quad (319) $$

The frequency of the signal at the entrance of the bandwidth is the seismic cut-off frequency $f_{seismic}$ of ground-based detectors; the terminal frequency is assumed for simplicity to be given

<footnote>$^{68}$ This procedure for computing analytically the orbital phase corresponds to what is called in the jargon the “Taylor T2 approximant”. We refer to Ref. [98] for discussions on the usefulness of defining several types of approximants for computing (in general numerically) the orbital phase.</footnote>
by the Schwarzschild innermost stable circular orbit: $f_{\text{ISCO}} = \frac{c^3}{6^{3/2}\pi Gm}$. Here we denote by $f = \Omega/\pi = 2/P$ the signal frequency of the dominant harmonics at twice the orbital frequency. As we see in Table 3, with the 3PN or 3.5PN approximations we reach an acceptable accuracy level of a few cycles say, that roughly corresponds to the demand made by data-analysts in the case of neutron-star binaries [139, 137, 138, 346, 105, 106]. Indeed, the above estimation suggests that the neglected 4PN terms will yield some systematic errors that are, at most, of the same order of magnitude, i.e., a few cycles, and perhaps much less.

<table><caption>**Table 3:** Post-Newtonian contributions to the accumulated number of gravitational-wave cycles $N_{\text{cycle}}$ for compact binaries detectable in the bandwidth of LIGO-VIRGO detectors. The entry frequency is $f_{\text{seismic}} = 10$ Hz and the terminal frequency is $f_{\text{ISCO}} = \frac{c^3}{6^{3/2}\pi Gm}$. The main origin of the approximation (instantaneous vs. tail) is indicated. See also Table 4 in Section 11 below for the contributions of spin-orbit effects.</caption><thead><tr><th>PN order</th><th></th><th>1.4 M<sub>☉</sub> + 1.4 M<sub>☉</sub></th><th>10 M<sub>☉</sub> + 1.4 M<sub>☉</sub></th><th>10 M<sub>☉</sub> + 10 M<sub>☉</sub></th></tr></thead><tbody><tr><td>N</td><td>(inst)</td><td>15952.6</td><td>3558.9</td><td>598.8</td></tr><tr><td>1PN</td><td>(inst)</td><td>439.5</td><td>212.4</td><td>59.1</td></tr><tr><td>1.5PN</td><td>(leading tail)</td><td>-210.3</td><td>-180.9</td><td>-51.2</td></tr><tr><td>2PN</td><td>(inst)</td><td>9.9</td><td>9.8</td><td>4.0</td></tr><tr><td>2.5PN</td><td>(1PN tail)</td><td>-11.7</td><td>-20.0</td><td>-7.1</td></tr><tr><td>3PN</td><td>(inst + tail-of-tail)</td><td>2.6</td><td>2.3</td><td>2.2</td></tr><tr><td>3.5PN</td><td>(2PN tail)</td><td>-0.9</td><td>-1.8</td><td>-0.8</td></tr></tbody></table>

## 9.4 Polarization waveforms for data analysis

The theoretical templates of the compact binary inspiral follow from insertion of the previous solutions for the 3.5PN-accurate orbital frequency and phase into the binary's two polarization waveforms $h_+$ and $h_-$ defined with respect to a choice of two polarization vectors $P = (P_i)$ and $Q = (Q_i)$ orthogonal to the direction $N$ of the observer; see Eqs. (69).

Our convention for the two polarization vectors is that they form with $N$ a right-handed triad, and that $P$ and $Q$ lie along the major and minor axis, respectively, of the projection onto the plane of the sky of the circular orbit. This means that $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 is also chosen for the origin of the orbital phase $\phi$. We denote by $i$ the inclination angle between the direction of the detector $N$ as seen from the binary's center-of-mass, and the normal to the orbital plane (we always suppose that the normal is right-handed with respect to the sense of motion, so that $0 \le i \le \pi$). We use the shorthands $c_i \equiv \cos i$ and $s_i \equiv \sin i$ for the cosine and sine of the inclination angle.

We shall include in $h_+$ and $h_-$ all the harmonics, besides the dominant one at twice the orbital frequency, consistent with the 3PN approximation [82, 11, 74]. In Section 9.5 we shall give all the modes $(\ell, m)$ in a spherical-harmonic decomposition of the waveform, and shall extend the dominant quadrupole mode $(2, 2)$ at 3.5PN order [197]. The post-Newtonian terms are ordered by means of the frequency-related variable $x = (\frac{Gm\Omega}{c^3})^{2/3}$; to ease the notation we pose

$$h_{+, \times} = \frac{2G\mu x}{c^2 R} \sum_{p=0}^{+\infty} x^{p/2} H_{+, \times}(\psi, c_i, s_i; \ln x) + \mathcal{O}\left(\frac{1}{R^2}\right). \qquad (320)$$

Note that the post-Newtonian coefficients will involve the logarithm $\ln x$ starting at 3PN order; see Eq. (127). They depend on the binary's phase $\phi$, explicitly given at 3.5PN order by Eq. (318),
through the very useful auxiliary phase variable $\psi$ that is “distorted by tails” [87, 11] and reads

$$ \psi \equiv \phi - \frac{2GM\Omega}{c^3} \ln \left( \frac{\Omega}{\Omega_0} \right). \qquad (321) $$

Here M denotes the binary’s ADM mass and it is very important to include all its relevant post-Newtonian contributions as given by Eq. (309). The constant frequency $\Omega_0$ can be chosen at will, for instance to be the entry frequency of some detector. For the plus polarization we have$^{69}$

$$ H_+ = -(1 + c_i^2) \cos 2\psi - \frac{1}{96} s_i^2 (17 + c_i^2), \qquad (322a) $$

$$ H_{1/2} = -s_i \Delta \left[ \cos \psi \left( \frac{5}{8} + \frac{1}{8} c_i^2 \right) - \cos 3\psi \left( \frac{9}{8} + \frac{9}{8} c_i^2 \right) \right], \qquad (322b) $$

$$ H_+ = \cos 2\psi \left[ \frac{19}{6} + \frac{3}{2} c_i^2 - \frac{1}{3} c_i^4 + \nu \left( -\frac{19}{6} + \frac{11}{6} c_i^2 + c_i^4 \right) \right] \\ - \cos 4\psi \left[ \frac{4}{3} s_i^2 (1 + c_i^2) (1 - 3\nu) \right], \qquad (322c) $$

$$ H_{3/2} = s_i \Delta \cos \psi \left[ \frac{19}{64} + \frac{5}{16} c_i^2 - \frac{1}{192} c_i^4 + \nu \left( -\frac{49}{96} + \frac{1}{8} c_i^2 + \frac{1}{96} c_i^4 \right) \right] \\ + \cos 2\psi \left[ -2\pi(1 + c_i^2) \right] \\ + s_i \Delta \cos 3\psi \left[ -\frac{657}{128} - \frac{45}{16} c_i^2 + \frac{81}{128} c_i^4 + \nu \left( \frac{225}{64} - \frac{9}{8} c_i^2 - \frac{81}{64} c_i^4 \right) \right] \\ + s_i \Delta \cos 5\psi \left[ \frac{625}{384} s_i^2 (1 + c_i^2) (1 - 2\nu) \right], \qquad (322d) $$

$$ H_+ = \pi s_i \Delta \cos \psi \left[ -\frac{5}{8} - \frac{1}{8} c_i^2 \right] \\ + \cos 2\psi \left[ \frac{11}{60} + \frac{33}{10} c_i^2 + \frac{29}{24} c_i^4 - \frac{1}{24} c_i^6 + \nu \left( \frac{353}{36} - 3 c_i^2 - \frac{251}{72} c_i^4 + \frac{5}{24} c_i^6 \right) \right. \\ \left. + \nu^2 \left( -\frac{49}{12} + \frac{9}{2} c_i^2 - \frac{7}{24} c_i^4 - \frac{5}{24} c_i^6 \right) \right] \\ + \pi s_i \Delta \cos 3\psi \left[ \frac{27}{8} (1 + c_i^2) \right] \\ + \frac{2}{15} s_i^2 \cos 4\psi \left[ 59 + 35 c_i^2 - 8 c_i^4 - \frac{5}{3} \nu (131 + 59 c_i^2 - 24 c_i^4) + 5 \nu^2 (21 - 3 c_i^2 - 8 c_i^4) \right] \\ + \cos 6\psi \left[ -\frac{81}{40} s_i^4 (1 + c_i^2) (1 - 5\nu + 5\nu^2) \right] \\ + s_i \Delta \sin \psi \left[ \frac{11}{40} + \frac{5 \ln 2}{4} + c_i^2 \left( \frac{7}{40} + \frac{\ln 2}{4} \right) \right] \\ + s_i \Delta \sin 3\psi \left[ \left( -\frac{189}{40} + \frac{27}{4} \ln(3/2) \right) (1 + c_i^2) \right], \qquad (322e) $$

<footnote>$^{69}$ Notice the obvious fact that the polarization waveforms remain invariant when we rotate by $\pi$ the separation direction between the particles and simultaneously exchange the labels of the two particles, i.e., when we apply the transformation $(\psi, \Delta) \to (\psi + \pi, -\Delta)$. Moreover, due to the parity invariance, the $H_+$’s are unchanged after the replacement $i \to \pi - i$, while the $H_\times$’s being the projection of $h_{ij}^{TT}$ on a tensorial product of two vectors of inverse parity types, is changed into its opposite.</footnote>
$$
\begin{equation}
\begin{aligned}
H_{+} = s_i \Delta \cos \psi & \left[ \frac{1771}{5120} - \frac{1667}{5120} c_i^2 + \frac{217}{9216} c_i^4 - \frac{1}{9216} c_i^6 \right. \\
& \left. \qquad + \nu \left( \frac{681}{256} + \frac{13}{768} c_i^2 - \frac{35}{768} c_i^4 + \frac{1}{2304} c_i^6 \right) \right. \\
& \left. \qquad + \nu^2 \left( -\frac{3451}{9216} + \frac{673}{3072} c_i^2 - \frac{5}{9216} c_i^4 - \frac{1}{3072} c_i^6 \right) \right] \\
& + \pi \cos 2\psi \left[ \frac{19}{3} + 3 c_i^2 - \frac{2}{3} c_i^4 + \nu \left( -\frac{16}{3} + \frac{14}{3} c_i^2 + 2 c_i^4 \right) \right] \\
& + s_i \Delta \cos 3\psi \left[ \frac{3537}{1024} - \frac{22977}{5120} c_i^2 - \frac{15309}{5120} c_i^4 + \frac{729}{5120} c_i^6 \right. \\
& \left. \qquad + \nu \left( -\frac{23829}{1280} + \frac{5529}{1280} c_i^2 + \frac{7749}{1280} c_i^4 - \frac{729}{1280} c_i^6 \right) \right. \\
& \left. \qquad + \nu^2 \left( \frac{29127}{5120} - \frac{27267}{5120} c_i^2 - \frac{1647}{5120} c_i^4 + \frac{2187}{5120} c_i^6 \right) \right] \\
& + \cos 4\psi \left[ -\frac{16\pi}{3} (1 + c_i^2) s_i^2 (1 - 3\nu) \right] \\
& + s_i \Delta \cos 5\psi \left[ -\frac{108125}{9216} + \frac{40625}{9216} c_i^2 + \frac{83125}{9216} c_i^4 - \frac{15625}{9216} c_i^6 \right. \\
& \left. \qquad + \nu \left( \frac{8125}{256} - \frac{40625}{2304} c_i^2 - \frac{48125}{2304} c_i^4 + \frac{15625}{2304} c_i^6 \right) \right. \\
& \left. \qquad + \nu^2 \left( -\frac{119375}{9216} + \frac{40625}{3072} c_i^2 + \frac{44375}{9216} c_i^4 - \frac{15625}{3072} c_i^6 \right) \right] \\
& + \Delta \cos 7\psi \left[ \frac{117649}{46080} s_i^5 (1 + c_i^2) (1 - 4\nu + 3\nu^2) \right] \\
& + \sin 2\psi \left[ -\frac{9}{5} + \frac{14}{5} c_i^2 + \frac{7}{5} c_i^4 + \nu \left( 32 + \frac{56}{5} c_i^2 - \frac{28}{5} c_i^4 \right) \right] \\
& + s_i^2 (1 + c_i^2) \sin 4\psi \left[ \frac{56}{5} - \frac{32 \ln 2}{3} + \nu \left( -\frac{1193}{30} + 32 \ln 2 \right) \right],
\end{aligned}
\tag{322f}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
H_3^+ ={}& \pi \Delta s_i \cos \psi \left[ \frac{19}{64} + \frac{5}{16} c_i^2 - \frac{1}{192} c_i^4 + \nu \left( -\frac{19}{96} + \frac{3}{16} c_i^2 + \frac{1}{96} c_i^4 \right) \right] \\
& + \cos 2\psi \left[ -\frac{465497}{11025} + \left( \frac{856}{105} \gamma_E - \frac{2\pi^2}{3} + \frac{428}{105} \ln(16x) \right) (1+c_i^2) \right. \\
& \qquad \left. - \frac{3561541}{88200} c_i^2 - \frac{943}{720} c_i^4 + \frac{169}{720} c_i^6 - \frac{1}{360} c_i^8 \right. \\
& \qquad \left. + \nu \left( \frac{2209}{360} - \frac{41\pi^2}{96} (1+c_i^2) + \frac{2039}{180} c_i^2 + \frac{3311}{720} c_i^4 - \frac{853}{720} c_i^6 + \frac{7}{360} c_i^8 \right) \right. \\
& \qquad \left. + \nu^2 \left( \frac{12871}{540} - \frac{1583}{60} c_i^2 - \frac{145}{108} c_i^4 + \frac{56}{45} c_i^6 - \frac{7}{180} c_i^8 \right) \right. \\
& \qquad \left. + \nu^3 \left( -\frac{3277}{810} + \frac{19661}{3240} c_i^2 - \frac{281}{144} c_i^4 - \frac{73}{720} c_i^6 + \frac{7}{360} c_i^8 \right) \right] \\
& + \pi \Delta s_i \cos 3\psi \left[ -\frac{1971}{128} - \frac{135}{16} c_i^2 + \frac{243}{128} c_i^4 + \nu \left( \frac{567}{64} - \frac{81}{16} c_i^2 - \frac{243}{64} c_i^4 \right) \right]
\end{aligned}
\end{equation}
$$
$$
\begin{align}
&+ s_i^2 \cos 4\psi \left[ -\frac{2189}{210} + \frac{1123}{210} c_i^2 + \frac{56}{9} c_i^4 - \frac{16}{45} c_i^6 \right. \nonumber \\
&\qquad \left. + \nu \left( \frac{6271}{90} - \frac{1969}{90} c_i^2 - \frac{1432}{45} c_i^4 + \frac{112}{45} c_i^6 \right) \right. \nonumber \\
&\qquad \left. + \nu^2 \left( -\frac{3007}{27} + \frac{3493}{135} c_i^2 + \frac{1568}{45} c_i^4 - \frac{224}{45} c_i^6 \right) + \nu^3 \left( \frac{161}{6} - \frac{1921}{90} c_i^2 - \frac{184}{45} c_i^4 + \frac{112}{45} c_i^6 \right) \right] \nonumber \\
&+ \Delta \cos 5\psi \left[ \frac{3125 \pi}{384} s_i^3 (1+c_i^2)(1-2\nu) \right] \nonumber \\
&+ s_i^4 \cos 6\psi \left[ \frac{1377}{80} + \frac{891}{80} c_i^2 - \frac{729}{280} c_i^4 + \nu \left( -\frac{7857}{80} - \frac{891}{16} c_i^2 + \frac{729}{40} c_i^4 \right) \right. \nonumber \\
&\qquad \left. + \nu^2 \left( \frac{567}{4} + \frac{567}{10} c_i^2 - \frac{729}{20} c_i^4 \right) + \nu^3 \left( -\frac{729}{16} - \frac{243}{80} c_i^2 + \frac{729}{40} c_i^4 \right) \right] \nonumber \\
&+ \cos 8\psi \left[ -\frac{1024}{315} s_i^6 (1+c_i^2)(1-7\nu+14\nu^2-7\nu^3) \right] \nonumber \\
&+ \Delta s_i \sin \psi \left[ -\frac{2159}{40320} - \frac{19 \ln 2}{32} + \left( -\frac{95}{224} - \frac{5 \ln 2}{8} \right) c_i^2 + \left( \frac{181}{13440} + \frac{\ln 2}{96} \right) c_i^4 \right. \nonumber \\
&\qquad \left. + \nu \left( \frac{1369}{160} + \frac{19 \ln 2}{48} + \left( -\frac{41}{48} - \frac{3 \ln 2}{8} \right) c_i^2 + \left( -\frac{313}{480} - \frac{\ln 2}{48} \right) c_i^4 \right) \right] \nonumber \\
&+ \sin 2\psi \left[ -\frac{428 \pi}{105} (1+c_i^2) \right] \nonumber \\
&+ \Delta s_i \sin 3\psi \left[ \frac{205119}{8960} - \frac{1971}{64} \ln(3/2) + \left( \frac{1917}{224} - \frac{135}{8} \ln(3/2) \right) c_i^2 \right. \nonumber \\
&\qquad \left. + \left( -\frac{43983}{8960} + \frac{243}{64} \ln(3/2) \right) c_i^4 \right. \nonumber \\
&\qquad \left. + \nu \left( -\frac{54869}{960} + \frac{567}{32} \ln(3/2) + \left( -\frac{923}{80} - \frac{81}{8} \ln(3/2) \right) c_i^2 \right. \right. \nonumber \\
&\qquad \left. \left. + \left( \frac{41851}{2880} - \frac{243}{32} \ln(3/2) \right) c_i^4 \right) \right] \nonumber \\
&+ \Delta s_i^3 (1+c_i^2) \sin 5\psi \left[ -\frac{113125}{5376} + \frac{3125}{192} \ln(5/2) + \nu \left( \frac{17639}{320} - \frac{3125}{96} \ln(5/2) \right) \right]. \tag{322g}
\end{align}
$$

For the cross polarizations we obtain

$$
H_{0 \times} = -2c_i \sin 2\psi, \tag{323a}
$$

$$
H_{1/2 \times} = s_i c_i \Delta \left[ -\frac{3}{4} \sin \psi + \frac{9}{4} \sin 3\psi \right], \qquad (323\text{b})
$$

$$
\begin{align}
H_{1 \times} &= c_i \sin 2\psi \left[ \frac{17}{3} - \frac{4}{3} c_i^2 + \nu \left( -\frac{13}{3} + 4 c_i^2 \right) \right] \nonumber \\
&\quad + c_i s_i^2 \sin 4\psi \left[ -\frac{8}{3} (1 - 3\nu) \right], \tag{323c}
\end{align}
$$

$$
\begin{align*}
H_{3/2 \times} = {}& s_i c_i \Delta \sin \psi \left[ \frac{21}{32} - \frac{5}{96} c_i^2 + \nu \left( -\frac{23}{48} + \frac{5}{48} c_i^2 \right) \right] \\
& - 4\pi c_i \sin 2\psi
\end{align*}
$$
$$ + s_i c_i \Delta \sin 3\psi \left[ -\frac{603}{64} + \frac{135}{64} c_i^2 + \nu \left( \frac{171}{32} - \frac{135}{32} c_i^2 \right) \right] \\ + s_i c_i \Delta \sin 5\psi \left[ \frac{625}{192} (1 - 2\nu) s_i^2 \right], \qquad (323d) $$

$$ H_{\frac{1}{2}} = s_i c_i \Delta \cos \psi \left[ -\frac{9}{20} - \frac{3}{2} \ln 2 \right] \\ + s_i c_i \Delta \cos 3\psi \left[ \frac{189}{20} - \frac{27}{2} \ln(3/2) \right] \\ - s_i c_i \Delta \left[ \frac{3\pi}{4} \right] \sin \psi \\ + c_i \sin 2\psi \left[ \frac{17}{15} + \frac{113}{30} c_i^2 - \frac{1}{4} c_i^4 + \nu \left( \frac{143}{9} - \frac{245}{18} c_i^2 + \frac{5}{4} c_i^4 \right) + \nu^2 \left( -\frac{14}{3} + \frac{35}{6} c_i^2 - \frac{5}{4} c_i^4 \right) \right] \\ + s_i c_i \Delta \sin 3\psi \left[ \frac{27\pi}{4} \right] \\ + \frac{4}{15} c_i s_i^2 \sin 4\psi \left[ 55 - 12 c_i^2 - \frac{5}{3} \nu (119 - 36 c_i^2) + 5 \nu^2 (17 - 12 c_i^2) \right] \\ + c_i \sin 6\psi \left[ -\frac{81}{20} s_i^4 (1 - 5\nu + 5\nu^2) \right], \qquad (323e) $$

$$ H_{\frac{1}{2}} = \frac{6}{5} s_i^2 c_i \nu \\ + c_i \cos 2\psi \left[ 2 - \frac{22}{5} c_i^2 + \nu \left( -\frac{282}{5} + \frac{94}{5} c_i^2 \right) \right] \\ + c_i s_i^2 \cos 4\psi \left[ -\frac{112}{5} + \frac{64}{3} \ln 2 + \nu \left( \frac{1193}{15} - 64 \ln 2 \right) \right] \\ + s_i c_i \Delta \sin \psi \left[ -\frac{913}{7680} + \frac{1891}{11520} c_i^2 - \frac{7}{4608} c_i^4 \right. \\ \left. + \nu \left( \frac{1165}{384} - \frac{235}{576} c_i^2 + \frac{7}{1152} c_i^4 \right) + \nu^2 \left( -\frac{1301}{4608} + \frac{301}{2304} c_i^2 - \frac{7}{1536} c_i^4 \right) \right] \\ + \pi c_i \sin 2\psi \left[ \frac{34}{3} - \frac{8}{3} c_i^2 + \nu \left( -\frac{20}{3} + 8 c_i^2 \right) \right] \\ + s_i c_i \Delta \sin 3\psi \left[ \frac{12501}{2560} - \frac{12069}{1280} c_i^2 + \frac{1701}{2560} c_i^4 \right. \\ \left. + \nu \left( -\frac{19581}{640} + \frac{7821}{320} c_i^2 - \frac{1701}{640} c_i^4 \right) + \nu^2 \left( \frac{18903}{2560} - \frac{11403}{1280} c_i^2 + \frac{5103}{2560} c_i^4 \right) \right] \\ + s_i^2 c_i \sin 4\psi \left[ -\frac{32\pi}{3} (1 - 3\nu) \right] \\ + \Delta s_i c_i \sin 5\psi \left[ -\frac{101875}{4608} + \frac{6875}{256} c_i^2 - \frac{21875}{4608} c_i^4 \right. \\ \left. + \nu \left( \frac{66875}{1152} - \frac{44375}{576} c_i^2 + \frac{21875}{1152} c_i^4 \right) + \nu^2 \left( -\frac{100625}{4608} + \frac{83125}{2304} c_i^2 - \frac{21875}{1536} c_i^4 \right) \right] \\ + \Delta s_i^5 c_i \sin 7\psi \left[ \frac{117649}{23040} (1 - 4\nu + 3\nu^2) \right], \qquad (323f) $$
$$
\begin{align}
H_3^\times ={}& \Delta s_i c_i \cos \psi \left[ \frac{11617}{20160} + \frac{21}{16} \ln 2 + \left( -\frac{251}{2240} - \frac{5}{48} \ln 2 \right) c_i^2 \right. \nonumber \\
& \left. + \nu \left( -\frac{2419}{240} - \frac{5}{24} \ln 2 + \left( \frac{727}{240} + \frac{5}{24} \ln 2 \right) c_i^2 \right) \right] \nonumber \\
& + c_i \cos 2\psi \left[ \frac{856 \pi}{105} \right] \nonumber \\
& + \Delta s_i c_i \cos 3\psi \left[ -\frac{36801}{896} + \frac{1809}{32} \ln(3/2) + \left( \frac{65097}{4480} - \frac{405}{32} \ln(3/2) \right) c_i^2 \right. \nonumber \\
& \left. + \nu \left( \frac{28445}{288} - \frac{405}{16} \ln(3/2) + \left( -\frac{7137}{160} + \frac{405}{16} \ln(3/2) \right) c_i^2 \right) \right] \nonumber \\
& + \Delta s_i^3 c_i \cos 5\psi \left[ \frac{113125}{2688} - \frac{3125}{96} \ln(5/2) + \nu \left( -\frac{17639}{160} + \frac{3125}{48} \ln(5/2) \right) \right] \nonumber \\
& + \pi \Delta s_i c_i \sin \psi \left[ \frac{21}{32} - \frac{5}{96} c_i^2 + \nu \left( -\frac{5}{48} + \frac{5}{48} c_i^2 \right) \right] \nonumber \\
& + c_i \sin 2\psi \left[ -\frac{3620761}{44100} + \frac{1712}{105} \gamma_E - \frac{4\pi^2}{3} + \frac{856}{105} \ln(16x) - \frac{3413}{1260} c_i^2 + \frac{2909}{2520} c_i^4 - \frac{1}{45} c_i^6 \right. \nonumber \\
& \left. + \nu \left( \frac{743}{90} - \frac{41\pi^2}{48} + \frac{3391}{180} c_i^2 - \frac{2287}{360} c_i^4 + \frac{7}{45} c_i^6 \right) \right. \nonumber \\
& \left. + \nu^2 \left( \frac{7919}{270} - \frac{5426}{135} c_i^2 + \frac{382}{45} c_i^4 - \frac{14}{45} c_i^6 \right) + \nu^3 \left( -\frac{6457}{1620} + \frac{1109}{180} c_i^2 - \frac{281}{120} c_i^4 + \frac{7}{45} c_i^6 \right) \right] \nonumber \\
& + \pi \Delta s_i c_i \sin 3\psi \left[ -\frac{1809}{64} + \frac{405}{64} c_i^2 + \nu \left( \frac{405}{32} - \frac{405}{32} c_i^2 \right) \right] \nonumber \\
& + s_i^2 c_i \sin 4\psi \left[ -\frac{1781}{105} + \frac{1208}{63} c_i^2 - \frac{64}{45} c_i^4 \right. \nonumber \\
& \left. + \nu \left( \frac{5207}{45} - \frac{536}{5} c_i^2 + \frac{448}{45} c_i^4 \right) + \nu^2 \left( -\frac{24838}{135} + \frac{2224}{15} c_i^2 - \frac{896}{45} c_i^4 \right) \right. \nonumber \\
& \left. + \nu^3 \left( \frac{1703}{45} - \frac{1976}{45} c_i^2 + \frac{448}{45} c_i^4 \right) \right] \nonumber \\
& + \Delta \sin 5\psi \left[ \frac{3125 \pi}{192} s_i^3 c_i (1 - 2\nu) \right] \nonumber \\
& + s_i^4 c_i \sin 6\psi \left[ \frac{9153}{280} - \frac{243}{35} c_i^2 + \nu \left( -\frac{7371}{40} + \frac{243}{5} c_i^2 \right) \right. \nonumber \\
& \left. + \nu^2 \left( \frac{1296}{5} - \frac{486}{5} c_i^2 \right) + \nu^3 \left( -\frac{3159}{40} + \frac{243}{5} c_i^2 \right) \right] \nonumber \\
& + \sin 8\psi \left[ -\frac{2048}{315} s_i^6 c_i (1 - 7\nu + 14\nu^2 - 7\nu^3) \right]. \tag{323g}
\end{align}
$$

Notice the non-linear memory zero-frequency (DC) term present in the Newtonian plus polarization $_0H_+$; see Refs. [427, 11, 189] for the computation of this term. Notice also that there is another DC term in the 2.5PN cross polarization $_{5/2}H_\times$, first term in Eq. (323f).

The practical implementation of the theoretical templates in the data analysis of detectors fol-
lows from the standard matched filtering technique. The raw output of the detector $o(t)$ consists
of the superposition of the real gravitational wave signal $h_{\text{real}}(t)$ and of noise $n(t)$. The noise is
assumed to be a stationary Gaussian random variable, with zero expectation value, and with (sup-
posedly known) frequency-dependent power spectral density $S_n(\omega)$. The experimenters construct
the correlation between $o(t)$ and a filter $q(t)$, i.e.,

$$c(t) = \int_{-\infty}^{+\infty} dt' o(t')q(t+t'), \qquad (324)$$

and divide $c(t)$ by the square root of its variance, or correlation noise. The expectation value of this ratio defines the filtered signal-to-noise ratio (SNR). Looking for the useful signal $h_{\text{real}}(t)$ in the detector's output $o(t)$, the data analysts adopt for the filter

$$\tilde{q}(\omega) = \frac{\tilde{h}(\omega)}{S_n(\omega)}, \qquad (325)$$

where $\tilde{q}(\omega)$ and $\tilde{h}(\omega)$ are the Fourier transforms of $q(t)$ and of the theoretically computed template $h(t)$. By the matched filtering theorem, the filter (325) maximizes the SNR if $h(t) = h_{\text{real}}(t)$. The maximum SNR is then the best achievable with a linear filter. In practice, because of systematic errors in the theoretical modelling, the template $h(t)$ will not exactly match the real signal $h_{\text{real}}(t)$; however if the template is to constitute a realistic representation of nature the errors will be small. This is of course the motivation for computing high order post-Newtonian templates, in order to reduce as much as possible the systematic errors due to the unknown post-Newtonian remainder.

To conclude, the use of theoretical templates based on the preceding 3PN/3.5PN waveforms, and having their frequency evolution built in via the 3.5PN phase evolution (318) [recall also the “tail-distorted” phase variable (321)], should yield some accurate detection and measurement of the binary signals, whose inspiral phase takes place in the detector’s bandwidth [105, 106, 159, 156, 3, 18, 111]. Interestingly, it should also permit some new tests of general relativity, because we have the possibility of checking that the observed signals do obey each of the terms of the phasing formula (318) – particularly interesting are those terms associated with non-linear tails – exactly as they are predicted by Einstein’s theory [84, 85, 15, 14]. Indeed, we don’t know of any other physical systems for which it would be possible to perform such tests.

## 9.5 Spherical harmonic modes for numerical relativity

The spin-weighted spherical harmonic modes of the polarization waveforms have been defined in Eq. (71). They can be evaluated either from applying the angular integration formula (72), or alternatively from using the relations (73)–(74) giving the individual modes directly in terms of separate contributions of the radiative moments $U_L$ and $V_L$. The latter route is actually more interesting [272] if some of the radiative moments are known to higher PN order than others. In this case the comparison with the numerical calculation for these particular modes can be made with higher post-Newtonian accuracy.

A useful fact to remember is that for non-spinning binaries, the mode $h^{\ell m}$ is entirely given by the mass multipole moment $U_L$ when $\ell + m$ is even, and by the current one $V_L$ when $\ell + 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. (73) splits to (see Ref. [197] for a proof)

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

$$h^{\ell m} = \frac{G}{\sqrt{2}Rc^{\ell+3}} i V^{\ell m} \quad (\text{when } \ell + m \text{ is odd}). \qquad (326b)$$
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 tail-distorted phase introduced in Eq. (321), and such that the dominant mode with $(\ell, m) = (2, 2)$ conventionally starts with one at the Newtonian order. We thus pose

$$h^{\ell m} = \frac{2G m \nu x}{R c^2} \sqrt{\frac{16\pi}{5}} \mathcal{H}^{\ell m} e^{-im\psi}. \quad (327)$$

We now list all the known results for $\mathcal{H}^{\ell m}$. We assume $m \ge 0$; the modes having $m < 0$ are easily deduced using $\mathcal{H}^{\ell, -m} = (-)^{\ell} \overline{\mathcal{H}}^{\ell m}$. The dominant mode $\mathcal{H}^{22}$, which is primarily important for numerical relativity comparisons, is known at 3.5PN order and reads [74, 197]

$$ \begin{aligned} \mathcal{H}^{22} = & 1 + x \left(-\frac{107}{42} + \frac{55}{42}\nu\right) + 2\pi x^{3/2} + x^2 \left(-\frac{2173}{1512} - \frac{1069}{216}\nu + \frac{2047}{1512}\nu^2\right) \\ & + x^{5/2} \left(-\frac{107\pi}{21} - 24i\nu + \frac{34\pi}{21}\nu\right) + x^3 \left(\frac{27027409}{646800} - \frac{856}{105}\gamma_E + \frac{428\pi}{105}i + \frac{2\pi^2}{3}\right) \\ & + \left(-\frac{278185}{33264} + \frac{41\pi^2}{96}\right)\nu - \frac{20261}{2772}\nu^2 + \frac{114635}{99792}\nu^3 - \frac{428}{105}\ln(16x) \\ & + x^{7/2} \left(-\frac{2173\pi}{756} + \left(-\frac{2495\pi}{378} + \frac{14333}{162}i\right)\nu + \left(\frac{40\pi}{27} - \frac{4066}{945}i\right)\nu^2\right) + \mathcal{O}\left(\frac{1}{c^8}\right). \end{aligned} \quad (328a) $$

Similarly, we report the subdominant modes $\mathcal{H}^{33}$ and $\mathcal{H}^{31}$ also known at 3.5PN order [195]

$$ \begin{aligned} \mathcal{H}^{33} = & -\frac{3}{4}i\sqrt{\frac{15}{14}}\Delta \left[ x^{1/2} + x^{3/2} \left(-4 + 2\nu\right) + x^2 \left(3\pi + i\left[-\frac{21}{5} + 6\ln(3/2)\right]\right) \right. \\ & \left. + x^{5/2} \left(\frac{123}{110} - \frac{1838\nu}{165} + \frac{887\nu^2}{330}\right) + x^3 \left(-12\pi + \frac{9\pi\nu}{2}\right. \right. \\ & \left. \left. + i\left[\frac{84}{5} - 24\ln(3/2) + \nu\left(-\frac{48103}{1215} + 9\ln(3/2)\right)\right]\right) \right. \\ & \left. + 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 + 6i\pi\left[-\frac{41}{35} + 3\ln(3/2)\right] \right. \right. \\ & \left. \left. + \frac{\nu}{8}\left[-\frac{7055}{429} + \frac{41}{8}\pi^2\right] - \frac{318841}{17160}\nu^2 + \frac{8237}{2860}\nu^3 - \frac{39}{7}\ln(16x)\right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \end{aligned} \quad (328b) $$

$$ \begin{aligned} \mathcal{H}^{31} = & \frac{i\Delta}{12\sqrt{14}} \left[ x^{1/2} + x^{3/2} \left(-\frac{8}{3} - \frac{2\nu}{3}\right) + x^2 \left(\pi + i\left[-\frac{7}{5} - 2\ln 2\right]\right) \right. \\ & \left. + x^{5/2} \left(\frac{607}{198} - \frac{136\nu}{99} - \frac{247\nu^2}{198}\right) + x^3 \left(-\frac{8\pi}{3} - \frac{7\pi\nu}{6}\right. \right. \\ & \left. \left. + i\left[\frac{56}{15} + \frac{16\ln 2}{3} + \nu\left(-\frac{1}{15} + \frac{7\ln 2}{3}\right)\right]\right) \right. \\ & \left. + x^{7/2} \left(\frac{10753397}{1513512} - 2\ln 2 \left[\frac{212}{105} + \ln 2\right] - \frac{26}{21}\gamma_E + \frac{\pi^2}{6} - 2i\pi \left[\frac{41}{105} + \ln 2\right] \right. \right. \\ & \left. \left. + \frac{\nu}{8}\left(-\frac{1738843}{19305} + \frac{41}{8}\pi^2\right) + \frac{327059}{30888}\nu^2 - \frac{17525}{15444}\nu^3 - \frac{13}{21}\ln x\right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right). \end{aligned} \quad (328c) $$

The other modes are known with a precision consistent with 3PN order in the full waveform [74]:

$$ \mathcal{H}^{21} = \frac{i}{3} \Delta \left[ x^{1/2} + x^{3/2} \left( -\frac{17}{28} + \frac{5\nu}{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{509\nu}{126} + \frac{79\nu^2}{168} \right) + x^3 \left( -\frac{17\pi}{28} + \frac{3\pi\nu}{14} + i \left( \frac{17}{56} + \nu \left( -\frac{353}{28} - \frac{3\ln 2}{7} \right) + \frac{17\ln 2}{14} \right) \right) \right] + \mathcal{O} \left( \frac{1}{c^7} \right), \quad (329a) $$

$$ \mathcal{H}^{20} = -\frac{5}{14\sqrt{6}} + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (329b) $$

$$ \mathcal{H}^{32} = \frac{1}{3}\sqrt{\frac{5}{7}} \left[ x(1-3\nu) + x^2 \left( -\frac{193}{90} + \frac{145\nu}{18} - \frac{73\nu^2}{18} \right) + x^{5/2} \left( 2\pi - 6\pi\nu + i \left( -3 + \frac{66\nu}{5} \right) \right) \right. \\ \left. + x^3 \left( -\frac{1451}{3960} - \frac{17387\nu}{3960} + \frac{5557\nu^2}{220} - \frac{5341\nu^3}{1320} \right) \right] + \mathcal{O} \left( \frac{1}{c^7} \right), \quad (329c) $$

$$ \mathcal{H}^{30} = -\frac{2}{5}i\sqrt{\frac{6}{7}}x^{5/2}\nu + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (329d) $$

$$ \mathcal{H}^{44} = -\frac{8}{9}\sqrt{\frac{5}{7}} \left[ x(1-3\nu) + x^2 \left( -\frac{593}{110} + \frac{1273\nu}{66} - \frac{175\nu^2}{22} \right) \right. \\ \left. + x^{5/2} \left( 4\pi - 12\pi\nu + i \left( -\frac{42}{5} + \nu \left( \frac{1193}{40} - 24\ln 2 \right) + 8\ln 2 \right) \right) \right. \\ \left. + x^3 \left( \frac{1068671}{200200} - \frac{1088119\nu}{28600} + \frac{146879\nu^2}{2340} - \frac{226097\nu^3}{17160} \right) \right] + \mathcal{O} \left( \frac{1}{c^7} \right), \quad (329e) $$

$$ \mathcal{H}^{43} = -\frac{9i}{4\sqrt{70}} \Delta \left[ x^{3/2}(1-2\nu) + x^{5/2} \left( -\frac{39}{11} + \frac{1267\nu}{132} - \frac{131\nu^2}{33} \right) \right. \\ \left. + x^3 \left( 3\pi - 6\pi\nu + i \left( -\frac{32}{5} + \nu \left( \frac{16301}{810} - 12\ln(3/2) \right) + 6\ln(3/2) \right) \right) \right] + \mathcal{O} \left( \frac{1}{c^7} \right), \quad (329f) $$

$$ \mathcal{H}^{42} = \frac{1}{63}\sqrt{5} \left[ x(1-3\nu) + x^2 \left( -\frac{437}{110} + \frac{805\nu}{66} - \frac{19\nu^2}{22} \right) + x^{5/2} \left( 2\pi - 6\pi\nu \right. \right. \\ \left. \left. + i \left( -\frac{21}{5} + \frac{84\nu}{5} \right) \right) + x^3 \left( \frac{1038039}{200200} - \frac{606751\nu}{28600} + \frac{400453\nu^2}{25740} + \frac{25783\nu^3}{17160} \right) \right] + \mathcal{O} \left( \frac{1}{c^7} \right), \quad (329g) $$

$$ \mathcal{H}^{41} = \frac{i}{84\sqrt{10}} \Delta \left[ x^{3/2}(1-2\nu) + x^{5/2} \left( -\frac{101}{33} + \frac{337\nu}{44} - \frac{83\nu^2}{33} \right) \right. \\ \left. + x^3 \left( \pi - 2\pi\nu + i \left( -\frac{32}{15} - 2\ln 2 + \nu \left( \frac{1661}{30} + 4\ln 2 \right) \right) \right) \right] + \mathcal{O} \left( \frac{1}{c^7} \right), \quad (329h) $$

$$ \mathcal{H}^{40} = -\frac{1}{504\sqrt{2}} + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (329i) $$

$$ \mathcal{H}^{55} = \frac{625i}{96\sqrt{66}} \Delta \left[ x^{3/2}(1-2\nu) + x^{5/2} \left( -\frac{263}{39} + \frac{688\nu}{39} - \frac{256\nu^2}{39} \right) \right. \\ \left. + x^3 \left( 5\pi - 10\pi\nu + i \left( -\frac{181}{14} + \nu \left( \frac{105834}{3125} - 20\ln(5/2) \right) + 10\ln(5/2) \right) \right) \right] + \mathcal{O} \left( \frac{1}{c^7} \right), \quad (329j) $$

$$ \mathcal{H}^{54} = -\frac{32}{9\sqrt{165}} \left[ x^2 (1 - 5\nu + 5\nu^2) + x^3 \left( -\frac{4451}{910} + \frac{3619\nu}{130} - \frac{521\nu^2}{13} + \frac{339\nu^3}{26} \right) \right] + \mathcal{O} \left( \frac{1}{c^7} \right), \quad (329k) $$

$$ \mathcal{H}^{53} = -\frac{9}{32}i\sqrt{\frac{3}{110}} \Delta \left[ x^{3/2}(1-2\nu) + x^{5/2} \left( -\frac{69}{13} + \frac{464\nu}{39} - \frac{88\nu^2}{39} \right) \right] $$
$$ + x^3 \left( 3\pi - 6\pi\nu + i \left( -\frac{543}{70} + \nu \left( \frac{83702}{3645} - 12 \ln(3/2) \right) + 6 \ln(3/2) \right) \right) \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (329l) $$

$$ \mathcal{H}^{52} = \frac{2}{27\sqrt{55}} \left[ x^2 (1 - 5\nu + 5\nu^2) + x^3 \left( -\frac{3911}{910} + \frac{3079\nu}{130} - \frac{413\nu^2}{13} + \frac{231\nu^3}{26} \right) \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (329m) $$

$$ \mathcal{H}^{51} = \frac{i}{288\sqrt{385}} \Delta \left[ x^{3/2}(1-2\nu) + x^{5/2} \left( -\frac{179}{39} + \frac{352\nu}{39} - \frac{4\nu^2}{39} \right) + x^3 \left( \pi - 2\pi\nu + i \left( -\frac{181}{70} - 2\ln 2 + \nu \left( \frac{626}{5} + 4\ln 2 \right) \right) \right) \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (329n) $$

$$ \mathcal{H}^{50} = \mathcal{O}\left(\frac{1}{c^7}\right), \quad (329o) $$

$$ \mathcal{H}^{66} = \frac{54}{5\sqrt{143}} \left[ x^2 (1 - 5\nu + 5\nu^2) + x^3 \left( -\frac{113}{14} + \frac{91\nu}{2} - 64\nu^2 + \frac{39\nu^3}{2} \right) \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (329p) $$

$$ \mathcal{H}^{65} = \frac{3125i x^{5/2}}{504\sqrt{429}} \Delta \left[ 1 - 4\nu + 3\nu^2 \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (329q) $$

$$ \mathcal{H}^{64} = -\frac{128}{495} \sqrt{\frac{2}{39}} \left[ x^2 (1 - 5\nu + 5\nu^2) + x^3 \left( -\frac{93}{14} + \frac{71\nu}{2} - 44\nu^2 + \frac{19\nu^3}{2} \right) \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (329r) $$

$$ \mathcal{H}^{63} = -\frac{81i x^{5/2}}{616\sqrt{65}} \Delta \left[ 1 - 4\nu + 3\nu^2 \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (329s) $$

$$ \mathcal{H}^{62} = \frac{2}{297\sqrt{65}} \left[ x^2 (1 - 5\nu + 5\nu^2) + x^3 \left( -\frac{81}{14} + \frac{59\nu}{2} - 32\nu^2 + \frac{7\nu^3}{2} \right) \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (329t) $$

$$ \mathcal{H}^{61} = \frac{i x^{5/2}}{8316\sqrt{26}} \Delta \left[ 1 - 4\nu + 3\nu^2 \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (329u) $$

$$ \mathcal{H}^{60} = \mathcal{O}\left(\frac{1}{c^7}\right). \quad (329v) $$

Notice that the modes with $m=0$ are zero except for the DC (zero-frequency) non-linear memory contributions. We already know that this effect arises at Newtonian order [see Eq. (322a)], hence the non zero values of the modes $\mathcal{H}^{20}$ and $\mathcal{H}^{40}$. See Ref. [189] for the DC memory contributions in the higher modes having $m=0$.

With the 3PN approximation all the modes with $\ell \ge 7$ can be considered as merely Newtonian. We give here the general Newtonian leading order expressions of any mode with arbitrary $\ell$ and non-zero $m$ (see the derivation in [272]):

$$ \mathcal{H}^{\ell m} = \frac{(-)^{(\ell-m+2)/2}}{2^{\ell+1} \left(\frac{\ell+m}{2}\right)! \left(\frac{\ell-m}{2}\right)! (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}(\nu) (im)^{\ell} x^{\ell/2-1} + \mathcal{O}\left(\frac{1}{c^{\ell-2}}\right) \quad (\text{for } \ell+m \text{ even}), \quad (330a) $$

$$ \mathcal{H}^{\ell m} = \frac{(-)^{(\ell-m-1)/2}}{2^{\ell-1} \left(\frac{\ell+m-1}{2}\right)! \left(\frac{\ell-m-1}{2}\right)! (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}(\nu) (im)^{\ell} x^{(\ell-1)/2} + \mathcal{O}\left(\frac{1}{c^{\ell-2}}\right) \quad (\text{for } \ell+m \text{ odd}), \quad (330b) $$

in which we employ the function $\sigma_{\ell}(\nu) = X_2^{\ell-1} + (-)^{\ell} X_1^{\ell-1}$, also given by Eqs. (305).
# 10 Eccentric Compact Binaries

Inspiralling compact binaries are usually modelled as moving in quasi-circular orbits since gravitational radiation reaction circularizes the orbit towards the late stages of inspiral [340, 339], as we discussed in Section 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 Kozai oscillations [283, 300] 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 [419]. On the other hand the population of stellar mass binaries in globular clusters is expected to have a thermal distribution of eccentricities [32]. In a study of the growth of intermediate black holes [235] in globular clusters it was found that the binaries have eccentricities between 0.1 and 0.2 in the eLISA bandwidth. Though, supermassive black hole binaries are powerful gravitational wave sources for eLISA, it is not known if they would be in quasi-circular or quasi-eccentric orbits. If a Kozai mechanism is at work, these supermassive black hole binaries could be in highly eccentric orbits and merge within the Hubble time [40]. Sources of the kind discussed above provide the prime motivation for investigating higher post-Newtonian order modelling for quasi-eccentric binaries.

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

In Section 7.3 we have given the equations of motion of non-spinning compact binary systems in the frame of the center-of-mass for general orbits at the 3PN and even 3.5PN order. We shall now investigate (in this section and the next one) the explicit solution to those equations. In particular, let us discuss the general “doubly-periodic” structure of the post-Newtonian solution, closely following Refs. [142, 143, 149].

The 3PN equations of motion admit, when neglecting the radiation reaction terms at 2.5PN order, 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 $\mathbf{J}$, given in harmonic coordinates at 3PN order by Eqs. (4.8)–(4.9) of Ref. [79].

The motion takes place in the plane orthogonal to $\mathbf{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 $\mathbf{J}$ are functions of $r$, $\dot{r}^2$, $v^2$ and $\mathbf{x} \times \mathbf{v}$. We adopt polar coordinates $(r, \phi)$ in the orbital plane, and express $E$ and the norm $J = |\mathbf{J}|$, thanks to $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, J], \qquad (331a) $$

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

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 [312]. 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 [143] 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="70" /> 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 (332)$$

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{S[r]}{\sqrt{\mathcal{R}[r]}}, \qquad (333)$$

which is such that the precession of the periastron per period is given by $\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 (334a)$$

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

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{\mathcal{R}[s]}}. \qquad (335)$$

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 dl S[r(l)], \qquad (336)$$

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 (337)$$

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 (336) by writing it in the form

$$\phi = \phi_p + \Omega (t - t_p) + W(\ell) \\ = \phi_p + K \ell + W(\ell). \qquad (338)$$

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. (336) 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 (339)$$

<footnote n="70">⁷⁰ The dependence on E and J will no longer be indicated but is always understood as implicit in what follows.</footnote>
Finally, the decomposition (338) 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$).

## 10.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 three types of eccentricities. It was introduced by Damour & Deruelle [149, 150] to study the problem of binary pulsar timing data including relativistic corrections at the 1PN order, where the relativistic periastron precession complicates the simpler Keplerian solution.

In the QK representation the radial motion is given in standard parametric form as

$$r = a_r (1 - e_r \cos u), \qquad (340)$$

where $u$ is the eccentric anomaly, 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 takes the usual form

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

where the mean anomaly is proportional to the time elapsed since the instant $t_p$ of passage at the periastron, $\ell = n(t - t_p)$ where $n = 2\pi/P$ is the mean motion and $P$ is the orbital period; see Eqs. (334). 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 + \mathcal{O}\left(\frac{1}{c^4}\right), \qquad (342)$$

where the true anomaly $v$ is defined by$^{71}$

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

The constant $K$ is the advance of periastron per orbital revolution defined by Eq. (333); it may be written as $K = \frac{\Phi}{2\pi}$ where $\Phi$ is the angle of return to the periastron.

Crucial to the 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

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

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

<footnote>$^{71}$ Comparing with Eqs. (338) we have also</footnote>
$$j = -\frac{2E h^2}{\mu^3}, \qquad (344b)$$

where $\mu = m\nu$ 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{J}{Gm}$. The equations to follow will then appear as expansions in powers of the small post-Newtonian parameter $\epsilon = O(1/c^2)$,<fnref n="72" /> with coefficients depending on $j$ and the dimensionless reduced mass ratio $\nu$; notice that the parameter $j$ is at Newtonian order, $j = O(1/c^0)$. We have [149]

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

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

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

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

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

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

The dependence of such relations on the coordinate system in use 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 will start only at 2PN order, see e.g., Eq. (347a).

The above QK representation of the compact binary motion at 1PN order has been generalized at the 2PN order in Refs. [170, 379, 420], and at the 3PN order by Memmesheimer, Gopakumar & Schäfer [312]. The construction of a generalized QK representation at 3PN order exploits the fact that the radial equation given by Eq. (331a) is a polynomial in $1/r$ (of seventh degree at 3PN order). 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 the 3PN order will obstruct the construction of the QK parametrization. Therefore Ref. [312] obtained it in the ADM coordinate system and also in the modified harmonic (MH) coordinates, obtained by applying the gauge transformation given in Eq. (204) on the SH coordinates. The equations of motion in the center-of-mass frame in MH coordinates have been given in Eqs. (222); see also Ref. [9] for details about the transformation between SH and MH coordinates.

At the 3PN order the radial equation in ADM or MH coordinates is still given by Eq. (340). However, the Kepler equation (341) and angular equation (342) acquire extra contributions and now become

$$\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), \qquad (346a)$$

$$\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), \qquad (346b)$$

in which the true anomaly $v$ is still given by Eq. (343). The new orbital elements $f_t, f_\phi, g_t, g_\phi, i_t, i_\phi, h_t$ and $h_\phi$ parametrize the 2PN and 3PN relativistic corrections.<fnref n="73" /> All the orbital elements

<footnote n="72">⁷² Note that this post-Newtonian parameter $\epsilon$ is precisely specified by Eq. (344a), while we only intended to define $\epsilon$ in Eq. (1) as representing a post-Newtonian estimate.</footnote>

<footnote n="73">⁷³ More precisely, $f_t, f_\phi, g_t, g_\phi$ are composed of 2PN and 3PN terms, but $i_t, i_\phi, h_t, h_\phi$ start only at 3PN order.</footnote>
are now to be related, similarly to Eqs. (345), 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. (340) and (346). 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 [170]. Their explicit expressions at 3PN order read

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

$$ K = 1 + \frac{3\epsilon}{j} + \frac{\epsilon^2}{4} \left[ \frac{3}{j}(-5 + 2\nu) + \frac{15}{j^2}(7 - 2\nu) \right] \\ + \frac{\epsilon^3}{128} \left[ \frac{24}{j}(5 - 5\nu + 4\nu^2) + \frac{1}{j^2}(-10080 + (13952 - 123\pi^2)\nu - 1440\nu^2) \right. \\ \left. + \frac{5}{j^3}(7392 + (-8000 + 123\pi^2)\nu + 336\nu^2) \right] + \mathcal{O}\left(\frac{1}{c^8}\right). \quad (347b) $$

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. (337). Moreover we can pose

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

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} \frac{x}{\epsilon} = & 1 + \epsilon \left[ -\frac{5}{4} + \frac{1}{12}\nu + \frac{2}{j} \right] \\ & + \epsilon^2 \left[ \frac{5}{2} + \frac{5}{24}\nu + \frac{1}{18}\nu^2 + \frac{1}{j^{1/2}}(-5 + 2\nu) + \frac{1}{j}(-5 + \frac{7}{6}\nu) + \frac{1}{j^2}\left(\frac{33}{2} - 5\nu\right) \right] \\ & + \epsilon^3 \left[ -\frac{235}{48} - \frac{25}{24}\nu - \frac{25}{576}\nu^2 + \frac{35}{1296}\nu^3 + \frac{1}{j}\left(\frac{35}{4} - \frac{5}{3}\nu + \frac{25}{36}\nu^2\right) \right. \\ & \left. + \frac{1}{j^{1/2}}\left(\frac{145}{8} - \frac{235}{24}\nu + \frac{29}{12}\nu^2\right) + \frac{1}{j^{3/2}}\left(-45 + \left(\frac{472}{9} - \frac{41}{96}\pi^2\right)\nu - 5\nu^2\right) \right. \\ & \left. + \frac{1}{j^2}\left(-\frac{565}{8} + \left(\frac{1903}{24} - \frac{41}{64}\pi^2\right)\nu - \frac{95}{12}\nu^2\right) \right. \\ & \left. + \frac{1}{j^3}\left(\frac{529}{3} + \left(-\frac{610}{3} + \frac{205}{64}\pi^2\right)\nu + \frac{35}{4}\nu^2\right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right). \end{aligned} \quad (349) $$
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. (340) and (346) are not gauge invariant; their expressions in terms of $\varepsilon$ and $j$ depend on the coordinate system in use. We refer to Refs. [312, 9] 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:

$$e_t^2 = 1 - j + \frac{\varepsilon}{4} [-8 + 8\nu + j(17 - 7\nu)] + \frac{\varepsilon^2}{8} [12 + 72\nu + 20\nu^2 + j(-112 + 47\nu - 16\nu^2) + 24j^{1/2}(5 - 2\nu) + \frac{16}{j}(4 - 7\nu) + \frac{24}{j^{1/2}}(-5 + 2\nu)] + \frac{\varepsilon^3}{6720} [23520 - 464800\nu + 179760\nu^2 + 16800\nu^3 + 525j(528 - 200\nu + 77\nu^2 - 24\nu^3) + 2520j^{1/2}(-265 + 193\nu - 46\nu^2) + \frac{6}{j}(-73920 + (260272 - 4305\pi^2)\nu - 61040\nu^2) + \frac{70}{j^{1/2}}(16380 + (-19964 + 123\pi^2)\nu + 3240\nu^2) + \frac{70}{j^{3/2}}(-10080 + (13952 - 123\pi^2)\nu - 1440\nu^2) + \frac{8}{j^2}(53760 + (-176024 + 4305\pi^2)\nu + 15120\nu^2)] + \mathcal{O}\left(\frac{1}{c^8}\right). \quad (350)$$

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

In the case of a circular orbit, the angular momentum variable, say $j_{\text{circ}}$, is related to the constant of energy $\varepsilon$ by the 3PN gauge-invariant expansion

$$j_{\text{circ}} = 1 + \left(\frac{9}{4} + \frac{\nu}{4}\right) \varepsilon + \left(\frac{81}{16} - 2\nu + \frac{\nu^2}{16}\right) \varepsilon^2 + \left(\frac{945}{64} + \left[-\frac{7699}{192} + \frac{41}{32}\pi^2\right]\nu + \frac{\nu^2}{2} + \frac{\nu^3}{64}\right) \varepsilon^3 + \mathcal{O}\left(\frac{1}{c^8}\right).$$

This permits to reduce various quantities to circular orbits, for instance, the periastron advance is found to be well defined in the limiting case of a circular orbit, and is given at 3PN order in terms of the PN parameter (230) [or (348)] by

$$K_{\text{circ}} = 1 + 3x + \left(\frac{27}{2} - 7\nu\right) x^2 + \left(\frac{135}{2} + \left[-\frac{649}{4} + \frac{123}{32}\pi^2\right] \nu + 7\nu^2\right) x^3 + \mathcal{O}\left(\frac{1}{c^8}\right).$$

See Ref. [291] for a comparison between the PN prediction for the periastron advance of circular orbits and numerical calculations based on self-force theory in the small mass ratio limit.

## 10.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 & Mathews [340, 339] at Newtonian level. The 1PN and 1.5PN corrections to the fluxes were provided in Refs. [416, 86, 267, 87, 366] 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 [149]. These results were extended to 2PN order in Refs. [224, 225] for the instantaneous terms (leaving aside the tails) using the
generalized QK representation [170, 379, 420]; the energy flux and waveform were in agreement with those of Ref. [424] obtained using a different method. Arun et al. [10, 9, 12] have fully generalized the results at 3PN order, including all tails and related hereditary contributions, by computing the averaged energy and angular momentum fluxes for quasi-elliptical orbits using the QK representation at 3PN order [312], and deriving the secular evolution of the orbital elements under 3PN gravitational radiation reaction.<fnref n="74" />

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 terms rapidly oscillating at the orbital period. An analytic approach, based on an improved method of variation of constants, has been discussed in Ref. [153] for dealing with this issue at the leading 2.5PN radiation reaction order.

The generalized QK representation of the motion discussed in Section 10.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="75" /> Actually the averaging procedure applies to the “instantaneous” parts of the fluxes, while the “hereditary” parts are treated separately for technical reasons [10, 9, 12]. Following the decomposition (308) 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. For the angular momentum flux one needs also to include a contribution from the memory effect [12]. 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 (351) $$

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 (348), and the “time” eccentricity $e_t$ given by (350). 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 (346a) 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 Section 7.3, 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 associated with Hadamard’s self-field regularization [see Eq. (221)]; 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 6 where the finite part $\mathcal{FP}$ contains the regularization factor (42). 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_{t}) + \beta_{l}(x, e_{t}) \sin u + \gamma_{l}(x, e_{t}) \ln(1 - e_{t} \cos u)}{(1 - e_{t} \cos u)^{l+1}}. \qquad (352) $$

The factor $du/d\ell$ has been inserted to prepare for the orbital average (351). 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_t$. 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

<footnote n="74">⁷⁴ On the other hand, for the computation of the gravitational waveform of eccentric binary orbits up to the 2PN order in the Fourier domain, see Refs. [401, 402].</footnote>

<footnote n="75">⁷⁵ Recall that the fluxes are defined in a general way, for any matter system, in terms of the radiative multipole moments by the expressions (68).</footnote>
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}$)$^{76}$

$$ \int_{0}^{2\pi} \frac{du}{2\pi} \frac{\sin u}{(1 - e_t \cos u)^{l+1}} = 0, \qquad (353a) $$

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

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

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

$$ Z(z, e_t) = \ln \left[ \frac{\sqrt{1 - e_t^2} + 1}{2} \right] + 2 \ln \left[ 1 + \frac{\sqrt{1 - e_t^2} - 1}{z + \sqrt{z^2 - e_t^2}} \right]. \qquad (354) $$

From Eq. (353a) 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 in the form [9]

$$ \langle \mathcal{F}_{\text{inst}} \rangle = \frac{32c^5}{5G} \nu^2 x^5 \left( \mathcal{I}_0 + x \mathcal{I}_1 + x^2 \mathcal{I}_2 + x^3 \mathcal{I}_3 \right), \qquad (355) $$

where we recall that the post-Newtonian parameter $x$ is defined by (348). The various instantaneous post-Newtonian pieces depend on the symmetric mass ratio $\nu$ and the time eccentricity $e_t$ in MH coordinates as

$$ \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 (356a) $$

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

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

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

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

The Newtonian coefficient $\mathcal{I}_0$ is nothing but the Peters & Mathews [340] enhancement function of eccentricity that enters in the orbital gravitational radiation decay of the binary pulsar; see Eq. (11). For ease of presentation we did not add a label on $e_t$ to indicate that it is the time eccentricity in MH coordinates; such MH-coordinates $e_t$ is given by Eq. (350). Recall that on the contrary $x$ is gauge invariant, so no such label is required on it.

The last term in the 3PN coefficient is proportional to some logarithm which directly arises from the integration formula (353c). Inside the logarithm we have posed

$$ x_0 \equiv \frac{Gm}{c^2 r_0}, \qquad (357) $$

exhibiting an explicit dependence upon the arbitrary length scale $r_0$; we recall that $r_0$ was introduced in the formalism through Eq. (42). Only after adding the hereditary contribution to the 3PN energy flux can we check the required cancellation of the constant $x_0$. The hereditary part is made of tails and tails-of-tails, and is of the form

$$ \langle \mathcal{F}_{\text{hered}} \rangle = \frac{32c^5}{5G} \nu^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 (358) $$

where the post-Newtonian pieces, only at the 1.5PN, 2.5PN and 3PN orders, read [10]

$$ \mathcal{K}_{3/2} = 4\pi \varphi(e_t), \qquad (359a) $$
$$ K_{5/2} = -\frac{8191}{672} \pi \psi(e_t) - \frac{583}{24} \nu \pi \zeta(e_t), \qquad (359b) $$

$$ 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}}{x_0}\right) \right] F(e_t), \qquad (359c) $$

where $\varphi(e_t)$, $\psi(e_t)$, $\zeta(e_t)$, $\kappa(e_t)$ and $F(e_t)$ are certain “enhancement” functions of the eccentricity.

Among them the four functions $\varphi(e_t)$, $\psi(e_t)$, $\zeta(e_t)$ and $\kappa(e_t)$ appearing in Eqs. (359) do not admit analytic closed-form expressions. They have been obtained in Refs. [10] (extending Ref. [87]) in the form of infinite series made out of quadratic products of Bessel functions. Numerical plots of these four enhancement factors as functions of eccentricity $e_t$ have been provided in Ref. [10]; we give in Figure 3 the graph of the function $\varphi(e_t)$ which enters the dominant 1.5PN tail term in Eq. (358).

<figure><img src="image_7.png" /><figcaption>**Figure 3:** Variation of the enhancement factor $\varphi(e)$ with the eccentricity $e$. This function agrees with the numerical calculation of Ref. [87] modulo a trivial rescaling with the Peters-Mathews function (356a). 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>

Furthermore their leading correction term $e_t^2$ in the limit of small eccentricity $e_t \ll 1$ can be obtained analytically as [10]

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

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

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

$$ \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 (360d) $$
On the other hand the function $F(e_t)$ in factor of the logarithm in the 3PN piece does admit some closed analytic form:

$$F(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 (361)$$

The latter analytical result is very important for checking that the arbitrary constant $x_0$ disappears from the final result. Indeed we immediately verify from comparing the last term in Eq. (356d) with Eq. (359c) that $x_0$ cancels out from the sum of the instantaneous and hereditary contributions in the 3PN energy flux. This fact was already observed for the circular orbit case in Ref. [81]; see also the discussions around Eqs. (93)–(94) and at the end of Section 4.2.

Finally we can check that the correct circular orbit limit, which is given by Eq. (314), is recovered from the sum $\langle \mathcal{F}_{\text{inst}} \rangle + \langle \mathcal{F}_{\text{hered}} \rangle$. The next correction of order $e_t^2$ when $e_t \to 0$ can be deduced from Eqs. (360)–(361) in analytic form; having the flux in analytic form may be useful for studying the gravitational waves from binary black hole systems with moderately high eccentricities, such as those formed in globular clusters [235].

Previously the averaged energy flux was represented using $x$ – the gauge invariant variable (348) – and the time eccentricity $e_t$ which however is gauge dependent. Of course it 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. (344). To this end it suffices to replace $e_t$ by its MH-coordinate expression (350) and to use Eq. (349) 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 $(x, \iota)$, where $x$ is given by (348) and we pose

$$\iota \equiv \frac{3x}{K - 1}. \qquad (362)$$

Such choice would be motivated by the fact that $\iota$ 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, \iota)$ 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}_i$. 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 (see Ref. [12]): The angular momentum flux is split into instantaneous $\mathcal{G}_{\text{inst}}$ and hereditary $\mathcal{G}_{\text{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_t$; finally these are obtained numerically as well as analytically to next-to-leading order $e_t^2$. 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 [12]. This extends earlier analyses at previous orders: Newtonian [339] as we have reviewed in Section 1.2; 1PN [86, 267]; 1.5PN [87, 366] and 2PN [224, 153]. Let us take the example of the mean motion $n$. From Eq. (347a) together with the definitions (344) 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 (363) $$

The usual balance equations for energy and angular momentum

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

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

have already been used at Newtonian order in Eqs. (9). Although heuristically assumed at 3PN order, they have been proved through 1.5PN order in Section 5.4. 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 (365) $$

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 [153, 279].

# 11 Spinning Compact Binaries

The post-Newtonian templates have been developed so far for compact binary systems which can be described with great precision by point masses without spins. Here by spin, we mean the intrinsic (classical) angular momentum S of the individual compact body. However, including the effects of spins is essential, as the astrophysical evidence indicates that stellar-mass black holes [2, 390, 311, 227, 323] and supermassive black holes [188, 101, 102] (see Ref. [364] for a review) can be generically close to maximally spinning. The presence of spins crucially affects the dynamics of the binary, in particular leading to orbital plane precession if they are not aligned with the orbital angular momentum (see for instance [138, 8]), and thereby to strong modulations in the observed signal frequency and phase.

In recent years an important effort has been undertaken to compute spin effects to high post-Newtonian order in the dynamics and gravitational 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 dynamics of particles with spins (either Lagrangian or Hamiltonian);

2. *Radiation.* The mass and current radiative multipole moments, including tails and all hereditary effects, are to be computed. 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 Part A.

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 = G m_a^2 \chi_a, \qquad (366) $$
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. (366)], 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 correct) 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 shall make a distinction between spin-orbit (SO) effects, which are linear in the spins, and spin-spin (SS) ones, which are quadratic. In this article we shall especially review the SO effects as they play the most important role in gravitational wave detection and parameter estimation. As we shall see a good deal is known on spin effects (both SO and SS), but still it will be important in the future to further improve our knowledge of the waveform and gravitational-wave phasing, by computing still higher post-Newtonian SO and SS terms, and to include at least the dominant spin-spin-spin (SSS) effect [305]. For the computations of SSS and even SSSS effects see Refs. [246, 245, 296, 305, 413].

The SO effects have been known at the leading level since the seminal works of Tulczyjew [411, 412], Barker & O'Connell [27, 28] and Kidder et al. [275, 271]. 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 in Refs. [394, 194, 165, 292, 352, 241] for the equations of motion or dynamics, and in Refs. [53, 54] for the gravitational radiation field. Note that Refs. [394, 194, 165, 241] employ traditional post-Newtonian methods (both harmonic-coordinates and Hamiltonian), but that Refs. [292, 352] 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 in Refs. [242, 244] using the Hamiltonian method for the equations of motion, in Ref. [297] using the EFT, and in Refs. [307, 90] using the harmonic-coordinates method. Here we shall focus on the harmonic-coordinates approach [307, 90, 89, 306] which is in fact best formulated using a Lagrangian, see Section 11.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) [89, 306]. An analytic solution for the SO precession effects will be presented in Section 11.2. Note that concerning the radiation field the highest known SO level actually contains specific tail-induced contributions at 3PN [54] and 4PN [306] orders, see Section 11.3.

The SS effects are known at the leading level corresponding to 2PN order from Barker & O'Connell [27, 28] in the equations of motion (see [271, 351, 110] for subsequent derivations), and from Refs. [275, 271] in the radiation field. Next-to-leading SS contributions are at 3PN order and have been obtained with Hamiltonian [387, 389, 388, 247, 241], EFT [354, 356, 355, 293, 299] and harmonic-coordinates [88] techniques (with [88] obtaining also the next-to-leading SS terms in the gravitational-wave flux). 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 [347]. They have been computed through 2PN order in the fluxes and orbital phase in Refs. [217, 218, 314]. The interaction terms $S_1 \times S_2$ can be computed using a simple pole-dipole formalism like the one we shall review in Section 11.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 in Refs. [294, 298] (EFT) and [243] (Hamiltonian). In this article we shall generally neglect the SS effects and refer for these to the literature quoted above.
## 11.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 [239, 19] and has also been developed in the context of the EFT approach [351]. Variants and alternatives (most importantly the associated Hamiltonian formalism) can be found in Refs. [389, 386, 25]. 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 [411, 412, 310, 331, 135, 409, 179].

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="77" /> 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 (367) $$

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 (368) $$

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 Refs. [239, 351] and the general spirit of effective field theories, 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) \rightarrow \Lambda_A^B e_B^\alpha(\tau)$ where $\Lambda_A^B$ is a constant Lorentz matrix;

3. It is reparametrization-invariant, i.e., its form is independent of the parameter $\tau$ used to follow the particle's worldline.

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

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

These assumptions confine the formalism to a “pole-dipole” model and to terms 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.

<footnote n="77">⁷⁷ 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^\alpha$, defined by $e_A^\beta e_A^\alpha = \delta^\beta_\alpha$, satisfies $\eta_{AB} e_A^\alpha e_B^\beta = g_{\alpha\beta}$. We have also the completeness relation $e_A^\beta e_B^\alpha = \delta^\beta_\alpha$.</footnote>
As it is written in (369), 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., Ref. [19])

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

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

$$p_{\alpha} \equiv \left. \frac{\partial L}{\partial u^{\alpha}} \right|_{\omega,g}, \qquad (371a)$$

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

Note that the right-hand side of Eq. (370) is necessarily symmetric by exchange of the indices $\alpha$ and $\beta$. Finally, imposing the invariance of the action (369) 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 (372)$$

where the functions $p_\alpha(u, \omega)$ and $S_{\alpha\beta}(u, \omega)$ must be reparametrization 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 Eq. (379) below.

We now investigate the unconstrained variations of the action (369) 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 (373)$$

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 (374)$$

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

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

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 (376)$$

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. (376) 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-Papapetrou [310, 331, 135] equation of motion

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

which involves the famous coupling of the spin tensor to the Riemann curvature.<fnref n="78" /> With a little more work, the equation of motion (377) can also be derived using an arbitrary coordinate system, making use of the parallel transport equation (376). 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. (370). We obtain the standard pole-dipole result [411, 412, 310, 331, 135, 409, 179]:

$$ 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 (378) $$

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 (377) and the equation of spin precession (375).

Up to now we have considered unconstrained variations of the action (369), 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 body with respect to which the spin angular momentum is defined (see Ref. [271] for a discussion). Here we adopt the Tulczyjew covariant SSC [411, 412]

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

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

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

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

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

where we denote $(pu) \equiv p_\mu u^\mu$. By further contracting Eq. (381) with $u^\alpha$ we obtain an explicit expression for $(pu)$, which can then be substituted back into (381) to provide the relation linking the four-momentum $p_\alpha$ to the four-velocity $u_\alpha$. It can be checked using (379) and (381) 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$.

<footnote n="78">⁷⁸ Our conventions for the Riemann tensor $R_{\alpha\beta\mu\nu}$ follow those of MTW [319].</footnote>

<footnote n="79">⁷⁹ The four-dimensional Levi-Civita tensor is defined by $\varepsilon_{\alpha\beta\mu\nu} \equiv \sqrt{-g} \varepsilon_{\alpha\beta\mu\nu}$ and $\varepsilon^{\alpha\beta\mu\nu} \equiv -\varepsilon^{\alpha\beta\mu\nu}/\sqrt{-g}$; here $\varepsilon_{\alpha\beta\mu\nu} = \varepsilon^{\alpha\beta\mu\nu}$ denotes the completely anti-symmetric Levi-Civita symbol such that $\varepsilon_{0123} = \varepsilon^{0123} = 1$. For convenience in this section we pose $c = 1$.</footnote>
Henceforth we shall restrict our attention to spin-orbit (SO) interactions, which are linear in the spins. 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 + \mathcal{O}(S^2)$. Hence, from Eq. (375) we deduce that $DS_{\alpha\beta}/d\tau = \mathcal{O}(S^2)$. The equation for the spin covariant vector $S_\alpha$ then reduces at linear order to

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

Thus the spin covector is parallel transported along the particle's trajectory at linear order in spin. We can also impose that the spin should be purely spatial for the comoving observer:

$$ S_\alpha u^\alpha = 0. \qquad (383) $$

From now on, we shall often omit writing the $\mathcal{O}(S^2)$ remainders.

In applications (e.g., the construction of gravitational wave templates for the compact 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 [90]. Indeed the evolution equations of such spin variables reduces, by construction, to ordinary precession equations, and these variables are secularly constant (see Ref. [423]).

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 (384) $$

If we choose for the time-like tetrad vector the four-velocity itself, $e_0^\alpha = u^\alpha$, $^{80}$ the time component tetrad projection $S_0$ vanishes because of the orthogonality condition (383). We have seen that $S_\alpha S^\alpha = s^2$ is conserved along the trajectory; because of (383) 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 (385) $$

which is indeed the relation defining a Euclidean conserved-norm spin variable $S_a$. $^{81}$ 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 [165, 90] 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 Section 11.3.

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

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

<footnote>$^{80}$ 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 (369).</footnote>

<footnote>$^{81}$ Beware that here we employ the usual slight ambiguity in the notation when using the same carrier letter $S$ to denote the tetrad components (384) 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>
where the precession tensor $\Omega_{ab}$ is related to the tetrad components $\omega_{AB}$ of the rotation tensor defined in (368) by $\Omega_{ab} = z\omega_{ab}$ where we pose $z \equiv d\tau/dt$, remembering the redshift variable (276). The antisymmetric character of the matrix $\Omega_{ab}$ guarantees that $S_a$ satisfies the Euclidean precession equation

$$ \frac{dS}{dt} = \Omega \times S, \qquad (387) $$

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

## 11.2 Equations of motion and precession for spin-orbit effects

The previous formalism can be generalized to self-gravitating systems consisting of two (or more generally $N$) spinning point particles. The metric generated by the system of particles, interacting only through gravitation, is solution of the Einstein field equations (18) with stress-energy tensor given by the sum of the individual stress-energy tensors (378) for each particles. The equations of motion of the particles are given by the Mathisson-Papapetrou equations (377) with “self-gravitating” metric evaluated at the location of the particles thanks to a regularization procedure (see Section 6). The precession equations of each of the spins are given by

$$ \frac{dS_a}{dt} = \Omega_a \times S_a, \qquad (388) $$

where $a = 1, 2$ labels the particles. The spin variables $S_a$ are the conserved-norm spins defined in Section 11.1. In the following it is convenient to introduce two combinations of the individuals spins defined by (with $X_a \equiv m_a/m$ and $\nu \equiv X_1X_2$)$^{82}$

$$ S \equiv S_1 + S_2, \qquad (389a) $$

$$ \Sigma \equiv \frac{S_2}{X_2} - \frac{S_1}{X_1}. \qquad (389b) $$

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.$^{83}$ We introduce an orthonormal moving triad $\{\mathbf{n}, \lambda, \ell\}$ defined by the unit separation vector $\mathbf{n} = \mathbf{x}/r$ (with $r = |\mathbf{x}|$) and the unit normal $\ell$ to the instantaneous orbital plane given by $\ell = \mathbf{n} \times \mathbf{v}/|\mathbf{n} \times \mathbf{v}|$; the orthonormal triad is then completed by $\lambda = \ell \times \mathbf{n}$. Those vectors are represented on Figure 4, which shows the geometry of the system. The orbital frequency $\Omega$ is defined for general orbits, not necessarily circular, by $\mathbf{v} = \dot{\mathbf{r}}\mathbf{n} + r\Omega\lambda$ where $\dot{\mathbf{r}} = \mathbf{n} \cdot \mathbf{v}$ represents

$$ S_1 = X_1 S - \nu \Sigma, $$

$$ S_2 = X_2 S + \nu \Sigma. $$

<footnote>$^{82}$ Notation adopted in Ref. [271]; the inverse formulas read</footnote>

<footnote>$^{83}$ 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 of motion: $G = 0$) are related to the relative position and velocity $\mathbf{x}$ and $\mathbf{v}$ by some expressions similar to Eqs. (224) but containing spin effects starting at order 1.5PN.</footnote>
the derivative of $r$ with respect to the coordinate time $t$. The general expression for the relative acceleration $\mathbf{a} \equiv d\mathbf{v}/dt$ decomposed in the moving basis $\{\mathbf{n}, \mathbf{\lambda}, \mathbf{\ell}\}$ is

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

Here we have introduced the orbital plane precession $\varpi$ of the orbit defined by $\varpi \equiv -\mathbf{\lambda} \cdot d\mathbf{\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. (227). Then $\mathbf{\lambda}$ represents the direction of the velocity, and the precession frequency $\varpi$ is proportional to the variation of $\mathbf{\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\varpi\Omega\mathbf{\ell} + \mathcal{O}\left(\frac{1}{c^5}\right). \qquad (391) $$

Since we neglect the radiation reaction damping there is no component of the acceleration along $\mathbf{\lambda}$. This equation represents the generalization of Eq. (226) 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 (228), while the precessional frequency $\varpi$ will entirely be due to spins.

Here we report the latest results for the spin-orbit (SO) contributions into these quantities at the next-to-next-to-leading level corresponding to 3.5PN order [307, 90]. 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

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

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

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

where this time the ordering post-Newtonian parameter is $x \equiv (Gm\Omega/c^3)^{2/3}$. The SO terms at the same level in the conserved energy associated with the equations of motion will be given in Eq. (415) below. In order to complete the evolution equations for quasi-circular orbits we need also the precession vectors $\mathbf{\Omega}_a$ of the two spins as defined by Eq. (388). These are given by

$$ \mathbf{\Omega}_1 = \frac{c^3 x^{5/2}}{Gm} \mathbf{\ell} \left\{ \begin{aligned} & \frac{3}{4} + \frac{1}{2}\nu - \frac{3}{4}\Delta + x \left[ \frac{9}{16} + \frac{5}{4}\nu - \frac{1}{24}\nu^2 + \Delta \left( -\frac{9}{16} + \frac{5}{8}\nu \right) \right] \\ & + x^2 \left[ \frac{27}{32} + \frac{3}{16}\nu - \frac{105}{32}\nu^2 - \frac{1}{48}\nu^3 + \Delta \left( -\frac{27}{32} + \frac{39}{8}\nu - \frac{5}{32}\nu^2 \right) \right] \end{aligned} \right\} + \mathcal{O}\left(\frac{1}{c^8}\right). \qquad (394) $$

We obtain $\mathbf{\Omega}_2$ from $\mathbf{\Omega}_1$ simply by exchanging the masses, $\Delta \to -\Delta$. At the linear SO level the precession vectors $\mathbf{\Omega}_a$ are independent of the spins.<fnref n="84" />

<footnote n="84">⁸⁴ 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; $\mathbf{\Omega}_a$ is the precession frequency of the a-th spin while $\varpi$ is the precession frequency of the orbital plane; and $\omega_a$ defined earlier in Eqs. (244) and (284) 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 Section 8.3; see Eq. (282) and the corotation condition (285).</footnote>
We now investigate an analytical solution for the dynamics of compact spinning binaries on quasi-circular orbits, including the effects of spin precession [54, 306]. This solution will be valid whenever the radiation reaction effects can be neglected, and is restricted to the linear SO level.

<figure><img src="image_8.png" /><figcaption>**Figure 4:** Geometric definitions for the precessional motion of spinning compact binaries [54, 306]. 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$; (v) 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>

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{dJ}{dt} = 0. \qquad (395) $$

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="85" />

$$ J = L + \frac{S}{c}. \qquad (396) $$

This split between $L$ and $S = S_1 + S_2$ is specified by our choice of spin variables, here the conserved-norm spins defined in Section 11.1. Note that although $L$ is called the “orbital” angular momentum, it actually includes both non-spin and spin contributions. We refer to Eq. (4.7) in [90] for the expression of $L$ at the next-to-next-to-leading SO level for quasi-circular orbits.

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 (397a) $$

<footnote n="85">⁸⁵ Recall from Eq. (366) that in our convention the spins have the dimension of an angular momentum times c.</footnote>
$$ \frac{d\boldsymbol{\lambda}}{dt} = -\Omega \boldsymbol{n} + \varpi \boldsymbol{\ell}, \qquad (397b) $$

$$ \frac{d\boldsymbol{\ell}}{dt} = -\varpi \boldsymbol{\lambda}. \qquad (397c) $$

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

$$ \frac{d\boldsymbol{n}}{dt} = \boldsymbol{\Omega} \times \boldsymbol{n}, \qquad (398a) $$

$$ \frac{d\boldsymbol{\lambda}}{dt} = (\boldsymbol{\Omega} + \varpi) \times \boldsymbol{\lambda}, \qquad (398b) $$

$$ \frac{d\boldsymbol{\ell}}{dt} = \varpi \times \boldsymbol{\ell}. \qquad (398c) $$

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

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

where $J \equiv |\mathbf{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 $\mathbf{n}$ with respect to some reference direction in the orbital plane denoted $\mathbf{x}_\ell$. Following Ref. [13], we pose

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

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

such that $\{\mathbf{x}_\ell, \mathbf{y}_\ell, \boldsymbol{\ell}\}$ forms an orthonormal basis. The motion takes place in the instantaneous orbital plane spanned by $\mathbf{n}$ and $\boldsymbol{\lambda}$, and the phase angle $\Phi$ is such that (see Figure 4):

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

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

from which we deduce

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

Combining Eqs. (402) with (399) we also get

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

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

$$ \frac{d\alpha}{dt} = \varpi \frac{\sin \Phi}{\sin \iota}, \qquad (404a) $$
$$ \frac{d\iota}{dt} = \varpi \cos \Phi, \qquad (404b) $$

$$ \frac{d\Phi}{dt} = \Omega - \varpi \frac{\sin \Phi}{\tan \iota}. \qquad (404c) $$

On the other hand, using the decomposition of the total angular momentum (396) together with the fact that the components of $L$ projected along $n$ and $\lambda$ are of the order $O(S)$, see e.g., Eq. (4.7) in Ref. [90], we deduce that $\sin \iota$ is itself a small quantity of order $O(S)$. Since we also have $\varpi = O(S)$, we conclude by direct integration of the sum of Eqs. (404a) and (404c) that

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

in which we have defined the “carrier” phase as

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

with $\phi_0$ the value of the carrier phase at some arbitrary initial time $t_0$. An important point we have used when integrating (406) is that the orbital frequency $\Omega$ is constant at linear order in the spins. Indeed, from Eq. (392) 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. (409c) 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. (405), 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. (399) and (403):

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

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

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. (407) together with the decomposition (396) and the explicit expression of $L$ in Ref. [90], we shall then be able to obtain the explicit time variation of the precession angles $(\alpha, \iota)$ and phase $\Phi$. Combining (388) and (397) we obtain

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

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

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

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

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

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

$$ \frac{dS_\ell^a}{dt} = \mathcal{O}(S^2). \qquad (409c) $$
We see that, as stated before, the spin components along $\ell$ are constant, and so is the orbital frequency $\Omega$ given by (392). At the linear SO level, the equations (409) can be decoupled and integrated as

$$S_n^a = S_{\perp}^a \cos \psi_a, \qquad (410a)$$

$$S_{\lambda}^{a} = -S_{\perp}^{a} \sin \psi_{a}, \qquad (410b)$$

$$S_{\ell}^{a} = S_{\parallel}^{a}. \qquad (410c)$$

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 (411)$$

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 triad $\{n(t), \lambda(t), \ell(t)\}$ in terms of some reference triad $\{n_0, \lambda_0, \ell_0\}$ at the reference time $t_0$ in Eqs. (411) and (406). The best way to express the result is to introduce the complex null vector $m \equiv \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 \iota e^{i\alpha} - \sin \iota_0 e^{i\alpha_0}) e^{-i\phi} \ell_0 + O(S^2), \qquad (412a)$$

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

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 \iota e^{i\alpha}$ and its complex conjugate $\sin \iota e^{-i\alpha}$, which follows from Eqs. (407) 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 (412) for the moving triad start at order $O(1/c^3)$.

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

Like in Section 9 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. Eqs. (66)–(68). 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; see e.g., Eqs. (76).$^{86}$

The source moments have been expressed in Eqs. (123) 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. (145). Now, computing spin effects, the matter tensor $T^{\alpha\beta}$ has been found to be given by (378) 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 spin contributions (necessarily SO) to the general mass and current source multipole moments ($\forall \ell \in \mathbb{N}$):

$$I_L^{\text{SO}} = \frac{2\ell \nu}{c^3(\ell + 1)} \left\{ \ell \left[ \sigma_{\ell}(\nu)(\mathbf{v} \times \mathbf{S})^{(i\ell)} - \sigma_{\ell+1}(\nu)(\mathbf{v} \times \mathbf{\Sigma})^{(i\ell)} \right] x^{L-1} \right\}$$

<footnote>$^{86}$ In this section we can neglect the gauge multipole moments $W_L, \dots, Z_L$.</footnote>
$$ -(\ell - 1) \left[ \sigma_{\ell}(\nu)(\mathbf{x} \times \mathbf{S})^{\langle i_{\ell} \rangle} - \sigma_{\ell+1}(\nu)(\mathbf{x} \times \mathbf{\Sigma})^{\langle i_{\ell} \rangle} \right] v^{i_{\ell-1}} x^{L-2} \right\} + \mathcal{O}\left(\frac{1}{c^5}\right), \quad (413a) $$

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

Paralleling the similar expressions (304) for the Newtonian approximation to the source moments in the non-spin case, we posed $\sigma_{\ell}(\nu) \equiv X_2^{\ell-1} + (-)^{\ell}X_1^{\ell-1}$ with $X_a = m_a/m$ [see also Eqs. (305)]. In Eqs. (413) we employ the notation (389) 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 [53, 89] for higher order post-Newtonian expressions of the source moments. If we insert the expressions (413) into tail integrals like (76), we find that some spin contributions originate from tails starting at the 3PN order [54].

Finally, skipping details, we are left with the following highest-order known result for the SO contributions to the gravitational wave energy flux, which is currently 4PN order [53, 89, 306]:<fnref n="87" />

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

We recall that $S_{\ell} \equiv \ell \cdot S$ and $\Sigma_{\ell} \equiv \ell \cdot \Sigma$, with $S$ and $\Sigma$ denoting the combinations (389), and the individual spins are the specific conserved-norm spins that have been introduced in Section 11.1. The result (414) superposes to the non-spin contributions given by Eq. (314). Satisfyingly it is in complete agreement in the test-mass limit where $\nu \to 0$ with the result of black-hole perturbation theory on a Kerr background obtained in Ref. [396].

Finally we can compute the spin effects in the time evolution of the binary's orbital frequency $\Omega$. We rely as in Section 9 on the equation (295) balancing the total emitted energy flux $\mathcal{F}$ with the variation of the binary's center-of-mass energy $E$. The non-spin contributions in $E$ have been provided for quasi-circular binaries in Eq. (232); the SO contributions to next-to-next-to-leading order are given by [307, 90]

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

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}}{dE/d\Omega} \right)_{\text{SO}}. \quad (416) $$

<footnote n="87">⁸⁷ Notice that the spin-orbit contributions due to the absorption by the black-hole horizons have to be added to this post-Newtonian result [349, 392, 5, 125].</footnote>
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. (415). Fortunately, this is the case of the conserved-norm spin variables, as proved in Ref. [423] 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 [54, 89]. In the center-of-mass frame, the only vectors at our disposal, except for the spins, are $\mathbf{n}$ and $\mathbf{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 $\mathbf{\Omega}_a$ in Eqs. (388): They must be pseudovectors, and, at linear order in spin, they must only depend on $\mathbf{n}$ and $\mathbf{v}$; so that they must be proportional to $\ell$, as can be explicitly seen for instance in Eq. (394). Now, the time derivative of the components along $\ell$ of the spins are given by $dS_\ell^a/dt = \mathbf{S}_a \cdot (d\ell/dt + \ell \times \mathbf{\Omega}_a)$. The second term vanishes because $\mathbf{\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. (409c). 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 Ref. [218] and references therein for related discussions).

<table><caption>**Table 4:** Spin-orbit contributions to the number of gravitational-wave cycles $\mathcal{N}_{\text{cycle}}$ [defined by Eq. (319)] for binaries detectable by ground-based detectors LIGO-VIRGO. The entry frequency is $f_{\text{seismic}} = 10$ Hz and the terminal frequency is $f_{\text{ISCO}} = \frac{c^3}{6^{3/2}\pi G m}$. For each compact object the magnitude $\chi_a$ and the orientation $\kappa_a$ of the spin are defined by $\mathbf{S}_a = G m_a^2 \chi_a \hat{\mathbf{S}}_a$ and $\kappa_a = \hat{\mathbf{S}}_a \cdot \ell$; remind Eq. (366). The spin-spin (SS) terms are neglected.</caption><thead><tr><th>PN order</th><th></th><th>1.4 M<sub>☉</sub> + 1.4 M<sub>☉</sub></th><th>10 M<sub>☉</sub> + 1.4 M<sub>☉</sub></th><th>10 M<sub>☉</sub> + 10 M<sub>☉</sub></th></tr></thead><tbody><tr><td>1.5PN</td><td>(leading SO)</td><td>65.6κ<sub>1</sub>χ<sub>1</sub> + 65.6κ<sub>2</sub>χ<sub>2</sub></td><td>114.0κ<sub>1</sub>χ<sub>1</sub> + 11.7κ<sub>2</sub>χ<sub>2</sub></td><td>16.0κ<sub>1</sub>χ<sub>1</sub> + 16.0κ<sub>2</sub>χ<sub>2</sub></td></tr><tr><td>2.5PN</td><td>(1PN SO)</td><td>9.3κ<sub>1</sub>χ<sub>1</sub> + 9.3κ<sub>2</sub>χ<sub>2</sub></td><td>33.8κ<sub>1</sub>χ<sub>1</sub> + 2.9κ<sub>2</sub>χ<sub>2</sub></td><td>5.7κ<sub>1</sub>χ<sub>1</sub> + 5.7κ<sub>2</sub>χ<sub>2</sub></td></tr><tr><td>3PN</td><td>(leading SO-tail)</td><td>-3.2κ<sub>1</sub>χ<sub>1</sub> - 3.2κ<sub>2</sub>χ<sub>2</sub></td><td>-13.2κ<sub>1</sub>χ<sub>1</sub> - 1.3κ<sub>2</sub>χ<sub>2</sub></td><td>-2.6κ<sub>1</sub>χ<sub>1</sub> - 2.6κ<sub>2</sub>χ<sub>2</sub></td></tr><tr><td>3.5PN</td><td>(2PN SO)</td><td>1.9κ<sub>1</sub>χ<sub>1</sub> + 1.9κ<sub>2</sub>χ<sub>2</sub></td><td>11.1κ<sub>1</sub>χ<sub>1</sub> + 0.8κ<sub>2</sub>χ<sub>2</sub></td><td>1.7κ<sub>1</sub>χ<sub>1</sub> + 1.7κ<sub>2</sub>χ<sub>2</sub></td></tr><tr><td>4PN</td><td>(1PN SO-tail)</td><td>-1.5κ<sub>1</sub>χ<sub>1</sub> - 1.5κ<sub>2</sub>χ<sub>2</sub></td><td>-8.0κ<sub>1</sub>χ<sub>1</sub> - 0.7κ<sub>2</sub>χ<sub>2</sub></td><td>-1.5κ<sub>1</sub>χ<sub>1</sub> - 1.5κ<sub>2</sub>χ<sub>2</sub></td></tr></tbody></table>

We are then allowed to apply Eq. (416) with conserved-norm spin variables at the SO level. We thus obtain the secular evolution of $\mathbf{\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 \nu} \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}\nu \right) S_{\ell} + \left( -\frac{41745}{448} + \frac{15}{8}\nu \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[ \begin{aligned} & \left( -\frac{8980424995}{6096384} + \frac{6586595}{6048}\nu - \frac{305}{288}\nu^2 \right) S_{\ell} \\ & + \left( -\frac{170978035}{387072} + \frac{2876425}{5376}\nu + \frac{4735}{1152}\nu^2 \right) \Delta\Sigma_{\ell} \end{aligned} \right] \end{aligned} \right\} \quad (417) $$
$$+x^{5/2} \left[ \left( \frac{2388425\pi}{3024} - \frac{9925\pi}{36}\nu \right) S_{\ell} + \left( \frac{3237995\pi}{12096} - \frac{258245\pi}{2016}\nu \right) \Delta\Sigma_{\ell} \right] + \mathcal{O}\left(\frac{1}{c^6}\right).$$

This expression, when added to the expression for the non-spin terms given by Eq. (318), 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 Figure 4, has been proved to be given by $\Phi = \phi - \alpha$ at linear order in spins, cf. Eq. (405). The precessional correction $\alpha$ can be computed at linear order in spins from the results of Section 11.2.

As we have done in Table 3 for the non-spin terms, we evaluate in Table 4 the SO contributions to the number of gravitational-wave cycles accumulated in the bandwidth of LIGO-VIRGO detectors, Eq. (319). The results show that the SO terms up to 4PN order can be numerically important, for spins close to maximal and for suitable orientations. They can even be larger than the corresponding non-spin contributions at 3.5PN order and perhaps at 4PN order (but recall that the non-spin terms at 4PN order are not known); compare with Table 3. We thus conclude that the SO terms are relevant to be included in the gravitational wave templates for an accurate extraction of the binary parameters. Although numerically smaller, the SS terms are also relevant; for these we refer to the literature quoted at the beginning of Section 11.

## Acknowledgments

The author's works reviewed in this article have been done in a series of collaborations with K.G. Arun, Alejandro Bohé, Alessandra Buonanno, Thibault Damour, Steven Detweiler, Gilles Esposito-Farèse, Guillaume Faye, Bala R. Iyer, Alexandre Le Tiec, Sylvain Marsat, Samaya Nissanke, Olivier Poujade, B.S. Sathyaprakash, Gerhard Schäfer, Bernard F. Whiting, Clifford M. Will, and Alan G. Wiseman.