<figure><img src="image_1.png" /></figure>

# Post-Newtonian theory for gravitational waves

Luc Blanchet¹

Received: 5 February 2024 / Accepted: 7 May 2024

© The Author(s) 2024

## Abstract

To be observed and analyzed by the network of current gravitational-wave detectors (LIGO, Virgo, KAGRA), and in anticipation of future third generation ground-based (Einstein Telescope, Cosmic Explorer) and space-borne (LISA) detectors, inspi-ralling compact binaries—binary star systems composed of neutron stars and/or black holes in their late stage of evolution prior the final coalescence—require high-accuracy predictions from general relativity. The orbital dynamics and emitted gravitational waves of these very relativistic systems can be accurately modelled using state-of-the-art post-Newtonian theory. In this article we review the multipolar-post-Minkowskian approximation scheme, merged to the standard post-Newtonian expansion into a single formalism valid for general isolated matter system. This cocktail of approximation methods (called MPM-PN) has been successfully applied to compact binary systems, producing equations of motion up to the fourth-post-Newtonian (4PN) level, and gravitational waveform and flux to 4.5PN order beyond the Einstein quadrupole formula. We describe the dimensional regularization at work in such high post-Newtonian calculations, for curing both ultra-violet and infra-red divergences. Several landmark results are detailed: the definition of multipole moments, the gravitational radiation reaction, the conservative dynamics of circular orbits, the first law of compact binary mechanics, and the non-linear effects in the gravitational-wave propagation (tails, iterated tails and non-linear memory). We also discuss the case of compact binaries moving on eccentric orbits, and the effects of spins (both spin-orbit and spin-spin) on the equations of motion and gravitational-wave energy flux and waveform.

This article is a revised version of https://doi.org/10.12942/lrr-2014-2

**Change summary** Major revision, updated and expanded.

**Change details** The main improvements, rewriting, and addition of new subsections have been done in Sects. 1, 2.4, 3.2, and 3.4. The number of total references has increased from 428 to 604.

Extended author information available on the last page of the article
**Keywords** Gravitational radiation · Post-Newtonian approximation · Post-Minkowskian approximation · Multipolar expansion · Dimensional regularization · Equations of motion · Inspiralling compact binary

Contents

<table><tbody><tr><td>1</td><td>Physical introduction</td><td>3</td></tr><tr><td>1.1</td><td>Analytic approximations and wave-generation formalism</td><td>5</td></tr><tr><td>1.1.1</td><td>Cocktail of approximation methods</td><td>5</td></tr><tr><td>1.1.2</td><td>The quadrupole moment formalism</td><td>10</td></tr><tr><td>1.1.3</td><td>Influence of the internal structure of compact bodies</td><td>16</td></tr><tr><td>1.1.4</td><td>The recoil by gravitational waves</td><td>21</td></tr><tr><td>1.2</td><td>Overview on motion and radiation</td><td>24</td></tr><tr><td>1.2.1</td><td>Post-Newtonian equations of motion</td><td>24</td></tr><tr><td>1.2.2</td><td>Post-Newtonian gravitational radiation</td><td>30</td></tr><tr><td>2</td><td>Post-Newtonian isolated sources</td><td>33</td></tr><tr><td>2.1</td><td>Non-linear iteration of the field equations outside the source</td><td>33</td></tr><tr><td>2.1.1</td><td>Einstein's field equations</td><td>33</td></tr><tr><td>2.1.2</td><td>Linearized vacuum equations</td><td>37</td></tr><tr><td>2.1.3</td><td>The multipolar post-Minkowskian solution in harmonic coordinates</td><td>41</td></tr><tr><td>2.1.4</td><td>Source multipole moments versus canonical moments</td><td>45</td></tr><tr><td>2.1.5</td><td>The multipolar post-Minkowskian solution in radiative coordinates</td><td>48</td></tr><tr><td>2.2</td><td>Matching to a post-Newtonian source</td><td>53</td></tr><tr><td>2.2.1</td><td>The matching equation</td><td>54</td></tr><tr><td>2.2.2</td><td>General expression of the multipole expansion</td><td>56</td></tr><tr><td>2.2.3</td><td>The source multipole moments</td><td>58</td></tr><tr><td>2.3</td><td>Interior field of a post-Newtonian source</td><td>63</td></tr><tr><td>2.3.1</td><td>Post-Newtonian iteration in the near zone</td><td>64</td></tr><tr><td>2.3.2</td><td>Post-Newtonian metric and radiation reaction effects</td><td>67</td></tr><tr><td>2.3.3</td><td>The 3.5PN metric for general matter systems</td><td>70</td></tr><tr><td>2.3.4</td><td>Radiation reaction potentials</td><td>75</td></tr><tr><td>2.4</td><td>Asymptotic gravitational waveform</td><td>79</td></tr><tr><td>2.4.1</td><td>The radiative multipole moments</td><td>79</td></tr><tr><td>2.4.2</td><td>Gravitational-wave tails and tails of tails</td><td>83</td></tr><tr><td>2.4.3</td><td>Formulas to integrate non-linearities</td><td>87</td></tr><tr><td>2.4.4</td><td>The non-linear memory effect</td><td>90</td></tr><tr><td>2.4.5</td><td>Radiative versus canonical moments</td><td>95</td></tr><tr><td>2.4.6</td><td>Canonical versus source moments</td><td>103</td></tr><tr><td>3</td><td>Compact binary systems</td><td>105</td></tr><tr><td>3.1</td><td>Regularization of the field of point particles</td><td>105</td></tr><tr><td>3.1.1</td><td>Hadamard's regularization</td><td>106</td></tr><tr><td>3.1.2</td><td>Dimensional regularization of the equations of motion</td><td>109</td></tr><tr><td>3.1.3</td><td>Dimensional regularization of the radiation field</td><td>115</td></tr><tr><td>3.1.4</td><td>Infra-red divergences and radiative modes</td><td>120</td></tr><tr><td>3.2</td><td>Newtonian-like equations of motion</td><td>125</td></tr><tr><td>3.2.1</td><td>The Fokker action for point particles</td><td>125</td></tr><tr><td>3.2.2</td><td>The 3.5PN acceleration and 3PN energy</td><td>129</td></tr><tr><td>3.2.3</td><td>The 4PN Lagrangian and Hamiltonian</td><td>137</td></tr><tr><td>3.2.4</td><td>Equations of motion in the frame of the center of mass</td><td>153</td></tr><tr><td>3.2.5</td><td>Equations of motion and energy for quasi-circular orbits</td><td>168</td></tr></tbody></table>
<table><tr><td>3.2.6</td><td>The 2.5PN metric in the near zone</td><td>173</td></tr><tr><td>3.3</td><td>Conservative dynamics of compact binaries</td><td>181</td></tr><tr><td>3.3.1</td><td>The innermost circular orbit</td><td>181</td></tr><tr><td>3.3.2</td><td>Dynamical stability of circular orbits</td><td>185</td></tr><tr><td>3.3.3</td><td>The first law of binary point-particle mechanics</td><td>190</td></tr><tr><td>3.3.4</td><td>Post-Newtonian approximation versus gravitational self-force</td><td>196</td></tr><tr><td>3.4</td><td>Gravitational waves from compact binaries</td><td>201</td></tr><tr><td>3.4.1</td><td>The binary's multipole moments for circular orbits</td><td>202</td></tr><tr><td>3.4.2</td><td>Post-adiabatic evaluation of hereditary integrals</td><td>209</td></tr><tr><td>3.4.3</td><td>The 4.5PN energy flux for circular orbits</td><td>214</td></tr><tr><td>3.4.4</td><td>The 4.5PN phase and frequency evolution</td><td>216</td></tr><tr><td>3.4.5</td><td>Spherical modes for data analysis and numerical relativity</td><td>220</td></tr><tr><td>3.5</td><td>Eccentric compact binaries</td><td>229</td></tr><tr><td>3.5.1</td><td>Doubly periodic structure of the motion of eccentric binaries</td><td>229</td></tr><tr><td>3.5.2</td><td>Quasi-Keplerian representation of the motion</td><td>232</td></tr><tr><td>3.5.3</td><td>Averaged energy and angular momentum fluxes</td><td>236</td></tr><tr><td>3.6</td><td>Spinning compact binaries</td><td>244</td></tr><tr><td>3.6.1</td><td>Lagrangian formalism for spinning point particles</td><td>246</td></tr><tr><td>3.6.2</td><td>Equations of motion and precession for circular binaries</td><td>252</td></tr><tr><td>3.6.3</td><td>Orbital precession at the spin-orbit level</td><td>255</td></tr><tr><td>3.6.4</td><td>Spin-orbit effects in the gravitational-wave flux and orbital phase</td><td>259</td></tr><tr><td>3.6.5</td><td>Spin-spin effects</td><td>262</td></tr><tr><td colspan="2">References</td><td>265</td></tr></table>

# 1 Physical introduction

The theory of gravitational radiation from isolated sources, is a fascinating science that can be explored by means of the analytical (i.e., mathematical) resolution of partial differential equations. Indeed, the Einstein (1915) field equations of general relativity (GR), when use is made of the harmonic-coordinate conditions, take the form of a quasi-linear hyperbolic differential system of equations, involving the D'Alembert (1743) wave operator. The resolution of that system of equations constitutes a *problème bien posé* in the sense of Hadamard (1932), Bruhat (1962), and which is amenable to an analytic solution using approximation methods.

Nowadays, the importance of the field lies in the exciting comparison of the theory with contemporary astrophysical observations, of binary pulsars like the historical Hulse and Taylor (1975) pulsar PSR 1913+16, and of gravitational waves produced by massive and rapidly evolving systems of neutron stars and black holes. These are now routinely analyzed on Earth by the world-wide network of large-scale laser interferometers, comprising the interferometers LIGO, Virgo and KAGRA (with the smaller detector GEO in support). The first direct detection of the coalescence of two black holes has been achieved on September 14, 2015 by the Advanced LIGO detectors (Abbott et al. 2016c). Since 2017 joint observations were made by both LIGO and Virgo (with thus a good sky localization of the events), including the historical event of the merger of two neutron stars detected on August 17, 2017 (Abbott et al. 2017b), which constituted the start of multi messenger astronomy with gravitational waves, with the almost simultaneous detection of a γ-ray burst by the satellites Fermi and Integral (Abbott et al. 2017a). Further ahead, the space-based
laser interferometer detector LISA—an outstanding mission of the European Space Agency to be launched in 2034—should be able to detect the collision and merger of supermassive black-hole binaries at cosmological distances, and learning about the formation of large scale structure, stellar evolution and the early Universe.

To prepare these experiments, the required theoretical work consists of carrying out a sufficiently general solution of the Einstein field equations, valid for a large class of matter systems, and describing the physical processes of the emission and propagation of the gravitational waves from the source to the distant detector, as well as their back-reaction onto the source. The solution should then be applied to specific sources like inspiralling compact binaries—i.e. binary sources losing energy by gravitational waves prior to merger.

For general sources it is hopeless to solve the problem via a rigorous deduction within the exact theory of general relativity, and we have to resort to approximation methods. Of course the ultimate aim of approximation methods is to extract from the theory some firm predictions to be compared with the outcome of experiments. However, we have to keep in mind that such methods are often missing a clear theoretical framework and are sometimes not related in a very precise mathematical way to the first principles of the theory (Futamase and Schutz 1983; Rendall 1992).

The flagship of approximation methods is the *post-Newtonian* approximation, which has been developed from the early days of general relativity by Droste (1917), Lorentz and Droste (1937). This approximation is at the origin of many of the great successes of general relativity, and it gives wonderful answers to the problems of motion and gravitational radiation of systems of compact objects. Three crucial applications are:

1. The motion of point-like objects at the first post-Newtonian approximation level (Einstein et al. 1938), is routinely taken into account to describe the solar system dynamics (motion of the centers of mass of planets);

2. The gravitational radiation-reaction force, which appears in the equations of motion at the second-and-a-half post-Newtonian (2.5PN) order (Damour and Deruelle 1981a, b; Damour 1983a, b), has been experimentally verified by the observation of the secular acceleration of the orbital motion of the Hulse–Taylor binary pulsar (Taylor et al. 1979; Taylor and Weisberg 1982; Taylor 1993) and for instance the double pulsar (Kramer and Wex 2009);

3. The analysis of gravitational waves emitted by inspiralling compact binaries—two neutron stars or black holes driven into coalescence by emission of gravitational radiation—necessitates the prior knowledge of the equations of motion and radiation field up to very high post-Newtonian order (Cutler et al. 1993a; Cutler and Flanagan 1994).

This article reviews the current status of the post-Newtonian approach to the problems of the motion of inspiralling compact binaries and their emission of gravitational waves. When the two compact objects approach each other toward merger, the post-Newtonian expansion will lose accuracy and should be taken over by numerical-relativity computations (Pretorius 2005; Campanelli et al. 2006; Baker
et al. 2006b). We refer to review articles (Grandclément and Novak 2009; Faber and Rasio 2012; Lehner and Pretorius 2014; Duez and Zlochower 2018) for numerical-relativity methods. Despite very intensive developments of numerical relativity, the post-Newtonian approximation remains indispensable for describing the inspiral phase of compact binaries to high accuracy, and for providing a powerful benchmark against which the numerical computations are tested.

Section 2 of the article deals with general post-Newtonian matter sources. The exterior field of the source is investigated by means of a combination of analytic multipolar and post-Minkowskian (MPM) approximations. The physical observables in the far-zone of the source are described by a specific set of radiative multipole moments. By matching the exterior solution to the metric of the post-Newtonian (PN) source in the near-zone the explicit expressions of the source multipole moments are obtained. The relationships between the radiative and source moments involve many non-linear multipole interactions, among them those associated with the tails (and tails-of-tails, etc.) of gravitational waves. This is the so-called MPM-PN formalism.

Section 3 is devoted to the application to compact binary systems, with particular emphasis on black hole binaries with spins. We present the equations of binary motion, and the associated Lagrangian and Hamiltonian, at the fourth post-Newtonian (4PN) order beyond the Newtonian acceleration. The gravitational-wave energy flux, taking consistently into account the relativistic corrections in the binary's moments as well as the various tail effects, is derived through 4PN order with respect to the quadrupole formalism. The binary's orbital phase, whose prior knowledge is crucial for searching and analyzing the signals from inspiralling compact binaries, is deduced from an energy balance argument at the 4.5PN order (in the simple case of circular orbits).

# 1.1 Analytic approximations and wave-generation formalism

## 1.1.1 Cocktail of approximation methods

The basic problem we face is to relate the asymptotic gravitational-wave form $h_{ij}$ generated by some isolated source (in a suitable asymptotic coordinate system), at the location of a detector in the wave zone of the source, to the material content of the source, i.e., its stress-energy tensor $T^{\alpha\beta}$, using approximation methods in general relativity.<fnref n="1" /> Therefore, a general wave-generation formalism must solve the field equations, and the non-linearity therein, by imposing some suitable approximation series in one or several small physical parameters. Some important approximations that we shall use in this article are the post-Newtonian method (or non-linear 1/c-expansion), the post-Minkowskian method or non-linear iteration (G-expansion), the multipole decomposition in irreducible representations of the rotation group (or equivalently a-expansion in the source radius), the far-zone expansion (1/R-expansion).

<footnote n="1">¹ In this article Greek indices $\alpha\beta\dots\mu\nu\dots$ take space-time values 0, 1, 2, 3 and Latin indices $ab\dots ij\dots$ spatial values 1, 2, 3. Cartesian coordinates are assumed throughout and boldface notation is often used for ordinary Euclidean vectors. Our signature is +2; hence the Minkowski metric reads $\eta_{\alpha\beta} = \text{diag}(-1, +1, +1, +1)$. As usual G and c are Newton's constant and the speed of light, consistently included unless otherwise specified.</footnote>
expansion in the distance to the source), and the perturbation in the small mass limit ($v$-expansion in the mass ratio of a binary system). In particular, the post-Newtonian expansion has provided us in the past with our best insights into the problems of motion and radiation. The most successful wave-generation formalisms make a *gourmet* cocktail of these approximation methods. For reviews on analytic approximations and applications to the motion and the gravitational-wave generation see Thorne (1983), Damour (1983a, 1987a, b), Thorne (1987), Will (1993), Blanchet (1997a, 2011b), Schäfer (2011), Maggiore (2008), Buonanno and Sathyaprakash (2015). For reviews on black-hole perturbations and the self-force approach see Poisson et al. (2011), Sasaki and Tagoshi (2003), Detweiler (2011), Barack (2011), Barack and Pound (2018).

The post-Newtonian approximation is valid under the assumptions of a weak gravitational field inside the source (we shall see later how to model neutron stars and black holes), and of slow internal motions.<fnref n="2" /> The main problem with this approximation, is its domain of validity, which is limited to the near zone of the source—the region surrounding the source that is of small extent with respect to the wavelength of the gravitational waves. A serious consequence is the *a priori* inability of the post-Newtonian expansion to incorporate the boundary conditions at infinity, which determine the radiation reaction force in the source's local equations of motion.

The post-Minkowskian expansion, by contrast, is uniformly valid, as soon as the source is weakly self-gravitating, over all space-time. In a sense, the post-Minkowskian method is more fundamental than the post-Newtonian one; it can be regarded as an “upstream” approximation with respect to the post-Newtonian expansion, because each coefficient of the post-Minkowskian series can in turn be re-expanded in a post-Newtonian fashion. Therefore, a way to take into account the boundary conditions at infinity in the post-Newtonian series is to control *first* the post-Minkowskian expansion. Notice that the post-Minkowskian method is also upstream (in the previous sense) with respect to the multipole expansion, when considered outside the source, and with respect to the far-zone expansion, when considered far from the source.

The most “downstream” approximation that we shall use in this article is the post-Newtonian one; therefore this is the approximation that dictates the allowed physical properties of our matter source. We assume mainly that the source is at once *slowly moving* and *weakly stressed*, and we abbreviate this by saying that the source is *post-Newtonian*. For post-Newtonian sources, the parameter defined from the components of the matter stress-energy tensor $T^{\alpha\beta}$ and the source's Newtonian potential $U$ by

$$ \epsilon \equiv \max \left\{ \left| \frac{T^{0i}}{T^{00}} \right|, \left| \frac{T^{ij}}{T^{00}} \right|^{1/2}, \left| \frac{U}{c^2} \right|^{1/2} \right\}, \quad (1) $$

is much less than one. This parameter represents essentially a slow motion estimate $\epsilon \sim v/c$, where $v$ denotes a typical internal velocity. By a slight abuse of notation,

<footnote n="2">² Establishing the post-Newtonian expansion rigorously has been the subject of numerous mathematically oriented works, see e.g., Rendall (1990, 1992, 1994).</footnote>
following, we shall henceforth write formally $\epsilon \equiv 1/c$, even though $\epsilon$ is dimensionless whereas $c$ has the dimension of a velocity. Thus, $1/c \ll 1$ in the case of post-Newtonian sources. The small post-Newtonian remainders will be denoted $O(1/c^n)$. Furthermore, following Chandrasekhar (1965), we shall refer to a small post-Newtonian term with formal order $O(1/c^n)$ relative to the Newtonian acceleration in the equations of motion, as being "$\frac{n}{2}$PN".

We have $|U/c^2|^{1/2} \ll 1/c$ for sources with negligible self-gravity, and whose dynamics are therefore driven by non-gravitational forces. However, we shall generally assume that the source is self-gravitating; in that case we see that it is necessarily weakly (but not negligibly) self-gravitating, i.e., $|U/c^2|^{1/2} = O(1/c)^3$. Note that the adjective “slow-motion” is a bit clumsy because we shall in fact consider very relativistic sources such as inspiralling compact binaries, for which $v/c$ can be as large as 50% in the last rotations, and whose description necessitates the control of high post-Newtonian approximations.

At the lowest-order in the Newtonian limit $1/c \to 0$, the gravitational waveform of a post-Newtonian matter source is generated by the time variations of the quadrupole moment of the source. We shall review in Sect. 1.1.2 the utterly important “Newtonian” quadrupole moment formalism of Einstein (1918), Landau and Lifshitz (1971). Taking into account higher post-Newtonian corrections in a wave generation formalism will mean including into the waveform the contributions of higher multipole moments, beyond the mass quadrupole moment. Post-Newtonian corrections of order $O(1/c^n)$ beyond the quadrupole formalism will still be denoted as $\frac{n}{2}$PN. But notice that the quadrupole formalism corresponds itself to a small radiation reaction effect of order $2.5$PN = $O(1/c^5)$ relative to the Newtonian acceleration. The lesson here is that building a post-Newtonian wave generation formalism must be concomitant to understanding the multipole expansion in general relativity.

The multipole expansion is one of the most useful tools of physics, but its use in general relativity is difficult because of the non-linearity of the theory and the tensorial character of the gravitational interaction. In the stationary case, the multipole moments are determined by the expansion of the metric at spatial infinity (Geroch 1970; Hansen 1974; Simon and Beig 1983), while, in the case of non-stationary fields, the moments, starting with the quadrupole, are defined at future null infinity. The multipole moments have been extensively studied in the linearized theory, which ignores the gravitational forces inside the source. Early studies have extended the Einstein quadrupole formula [given by Eq. (5) below] to include the current-quadrupole and mass-octupole moments (Papapetrou 1962, 1971), and obtained the corresponding formulas for linear momentum (Papapetrou 1962, 1971; Bekenstein 1973; Press 1977) and angular momentum (Peters 1964; Cooperstock and Booth 1969). The general structure of the infinite multipole series in the linearized theory was investigated by several works (Sachs and Bergmann 1958; Sachs 1961; Pirani 1965; Thorne 1980), from which it emerged that the expansion is

<footnote n="3">³ Note that in the case of Newtonian binary systems, for high eccentricities (with say $e \to 1^{-}$) the Newtonian potential $U$ can be numerically much larger than the estimate $O(1/c^2) \sim v^2/c^2$ around the apastron of the orbit.</footnote>
characterized by two and only two sets of moments: Mass-type and current-type moments. Below we shall use a particular multipole decomposition of the linearized (vacuum) metric, parametrized by symmetric and trace-free (STF) mass and current moments, as given by Thorne (1980). The expressions of the multipole moments, valid in the linearized theory but irrespective of a slow motion hypothesis, have been worked out by Mathews (1962), Campbell and Morgan (1971), Campbell et al. (1977) and culminated with Damour and Iyer (1991a) obtaining the closed-form expressions of the time-dependent STF mass and current multipole moments as integrals over the source in linearized gravity.

In the full non-linear theory, the (radiative) multipole moments can be read off the coefficient of $1/R$ in the expansion of the metric when $R \to +\infty$, with a null coordinate $U = \text{const}$. The solutions of the field equations in the form of a far-field expansion (power series in $1/R$) have been constructed, and their properties elucidated, by Bondi et al. (1962) and Sachs (1962). The precise way under which such radiative space-times fall off asymptotically has been formulated geometrically by Penrose (1963), Penrose (1965), Geroch and Horowitz (1978) in the concept of an asymptotically simple space-time. The interest of this structure of the asymptotic field arises also from the fact that it is preserved under an infinite set of residual symmetries, the Bondi-Metzner-Sachs (BMS) group, generated by supertranslations and arbitrary diffeomorphisms on the two-sphere at infinity (see Ashtekar et al. 2018 for a review). The Bondi-Sachs-Penrose approach is very powerful, but it can answer a priori only a part of the problem, because it gives information on the field only in the limit where $R \to +\infty$, which cannot be connected in a direct way to the actual matter content and dynamics of the source. In particular the multipole moments that one considers in this approach are those measured at infinity—we call them the *radiative* multipole moments. These moments are distinct, because of non-linearities, from some more natural *source* multipole moments, which are operationally given by explicit integrals extending over the matter and gravitational field of the source.

The problem of non-linearities can be attacked by the post-Minkowskian expansion, valid at once inside and outside the source (Thorne and Kovács 1975; Crowley and Thorne 1977). This approximation has been developed in the pioneering work of Bertotti and Plebański (1960). Ledvinka et al. (2008) obtained a closed form expression for the Hamiltonian of $N$ particles in the first post-Minkowskian (1PM) approximation (see also Blanchet and Fokas 2018). The equations of motion of systems of particles were derived up to 2PM order by Westpfahl and Goller (1979), Westpfahl and Hoyler (1980), Bel et al. (1981), Westpfahl (1985). An important revival of the post-Minkowskian approximation took place recently using scattering amplitudes and the effective field theory (EFT): Bern et al. (2019a, b) derived the 3PM corrections to the conservative two-body Hamiltonian for spinless compact objects; this was extended to 4PM order by Bern et al. (2021), Dlapa et al. (2022) (but for the potential modes only). The problem of gravitational radiation has also progressed using the EFT post-Minkowskian method (Mougiakakos et al. 2021; Brandhuber et al. 2023; Herderschee et al. 2023; Elkhidir et al. 2023; Georgoudis et al. 2023). The amplitude-based derivation of the gravitational waveform generated by the scattering of two scalar black holes has been
successfully compared to the result of the MPM-PN formalism up to 2.5PN order (Bini et al. 2023, 2024; Georgoudis et al. 2024a, b).

A different way of defining the multipole expansion within the complete non-linear theory was proposed by Blanchet and Damour (1986), Blanchet (1987), following works by Bonnor (1959), Bonnor and Rotenberg (1961), Bonnor and Rotenberg (1966), Hunter and Rotenberg (1969) and Thorne (1980). In this approach the basic multipole moments are the *source* moments, rather than the radiative ones.<fnref n="4" /> In a first stage, the moments are left unspecified, as being some arbitrary functions of time, supposed to describe an actual physical source. They are iterated by means of a post-Minkowskian expansion of the vacuum field equations (valid in the source's exterior). Technically, the post-Minkowskian approximation scheme is greatly simplified by the assumption of a multipolar expansion, as one can consider separately the iteration of the different multipole pieces composing the exterior field. In this multipolar-post-Minkowskian (MPM) formalism, which is physically valid over the entire weak-field region outside the source, and in particular in the wave zone (going up to future null infinity), the radiative multipole moments are obtained in the form of some non-linear functionals of the more basic source moments. *A priori*, the method is not limited to post-Newtonian sources; however, in practice, the *closed-form* expressions of the source multipole moments can be established only in the case where the source is post-Newtonian (Blanchet and Damour 1989; Blanchet 1995, 1998b). The reason is that in this case the domain of validity of the post-Newtonian iteration (viz. the near zone) overlaps the exterior weak-field region, so that there exists an intermediate zone in which the post-Newtonian and multipolar expansions can be matched together. This is an application of the method of matched asymptotic expansions in general relativity (Burke and Thorne 1970; Burke 1971; Anderson et al. 1982; Poujade and Blanchet 2002; Blanchet et al. 2005c). The previous MPM-PN construction, i.e., the MPM expansion and its asymptotic matching to the near zone of the PN source, constitutes a powerful gravitational-wave generation formalism for arbitrary isolated sources, able to take into account, in principle, any post-Newtonian correction in the waveform generated by the source. The relationships between the radiative moments and the source moments include many non-linear multipole interactions, as the source moments mix with each other as they “propagate” from the source to the detector. Such multipole interactions include the famous tail effect, corresponding to the coupling between the non-static moments with the mass of the source, and the non-linear memory effect (Blanchet and Damour 1988, 1992). Furthermore Blanchet (1987) proved that this construction, when considered in a neighbourhood of future null infinity, recovers the findings of the Bondi et al. (1962), Sachs (1962) approach to gravitational radiation.

A different wave-generation formalism has been devised by Will and Wiseman (1996) (see also Will 1999; Pati and Will 2000, 2002). This formalism, called “Direct Integration of the Relaxed Equations—DIRE”, has exactly the same scope as the MPM-PN method, i.e., it applies to any isolated post-Newtonian sources, in principle

<footnote n="4">⁴ An alternative definition of the source multipole moments, in terms of canonical Noether charges associated with residual gauge transformations in harmonic coordinates, has been proposed by Compère et al. (2018).</footnote>
up to any PN order. However it differs in the definition of the source multipole moments, that are some well-defined (compact-support) versions of the Epstein and Wagoner (1975), Thorne (1980) moments rather than the Blanchet and Damour (1989) moments, and in the calculation of tails and related non-linear effects. In both formalisms the source multipole moments, which involve a whole series of relativistic PN corrections, must be coupled together in a complicated way in the true non-linear solution; such non-linear couplings form an integral part of the radiative moments at infinity and thereby of the observed signal. The equivalence, at the most general level, between the MPM-PN framework and the DIRE formalism, has been proved in Sect. 4.3 of Blanchet (2014).

### 1.1.2 The quadrupole moment formalism

The lowest-order wave generation formalism is the famous Einstein (1918) quadrupole moment formalism. In his original paper Einstein derived it in the case where the dynamics of the source is driven by non gravitational forces. However a simple derivation of the quadrupole formalism due to Landau and Lifshitz (1971) showed that it is actually also valid for self-gravitating sources, whose dynamics is due to the gravitational force—for instance a Newtonian binary star system.<fnref n="5" />

The quadrupole formalism applies to a general isolated matter source which is post-Newtonian in the sense of existence of the small post-Newtonian parameter $\epsilon$ defined by Eq. (1). However, the quadrupole formalism is valid in the Newtonian limit $\epsilon \to 0$; it can rightly be qualified as “Newtonian” because the quadrupole moment of the matter source is Newtonian and its evolution obeys Newton’s laws of gravity. In this formalism the gravitational field $h_{ij}^{TT}$ (written here as a perturbation of the “gothic” spatial metric $g^{ij} \equiv \sqrt{-g}g^{ij}$) is expressed in a transverse and traceless (TT) coordinate system covering the far zone of the source at retarded times,<fnref n="6" /> as

$$h_{ij}^{TT} = -\frac{2G}{c^4 R} \perp_{ijab}(N) \left\{ \frac{d^2 Q_{ab}}{dT^2} (T - R/c) + \mathcal{O}\left(\frac{1}{c}\right) \right\} + \mathcal{O}\left(\frac{1}{R^2}\right), \quad (2)$$

where $R = |X|$ is the distance to the source, $T - R/c$ is the retarded time, $N = X/R$ is the unit direction from the source to the far away observer, and

$$\perp_{ijab} \equiv \frac{1}{2} [\perp_{ai}\perp_{jb} + \perp_{aj}\perp_{ib} - \perp_{ij}\perp_{ab}], \quad (3)$$

is the TT projection operator, with $\perp_{ij} \equiv \delta_{ij} - N_i N_j$ being the projector onto the plane orthogonal to $N$. The source’s quadrupole moment takes the familiar Newtonian form (see the footnote 23 for further notation)

<footnote n="5">⁵ The fact that the Newtonian conservation laws forbid the emission of dipole gravitational waves was pointed out by Abraham in 1914 (see Kennefick 2007).</footnote>

<footnote n="6">⁶ The TT coordinate system can be extended to the near zone of the source as well; see, e.g., Kopeikin et al. (1999).</footnote>
$$Q_{ij}(t) = \int_{\text{source}} \mathrm{d}^3\mathbf{x} \rho(\mathbf{x}, t) \left(x_i x_j - \frac{1}{3} \delta_{ij} \mathbf{x}^2\right), \quad (4)$$

where $\rho$ is the Newtonian (essentially baryonic) mass density, obeying the usual continuity equation $\partial_t \rho + \partial_i(\rho v^i) = 0$. The total gravitational power emitted by the source in all directions around the source is given by the Einstein quadrupole formula

$$\mathcal{F} \equiv \left(\frac{\mathrm{d}E}{\mathrm{d}T}\right)^{\mathrm{GW}} = \frac{G}{c^5} \left\{ \frac{1}{5} \frac{\mathrm{d}^3 Q_{ab}}{\mathrm{d}T^3} \frac{\mathrm{d}^3 Q_{bc}}{\mathrm{d}T^3} + \mathcal{O}\left(\frac{1}{c^2}\right) \right\}. \quad (5)$$

Our notation $\mathcal{F}$ stands for the total gravitational energy flux or gravitational “lu- minosity” of the source. Similarly, the total angular momentum flux is (Papapetrou 1971; Thorne 1980)

$$\mathcal{G}_i \equiv \left(\frac{\mathrm{d}J_i}{\mathrm{d}T}\right)^{\mathrm{GW}} = \frac{G}{c^5} \left\{ \frac{2}{5} \epsilon_{iab} \frac{\mathrm{d}^2 Q_{ac}}{\mathrm{d}T^2} \frac{\mathrm{d}^3 Q_{bc}}{\mathrm{d}T^3} + \mathcal{O}\left(\frac{1}{c^2}\right) \right\}, \quad (6)$$

where $\epsilon_{abc}$ denotes the standard Levi-Civita symbol with $\epsilon_{123} = 1$. We shall give in Eqs. (35)–(36) the related fluxes of linear momentum and position of the center of mass.

Associated with the latter energy and angular momentum fluxes, there is also a quadrupole formula for the radiation reaction force, which reacts on the source’s dynamics in consequence of the emission of waves. This force will inflect the time evolution of the orbital phase of the binary pulsar and is responsible for the inspiral of compact binaries. At the position $(\mathbf{x}, t)$ in a particular coordinate system covering the source, the reaction force density can be written as (Burke and Thorne 1970; Burke 1971; Misner et al. 1973)

$$F_i^{\text{reac}} = \frac{G}{c^5} \rho \left\{ -\frac{2}{5} x^a \frac{\mathrm{d}^5 Q_{ia}}{\mathrm{d}t^5} + \mathcal{O}\left(\frac{1}{c^2}\right) \right\}. \quad (7)$$

This is the gravitational analogue of the damping force of electromagnetism. How- ever, notice that gravitational radiation reaction is inherently coordinate dependent, so the expression of the force depends on the coordinate system which is used. Consider now the energy and angular momentum of an isolated matter system made of some perfect fluid, say

$$E = \int \mathrm{d}^3 \mathbf{x} \rho \left[ \frac{\mathbf{v}^2}{2} + \Pi - \frac{U}{2} \right] + \mathcal{O}\left(\frac{1}{c^2}\right), \quad (8a)$$

$$J_i = \int \mathrm{d}^3 \mathbf{x} \rho \epsilon_{iab} x_a v_b + \mathcal{O}\left(\frac{1}{c^2}\right). \quad (8b)$$

The specific internal energy of the fluid is denoted $\Pi$, and obeys the usual ther- modynamic relation $d\Pi = -Pd(1/\rho)$ where $P$ is the pressure; the gravitational potential obeys the Poisson equation $\Delta U = -4\pi G\rho$. We compute the mechanical losses of energy and angular momentum from the time derivatives of $E$ and $J_i$. We employ the Eulerian equation of motion $\rho \mathrm{d}v^i/\mathrm{d}t = -\partial_t P + \rho\partial_i U + F_i^{\text{reac}}$ and
continuity equation $\partial_t \rho + \partial_i (\rho v^i) = 0$. Note that we add the small dissipative radiation-reaction contribution $F_i^{\text{reac}}$ in the equation of motion but neglect all conservative post-Newtonian corrections. The result is

$$ \frac{dE}{dt} = \int d^3\mathbf{x} v^i F_i^{\text{reac}} = -\mathcal{F} + \frac{df}{dt}, \qquad (9a) $$

$$ \frac{dJ_i}{dt} = \int d^3\mathbf{x} \, \epsilon_{iab} x_a F_b^{\text{reac}} = -\mathcal{G}_i + \frac{dg_i}{dt}, \qquad (9b) $$

where one recognizes the fluxes at infinity given by Eqs. (5)–(6), and where the second terms denote some total time derivatives made of quadratic products of derivatives of the quadrupole moment. Looking only for secular effects, we apply an average over time on a typical period of variation of the system; the latter time derivatives will be in average numerically small in the case of quasi-periodic motion (Breuer and Rudolph 1981). Hence we obtain

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

$$ \left\langle \frac{dJ_i}{dt} \right\rangle = -\langle \mathcal{G}_i \rangle, \qquad (10b) $$

where the brackets denote the time averaging over an orbit.<fnref n="7" /> These balance equations encode the secular decreases of energy and angular momentum by gravitational radiation emission.

The cardinal virtues of the Einstein (1918), Landau and Lifshitz (1971) quadrupole formalism are: Its generality—the only restrictions are that the source be Newtonian and bounded; its simplicity, as it necessitates only the computation of the time derivatives of the Newtonian quadrupole moment (using the Newtonian laws of motion); and, most importantly, its agreement with the observation of the dynamics of the binary pulsar PSR 1913+16 by Taylor et al. (1979), Taylor and Weisberg (1982), Taylor (1993). Indeed let us apply the balance equations (10) to a system of two point masses moving on an eccentric orbit modelling the binary pulsar —the classic references are Peters and Mathews (1963), Peters (1964); see also Dyson (1963), Esposito and Harrison (1975), Wagoner (1975). We use the binary’s Newtonian energy and angular momentum,

$$ \langle f \rangle \equiv \lim_{T \to +\infty} \frac{1}{2T} \int_{-T}^{T} dt f(t). $$

For periodic functions this reduces to the usual average over the period $P$,

$$ \langle f \rangle = \frac{1}{P} \int_{0}^{P} dt f(t). $$

<footnote n="7">⁷ In the most general case, the time average $\langle f \rangle$ of a function $f(t)$ is defined as</footnote>
$$E = -\frac{Gm_1m_2}{2a}, \qquad (11a)$$

$$J = m_1m_2 \sqrt{\frac{Ga(1 - e^2)}{m_1 + m_2}}, \qquad (11b)$$

where $a$ and $e$ are the semi-major axis and eccentricity of the orbit and $m_1$ and $m_2$ are the two masses. From the energy balance equation (10a) (averaged over the orbital period $P$) we obtain first the secular evolution of $a$; next changing from $a$ to the orbital period $P$ using Kepler's third law,<fnref n="8" /> we get the secular evolution of the orbital period $P$ as

$$\left\langle \frac{dP}{dt} \right\rangle = -\frac{192\pi}{5c^5} \left(\frac{2\pi G}{P}\right)^{5/3} \frac{m_1 m_2}{(m_1 + m_2)^{1/3}} \frac{1 + \frac{73}{24}e^2 + \frac{37}{96}e^4}{(1 - e^2)^{7/2}}. \qquad (12)$$

The last factor, depending on the eccentricity, comes out from the orbital average and is known as the Peters and Mathews (1963) “enhancement” factor,<fnref n="9" /> so designated because in the case of the binary pulsar, which has a rather large eccentricity $e \simeq 0.617$, it enhances the effect by a factor $\sim 12$. Numerically, one finds $\langle dP/dt \rangle = -2.4 \times 10^{-12}$, a dimensionless number in excellent agreement with the observations of the binary pulsar (Taylor et al. 1979; Taylor and Weisberg 1982; Taylor 1993). More recently even more impressive tests were performed with the double pulsar (Kramer and Wex 2009; Kramer et al. 2021). On the other hand the secular evolution of the eccentricity $e$ is deduced from the angular momentum balance equation (10b) [together with the previous result (12)], as

$$\left\langle \frac{de}{dt} \right\rangle = -\frac{608\pi}{15c^5} \frac{e}{P} \left(\frac{2\pi G}{P}\right)^{5/3} \frac{m_1 m_2}{(m_1 + m_2)^{1/3}} \frac{1 + \frac{121}{304}e^2}{(1 - e^2)^{5/2}}. \qquad (13)$$

Interestingly, the system of equations (12)–(13) can be thoroughly integrated in closed analytic form. This yields the law of evolution of the eccentricity $e(t)$ in terms of the period $P(t)$ at any instant (Peters 1964):

$$\frac{e^2}{(1 - e^2)^{19/6}} \left(1 + \frac{121}{304} e^2\right)^{145/121} = c_0 P^{19/9}, \qquad (14)$$

where $c_0$ denotes an integration constant to be determined by the initial conditions at the start of the binary evolution. When $e \ll 1$ the latter relation gives approximately $e^2 \simeq c_0 P^{19/9}$.

Inspiring compact binaries, containing neutron stars and/or black holes, constitute the bread-and-butter sources of gravitational waves for the detectors LIGO, Virgo and KAGRA on ground. The two compact objects steadily lose their orbital binding energy by emission of gravitational radiation; as a result, the orbital

<footnote n="8">⁸ Namely $(Gm)^1 = \Omega^2 a^3$, where $m = m_1 + m_2$ is the total mass and $\Omega = 2\pi/P$ is the orbital frequency. This law is also appropriately called the 1–2–3 law (Misner et al. 1973).</footnote>

<footnote n="9">⁹ See also Eqs. (549)–(551) below for more details on its derivation.</footnote>
separation between them decreases, and the orbital frequency increases. Thus, the frequency of the gravitational-wave signal, which equals twice the orbital frequency for the dominant harmonics, “chirps” in time (i.e., the signal becomes higher and higher pitched) until the two objects collide and merge.

The orbit of most inspiralling compact binaries can be considered to be circular, apart from the gradual inspiral, because the gravitational radiation reaction forces tend to rapidly circularize the motion. This effect is due to the emission of angular momentum by gravitational waves, resulting in a secular decrease of the eccentricity of the orbit, which has been computed within the quadrupole formalism in Eq. (13). For instance, suppose that the inspiralling compact binary was long ago (a few hundred million years ago) a system similar to the binary pulsar system, with an orbital frequency $\Omega_0 \equiv 2\pi/P_0 \sim 10^{-4}$ rad/s and a rather large orbital eccentricity $e_0 \sim 0.6$. When it becomes visible by the detectors on ground, i.e., when the gravitational-wave signal frequency reaches about $f \equiv \Omega/\pi \sim 10$ Hz, the eccentricity of the orbit should be $e \sim 10^{-6}$ according to the formula (14). This is a very small eccentricity, even when compared to high-order relativistic corrections. Only non-isolated binary systems could have a non negligible eccentricity. For instance, the von Zeipel (1910), Lidov (1962), Kozai (1962) (ZLK) mechanism<fnref n="10" /> is one important scenario that produces eccentric binaries and involves the interaction between the binary and the environment of the dense cores of globular clusters (Miller and Hamilton 2002). If the mutual inclination angle of the inner compact binary is strongly tilted with respect to the outer compact star (in a “hierarchical” triplet), then a resonance occurs and can increase the eccentricity of the inner binary to large values. This is one motivation for looking at the waves emitted by inspiralling binaries in non-circular, quasi-elliptical orbits (see Sect. 3.5).

For a long while, it was thought that the quadrupole formula would be sufficient for sources of gravitational radiation to be observed directly on Earth—as it had proved to be amply sufficient in the case of the binary pulsar. However, the work of Cutler et al. (1993a)<fnref n="11" /> showed that this is not true, as one has to include post-Newtonian corrections up to high order beyond the quadrupole formula in order to prepare for the data analysis of future detectors. This expectation has been confirmed by many measurement-analyses (Cutler et al. 1993b; Finn and Chernoff 1993; Cutler and Flanagan 1994; Tagoshi and Nakamura 1994; Poisson 1995; Poisson and Will 1995; Królak et al. 1995; Damour et al. 1998, 2000a, 2002, 2003; Buonanno et al. 2003a, b; Ajith et al. 2005; Arun et al. 2005; Buonanno et al. 2009a) which have demonstrated that the post-Newtonian precision needed for the LIGO/Virgo detectors corresponds at least, in the case of neutron-star binaries, to the $3\text{PN} \sim 1/c^6$ approximation beyond the quadrupole formula.

In practice, the post-Newtonian prediction for the inspiral phase has to be matched to numerical-relativity results for the subsequent merger and ringdown phases. The match proceeds essentially through two routes: Either the effective-one-body (EOB)

<footnote n="10">¹⁰ See Ito and Ohtsuka (2019) for an historical discussion, and justification of the appellation ZLK for designating this effect.</footnote>

<footnote n="11">¹¹ This work entitled: “The last three minutes: Issues in gravitational-wave measurements of coalescing compact binaries” is sometimes coined the “3mn Caltech paper”.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Observational constraints on the post-Newtonian parameters, i.e. the coefficients in the phasing formula (487), from measurements of the black hole events GW150914 and GW151226 (top panel) and from the neutron star event GW170817 (bottom panel) (Abbott et al. 2016a, b, 2019). The limits are obtained by assuming the values predicted by general relativity for all the PN parameters but for one. This one is allowed to vary and its deviation with respect to GR is measured by the technique of matched filtering. The 1.5PN parameter agrees with the GR prediction within a fractional accuracy of the order of 10%, which constitutes an interesting test of the tail effect. Images reproduced with permission from [top] Abbott et al. (2016a), copyright by the author(s) and [bottom] from Abbott et al. (2019), copyright by APS</figcaption></figure>
waveform (Buonanno and Damour 1999, 2000; Damour et al. 2000c; Damour and Nagar 2011) that builds on post-Newtonian results and extend their domain of validity to facilitate the comparison with numerical relativity (Buonanno et al. 2009b; Pan et al. 2010), or the Hybrid waveform obtained by direct matching between the PN expanded prediction and the numerical computations (Ajith et al. 2008; Santamaría et al. 2010).

Strategies to detect and analyze the very weak signals from compact binary inspiral are based on 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 (supposedly 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 (15)$$

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 the Wiener (1942) filter

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

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, (16) 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 prediction 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. As we shall see, the signal at such high PN order contains the signature of several non-linear effects which are specific to general relativity. We thus have the possibility of probing, via matched filtering, some aspects of the non-linear structure of Einstein's theory (Blanchet and Sathyaprakash 1994, 1995; Arun et al. 2006a, b). See Fig. 1.

### 1.1.3 Influence of the internal structure of compact bodies

The main point about modelling the inspiralling compact binary is that a model made of two structureless point particles, characterized solely by mass parameters $m_a$ and possibly the spins $S_a$ (with a labelling the particles), is sufficient in first approximation. Indeed, most of the non-gravitational effects usually plaguing the dynamics of binary star systems, such as the effects of a magnetic field, of an interstellar medium, of the internal structure of extended bodies, are dominated by
gravitational effects. In particular the effects due to the finite size of the compact bodies are small, although not negligibly small.

Consider the influence of the Newtonian quadrupole moments $q_a^{ij}$ through the tidal interaction between extended neutron stars without spins during the inspiral phase. It is known that for inspiralling compact binaries the neutron stars are not co-rotating because the tidal synchronization time is much larger than the time left till the coalescence. As shown by Kochanek (1992) the best models for the fluid motion inside the two neutron stars are the so-called Roche-Riemann ellipsoids, which have tidally locked figures (the quadrupole moments facing each other at any instant during the inspiral), but for which the fluid motion has zero circulation in the inertial frame.

Here we perform a simple calculation to Newtonian order within such a model. The equations of motion of N extended spinless bodies, to linear order in the quadrupole moments, are (with a, b = 1, ..., N)

$$m_a \frac{dv_a^i}{dt} = G \sum_{b \neq a} \left[ m_a m_b \frac{\partial}{\partial y_a^i} \left( \frac{1}{r_{ab}} \right) + \frac{1}{2} \left( m_a q_b^{jk} + m_b q_a^{jk} \right) \frac{\partial^3}{\partial y_a^i \partial y_a^j \partial y_a^k} \left( \frac{1}{r_{ab}} \right) \right], \quad (17)$$

where $m_a$ are the masses, and we denote the position and velocity of the center of mass of the bodies by $\mathbf{y}_a(t)$ and $\mathbf{v}_a(t) = \frac{dy_a}{dt}$, with the Euclidean separation between centers of mass being $r_{ab} = |\mathbf{y}_a - \mathbf{y}_b|$. The quadrupole moments of the bodies, supposed to be made of a perfect fluid, read

$$q_a^{ij} = \int_{V_a} d^3\mathbf{z}_a \rho_a \left( z_a^j z_a^i - \frac{1}{3} \delta^{ij} z_a^2 \right), \quad (18)$$

with $V_a$ the volume of the body, $z_a = \mathbf{x} - \mathbf{y}_a(t)$ the distance between a generic point $\mathbf{x}$ inside the body and the center of mass, $\rho_a = \rho(\mathbf{y}_a + z_a, t)$ the Newtonian mass density of the body, $\rho(\mathbf{x}, t)$ being the usual Eulerian density. The mass-centred condition reads

$$\int_{V_a} d^3\mathbf{z}_a \rho_a z_a^i = 0. \quad (19)$$

The conserved energy of the N-body system is the sum of the internal (Newtonian) energies $e_a$ and of the orbital contributions, including the effect of the quadrupoles:

$$E = \sum_a \left\{ e_a + \frac{1}{2} m_a \mathbf{v}_a^2 - \frac{G}{2} \sum_{b \neq a} \frac{m_a m_b}{r_{ab}} - \frac{1}{2} q_a^{ij} \mathcal{E}_a^{ij} \right\}, \quad (20)$$

where we have introduced the tidal field acting on body a and due to the other bodies $b \neq a$:
$$ \mathcal{E}_a^{ij} = G \sum_{b \neq a} m_b \frac{\partial^2}{\partial y_a^j \partial y_b^i} \left( \frac{1}{r_{ab}} \right). \qquad (21) $$

Posing $w_a = dz_a/dt$ for the internal velocity field of body a, $\Pi_a = \Pi(y_a + z_a, t)$ for the specific internal energy satisfying the thermodynamical relation $d\Pi = -P d(1/\rho)$ (with P the pressure), and $u_a$ for the internal self-gravity given by the Poisson integral over the volume of the body, we have

$$ e_a = \int_{\mathcal{V}_a} d^3 z_a \rho_a \left( \frac{1}{2} w_a^2 + \Pi_a - \frac{u_a}{2} \right). \qquad (22) $$

The coupling of the quadrupole moment $q_a^{ij}$ with the external tidal field $\mathcal{E}_a^{ij}$ of the other bodies implies a variation of the internal energy given by

$$ \frac{\mathrm{d}e_a}{\mathrm{d}t} = \frac{1}{2} \frac{\mathrm{d}q_a^{ij}}{\mathrm{d}t} \mathcal{E}_a^{ij}. \qquad (23) $$

We consider the case where the quadrupole moment is induced by the tidal field of the other bodies. To linear order, we can introduce a coefficient $\lambda_a$ characterizing the response of the body under the influence of the external field, such that

$$ q_a^{ij} = \lambda_a \mathcal{E}_a^{ij}. \qquad (24) $$

In the case of the induced quadrupole moments (24), the total energy of the system becomes

$$ E = \sum_a \left\{ \frac{1}{2} m_a v_a^2 - \frac{G}{2} \sum_{b \neq a} \frac{m_a m_b}{r_{ab}} - \frac{1}{4} \lambda_a \mathcal{E}_a^{ij} \mathcal{E}_a^{ij} \right\}. \qquad (25) $$

The response coefficient $\lambda_a$ is conventionally written as $\lambda_a = \frac{2}{3G} k_a r_a^5$, where $r_a$ is the radius of the body and $k_a$ is the body's dimensionless second Love (1911) number which depends on the internal structure of the body, and is, typically, of the order unity. Next we characterize the internal structure by the dimensionless “deformability” or “polarizability” parameter

$$ \Lambda_a = \frac{c^{10}}{G^4 m_a^5} \lambda_a = \frac{2}{3} k_a \left( \frac{c^2 r_a}{G m_a} \right)^5 = \frac{2 k_a}{3 K_a^5}. \qquad (26) $$

where the dimensionless ratio $K_a \equiv \frac{Gm_a}{r_a c^2}$ is the “compactness” parameter and typically equals $\sim 0.2$ for neutron stars (depending on their equation of state).

Consider a compact binary system ($N=2$) moving on an exact circular orbit. From Eq. (24) we see that the two quadrupole moments face each other, and remain constant along the circular orbit. The equation of the relative motion reduces to $d\mathbf{v}/dt = -\Omega^2\mathbf{x}$, where $\mathbf{x} = \mathbf{y}_1 - \mathbf{y}_2$ and $\mathbf{v} = d\mathbf{x}/dt$ (with $r \equiv r_{12}$). We find from (17) the orbital frequency
$$ \Omega^2 = \frac{Gm}{r^3} \left[ 1 + 9\nu(X_1^3 \Lambda_1 + X_2^3 \Lambda_2)\gamma^5 \right]. \qquad (27) $$

We pose $X_a = m_a/m$ with $m = m_1 + m_2$ so that $\nu = X_1X_2$ is the symmetric mass ratio (353), denote $\gamma = \frac{Gm}{rc^2}$ and employ the notation (26). In turn the conserved energy (25) for circular orbits reduces to

$$ E = -\frac{Gm^2\nu}{2r} \left[ 1 - 6\nu(X_1^3 \Lambda_1 + X_2^3 \Lambda_2)\gamma^5 \right]. \qquad (28) $$

As the effect of the deformation of the bodies computed here is purely “Newtonian”, we see that the c’s we have introduced into our definitions naturally cancel out in Eqs. (27) and (28).

The above dynamics is *conservative*, as we have neglected the dissipative radiation reaction effect on the orbit. This effect is taken into account when we impose the flux-balance equation (10a); there is no need to impose the angular momentum balance equation (10b) for circular orbits. The total quadrupole moment of the system is the sum of the orbital one and of the intrinsic moments of the bodies, given by (18), hence

$$ Q^{ij} = m\nu(x^i x^j - \frac{1}{3}\delta^{ij}r^2) + q_1^{ij} + q_2^{ij}. \qquad (29) $$

Plugging this into the flux formula (5), computing the time derivatives using the equations of motion including the contributions from the quadrupole moments, see Eq. (27), and keeping only the terms linear in these quadrupoles, yields the flux (still for exact circular orbits) as

$$ \mathcal{F} = \frac{32G}{5c^5} r^4 \omega^6 m^2 \nu^2 \left[ 1 + 6(X_1^4 \Lambda_1 + X_2^4 \Lambda_2)\gamma^5 \right]. \qquad (30) $$

Next we reexpress the invariants E and $\mathcal{F}$ in terms of the orbital frequency $\omega$ instead of the coordinate-dependent separation distance $r$ using Eq. (27). Posing $x = (Gm\Omega/c^3)^{2/3}$ we obtain

$$ E = -\frac{1}{2}mvc^2x \left[ 1 - 9\nu(X_1^3 \Lambda_1 + X_2^3 \Lambda_2)\gamma^5 \right], \qquad (31a) $$

$$ \mathcal{F} = \frac{32c^5}{5G} x^5 \nu^2 \left\{ 1 + 6 \left[ (X_1 + 2\nu)X_1^3 \Lambda_1 + (X_2 + 2\nu)X_2^3 \Lambda_2 \right] x^5 \right\}. \qquad (31b) $$

At this stage we can draw a firm conclusion: The effect of the internal structure of non-spinning bodies is proportional to $x^5$, and is thus comparable to a relativistic effect occuring at the 5PN order. This is a higher post-Newtonian order than the one at which we shall obtain the phasing formula in Sect. 3.4.4. Therefore, in that sense the internal structure of the two compact bodies is “effaced” and their dynamics and radiation depend only, in first approximation, on the masses and possibly the spins (Damour 1983a).
But the numerical value of the coefficient involving the $Λ_a$'s is to be taken into account. For instance $\frac{Gm_a}{c^2r_a} \sim 0.15$ for neutron stars (hence $x^5 \sim 8 \times 10^{-5}$ at the merger), and the numerical estimates of the Love numbers for neutron stars are of the order of one or say, a tenth (Hinderer 2008; Damour and Nagar 2009; Binnington and Poisson 2009). Therefore the deformability parameters (325) for compact bodies should be of the order of $Λ_a \sim 1000$, depending on the equation of state, see Fig. 2.

The phase and frequency evolution follow from the balance equation (447) (without need to average), where both E and $\mathcal{F}$ have been computed for the conservative dynamics in Eqs. (31). This approximation is justified as we are interested in the secular, adiabatic evolution of the orbit over a radiation reaction time scale. The secular variation of the energy E is immediately deduced from (31a) as

$$ \frac{dE}{dt} = -\frac{1}{2}mvc^2 \left[ 1 - 54v \left( X_1^3 \Lambda_1 + X_2^3 \Lambda_2 \right) x^5 \right] \dot{x}. \qquad (32) $$

Combining this with (31b), we get an ordinary differential equation for $x$. Solving for $x(t)$ and the phase evolution $\phi = \int dt \Omega(t)$, we thereby obtain the quadrupole finite size effect due to the internal structure as (Flanagan and Hinderer 2008; Damour et al. 2012; Favata 2014)$^{12}$

$$ x = \frac{1}{4} \tau^{-1/4} \left[ 1 + \frac{39}{8192} \tilde{\Lambda} \tau^{-5/4} \right], \qquad (33a) $$

$$ \phi = \phi_0 - \frac{x^{-5/2}}{32v} \left[ 1 + \frac{39}{8} \tilde{\Lambda} x^5 \right], \qquad (33b) $$

where $\tau = \frac{vc^3}{5Gm}$ ($t_0 - t$) and $t_0$ is the instant of coalescence. The effect depends on the following combination of the two deformability parameters:

$$ \tilde{\Lambda} = \frac{16}{13} \left[ (X_1 + 11v)X_1^3 \Lambda_1 + (X_2 + 11v)X_2^3 \Lambda_2 \right], \qquad (34) $$

so normalized that in the case of two identical neutron stars (with the same mass, $X_1 = X_2 = \frac{1}{2}$, and the same equation of state) it reduces to $\tilde{\Lambda} = \Lambda_1 = \Lambda_2$. The 5PN finite size effect (33)–(34) is to be contrasted with the point-mass result, ignoring the internal structure, which will be obtained in Sect. 3.4.4. Remarkably, it has been possible to put a bound on the tidal deformability of neutron stars, and a constraint on several possible equations of state, with the binary neutron star event GW170817 (Abbott et al. 2017b), see Fig. 2.

<footnote>$^{12}$ The leading-order 5PN tidal phase was extended to higher order by Vines et al. (2011), Bini et al. (2012), Vines and Flanagan (2013), Abdelsalhin et al. (2018), Banihashemi and Vines (2020), Landry (2018), Henry et al. (2020a, b) and is now complete up to the 7.5PN order, including mass quadrupole, current quadrupole and mass octupole tidal contributions, with tail terms arising at 6.5PN and 7.5PN orders.</footnote>
### 1.1.4 The recoil by gravitational waves

A related topic is the loss of linear momentum by gravitational radiation and the resulting gravitational recoil (or “kick”) of black-hole binary systems. This phenomenon has potentially important astrophysical consequences (Merritt et al. 2004). In models of formation of massive black holes involving successive mergers of smaller “seed” black holes, a recoil with sufficient velocity could eject the system from dwarf galaxies or globular clusters and terminate the process. Even in galaxies whose potential wells are deep enough to confine the recoiling system, displacement of the system from the center could have important dynamical consequences for the galactic core.

To compute the recoil we need the flux formula for the linear momentum. Integrating the leading 2.5PN radiation-reaction force (7) over the source results in a total time derivative, so there is no net recoiling force at that order, and the effect is a subdominant 3.5PN one. We have (Bonnor and Rotenberg 1961; Peres 1962; Papapetrou 1962, 1971; Bekenstein 1973; Thorne 1980)

$$ \mathcal{F}_{\text{P}}^{i} \equiv \left(\frac{\mathrm{d}P_{i}}{\mathrm{d}T}\right)^{\text{GW}} = \frac{G}{c^7} \left[ \frac{2}{63} \mathcal{O}_{iab}^{(4)} \mathcal{Q}_{ab}^{(3)} + \frac{16}{45} \epsilon_{iab} \mathcal{Q}_{ac}^{(3)} \mathcal{C}_{bc}^{(3)} \right] + \mathcal{O}\left(\frac{1}{c^9}\right), \quad (35) $$

where the Newtonian source moments are $Q_{ij} = \int d^3\mathbf{x} \rho \hat{x}_{ij}$ (mass quadrupole), $O_{ijk} = \int d^3\mathbf{x} \rho \hat{x}_{ijk}$ (mass octupole) and $C_{ij} = \int d^3\mathbf{x} \rho \epsilon_{ab(i} \hat{x}_{j)a} v_b$ (current quadrupole); see the footnote 23 for notation. Furthermore, besides the fluxes associated to energy E, angular momentum $J_i$ and linear momentum $P_i$, there is also a flux associated to the center of mass $G_i = P_i t + Z_i$, where $Z_i$ denotes the initial “position” of the center of mass. The effect is also subdominant at 3.5PN order and reads (Kozameh and Quirega 2016; Kozameh et al. 2018; Nichols 2018; Blanchet and Faye 2019; Compère et al. 2020)<fnref n="13" />

$$ \mathcal{F}_{\text{G}}^{i} \equiv \left(\frac{\mathrm{d}G_{i}}{\mathrm{d}T}\right)^{\text{GW}} = \frac{2G}{21c^7} \mathcal{O}_{iab}^{(3)} \mathcal{Q}_{ab}^{(3)} + \mathcal{O}\left(\frac{1}{c^9}\right), \quad (36) $$

The meaning of the latter flux formulas is that the secular evolution by gravitational wave emission of the linear momentum **P** and of the center-of-mass position **G** attributable to an isolated (post-Newtonian) system obey

$$ \frac{\mathrm{d}\mathbf{P}}{\mathrm{d}t} = -\mathcal{F}_{\text{P}}, \quad (37a) $$

$$ \frac{\mathrm{d}\mathbf{G}}{\mathrm{d}t} = \mathbf{P} - \mathcal{F}_{\text{G}}, \quad (37b) $$

where $\mathcal{F}_P$ and $\mathcal{F}_G$ are the fluxes given by (35) and (36). Consider the case where the source is stationary before some instant $t_0$, then emits a pulse of gravitational waves with finite duration between times $t_0$ and $t_1$, and finally comes back to a stationary

<footnote n="13">¹³ We provide the full multipole expansions of the various fluxes in terms of radiative-type moments in Eqs. (179) and (181) below.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2 Observational constraints on the tidal deformability and the inner equation of state of neutron stars obtained with the GW170817 event. The parameters $\Lambda_a$ are defined by (26). Contours enclosing 90% and 50% of the probability density are shown with dashed lines. The predictions for tidal deformability given by a set of representative equations of state are given with grey lines. For a stiff equation of state the pressure increases a lot for a given increase in density (for instance $P \propto \rho''$ with a large polytropic index $\gamma$), thus it gives more resistance to the gravitational force and the neutron star is less compact. The stiffest equations of state are excluded, while the softest (which predict more compact neutron stars) are still allowed; they appear in the dark blue region. The constraints are shown for a low-spin scenario, with dimensionless spin parameter $|\chi| \le 0.05$, probably favored for neutron stars. Image reproduced with permission from Abbott et al. (2017b), copyright by the author(s)</figcaption></figure>

state at later times $t > t_1$. In this situation, it is straightforward to find the form of the solution to Eqs. (37). Initially, the linear momentum is constant, so that, by applying a Lorentz boost, we can put ourselves in the rest frame of the system, thus achieving $\mathbf{P}_0 = 0$ (for $t < t_0$). Furthermore, we can translate the origin of our coordinate system in such a way that it coincides with the center of mass of the system, hence $\mathbf{G}_0 = 0$ initially. Then, by integrating (37), we get (for $t_0 < t < t_1$)

$$ \mathbf{P}(t) = - \int_{t_0}^{t} dt' \mathcal{F}_{\mathrm{P}}(t'), \qquad (38a) $$

$$ \mathbf{G}(t) = - \int_{t_0}^{t} dt' (t - t') \mathcal{F}_{\mathrm{P}}(t') - \int_{t_0}^{t} dt' \mathcal{F}_{\mathrm{G}}(t'). \qquad (38b) $$

After the period of emission (for $t > t_1$), the source is again stationary but has acquired a net constant linear momentum $\mathbf{P}_1$, and the motion of its center of mass is uniform, $\mathbf{G}_1 = \mathbf{P}_1 t + \mathbf{Z}_1$. We find
$$ \mathbf{P}_1 = - \int_{t_0}^{t_1} dt' \mathcal{F}_P(t'), \qquad (39a) $$

$$ \mathbf{Z}_1 = \int_{t_0}^{t_1} dt' \left[ t' \mathcal{F}_P(t') - \mathcal{F}_G(t') \right]. \qquad (39b) $$

The final value of the linear momentum $\mathbf{P}_1$ yields the total gravitational recoil velocity of the source as measured in the asymptotic Minkowskian frame. On the other hand, the cumulative effect of the center-of-mass flux $\mathcal{F}_G$ results in a modification of the position $\mathbf{Z}_1$ of the center of mass after the GW emission.

Consider the case of a Newtonian binary system with no spins and moving on a quasi-circular orbit. The right-side of the linear momentum flux equation (37a) is straightforward to evaluate with result (Fitchett 1983)

$$ \frac{d\mathbf{P}}{dt} = \frac{464}{105} \frac{G^4 m^5 \Omega}{c^7 r^4} \Delta v^2 \boldsymbol{\lambda} + \mathcal{O}\left(\frac{1}{c^9}\right). \qquad (40) $$

Here $m = m_1 + m_2$ and $v = m_1 m_2 / m^2$, assuming $m_1 \ge m_2$ (thus $\Delta = \frac{m_1 - m_2}{m} = \sqrt{1 - 4v}$), the orbital frequency of the circular orbit is $\Omega = \sqrt{Gm/r^3}$, and the unit vector $\boldsymbol{\lambda}$ is orthogonal to the orbital separation $\mathbf{n} = (\mathbf{y}_1 - \mathbf{y}_2)/r$ in the orbital plane, such that $\dot{\mathbf{n}} = \Omega \boldsymbol{\lambda}$ and $\dot{\boldsymbol{\lambda}} = -\Omega \mathbf{n}$. The relation (40) holds at any time along the orbit and can be integrated, yielding

$$ \mathbf{P} = \frac{464}{105} \frac{G^4 m^5}{c^7 r^4} \Delta v^2 \mathbf{n} + \mathcal{O}\left(\frac{1}{c^9}\right). \qquad (41) $$

We assume from now on that a Lorentz boost and a shift of the origin of the coordinate system have been applied to set $\mathbf{P}$ and $\mathbf{G}$ to zero when averaged over an orbit (neglecting the radiation-reaction decay). Then, we use the center-of-mass balance equation (37b), which leads to

$$ \frac{d\mathbf{G}}{dt} = \mathbf{P} + \frac{544}{105} \frac{G^4 m^5}{c^7 r^4} \Delta v^2 \mathbf{n} + \mathcal{O}\left(\frac{1}{c^9}\right). \qquad (42) $$

Combining the results (41)–(42) and integrating (Blanchet and Faye 2019),

$$ \mathbf{G} = -\frac{48}{5} \frac{G^4 m^5}{c^7 r^4 \Omega} \Delta v^2 \boldsymbol{\lambda} + \mathcal{O}\left(\frac{1}{c^9}\right). \qquad (43) $$

The results give the instantaneous values of the linear momentum and center-of-mass position of a circular binary (neglecting the orbital decay); these are equal to minus those that can be attributed to the gravitational wave field.

The recoil of binary black hole systems has been estimated by Fitchett (1983), Wiseman (1992), Kidder (1995), Blanchet et al. (2005d), Racine et al. (2009), Le Tiec et al. (2010) using post-Newtonian methods,<fnref n="14" /> and by Oohara and Nakamura (1983), Fitchett and Detweiler (1984), Nakamura et al. (1987), Favata et al. (2004)

<footnote n="14">¹⁴ Notably Kidder (1995) gives the dominant PN contribution of the spins to the recoil.</footnote>
using perturbation methods. In parallel the problem of gravitational recoil of coalescing binaries has attracted considerable attention from the numerical relativity community. The numerical computations led to increasingly accurate estimates of the kick velocity from the merger along quasicircular orbits of binary black holes without spins (Campanelli 2005; Baker et al. 2006a) and with spins (Campanelli et al. 2007). In particular they revealed the interesting result that very large kick velocities can be generated in the case of spinning black holes for particular spin configurations.

Unfortunately the post-Newtonian approach is not ideally suited for this problem because most of the recoil is generated in the strong field regime close to the merger. As a result the contribution of the “plunge” (difficult to model by post-Newtonian theory) dominates over that of the inspiral phase, as shown by the curve “up to the ISCO” in Fig. 3. The post-Newtonian corrections in the recoil formula, up to 2PN order beyond the dominant 3.5PN effect given by Eq. (40), have been computed using the gravitational-wave generation formalism of Sect. 2, with result [where $x \equiv (\frac{Gm\Omega}{c^3})^{2/3}$] (Blanchet et al. 2005d)

$$ \frac{d\mathbf{P}}{dt} = \frac{464}{105} \Delta v^2 x^{11/2} \left[ 1 + \left( -\frac{452}{87} - \frac{1139}{522}v \right) x + \frac{309}{58} \pi x^{3/2} + \left( -\frac{71345}{22968} + \frac{36761}{2088}v + \frac{147101}{68904}v^2 \right) x^2 \right] \lambda + \mathcal{O}\left(\frac{1}{c^{12}}\right). \quad (44) $$

This result has been extended to 2.5PN order by Mishra et al. (2012).

The recoil of non-spinning black-hole binaries generated by the inspiral and merger phases to 2PN relative order is shown in Fig. 3. The final result, however, must include also the recoil generated during the ringdown phase, and is presented in Fig. 4. The result is in good agreement with numerical computations for non-spinning black hole binaries (Campanelli 2005; Baker et al. 2006a; Herrmann et al. 2007; González et al. 2007, 2009). See notably the sequence of dots and accompanying dashed lines in Fig. 4 (in blue), which are from an exhaustive series of numerical simulations by González et al. (2007), González et al. (2009). Also shown in Fig. 4 are the comparisons with other analytic or semi-analytic methods: an application of the effective-one-body formalism (Damour and Gopakumar 2006); a close-limit calculation with Bowen-York type initial conditions (Sopuerta et al. 2006). An empirical fit to the final result from Le Tiec et al. (2010) in Fig. 4 is

$$ V_{\text{recoil}} \approx 9.5 \Delta v^2 (1 + 0.3 v) \times 10^3 \text{ km/s}, \quad (45) $$

where the leading $\Delta v^2$ dependence on the mass ratio comes from the lowest-order “Newtonian” expression, see Eq. (44).

## 1.2 Overview on motion and radiation

### 1.2.1 Post-Newtonian equations of motion

By equations of motion we mean the explicit expression of the accelerations of the bodies in terms of the positions and velocities. In Newtonian gravity, writing the
<figure><img src="image_4.png" /><figcaption>**Fig. 3** The gravitational recoil of non-spinning black-hole binaries generated by the inspiral + merger phases (up to the horizon), as a function of the symmetric mass ratio $\eta \equiv v$. The maximum recoil due to the inspiral phase up to the innermost stable circular orbit (ISCO) is of the order of 22 km s⁻¹. The recoil accumulated during the plunge phase, from the ISCO up to the horizon, is obtained by integrating the 2PN momentum flux formula (44) along a plunge geodesic of the Schwarzschild metric within an effective one-body approach. The maximum recoil due to the inspiral + merger phases (ignoring the ringdown) amounts to about 250 km s⁻¹ (red curve). Image reproduced with permission from Blanchet et al. (2005d), copyright by AAS</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4** The total recoil of non-spinning black hole binaries generated by the inspiral + merger + ringdown phases (green curve). The recoil up to the merger is reproduced from Fig. 3 (red curve). The ringdown contribution is computed using a “close-limit approximation” for black hole binaries that uses 2PN-accurate initial data (Le Tiec and Blanchet 2010). The effect of the ringdown phase on the recoil velocity is to produce an “anti-kick”, i.e. to reduce the value of the total recoil with respect to that computed up to the horizon. Thus the maximum recoil of non-spinning black-hole binaries is around 180 km s⁻¹ at a mass ratio of $v_{max} \approx 0.2$. Image reproduced with permission from Le Tiec et al. (2010), copyright by IOP</figcaption></figure>
equations of motion for a system of $N$ particles is trivial; in general relativity, even writing the equations in the case $N = 2$ is difficult. The first relativistic terms, at the 1PN order, were derived by Droste (1917), Lorentz and Droste (1937). Subsequently, Einstein et al. (1938) obtained the 1PN corrections by means of their famous “surface-integral” method, in which the equations of motion are deduced from the *vacuum* field equations, and are therefore applicable to any compact objects (be they neutron stars, black holes, or, perhaps, naked singularities). The 1PN-accurate equations were also obtained, for the motion of the centers of mass of compact bodies, by Fock (1939) (see also Petrova 1949; Papapetrou 1951a).

The 2PN approximation was tackled by Ohta et al. (1973), Okamura et al. (1973), Ohta et al. (1974a, b), who considered the post-Newtonian iteration of the Hamiltonian of point-particles. We refer here to the Hamiltonian as a “Fokker-type” Hamiltonian, which is obtained from the matter-plus-field Arnowitt–Deser–Misner (ADM) Hamiltonian by eliminating the field degrees of freedom. The 2.5PN equations of motion were obtained in harmonic coordinates by Damour and Deruelle (1981a, b), Damour (1982), Damour (1983a), building on a non-linear (post-Minkowskian) iteration of the metric of two particles initiated by Bel et al. (1981). The corresponding result for the ADM-Hamiltonian of two particles at the 2PN order was given by Damour and Schäfer (1985) (see also Schäfer 1985, 1986). The 2.5PN equations of motion have also been derived in the case of two *extended* compact objects by Kopeikin (1985), Grishchuk and Kopeikin (1986). The 2.5PN equations of two point masses as well as the near zone gravitational field in harmonic-coordinate were computed by Blanchet et al. (1998).

Up to the 2PN level the equations of motion are conservative. Only at the 2.5PN order does the first non-conservative effect appear, associated with the gravitational radiation emission. The equations of motion up to that level have been used for the study of the radiation damping of the binary pulsar—its orbital $\dot{P}$ (Damour 1983a, b; Damour and Taylor 1991). The result is in agreement with the prediction of the quadrupole formalism given by Eq. (12). An important point is that the 2.5PN equations of motion have been proved to hold in the case of binary systems of strongly self-gravitating bodies (Damour 1983a). This is via the effacing principle for the internal structure of the compact bodies. As a result, the equations depend essentially only on the “Schwarzschild” masses $m_a$ of the neutron stars (and their spins). Notably the compactness parameters $K_a = \frac{Gm_a}{r_a c^2}$ do not enter the equations of motion before the 5PN level (see the discussion in Sect. 1.1.3). This fact has been explicitly verified up to the 2.5PN order by Kopeikin (1985), Grishchuk and Kopeikin (1986), who made a physical computation *à la* Fock, taking into account the internal structure of two self-gravitating extended compact bodies. The 2.5PN equations of motion have also been obtained by Itoh et al. (2000, 2001) in harmonic coordinates, using a variant (but, much simpler and more developed) of the surface-integral approach of Einstein et al. (1938), that is valid for compact bodies, independently of the strength of the internal gravity.

At the 3PN order the equations of motion have been worked out by several groups, using different methods, and with equivalent results:
1. Jaranowski and Schäfer (1998, 1999, 2000), Damour et al. (2000d, 2001a, b), employ the ADM-Hamiltonian canonical formalism of general relativity;

2. Blanchet and Faye (2000a, b, 2001a, b); de Andrade et al. (2001); Blanchet and Iyer (2003); Blanchet et al. (2004a) compute directly the equations of motion (instead of a Hamiltonian) in harmonic coordinates;

3. Itoh and Futamase (2003), Itoh (2004), Futamase and Itoh (2007) obtain the complete 3PN equations of motion in harmonic coordinates, without need of a self-field regularization;

4. Foffa and Sturani (2011) derive the 3PN Lagrangian in harmonic coordinates within the effective field theory.

The 3PN equations of motion contained, at some point, some unspecified numerical coefficient, called an “ambiguity”, which is due to an incompleteness of the Hadamard self-field (ultra-violet, UV) regularization method, see Sect. 3.1.1. This coefficients has been fixed by means of dimensional regularization, both within the ADM-Hamiltonian formalism (Damour et al. 2001a), and the harmonic-coordinates equations of motion (Blanchet et al. 2004a), while dimensional regularization is routinely employed by the effective field theory (Foffa and Sturani 2011). All these works have demonstrated the power of dimensional regularization and its adequateness to the problem of interacting point masses in classical general relativity. By contrast, notice, interestingly, that the surface-integral method (Itoh and Futamase 2003; Itoh 2004) by-passes the need of a UV regularization. The regularizations—both UV and infra-red (IR)—are reviewed in Sect. 3.1.

The effective field theory (EFT) approach—sometimes coined the non-relativistic general relativity (NRGR)—to the problems of motion and radiation of compact binaries, has been pioneered by Goldberger and Rothstein (2006) and extensively developed since then (see Foffa and Sturani 2014; Porto 2016; Levi 2020 for reviews). It borrows techniques from quantum field theory and treats the gravitational interaction between point particles as a classical limit of a quantum field theory, i.e., in the “tree level” approximation. The theory is based on the effective action, defined from a Feynman path integral that integrates over the degrees of freedom that mediate the gravitational interaction. The phase factor in the path integral is built from the standard Einstein-Hilbert action for gravity, augmented by a harmonic gauge fixing term and by the action of particles. The Feynman diagrams naturally show up as a perturbative technique for solving iteratively the Green’s functions.

By itself, computing the equations of motion and radiation field using Feynman diagrams in classical general relativity is not a new idea: Bertotti and Plebański (1960) defined the diagrammatic tree-level perturbative expansion of the Green's functions in classical general relativity; Hari Dass and Soni (1982)<fnref n="15" /> showed how to derive the classical energy-loss formula at Newtonian approximation using tree-level propagating gravitons; Feynman diagrams have been used for the equations of motion up to 2PN order in general relativity (Ohta et al. 1973; Okamura et al. 1973; Ohta et al. 1974a, b) and in scalar-tensor theories (Damour and Esposito-Farèse 1996). Nevertheless, the systematic EFT treatment has proved to be powerful and

<footnote n="15">¹⁵ This reference has an eloquent title: “Feynman graph derivation of the Einstein quadrupole formula”.</footnote>
innovative for the field, e.g., with the introduction of a decomposition of the metric into “Kaluza-Klein type” potentials (Kol and Smolkin 2008), the interesting link with the renormalization group equation (Goldberger and Ross 2010; Goldberger et al. 2014), and the systematization of the generation of diagrams in high PN approximations (Foffa and Sturani 2011).

The 3.5PN terms in the equations of motion correspond to the 1PN relative corrections in the radiation reaction force. They were derived by Iyer and Will (1993, 1995) for point-particle binaries in a general gauge and center-of-mass frame, relying on energy and angular momentum balance equations and the known expressions of the 1PN fluxes. These works have been extended to 2PN order (Gopakumar et al. 1997) and to include the leading spin-orbit effects (Zeng and Will 2007). The result has been established from “first principles” (i.e., not relying on balance equations) in various works at 1PN order (Jaranowski and Schäfer 1997; Pati and Will 2002; Königsdörffer et al. 2003; Nissanke and Blanchet 2005; Itoh 2009) and 2PN order (Leibovich et al. 2023). The 1PN radiation reaction force has also been obtained for general extended fluid systems in the extended Burke and Thorne (1970), Burke (1971) gauge by Blanchet (1993, 1997b). Known also is the contribution of gravitational-wave tails in the equations of motion, which arises at the 4PN order and contains both a conservative part and a 1.5PN relative modification of the radiation damping force (Blanchet and Damour 1988). This 4PN tail-induced correction to the equations of motion has also been derived within the EFT approach (Foffa and Sturani 2013; Galley et al. 2016).

The state-of-the-art on equations of motion is the 4PN approximation, with result in complete mutual agreement between the different derivations:

1. After partial results have been reported by Jaranowski and Schäfer (2012, 2013, 2015) using the ADM Hamiltonian formalism, and by Foffa and Sturani (2012) using the EFT, the first derivation of the complete 4PN dynamics was obtained by Damour et al. (2014), combining the local contributions with the non-local (in time) term due to gravitational-wave tails (Blanchet and Damour 1988). However this calculation contains one “ambiguity” parameter, due to the absence of control of the IR divergences, which has been fixed by matching the near-zone computation to auxiliary results imported from the gravitational self-force (GSF) theory (Bini and Damour 2013, see also Bini et al. 2016; Kavanagh et al. 2015; Hopper et al. 2016). The non-local 4PN Hamiltonian has been transformed into a local one using an order reduction and the dynamics has been transcribed into the effective one-body (EOB) formalism (Damour et al. 2015);

2. The 4PN equations of motion were obtained independently by Bernard et al. (2016, 2017a, b), Marchand et al. (2018), Bernard et al. (2018) using the Fokker (1929) Lagrangian approach in harmonic coordinates. This second derivation was actually the first one to derive the result from scratch without ambiguity parameter and without resorting to an auxiliary GSF calculation. This was
achieved by using the dimensional regularization for both the UV and IR type divergences.<fnref n="16" /> The invariants of the motion (including periastron advance) and radiation-reaction terms were given by Bernard et al. (2018);

3. The effective field theory (EFT) approach also obtained the complete, first-principle and ambiguity-free results at 4PN order (Foffa and Sturani 2012, 2013; Galley et al. 2016; Foffa et al. 2017; Porto and Rothstein 2017; Foffa and Sturani 2019; Foffa et al. 2019; Blümlein et al. 2020a).

Much efforts have been spent trying to extend the validity of the equations of motion beyond the 4PN order. The EFT approach in harmonic coordinates has been pushed to 5PN and even 6PN orders (including tails and hereditary effects) using *ab initio* brute force calculations (Blümlein et al. 2020b, 2021, 2022), and, in particular, by exploiting some property of factorisation of Feynman diagrams (Foffa and Sturani 2020, 2021; Foffa et al. 2021; Almeida et al. 2023). On the other hand a new methodology for deriving the equations of motion was introduced by Bini et al. (2019), combining several formalisms: post-Newtonian, post-Minkowskian, multi-polar-post-Minkowskian, gravitational self-force, and effective one-body. This method (coined as “Tutti-Frutti”) has been applied to the rederivation of the 3PN conservative Hamiltonian (Bern et al. 2019b) and to the computation of new coefficients in the equations of motion at 5PN and 6PN orders (Bini et al. 2020a, b, c, 2021). Unfortunately, despite these efforts the results are still incomplete: The 5PN Hamiltonian is determined except for two unknown numerical coefficients, in front of terms proportional to the square of the symmetric mass ratio $v$ and to the fifth and sixth power of the gravitational constant $G$. Such undetermined terms $\propto G^6 v^2$ and $G^5 v^2$ thus belong to a high post-Minkowskian approximation, and at the same time cannot be determined from the comparison with linear GSF; one would need at least the second-order GSF which is the current challenge in self-force calculations (Pound et al. 2020; Warburton et al. 2021; Wardell et al. 2023). Similarly, at the 6PN order the Hamiltonian contains four unknown numerical coefficients in factor of terms $\propto G^7 v^3$, $G^7 v^2$, $G^6 v^2$ and $G^5 v^2$. More work should be done to compute these coefficients.

Most of the works reviewed in this article concern general relativity. However, let us mention that the equations of motion of compact binaries in scalar-tensor theories have been developed up to 3PN order (Mirshekari and Will 2013; Bernard 2018, 2019).

An important body of works in GR concerns the effects of spins on the equations of motion of compact bodies. The leading spin-orbit (SO) effect arises at the 1.5PN order while the leading spin-spin (SS) effect appears at 2PN order (see Sect. 3.6). They have been computed using the traditional PN approach by Kidder et al. (1993b), Kidder (1995). The next-to-leading SO effect, i.e., 1PN relative order corresponding to 2.5PN order, was obtained by Tagoshi et al. (2001), Faye et al. (2006). The results were also retrieved by two subsequent calculations, one using the

<footnote n="16">¹⁶ In fact a specific regularization procedure (called the “*Bë*” regularization) was used for the IR divergences, proven to be equivalent to dimensional regularization (Marchand et al. 2018). We comment more on this procedure in Sects. 3.1.3 and 3.1.4.</footnote>
ADM Hamiltonian (Damour et al. 2008) and the other with the EFT method (Levi 2010a; Porto 2010). The ADM calculation was generalized by Hartung and Steinhoff (2011a) to the N-body case and extended by Steinhoff et al. (2008a, b, c, ), Hergt et al. (2010), Hartung and Steinhoff (2011c) to the next-to-leading spin-spin effects, including both the coupling between different spins and spin square terms, and the next-to-next-to-leading SS interactions between different spins at the 4PN order was obtained by Hartung and Steinhoff (2011c). Using the EFT method the next-to-leading 3PN SS and spin-squared contributions were derived by Porto and Rothstein (2006, 2008a, b), Levi (2010b), Levi and Steinhoff (2015a), and the next-to-next-to-leading 4PN SS interactions for different spins by Levi (2012) and for spin-squared terms by Levi and Steinhoff (2016). The next-to-next-to-leading order SO effects, corresponding to 3.5PN order equivalent to 2PN relative order, were obtained in the ADM-coordinates Hamiltonian by Hartung and Steinhoff (2011b), Hartung et al. (2013) and in the harmonic-coordinates equations of motion by Marsat et al. (2013b), Bohé et al. (2013b), with complete equivalence between the results. Comparisons between the EFT and ADM Hamiltonian schemes for high-order SO and SS couplings can be found in Levi and Steinhoff (2014, 2015a, b). The next-to-next-to-next-to-leading order (NNNLO) has been worked out by Antonelli et al. (2020b) (SO) and Antonelli et al. (2020a) (SS for different spins) with traditional method, and by Kim et al. (2023b), Mandal et al. (2023b) (SO) and Kim et al. (2023a), Mandal et al. (2023a) (SS) from the EFT.

### 1.2.2 Post-Newtonian gravitational radiation

The problem of the computation of the gravitational waveform (and the energy flux) is solved by application of a wave generation formalism valid for general isolated matter systems (see Sect. 1.1.1). The earliest computation at the 1PN level beyond the quadrupole moment formalism was achieved by Wagoner and Will (1976),<fnref n="17" /> based on the Epstein and Wagoner (1975) multipole moments, applied to compact binaries moving on eccentric orbits. This 1PN level calculation was redone and confirmed by Blanchet and Schäfer (1989) using the Blanchet and Damour (1986) multipole moments.

At the 1.5PN order in the radiation field (beyond the quadrupole formula), appears the first “hereditary” contribution, which is a priori sensitive to the entire past history of the source, i.e., depends on the source at all previous times up to $-\infty$ in the past up to current time (Blanchet and Damour 1988, 1992). At 1.5PN order this hereditary term represents the dominant contribution of gravitational-wave tails in the wave zone—a quadratic coupling between the quadrupole moment and the static mass of the source. It has been evaluated for compact binaries by Poisson (1993) in the small mass ratio limit, and in the general case by Wiseman (1993), Blanchet and Schäfer (1993). The 1.5PN tail term (and the radiation reaction contribution therein) is also known from Foffa and Sturani (2013), Galley et al. (2016).

<footnote n="17">¹⁷ At a time when PN corrections to the emission of gravitational waves had only a purely academic interest.</footnote>
Applying the general multipole moments at 2PN level (Blanchet 1995), the energy flux of compact binaries was completed to the 2PN order by Blanchet et al. (1995a), Gopakumar and Iyer (1997), and, independently, by Will and Wiseman (1996), Will (1999), using the DIRE formalism; see Blanchet et al. (1995b) for the joint report of these calculations. The waveform to 2PN order has been computed by Blanchet et al. (1996). The energy flux and waveform have also been computed to 2PN order using the EFT approach, by Goldberger and Ross (2010), Leibovich et al. (2020) with results agreeing with those of “traditional” methods.

Higher-order tail effects arise at the 2.5PN and have been incorporated in the 2.5PN wave generation formalism (Blanchet 1996). In the waveform at this order appears the leading contribution of a different hereditary effect called the “non-linear memory” (Blanchet 1990; Christodoulou 1991; Wiseman and Will 1991; Thorne 1992; Blanchet and Damour 1992; Blanchet 1998c). The non-linearity is interpreted as due to the gravitational wave emission of gravitons themselves, following formulas for the emission by massless particles (Braginsky and Thorne 1987a; Thorne 1992). The non-linear memory takes the form of a simple anti-derivative of an “instantaneous” term, and therefore becomes instantaneous in the energy flux, involving the time derivative of the waveform. In principle the memory contribution must be computed using some model for the evolution of the binary system in the past. Because of the cumulative effect of integration over the whole past, the memory term, though originating from 2.5PN order, finally contributes in the waveform at Newtonian order (Wiseman and Will 1991; Arun et al. 2004). It represents a part of the waveform whose amplitude grows on a radiation reaction time scale, but is nearly constant over one orbital period; it is therefore essentially a *zero-frequency* effect (or DC effect). The non-linear memory in the waveform of inspiralling compact binaries has been computed to high post-Newtonian order by Favata (2009, 2011c). We give more details on this effect in Sect. 2.4.4.

At the 3PN order appears the first cubic process, which is due to tails generated by tails themselves—the so-called “tails of tails” (Blanchet 1998a; Almeida et al. 2021). The 3PN approximation also involves, besides the tails of tails, many non-tail contributions coming from the relativistic corrections in the (source type) multipole moments of the compact binary. The mass quadrupole moment at 3PN order has been obtained by Blanchet et al. (2002b), Blanchet and Iyer (2005), Blanchet et al. (2005a). This calculation crucially requires dimensional regularization for treating UV divergences, see Sect. 3.1.3. The wave generation from compact binaries is then complete to 3.5PN order (Blanchet et al. 2002a, 2004b). The formalism is valid for general orbits, but often for applications to compact binaries the results are presented for quasi-circular orbits. In the case of non-circular orbits the 3PN energy and angular momentum fluxes, and the associated balance equations, are known (Arun et al. 2008a, b, 2009a; Loutrel and Yunes 2017), as well as the gravitational-wave modes, taking into account memory and tail contributions (Ebersold et al. 2019; Boetzel et al. 2019). We discuss the problem of eccentric orbits in Sect. 3.5.

In recent years the gravitational wave generation by compact binaries was developed up to still higher 4.5PN order. A central part of this program is the control of the mass quadrupole moment of compact binaries (without spins) with 4PN precision. In a preliminary calculation at 4PN order by Marchand et al. (2020), the
UV divergences were properly treated by means of dimensional regularization, but the IR ones were regularized with the Hadamard partie finie regularization. This calculation was then completed by Larrouturou et al. (2022a, b) with: (i) the contributions from the dimensional regularization of the IR divergences; (ii) the non-local (in time) effect in the source quadrupole moment, due to the radiation modes associated with propagating tails at infinity. The latter effect is the analogue of the 4PN tail term in the conservative equations of motion and Lagrangian, see Eq. (348). The 3PN mass octupole moment and 3PN current quadrupole moment were also needed and derived by Faye et al. (2015), Henry et al. (2021). By contrast, in scalar-tensor theory the waveform and flux have been so far developed to much lower order, see Lang (2014, 2015), Sennett et al. (2016), Bernard et al. (2022).

Interestingly, at the 4PN order in GR there is a new type of cubic effect, corresponding to the coupling between the static mass of the source and two varying quadrupole moments. Physically it can be seen as a combination between the tails produced by the non-linear memory and the memory associated with the tail; thus we coin this effect “tails-of-memory”. A thorough computation of the cubic tails of memory was performed by Trestini and Blanchet (2023) by means of the MPM formalism, together with a simpler cubic interaction between mass, spin and quadrupole moment, also arising at 4PN order. In particular it was shown that the tail of memory contains a DC contribution in perfect agreement with the general formalism for memory in Sect. 2.4.4. An essential feature of this calculation is that it uses the MPM algorithm in radiative, non-harmonic coordinates as described in Sect. 2.1.5. A minor price to be paid is that the multipole moments in the radiative MPM construction are to be properly adjusted with respect to those in the corresponding harmonic MPM algorithm (Trestini et al. 2023). Finally at the 4.5PN order there is the first quartic non-linear multipole interaction, between three static masses and the quadrupole; this effect, rightly called “tail-of-tail-of-tail” has been computed by Marchand et al. (2016), Messina and Nagar (2017).

These calculations lead to the complete 4.5PN energy flux for circular orbits, and the corresponding 4.5PN gravitational-wave frequency and phase evolution, reported in Sects. 3.4.3 and 3.4.4 following Blanchet et al. (2023b, c). Concerning the gravitational-wave modes (reported in Sect. 3.4.5), the dominant quadrupolar mode $(\ell, m) = (2, 2)$ is derived at 4PN order while, all the other modes are computed to 3.5PN order in the waveform (Blanchet et al. 2008; Henry 2023).

The post-Newtonian results for the waveform and energy flux are in complete agreement (up to the 4.5PN order) with the results given by the very different technique of linear black-hole perturbations, valid when the mass of one of the bodies is small compared to the other. Linear black-hole perturbations in the radiation field, triggered by the geodesic motion of the small particle around the black hole, have been applied to this problem by Poisson (1993) at the 1.5PN order (following pioneering work by Gal'tsov et al. 1980), by Tagoshi and Nakamura (1994) using a numerical code up to the 4PN order, and by Sasaki (1994), Tagoshi and Sasaki (1994), Tanaka et al. (1996) (see also Mino et al. 1997a), analytically up to the 5.5PN order. They have been generalized to a rotating black hole by Tagoshi et al. (1996). More recently the method has been improved and extended up to extremely high 22PN order by Fujita (2012a, b) using the Mano et al. (1996a, b) method for
analytical solution of the homogeneous Teukolsky equation—but still for linear black-hole perturbations. The polarization modes and associated factorized resummed waveforms for circular orbits are known to 5.5PN order (Fujita and Iyer 2010). For eccentric orbits (still in the small mass ratio limit) the energy and angular momentum fluxes are known to 9PN order (Forseth et al. 2016; Munna et al. 2020).

The spins (SO and SS) affect the gravitational waves through a modulation of their amplitude, phase and frequency. The orbital plane precesses in the case where the spins are not aligned or anti-aligned with the orbital angular momentum (Apostolatos et al. 1994). The leading SO and SS contributions in the waveform and flux of compact binaries are known from Kidder et al. (1993b), Kidder (1995), Mikóczi et al. (2005); the next-to-leading SO terms at order 2.5PN were obtained by Blanchet et al. (2006) after a previous attempt by Owen et al. (1998); the 3PN SO contribution is due to tails and was computed by Blanchet et al. (2011), after intermediate results at the same order (but including SS terms) given by Porto et al. (2011). The next-to-next-to-leading SO contributions in the multipole moments and the energy flux, corresponding to 3.5PN order, and the next-to-leading SO tail corresponding to 4PN order, have been obtained by Bohé et al. (2013a), Marsat et al. (2013a), Cho et al. (2022a). The next-to-leading 3PN SS and spin-squared contributions in the radiation field were derived by Bohé et al. (2015), Cho et al. (2021). The SO and SS effects in the spherical harmonic modes for non-circular orbits with non-precessing components have been computed by Paul and Mishra (2023) through 2PN order. In Sect. 3.6 we shall give the results for the SO contributions up to the 4PN order, and for the SS contributions up to the 3PN order, in the energy flux and gravitational wave phase evolution.

# 2 Post-Newtonian isolated sources

## 2.1 Non-linear iteration of the field equations outside the source

### 2.1.1 Einstein's field equations

The field equations of general relativity are obtained by varying the action

$$S = \frac{c^3}{16\pi G} \int d^4 x \sqrt{-g} R + S_m[\Psi, g_{\alpha\beta}], \qquad (46)$$

with respect to the space-time metric $g_{\alpha\beta}$, where the first term is the famous Einstein-Hilbert action for the gravitational field, and the second term is the matter action. The field equations form a system of ten second-order partial differential equations obeyed by the metric,

$$G^{\alpha\beta}[g, \partial_g, \partial^2 g] = \frac{8\pi G}{c^4} T^{\alpha\beta}[\Psi, g], \qquad (47)$$

where the Einstein curvature tensor $G^{\alpha\beta} \equiv R^{\alpha\beta} - \frac{1}{2}Rg^{\alpha\beta}$ is generated, through the gravitational coupling constant $\kappa = 8\pi G/c^4$, by the stress-energy tensor $T^{\alpha\beta} \equiv$
$\frac{2}{\sqrt{-g}} \delta S_m / \delta g_{\alpha\beta}$ of the matter fields $\Psi$. Among these ten equations, four govern, via the contracted Bianchi identity, the evolution of the matter system, i.e. the covariant conservation of the matter stress-energy tensor,

$$ \nabla_{\mu} G^{\alpha\mu} = 0 \quad \Rightarrow \quad \nabla_{\mu} T^{\alpha\mu} = 0. \qquad (48) $$

The matter field equations can be obtained by varying the matter action $S_m$ in (46) with respect to each of the matter fields $\Psi$. Satisfying the matter field equations implies the conservation law (48) for the matter stress-energy tensor. The space-time geometry is constrained by the six remaining equations, which place six independent constraints on the ten components of the metric $g_{\alpha\beta}$, leaving four of them to be fixed by a choice of the coordinate system.

A widely adopted choice of coordinates are the so-called *harmonic coordinates*, also known as De Donder (1921) coordinates. We define, as a basic variable, the gravitational-field amplitude (or “gothic” metric deviation from Minkowski space-time) as

$$ h^{\alpha\beta} \equiv \sqrt{-g} g^{\alpha\beta} - \eta^{\alpha\beta}, \qquad (49) $$

where $g^{\alpha\beta}$ denotes the contravariant metric (satisfying $g^{\alpha\mu}g_{\mu\beta} = \delta^2_{\beta}$), where $g$ is the determinant of the covariant metric, $g \equiv \det(g_{\alpha\beta})$, and where $\eta^{\alpha\beta}$ represents an auxiliary Minkowskian metric $\eta^{\alpha\beta} \equiv \text{diag}(-1, 1, 1, 1)$. The harmonic-coordinate condition, which accounts exactly for the four equations (48) corresponding to the conservation of the matter tensor, reads,<fnref n="18" />

$$ H^{\alpha} \equiv \partial_{\mu} h^{\alpha\mu} = 0. \qquad (50) $$

Equation (50) introduces into the definition of our coordinate system a preferred Minkowskian structure, with Minkowski metric $\eta_{\alpha\beta}$. Of course, this is not contrary to the spirit of general relativity, where there is only one physical metric $g_{\alpha\beta}$ without any flat prior geometry, because the coordinates are not governed by geometry (so to speak), but rather can be chosen at convenience, depending on physical phenomena under study. The coordinate condition (50) is especially useful when studying gravitational waves as perturbations of space-time propagating on the fixed background metric $\eta_{\alpha\beta}$. This view is perfectly legitimate and represents a fruitful and rigorous way to think of the problem using approximation methods. Indeed, the metric $\eta_{\alpha\beta}$, originally introduced in the coordinate condition (50), does exist at any finite order of approximation (neglecting higher-order terms), and plays the role of some physical “prior” flat geometry at any order of approximation. The harmonic

$$ H^\alpha = \sqrt{-g} \Box_g x^\alpha, $$

<footnote n="18">¹⁸ Considering the coordinates $x^\alpha$ as a set of four scalars a simple calculation shows that</footnote>

<footnote>where $\Box_g \equiv g^{\mu\nu} \nabla_\mu \nabla_\nu$ denotes the curved d'Alembertian operator. Hence the harmonic-coordinate condition $H^\alpha = 0$ tells that the coordinates $x^\alpha$ themselves, considered as scalars, are harmonic, i.e., obey the vacuum (curved) d'Alembertian equation.</footnote>
gauge condition will be used for many investigations in this article<fnref n="19" />; however for some applications it should be avoided, see in particular Sect. 2.1.5.

The Einstein field equations for the gothic metric deviation (49) in a general coordinate system can be written in the following form,

$$ \Box h^{\alpha\beta} - \partial_{\alpha}H^{\beta} = \frac{16\pi G}{c^4} \tau^{\alpha\beta}, \qquad (51) $$

where $\Box \equiv \Box_{\eta} = \eta^{\mu\nu}\partial_{\mu}\partial_{\nu}$ is the flat d'Alembertian operator, and we have posed as useful shorthand (where $H^{\alpha} = \partial_{\alpha}h^{\alpha\mu}$ is non-zero in general)

$$ \partial H^{\alpha\beta} \equiv \partial^{\alpha} H^{\beta} + \partial^{\beta} H^{\alpha} - \eta^{\alpha\beta} \partial_{\mu} H^{\mu}. \qquad (52) $$

The source term $\tau^{\alpha\beta}$ can rightly be interpreted as the stress-energy pseudo-tensor (actually, $\tau^{\alpha\beta}$ is a Lorentz-covariant tensor) of the matter fields, described by $T^{\alpha\beta}$, and the gravitational field, given by the gravitational source term $\Lambda^{\alpha\beta}$, i.e.,

$$ \tau^{\alpha\beta} = |g| T^{\alpha\beta} + \frac{c^4}{16\pi G} \Lambda^{\alpha\beta}. \qquad (53) $$

The pseudo-tensor is conserved in an ordinary sense as a consequence of the Einstein field equations (51), and this in turn is equivalent to the matter equations of motion (48). Applying an ordinary divergence on (51), and using the fact that $\partial_{\mu}\partial^{\mu}H^{\alpha\mu} = \Box H^{\alpha}$, we obtain

$$ \partial_{\mu} \tau^{\alpha\mu} = 0 \iff \nabla_{\mu} T^{\alpha\mu} = 0. \qquad (54) $$

The exact expression of $\Lambda^{\alpha\beta}$ reads

$$ \Lambda^{\alpha\beta} = h^{\mu\alpha}\partial_{\mu}H^{\beta} + h^{\mu\beta}\partial_{\mu}H^{\alpha} - \partial_{\mu}(h^{\alpha\beta}H^{\mu}) + \Lambda_{\text{harm}}^{\alpha\beta}, \qquad (55) $$

where $\Lambda_{\text{harm}}^{\alpha\beta}$ is the source term when assuming the harmonic gauge condition (50), and given, including all non-linearities, by<fnref n="20" />

$$ \begin{aligned} \Lambda_{\text{harm}}^{\alpha\beta} = & -h^{\mu\nu}\partial_{\mu}^{2}h^{\alpha\beta} + \partial_{\mu}h^{\alpha\nu}\partial_{\nu}h^{\beta\mu} + \frac{1}{2}g^{\alpha\beta}g_{\mu\nu}\partial_{\lambda}h^{\mu\tau}\partial_{\tau}h^{\nu\lambda} \\ & -g^{\alpha\mu}g_{\nu\tau}\partial_{\lambda}h^{\beta\tau}\partial_{\mu}h^{\nu\lambda} - g^{\beta\mu}g_{\nu\tau}\partial_{\lambda}h^{\tau\kappa}\partial_{\mu}h^{\nu\lambda} + g_{\mu\nu}g^{\lambda\tau}\partial_{\lambda}h^{\kappa\mu}\partial_{\tau}h^{\beta\nu} \\ & + \frac{1}{8}(2g^{\alpha\mu}g^{\beta\nu} - g^{\alpha\beta}g^{\mu\nu})(2g_{\lambda\tau}g_{\epsilon\pi} - g_{\tau\epsilon}g_{\lambda\pi})\partial_{\mu}h^{\lambda\pi}\partial_{\nu}h^{\tau\epsilon}. \end{aligned} \qquad (56) $$

As is clear from its expression, $\Lambda^{\alpha\beta}$ is made of terms at least quadratic in the gravitational-field strength $h$ and its first and second space-time derivatives, while the left-hand side of (51) is linear in $h$. It looks hawkward, from a mathematical point of view, notably for setting the initial value problem, to include second order space-time (and in particular time) derivatives in the right-hand side of the differential equation

<footnote n="19">¹⁹ It is conveniently imposed by adding a gauge-fixing term in the Einstein-Hilbert action, see e.g. Damour and Esposito-Farèse (1996), Foffa and Sturani (2014) and Sect. 3.2.1 below.</footnote>

<footnote n="20">²⁰ In $D = d + 1$ space-time dimensions, only one coefficient in this expression is modified; see Eq. (263) below.</footnote>
(51). But of course this is because these second time derivatives appear in terms quadratic in $h$ and can be treated perturbatively in the approximation scheme. In the following, for the highest post-Newtonian order that we shall consider, we will need the quadratic, cubic and quartic pieces of $\Lambda^{\alpha\beta}$; with obvious notation, we can write them as

$$ \Lambda_{\text{harm}}^{\alpha\beta} = N^{\alpha\beta}[h, h] + M^{\alpha\beta}[h, h, h] + \mathcal{O}(h^4). \quad (57) $$

These various terms can be straightforwardly computed from expanding Eq. (56); The leading quadratic and cubic pieces are explicitly given by<fnref n="21" />

$$ \begin{aligned} N^{\alpha\beta} = & -h^{\mu\nu} \partial_{\mu\nu}^2 h^{\alpha\beta} + \frac{1}{2} \partial^{\alpha} h_{\mu\nu} \partial^{\beta} h^{\mu\nu} - \frac{1}{4} \partial^{\alpha} h \partial^{\beta} h + \partial_{\nu} h^{\alpha\mu} (\partial^{\nu} h_{\mu}^{\beta} + \partial_{\mu} h^{\beta\nu}) \\ & - 2 \partial^{(\alpha} h_{\mu\nu} \partial^{\mu} h^{\beta)\nu} + \eta^{\alpha\beta} \left[ -\frac{1}{4} \partial_{\tau} h_{\mu\nu} \partial^{\tau} h^{\mu\nu} + \frac{1}{8} \partial_{\mu} h \partial^{\mu} h + \frac{1}{2} \partial_{\mu} h_{\nu\tau} \partial^{\nu} h^{\mu\tau} \right]. \end{aligned} \quad (58a) $$

$$ \begin{aligned} M^{\alpha\beta} = & -h^{\mu\nu} \left( \partial^{\alpha} h_{\mu\tau} \partial^{\beta} h_{\nu}^{\tau} + \partial_{\tau} h_{\mu\tau}^{\alpha} \partial^{\tau} h_{\nu}^{\beta} - \partial_{\mu} h_{\tau}^{\alpha} \partial_{\nu} h^{\beta\tau} \right) \\ & + h^{\alpha\beta} \left[ -\frac{1}{4} \partial_{\tau} h_{\mu\nu} \partial^{\tau} h^{\mu\nu} + \frac{1}{8} \partial_{\mu} h \partial^{\mu} h + \frac{1}{2} \partial_{\mu} h_{\nu\tau} \partial^{\nu} h^{\mu\tau} \right] + \frac{1}{2} h^{\mu\nu} \partial^{(\alpha} h_{\mu\nu} \partial^{\beta)} h \\ & + 2 h^{\mu\nu} \partial_{\tau} h_{\mu}^{(\alpha} \partial^{\beta)} h_{\nu}^{\tau} + h^{\mu(\alpha} \left[ \partial^{\beta)} h_{\nu\tau} \partial_{\mu} h^{\nu\tau} - 2 \partial_{\nu} h_{\tau}^{\beta)} \partial_{\mu} h^{\nu\tau} - \frac{1}{2} \partial^{\beta)} h \partial_{\mu} h \right] \\ & + \eta^{\alpha\beta} \left[ \frac{1}{8} h^{\mu\nu} \partial_{\mu} h \partial_{\nu} h - \frac{1}{4} h^{\mu\nu} \partial_{\tau} h_{\mu\nu} \partial^{\tau} h - \frac{1}{4} h^{\tau\lambda} \partial_{\tau} h_{\mu\nu} \partial_{\lambda} h^{\mu\nu} \right. \\ & \left. - \frac{1}{2} h^{\tau\lambda} \partial_{\mu} h_{\tau\nu} \partial^{\nu} h_{\lambda}^{\mu} + \frac{1}{2} h^{\tau\lambda} \partial_{\mu} h_{\tau}^{\nu} \partial^{\mu} h_{\lambda\nu} \right], \end{aligned} \quad (58b) $$

In this article, we shall look for approximate solutions of the field equations (50)–(51) under the following four hypotheses:

1. The matter stress-energy tensor $T^{\alpha\beta}$ is of spatially compact support, i.e., can be enclosed into some time-like world tube, say $r \le a$, where $r = |x|$ is the coordinate radial distance. Outside the domain of the source, when $r > a$, the gravitational source term, according to (54), is divergence-free,

$$ \partial_{\mu} \Lambda^{\alpha\mu} = 0 \quad (\text{when } r > a); \quad (59) $$

2. The matter distribution inside the source is smooth: $T^{\alpha\beta} \in C^{\infty}(\mathbb{R}^3)$.<fnref n="22" /> We have in mind a smooth hydrodynamical fluid system, without any singularities nor

<footnote n="21">²¹ See Eqs. (3.8c) in Blanchet and Faye (2001a) for the quartic terms. We denote e.g., $h_{\mu}^{\alpha} = \eta_{\mu\nu} h^{\nu\alpha}$, $h = \eta_{\mu\nu} h^{\mu\nu}$, and $\partial^{\alpha} = \eta^{\alpha\mu} \partial_{\mu}$. A parenthesis around a pair of indices denotes the usual symmetrization: $T^{(\alpha\beta)} = \frac{1}{2} (T^{\alpha\beta} + T^{\beta\alpha})$.</footnote>

<footnote n="22">²² $\mathbb{N}$, $\mathbb{Z}$, $\mathbb{R}$, and $\mathbb{C}$ are the usual sets of non-negative integers, integers, real numbers, and complex numbers; $C^p(\Omega)$ is the set of $p$-times continuously differentiable functions on the open domain $\Omega$ ($p \le +\infty$).</footnote>
shocks (*a priori*), that is described by some Euler-type equations including high relativistic corrections. In particular, we *a priori* exclude the presence of any black hole singularities; however, we shall return to this question in Sect. 3 when we look for a model describing compact objects;

3. The source is post-Newtonian in the sense of the existence of the small parameter defined by Eq. (1). For such a source we assume the legitimacy of the method of matched asymptotic expansions for identifying the inner post-Newtonian field and the outer multipolar decomposition in the source's exterior near zone;

4. The gravitational field has been independent of time (stationary) in some remote past, i.e., before some finite instant $-T$ in the past, namely

$$ \frac{\partial}{\partial t} [h^{\alpha\beta}(\mathbf{x}, t)] = 0 \quad \text{when} \quad t \le -T. \qquad (60) $$

The latter condition is a means to impose, by “brute force”, the *no-incoming* radiation condition (Damour 1983a), ensuring that the matter source is isolated from the rest of the Universe and does not receive any radiation from infinity. Ideally, the no-incoming radiation condition should be imposed at past null infinity. As we shall see, this condition entirely fixes the radiation reaction forces inside the isolated source. We shall later argue (in Sect. 2.4.2) that the condition of stationarity in the past (60), although weaker than the ideal no-incoming radiation condition, does not entail essential physical restriction on the general validity of the formulas we derive. Furthermore, the condition (60) is actually better suited in the case of real astrophysical sources like inspiralling compact binaries, for which we do not know the details of the initial formation and remote past evolution. In practice the initial instant $-T$ can be set when the two supernovæ of the progenitor stars yielded two compact objects, perhaps right after the “common enveloppe” (CE) phase starting the close evolution of the compact binary system.

### 2.1.2 Linearized vacuum equations

In what follows we shall solve the field equations (51), in the *vacuum* region outside the compact-support source, in the form of a formal non-linearity or *post-Minkowskian* expansion, considering the field variable $h^{\alpha\beta}$ as a non-linear metric perturbation of Minkowski space-time, and further assuming formal multipole expansions. At the linearized level (or first-post-Minkowskian approximation), we write:

$$ h_{\text{MPM}}^{\alpha\beta} = G h_{(1)}^{\alpha\beta} + \mathcal{O}(G^2), \qquad (61) $$

where the subscript MPM refers to the formal multipolar-post-Minkowskian expansion, which is physically meaningful only in the exterior of the source, and where we have introduced Newton's constant $G$ as a book-keeping parameter, enabling one to label conveniently the successive post-Minkowskian approximations. Since $h^{\alpha\beta}$ is a dimensionless variable, with our convention the linear coefficient $h_{(1)}^{\alpha\beta}$
in Eq. (61) has the dimension of the inverse of $G$ (which should be a mass squared in a system of units where $\hbar = c = 1$). In vacuum, further assuming the harmonic coordinate condition (50), the linearized metric coefficient $h_{(1)}^{\alpha\beta}$ satisfies

$$ \Box h_{(1)}^{\alpha\beta} = 0, \qquad (62a) $$

$$ \partial_{\mu} h_{(1)}^{\alpha\mu} = 0. \qquad (62b) $$

We want to solve those equations by means of an infinite multipolar series valid outside a time-like world tube containing the source. Indeed the multipole expansion is the appropriate method for describing the physics of the source as seen from its exterior ($r > a$). On the other hand, the post-Minkowskian series is physically valid in the weak-field region, which surely includes the exterior of any source, starting at a sufficiently large distance. For post-Newtonian sources the exterior weak-field region, where both multipole and post-Minkowskian expansions are valid, simply coincides with the exterior region $r > a$. It is therefore quite natural, and even, one would say inescapable when considering general sources, to combine the post-Minkowskian approximation with the multipole decomposition. This is the original idea of the “double-expansion” series of Bonnor (1959), Bonnor and Rotenberg (1961), Bonnor and Rotenberg (1966), Hunter and Rotenberg (1969), which combines the $G$-expansion (or $m$-expansion in their notation) with the $a$-expansion (equivalent to the multipole expansion, since the $\ell$-th order multipole moment scales with the source radius like $a^\ell$).

The multipolar-post-Minkowskian (MPM) method will be implemented systematically, using symmetric-trace-free (STF) harmonics to describe the multipole expansion (Thorne 1980), and looking for a definite “algorithm” for implementing the approximation scheme (Blanchet and Damour 1986). The solution of the system of equations (62) takes the form of a series of retarded multipolar waves<fnref n="23" />

<footnote n="23">²³ Our notation is the following: $L = i_1 i_2 \cdots i_\ell$ denotes a multi-index, made of $\ell$ (spatial) indices. Similarly, we write for instance $P = j_1 \cdots j_p$ (in practice, we generally do not need to write explicitly the “carrier” letter $i$ or $j$), or $aL - 1 = a i_1 \cdots i_{\ell-1}$. Always understood in expressions such as Eq. (63) are $\ell$ summations over the indices $i_1, \dots, i_\ell$ ranging from 1 to 3. The derivative operator $\partial_L$ is a short-hand for $\partial_j \cdots \partial_\ell$. The function $K_L$ (for any space-time indices $\alpha\beta$) is symmetric and trace-free (STF) with respect to the $\ell$ indices composing $L$. This means that for any pair of indices $i_p, i_q \in L$, we have $K_{\cdots i_p \cdots i_q} = K_{\cdots i_q \cdots i_p}$ and that $\delta_{i_p i_q} K_{\cdots i_p \cdots i_q} = 0$ (see Thorne 1980 and Appendices A and B in Blanchet and Damour 1986 for reviews about the STF formalism). The STF projection is denoted with a hat, so $K_L \equiv \hat{K}_L$, or sometimes with carets around the indices, $K_L \equiv K_{(L)}$. In particular, $\hat{n}_L = n_{(L)}$ is the STF projection of the product of unit vectors $n_L = n_{i_1} \cdots n_{i_\ell}$, for instance $\hat{n}_{ij} = n_{(ij)} = n_{ij} - \frac{1}{3}\delta_{ij}$ and $\hat{n}_{ijk} = n_{(ijk)} = n_{ijk} - \frac{1}{3}(\delta_{ij}n_k + \delta_{ik}n_j + \delta_{jk}n_i)$; an expansion into STF tensors $\hat{n}_L = \hat{n}_L(0, \phi)$ is equivalent to the usual expansion in spherical harmonics $Y_{lm} = Y_{lm}(0, \phi)$, see Eqs. (189) below. Similarly, we denote $x_L = x_{i_1} \cdots x_{i_\ell} = r^\ell n_L$ where $r = |x|$, and $\hat{x}_L = x_{(L)} = \text{STF}[x_L]$. The Levi-Civita antisymmetric symbol is denoted $\epsilon_{ijk}$ (with $\epsilon_{123} = 1$). Parenthesis refer to symmetrization, $T_{(ij)} = \frac{1}{2}(T_{ij} + T_{ji})$. Superscripts (q) indicate q successive time derivations.</footnote>
$$h_{(1)}^{\alpha\beta} = \sum_{\ell=0}^{+\infty} \hat{\partial}_L \left( \frac{K_L^{\alpha\beta}(t-r/c)}{r} \right), \qquad (63)$$

where $r = |\mathbf{x}|$, and where the functions $K_L^{\alpha\beta} \equiv K_{i_1 \dots i_i}^{\alpha\beta}$ are smooth functions of the retarded time $u \equiv t - r/c$ [i.e., $K_L(u) \in C^\infty(\mathbb{R})$], which become constant in the past, when $t \le -T$, see Eq. (60). Since a monopolar wave satisfies $\Box(K_L(u)/r) = 0$ and the d'Alembertian commutes with the multi-derivative $\hat{\partial}_L$, it is evident that Eq. (63) represents the most general solution of the wave equation (62a); but see Sect. 2 of Blanchet and Damour (1986) for a rigorous proof based on the Euler-Poisson-Darboux equation. The gauge condition (62b), however, is not fulfilled in general, and to satisfy it we must algebraically decompose the set of functions $K_L^{00}$, $K_L^{0i}$, $K_L^{ij}$ into ten tensors which are STF with respect to all their indices, including the spatial indices $i$, $ij$. Imposing the condition (62b) reduces the number of independent tensors to six, and we find that the solution takes an especially simple “canonical” form, parametrized by only two moments, plus some arbitrary linearized gauge transformation.

**Theorem 1** (Sachs and Bergmann 1958; Pirani 1965; Thorne 1980; Blanchet and Damour 1986) The most general solution of the linearized field equations (62) outside some time-like world tube enclosing the source ($r > a$), and stationary in the past [see Eq. (60)], reads

$$h_{(1)}^{\alpha\beta} = h_{\text{can}(1)}^{\alpha\beta} + \hat{\partial}_{(1)} \varphi_{(1)}^{\alpha\beta}. \qquad (64)$$

The first term is referred to as a “canonical” solution and depends on two STF-tensorial multipole moments, $I_L(u)$ and $J_L(u)$, which are arbitrary functions of time except for the laws of conservation of the monopole: $I = \text{const}$, and dipoles: $I_i = \text{const}$, $J_i = \text{const}$. It is given by

$$h_{\text{can}(1)}^{00} = -\frac{4}{c^2} \sum_{\ell \ge 0} \frac{(-)^\ell}{\ell!} \hat{\partial}_L \left( \frac{1}{r} I_L(u) \right), \qquad (65a)$$

$$h_{\text{can}(1)}^{0i} = \frac{4}{c^3} \sum_{\ell \ge 1} \frac{(-)^\ell}{\ell!} \left\{ \hat{\partial}_{L-1} \left( \frac{1}{r} I_{i(L-1)}^{(1)}(u) \right) + \frac{\ell}{\ell+1} \epsilon_{iab} \hat{\partial}_{aL-1} \left( \frac{1}{r} J_{bL-1}(u) \right) \right\}, \qquad (65b)$$

$$h_{\text{can}(1)}^{ij} = -\frac{4}{c^4} \sum_{\ell \ge 2} \frac{(-)^\ell}{\ell!} \left\{ \hat{\partial}_{L-2} \left( \frac{1}{r} I_{ijL-2}^{(2)}(u) \right) + \frac{2\ell}{\ell+1} \hat{\partial}_{aL-2} \left( \frac{1}{r} \epsilon_{ab(i} J_{j)bL-2}^{(1)}(u) \right) \right\}. \qquad (65c)$$

The second term, $\hat{\partial}_{(1)} \varphi_{(1)}^{\alpha\beta} \equiv \hat{\partial}_{(1)} \varphi_{(1)}^{\alpha\beta} + \hat{\partial}_{(1)} \varphi_{(1)}^{\alpha} - \eta^{\alpha\beta} \hat{\partial}_{(1)} \varphi_{(1)}^{\mu}$, represents a linearized gauge transformation, with gauge vector $\varphi_{(1)}^{\alpha}$ parametrized by four other multipole moments, say $W_L(u)$, $X_L(u)$, $Y_L(u)$ and $Z_L(u)$ [see Eqs. (66)].

The conservation of the lowest-order moments gives the constancy of the total mass of the source, $M \equiv I = \text{const}$, center-of-mass position, $X_i \equiv I_i/I = \text{const}$, total
linear momentum $P_i \equiv I_i^{(1)} = 0,^{24}$ and total angular momentum, $J_i = \text{const.}$ It is always possible to achieve $X_i = 0$ by translating the origin of our coordinates to the center of mass. The total mass M is the ADM mass of the Hamiltonian formulation of general relativity. Note that the quantities M, $X_i$, $P_i$ and $J_i$ include the contributions due to the waves emitted by the source. They describe the initial state of the source, before the emission of gravitational waves.

The multipole functions $I_L(u)$ and $J_L(u)$, which thoroughly encode the physical properties of the source at the linearized level (because the other moments $W_L, \dots, Z_L$ parametrize a gauge transformation), will be referred to as the *mass-type* and *current-type* source multipole moments. Beware, however, that at this stage the moments are not specified in terms of the stress-energy tensor $T^{\alpha\beta}$ of the source: Theorem 1 follows merely from the algebraic and differential properties of the vacuum field equations outside the source.

For completeness, we give the expressions of the components of the gauge-vector $\varphi^{\alpha}_{(1)}$ entering Eq. (64):

$$ \varphi_{(1)}^{0} = \frac{4}{c^3} \sum_{\ell \ge 0} \frac{(-)^{\ell}}{\ell!} \partial_L \left( \frac{1}{r} W_L(u) \right), \qquad (66a) $$

$$ \varphi_{(1)}^{i} = -\frac{4}{c^4} \sum_{\ell \ge 0} \frac{(-)^{\ell}}{\ell!} \partial_{iL} \left( \frac{1}{r} X_L(u) \right) \\ - \frac{4}{c^4} \sum_{\ell \ge 1} \frac{(-)^{\ell}}{\ell!} \left\{ \partial_{L-1} \left( \frac{1}{r} Y_{iL-1}(u) \right) + \frac{\ell}{\ell+1} \epsilon_{iab} \partial_{aL-1} \left( \frac{1}{r} Z_{bL-1}(u) \right) \right\}. \qquad (66b) $$

Because the theory is covariant with respect to non-linear diffeomorphisms and not merely with respect to linear gauge transformations, the moments $W_L, \dots, Z_L$ do play a physical role starting at the non-linear level, in the following sense. If one takes these moments equal to zero and continues the post-Minkowskian iteration [see Sect. 2.1.3] one ends up with a metric depending on $I_L$ and $J_L$ only, but that metric will not describe the same physical source as the one which would have been constructed starting from the six moments $I_L, J_L, \dots, Z_L$ altogether. In other words, the two non-linear metrics associated with the sets of multipole moments $\{I_L, J_L, 0, \dots, 0\}$ and $\{I_L, J_L, W_L, \dots, Z_L\}$ are not physically equivalent—they are not isometric. We shall point out in Sect. 2.1.4 below that the full set of moments $\{I_L, J_L, W_L, \dots, Z_L\}$ is in fact physically equivalent to some other reduced set of moments $\{M_L, S_L, 0, \dots, 0\}$, but with some moments $M_L, S_L$ that differ from $I_L, J_L$ by non-linear corrections, see e.g. Equation (251). The moments $M_L, S_L$ are called “canonical” moments; they play an important role when describing the radiation emitted by the isolated source. All the multipole moments $I_L, J_L, W_L, X_L, Y_L, Z_L$ will be computed in Sect. 2.2.3.

<footnote n="24">24 The constancy of the center of mass $X_i$—rather than a linear variation with time—results from our assumption of stationarity before the date $-T$, see Eq. (60). Hence, $P_i = 0$.</footnote>
### 2.1.3 The multipolar post-Minkowskian solution in harmonic coordinates

By Theorem 1 we know the most general solution of the linearized equations in the exterior of the source. We then tackle the problem of the post-Minkowskian iteration of that solution. We consider the full post-Minkowskian series

$$h_{\text{MPM}}^{\alpha\beta} = \sum_{n=1}^{+\infty} G^n h_{(n)}^{\alpha\beta}, \qquad (67)$$

where the first term is composed of the result given by Eqs. (64)–(66). In this article, we shall always understand the infinite sums such as the one in Eq. (67) in the sense of formal power series, i.e., as an ordered collection of coefficients, $(h_{(n)}^{\alpha\beta})_{n \in \mathbb{N}}$. We do not attempt to control the mathematical nature of the series and refer to the mathematical-physics literature for discussion in the present context (Christodoulou and Schmidt 1979; Damour and Schmidt 1990; Rendall 1990, 1992, 1994).

We substitute the post-Minkowski ansatz (67) into the vacuum harmonic coordinate Einstein field equations (50)–(51), i.e., with $\tau^{\alpha\beta}$ simply given by the gravitational source term $\Lambda^{\alpha\beta} = \Lambda_{\text{harm}}^{\alpha\beta}$, and we equate term by term the factors of the successive powers of our book-keeping parameter $G$. We get an infinite set of equations for each of the $h_{(n)}^{\alpha\beta}$'s: namely, $\forall n \ge 2$,

$$\Box h_{(n)}^{\alpha\beta} = \Lambda_{(n)}^{\alpha\beta} [h_{(1)}, h_{(2)}, \dots, h_{(n-1)}], \qquad (68a)$$

$$\partial_{\mu} h_{(n)}^{\alpha\mu} = 0. \qquad (68b)$$

The right-hand side of the wave equation (68a) is obtained from inserting the previous iterations, known up to the order $n-1$, into the gravitational source term. In more details, the series of equations (68a) reads

$$\Box h_{(2)}^{\alpha\beta} = N^{\alpha\beta} [h_{(1)}, h_{(1)}], \qquad (69a)$$

$$\Box h_{(3)}^{\alpha\beta} = M^{\alpha\beta} [h_{(1)}, h_{(1)}, h_{(1)}] + N^{\alpha\beta} [h_{(1)}, h_{(2)}] + N^{\alpha\beta} [h_{(2)}, h_{(1)}]. \qquad (69b)$$

∴

The quadratic and cubic pieces of $\Lambda_{\text{harm}}^{\alpha\beta}$ are defined by Eqs. (58).

Let us now proceed by induction. Some $n \in \mathbb{N}$ being given, we assume that we succeeded in constructing, starting from the linearized solution $h_{(1)}$, the sequence of post-Minkowskian solutions $h_{(2)}$, $h_{(3)}$, ..., $h_{(n-1)}$, and from this we want to infer the next solution $h_{(n)}$. The right-hand side of Eq. (68a), $\Lambda_{(n)}^{\alpha\beta}$, is known by induction hypothesis. Thus the problem is that of solving a flat wave equation whose source is given. The point is that this wave equation, instead of being valid everywhere in $\mathbb{R}^3$, is physically correct only outside the matter source ($r > a$), and it makes no sense to solve it by means of the usual retarded integral. Technically speaking, the right-hand
side of Eq. (68a) is composed of the product of many multipole expansions, which are singular at the origin of the spatial coordinates $r = 0$, and which make the retarded integral divergent at that point. This does not mean that there are no solutions to the wave equation, but simply that the retarded integral does not constitute the appropriate solution in that context.

What we need is a solution which takes the same structure as the source term $\Lambda_{(n)}^{\alpha\beta}$, i.e., is expanded into multipole contributions, with a singularity at $r = 0$, and satisfies the d'Alembertian equation as soon as $r > 0$. Such a particular solution can be obtained by means of a mathematical trick, in which one first “regularizes” the source term $\Lambda_{(n)}^{\alpha\beta}$ by multiplying it by the factor $r^B$, where $r = |\mathbf{x}|$ is the spatial radial distance and $B$ is a complex number, $B \in \mathbb{C}$ (Blanchet and Damour 1986).

Let us assume, for definiteness, that $\Lambda_{(n)}^{\alpha\beta}$ is composed of multipolar pieces with maximal multipolarity $\ell_{\max}$. This means that we start the iteration from the linearized metric (64)–(66) in which the multipolar sums are actually finite.<fnref n="25" /> The divergences when $r \to 0$ of the source term are typically power-like, say $1/r^k$ (there are also powers of the logarithm of $r$), and with the previous assumption there will exist a maximal order of divergency, say $k_{\max}$. Thus, when the real part of $B$ is large enough, i.e., $\Re(B) > k_{\max} - 3$, the “regularized” source term $r^B \Lambda_{(n)}^{\alpha\beta}$ is regular enough when $r \to 0$ so that one can perfectly apply the retarded integral operator. This defines the $B$-dependent retarded integral, when $\Re(B)$ is large enough,

$$ \mathcal{I}^{\alpha\beta}(B) \equiv \Box_{\text{ret}}^{-1} [\tilde{r}^B \Lambda_{(n)}^{\alpha\beta}], \quad (70) $$

where the symbol $\Box_{\text{ret}}^{-1}$ stands for the usual retarded d'Alembertian integral, extending over all space $\mathbf{x}' \in \mathbb{R}^3$ and given by

$$ (\Box_{\text{ret}}^{-1} f)(\mathbf{x}, t) = -\frac{1}{4\pi} \iiint \frac{d^3\mathbf{x}'}{|\mathbf{x} - \mathbf{x}'|} f(\mathbf{x}', t - |\mathbf{x} - \mathbf{x}'|/c). \quad (71) $$

Subject to the past-stationarity condition (60), it is clear that the retarded integral will always be convergent at infinity.

It is convenient to introduce together with the regularizing factor $r^B$ some arbitrary constant length scale $r_0$ in order to make it dimensionless. Everywhere in this article we pose

$$ \tilde{r} \equiv \frac{r}{r_0}. \quad (72) $$

The fate of the constant $r_0$ in a detailed calculation will be interesting to follow, as we shall see. Now the point for our purpose is that the function $\mathcal{I}^{\alpha\beta}(B)$ on the complex plane, which was originally defined only when $\Re(B) > k_{\max} - 3$, admits a unique

<footnote n="25">²⁵ This assumption is justified because we are ultimately interested in the radiation field at some given finite post-Newtonian precision like 4PN, and because only a finite number of multipole moments can contribute at any finite order of approximation. With a finite number of multipoles in the linearized metric (64)–(66), there is a maximal multipolarity $\ell_{\max}(n)$ at any post-Minkowskian order $n$, which grows linearly with $n$.</footnote>
analytic continuation to all values of $B \in \mathbb{C}$ except at some integer values. Furthermore, the analytic continuation of $\mathcal{I}^{\alpha\beta}(B)$ can be expanded, when $B \to 0$ (namely the limit of interest to us) into a Laurent expansion involving in general some multiple poles. The key idea, as we shall prove, is that the finite part, or the coefficient of the zeroth power of $B$ in that expansion, represents the particular solution we are looking for. We write the Laurent expansion of $\mathcal{I}^{\alpha\beta}(B)$, when $B \to 0$, in the form

$$ \mathcal{I}^{\alpha\beta}(B) = \sum_{p=p_0}^{+\infty} I_p^{\alpha\beta} B^p, \qquad (73) $$

where $p \in \mathbb{Z}$, and the various coefficients $I_p^{\alpha\beta}$ are functions of the field point $(\mathbf{x}, t)$. When $p_0 \le -1$ there are poles; and $-p_0$, which depends on the post-Minkowskian order $n$, refers to the maximal order of these poles (we have $p_0 \to -\infty$ when $n \to +\infty$). By applying the d'Alembertian operator onto both sides of Eq. (73), and equating the different powers of $B$, we arrive at

$$ p_0 \le p \le -1 \quad \Rightarrow \quad \Box I_p^{\alpha\beta} = 0, \qquad (74a) $$

$$ p \ge 0 \quad \Rightarrow \quad \Box I_p^{\alpha\beta} = \frac{(\ln r)^p}{p!} \Lambda_{(n)}^{\alpha\beta}. \qquad (74b) $$

As we see, the case $p = 0$ shows that the finite-part coefficient in Eq. (73), namely $I_0^{\alpha\beta}$, is a particular solution of the requested equation: $\Box I_0^{\alpha\beta} = \Lambda_{(n)}^{\alpha\beta}$. Furthermore, we can prove that this solution, by its very construction, owns the same structure made of a multipolar expansion singular at $r = 0$ as the corresponding source.

Let us forget about the intermediate name $I_0^{\alpha\beta}$, and denote, from now on, the latter solution by $u_{(n)}^{\alpha\beta} \equiv I_0^{\alpha\beta}$, or, in more explicit terms,

$$ u_{(n)}^{\alpha\beta} = \mathrm{FP}_{B=0} \Box_{\mathrm{ret}}^{-1} [\tilde{r}^B \Lambda_{(n)}^{\alpha\beta}], \qquad (75) $$

where the finite-part symbol $\mathrm{FP}_{B=0}$ means the previously detailed operations of considering the analytic continuation, taking the Laurent expansion, and picking up the finite-part coefficient when $B \to 0$. The story is not complete, however, because $u_{(n)}^{\alpha\beta}$ does not fulfill the constraint of harmonic coordinates (68b); its divergence, say $w_{(n)}^{\alpha} = \partial_{\mu}u_{(n)}^{\alpha\mu}$, is different from zero in general. From the fact that the source term is divergence-free in vacuum, $\partial_{\mu}\Lambda_{(n)}^{\alpha\mu} = 0$ [see Eq. (59)], we find instead

$$ w_{(n)}^{\alpha} = \mathrm{FP}_{B=0} \Box_{\mathrm{ret}}^{-1} \left[ B \tilde{r}^B \frac{n_i}{r} \Lambda_{(n)}^{\alpha i} \right]. \qquad (76) $$

The factor $B$ comes from the differentiation of the regularization factor $\tilde{r}^B$. So, $w_{(n)}^{\alpha}$ is zero only in the special case where the Laurent expansion of the retarded integral in Eq. (76) does not develop any simple pole when $B \to 0$. Fortunately, when it does, the structure of the pole is quite easy to control. We find that it necessarily consists of
an homogeneous solution of the source-free d'Alembertian equation, and, what is more (as consequence of its stationarity in the past), that solution is a retarded one. Hence, taking into account the index structure of $w_{(n)}^{\alpha}$, there must exist four STF-tensorial functions of $u = t - r/c$, say $N_L$, $P_L$, $Q_L$ and $R_L$, such that

$$w_{(n)}^{0} = \sum_{\ell=0}^{+\infty} \partial_{L} [r^{-1} N_{L}(u)], \qquad (77a)$$

$$w_{(n)}^{j} = \sum_{\ell=0}^{+\infty} \partial_{iL} [r^{-1} P_{L}(u)] + \sum_{\ell=1}^{+\infty} \left\{ \partial_{L-1} [r^{-1} Q_{iL-1}(u)] + \epsilon_{iab} \partial_{aL-1} [r^{-1} R_{bL-1}(u)] \right\}. \qquad (77b)$$

From that expression we are able to find a new object, say $v_{(n)}^{\alpha\beta}$, which takes the same structure as $w_{(n)}^{\alpha}$ (a retarded solution of the source-free wave equation) and, furthermore, whose divergence is exactly the opposite of the divergence of $u_{(n)}^{\alpha\beta}$, i.e. $\partial_{\mu}v_{(n)}^{\alpha\mu} = -w_{(n)}^{\alpha}$. Such a $v_{(n)}^{\alpha\beta}$ is not unique, but we shall see that it is just necessary to make a choice for $v_{(n)}^{\alpha\beta}$ (the simplest one) in order to obtain the general solution. The formulas that we adopt are

$$v_{(n)}^{00} = -r^{-1}N^{(-1)} + \partial_a \left[ r^{-1} \left( -N_a^{(-1)} + Q_a^{(-2)} - 3P_a \right) \right], \qquad (78a)$$

$$v_{(n)}^{0i} = r^{-1}(-Q_i^{(-1)} + 3P_i^{(1)}) - \epsilon_{iab}\partial_a[r^{-1}R_b^{(-1)}] - \sum_{\ell=2}^{+\infty}\partial_{L-1}[r^{-1}N_{iL-1}], \qquad (78b)$$

$$v_{(n)}^{ij} = -\delta_{ij}r^{-1}P + \sum_{\ell=2}^{+\infty} \left\{ \begin{aligned} & 2\delta_{ij}\partial_{L-1}[r^{-1}P_{L-1}] - 6\partial_{L-2}(i[r^{-1}P_{jL-2}]) \\ & + \partial_{L-2}[r^{-1}(N_{ijL-2}^{(1)} + 3P_{ijL-2}^{(2)} - Q_{ijL-2})] - 2\partial_{aL-2}[r^{-1}\epsilon_{ab}(iR_{j}bL-2)] \end{aligned} \right\}. \qquad (78c)$$

Notice the presence of anti-derivatives, denoted e.g. $N^{(-1)}(u) = \int_{-\infty}^{u} dv N(v)$; there is no problem with the limit $v \to -\infty$ since all the corresponding functions are zero when $t \le -T$. The choice made in Eqs. (78) is dictated by the fact that the 00 component involves only some monopolar and dipolar terms, and that the spatial trace $ii$ is monopolar: $v_{(n)}^{ii} = -3r^{-1}P$. In the following we shall refer to the previous operation $w_{(n)}^{\alpha} \to v_{(n)}^{\alpha\beta}$ as the “harmonicity” algorithm and denote it by $V^{\alpha\beta}$ so that

$$v_{(n)}^{\alpha\beta} \equiv V^{\alpha\beta}[w_{(n)}]. \qquad (79)$$

Finally we see that we solve at once the d'Alembertian equation (68a) and the coordinate condition (68b) if we pose
$$h_{(n)}^{\alpha\beta} = u_{(n)}^{\alpha\beta} + v_{(n)}^{\alpha\beta}. \qquad (80)$$

That is, we have succeeded in finding a solution of the field equations at the n-th post-Minkowskian order. By induction the same method applies to any order n, and, therefore, we have constructed a complete post-Minkowskian series (67) based on the linearized approximation $h_{(1)}^{\alpha\beta}$ given by Eqs. (64)–(66). The previous procedure constitutes an algorithm, which has been implemented on a computer programme by Blanchet et al. (2008), Faye et al. (2012).

### 2.1.4 Source multipole moments versus canonical moments

The six sets of multipole moments $I_L(u), \dots, Z_L(u)$ contain the physical information about any isolated source as seen in its exterior. Indeed we have the following result.

**Theorem 2 (Blanchet and Damour 1986)** The most general solution of the harmonic-coordinates Einstein field equations in the vacuum region outside an isolated source, admitting some post-Minkowskian and multipolar expansions, is given by the previous construction as

$$h_{\text{MPM}}^{\alpha\beta} = \sum_{n=1}^{+\infty} G^n h_{(n)}^{\alpha\beta} [I_L, J_L, \dots, Z_L]. \qquad (81)$$

It depends on two sets of arbitrary STF-tensorial functions of time $I_L(u)$ and $J_L(u)$ (satisfying the conservation laws) defined by Eqs. (1), and on four supplementary functions $W_L(u), \dots, Z_L(u)$ parametrizing the gauge vector (66).

The proof is quite easy. With Eq. (80) we obtained a particular solution of the system of equations (68). To it we should add the most general solution of the corresponding homogeneous system of equations, which is obtained by setting $A_{(n)}^{\alpha\beta} = 0$ into Eqs. (68). But this homogeneous system of equations is nothing but the linearized vacuum field equations (62), to which we know the most general solution $h_{(1)}^{\alpha\beta}$ given by Eqs. (64)–(66). Thus, we must add to our particular solution $h_{(n)}^{\alpha\beta}$ a general homogeneous solution that is necessarily of the type $h_{(1)}^{\alpha\beta}$ [$\delta I_L, \dots, \delta Z_L$], where $\delta I_L, \dots, \delta Z_L$ denote some corrections to the multipole moments at the n-th post-Minkowskian order (with the monopole $\delta I$ and dipoles $\delta I_i$, $\delta J_i$ being constant). It is then clear, since precisely the linearized metric is a linear functional of all these moments, that the previous corrections to the moments can be absorbed into a re-definition of the original ones $I_L, \dots, Z_L$ by posing

$$I_L^{\text{new}} = I_L + G^{n-1} \delta I_L, \qquad (82a)$$

$$\vdots \qquad (82b)$$

$$Z_L^{\text{new}} = Z_L + G^{n-1} \delta Z_L.$$

After re-arranging the metric in terms of these new moments, taking into account the
fact that the precision of the metric is limited to the n-th post-Minkowskian order, and dropping the superscript “new”, we find exactly the same solution as the one we had before (indeed, the moments are arbitrary functions of time)—hence the proof.

However, it is always possible to find two, and only two, sets of STF multipole moments, $M_L(u)$ and $S_L(u)$, for parametrizing the most general isolated source as well. The route for constructing such a general solution is to get rid of the moments $W_L, X_L, Y_L, Z_L$ at the linearized level by performing the linearized gauge transformation $\delta x^\alpha = \phi_{(1)}^\alpha$, where $\phi_{(1)}^\alpha$ is the gauge vector given by Eqs. (66). So, at the linearized level, we have only the two types of moments $M_L$ and $S_L$, parametrizing $h_{\text{can}(1)}^{\alpha\beta}$ by the same formulas as in Eqs. (1), hence the new “canonical” algorithm starts with

$$h_{\text{can MPM}}^{\alpha\beta} = G h_{\text{can}(1)}^{\alpha\beta} [M_L, S_L] + \mathcal{O}(G^2). \quad (83)$$

We must be careful to denote these moments with names different from $I_L$ and $J_L$ because $M_L$ and $S_L$ will ultimately differ from $I_L$ and $J_L$ for a given physical source (see Theorem 3). Then we apply exactly the same post-Minkowskian algorithm, following the formulas (75)–(80) as we did above, but starting from the gauge-transformed linear metric $h_{\text{can}(1)}^{\alpha\beta}$ instead of $h_{(1)}^{\alpha\beta}$. That is, we successively define

$$ \left. \begin{aligned} u_{\text{can}(n)}^{\alpha\beta} &= \text{FP} \Box_{\text{ret}}^{-1} [\tilde{r}^B \Lambda_{\text{can}(n)}^{\alpha\beta}] \\ w_{\text{can}(n)}^{\alpha} &= \partial_{\mu} u_{\text{can}(n)}^{\alpha\mu} \\ v_{\text{can}(n)}^{\alpha\beta} &= \mathcal{V}^{\alpha\beta}[w_{\text{can}(n)}] \end{aligned} \right\} \Rightarrow h_{\text{can}(n)}^{\alpha\beta} = u_{\text{can}(n)}^{\alpha\beta} + v_{\text{can}(n)}^{\alpha\beta}, \quad (84) $$

where the source term $\Lambda_{\text{can}(n)}^{\alpha\beta}$ depends on previous iterations $h_{\text{can}(m)}^{\alpha\beta}$ (with $m \le n - 1$), and where the harmonicity algorithm $\mathcal{V}^{\alpha\beta}$ is defined in Eq. (79).

**Theorem 3** The most general solution of the Einstein field equations in the vacuum region outside an isolated source is given by

$$h_{\text{can MPM}}^{\alpha\beta} = \sum_{n=1}^{+\infty} G^n h_{\text{can}(n)}^{\alpha\beta} [M_L, S_L]. \quad (85)$$

The metric (85) is isometric to the original one (81) (i.e., describe the same physical source) if and only if the canonical moments $M_L$ and $S_L$ are related to the source moments $I_P, J_P, \dots, Z_P$ by non-linear (post-Minkowskian) relations

$$M_L = I_L + \sum_{n=2}^{+\infty} G^{n-1} \mathcal{M}_{nL} [I_P, J_P, \dots, Z_P], \quad (86a)$$

$$S_L = J_L + \sum_{n=2}^{+\infty} G^{n-1} S_{nL} [I_P, J_P, \dots, Z_P]. \quad (86b)$$
Obviously the canonical post-Minkowskian algorithm (85) yields some simpler calculations than the general algorithm (81) as we have only two multipole moments to iterate. We shall give in Eqs. (251)–(252) the most up to date relations between the canonical moments and source moments.

So why not consider from the start that the best description of the isolated source is provided by only the two types of multipole moments, $M_L$ and $S_L$, instead of the six types, $I_L, J_L, \dots, Z_L$? The reason is that we shall determine in Theorem 7 the explicit closed-form expressions of the six source moments $I_L, J_L, \dots, Z_L$, but that, by contrast, it seems to be impossible to obtain some similar closed-form expressions for the canonical moments $M_L$ and $S_L$. The only thing we can do is to write down the explicit non-linear algorithm that computes $M_L, S_L$ starting from $I_L, J_L, \dots, Z_L$, see Eqs. (3). In consequence, it is better to view the moments $I_L, J_L, \dots, Z_L$ as more “fundamental” than $M_L$ and $S_L$, in the sense that they appear to be more tightly related to the description of the source, since they admit closed-form expressions as some explicit integrals over the source. Hence, we choose to refer collectively to the six moments $I_L, J_L, \dots, Z_L$ as the multipole moments of the source. This being said, the moments $M_L$ and $S_L$ are generally very useful in practical computations, not only because they yield a simpler post-Minkowskian iteration but also because they are closely related to the two independent degrees of freedom of gravitational waves propagating at infinity from the source. Then, one can come back to the more fundamental source-rooted moments by using the relations (3), in particular the fact that $M_L$ and $S_L$ differ from the corresponding $I_L$ and $J_L$ only by high-order post-Newtonian terms starting at 2.5PN order; see Eqs. (251)–(252) below. Indeed, this is to be expected because the physical difference between both types of moments stems only from non-linearities.

Finally, let us point out that in our presentation of the post-Minkowskian algorithms, either the general one (81) or the canonical one (85), we have for the moment omitted a crucial recursive hypothesis, which is required in order to prove that at each post-Minkowskian order $n$, the inverse d'Alembertian operator can be applied in the way we did—notably the fact that the B-dependent retarded integral can be analytically continued down to a neighbourhood of $B = 0$. This hypothesis is that the “near-zone” expansion, i.e., when $r \to 0$, of each one of the post-Minkowskian coefficients $h_{(n)}^{\alpha\beta}$ [or $h_{\text{can}(n)}^{\alpha\beta}$] has a certain structure; this hypothesis is established as a theorem once the mathematical induction succeeds.

**Theorem 4** (Blanchet and Damour 1986) *The general structure of the expansion of the post-Minkowskian exterior metric in the near-zone (when $r \to 0$) is of the type: $\forall N \in \mathbb{N},^{26}$*

$$h_{(n)}^{\alpha\beta}(\mathbf{x}, t) = \sum \hat{n}_L r^m (\ln r)^p F_{L,m,p,n}^{\alpha\beta}(t) + o(r^N), \quad (87)$$

where $m \in \mathbb{Z}$, with $m_0 \le m \le N$ (and $m_0$ becoming more and more negative as $n$ grows), $p \in \mathbb{N}$ with $p \le n - 1$. The functions $F_{L,m,p,n}$ are multilinear functionals of

<footnote n="26">26 We employ the Landau symbol $o$ for remainders with its standard meaning. Thus, $f(r) = o[g(r)]$ when $r \to 0$ means that $f(r)/g(r) \to 0$ when $r \to 0$. Furthermore, we generally assume some differentiability properties such as $d^n f(r)/dr^n = o[g(r)/r^n]$.</footnote>
the source multipole moments $I_L, \dots, Z_L$ or the canonical moments $M_L, S_L$ (depending on the construction under consideration, $h_{(n)}^{\alpha\beta}$ or $h_{\text{can}(n)}^{\alpha\beta}$).

As we see, the near-zone expansion involves, besides the simple powers of $r$, some integer powers of the logarithm of $r$, with a maximal power of $n-1$. As a corollary of that theorem, we find, by restoring all the powers of $c$ in Eq. (87) and using the fact that each $r$ goes into the combination $r/c$, that the general structure of the post-Newtonian expansion ($c \to +\infty$) is of the type

$$h_{(n)}(c) \approx \sum_{p,q \in \mathbb{N}} \frac{(\ln c)^p}{c^q}, \qquad (88)$$

where $p \le n-1$ (and $q \ge 2$). The post-Newtonian expansion proceeds not only with the normal powers of $1/c$ but also with powers of the logarithm of $c$. It is remarkable that the structure of the post-Newtonian expansion (at high orders) is not more complicated than (88); for instance there are no terms behaving like $\sim \ln(\ln c)$.

### 2.1.5 The multipolar post-Minkowskian solution in radiative coordinates

We now tackle the important problem of the far-zone structure of this MPM construction and the link with observables quantities at infinity. Paralleling the structure of the near-zone expansion given by (87), we have a similar result concerning the *far-zone* expansion at Minkowskian future null infinity, i.e., when $r \to +\infty$ with $u = t - r/c = \text{const}: \forall N \in \mathbb{N}$,

$$h_{(n)}^{\alpha\beta}(\mathbf{x}, t) = \sum \frac{\hat{n}_L (\ln r)^p}{r^k} G_{L,k,p,n}^{\alpha\beta}(u) + o\left(\frac{1}{r^N}\right), \qquad (89)$$

where $k, p \in \mathbb{N}$, with $1 \le k \le N$, and where, likewise in the near-zone expansion (87), some integer powers of logarithms, such that $p \le n-1$, appear.

The appearance of logarithms in the far-zone expansion of the harmonic-coordinates metric has been known since the work of Fock (1959). One knows also that this is a coordinate effect, because the study of the “asymptotic” structure of space-time at future null infinity has been consistently performed in a coordinate system that avoids the appearance of far zone logarithms (Bondi et al. 1962; Sachs 1962; Newman and Unti 1963). This asymptotic structure was further elucidated thanks to the tools of the Newman and Penrose (1962) formalism and conformal compactifications leading to the concept of asymptotically “simple” spacetimes in the sense of Penrose (1963, 1965). More generally, a large class of so-called *radiative* coordinate systems exist (Papapetrou 1969; Madore 1970a, b), in which the far-zone expansion of the metric proceeds with simple powers of the inverse radial distance. Hence, the logarithms are simply an artifact of the use of harmonic coordinates (Isaacson and Winicour 1968; Blanchet 1987). An algorithm to go from harmonic coordinates to Newman and Unti (1963) coordinates has been implemented by Blanchet et al. (2021, 2023a). The following theorem shows that general MPM construction of the metric in the exterior of the source, when developed at future null
infinity, is consistent with the Bondi et al. (1962), Sachs (1962), Penrose (1963, 1965) approach to gravitational radiation.

**Theorem 5** (Blanchet 1987) The most general multipolar-post-Minkowskian solution, stationary in the past [see Eq.(60)], admits some radiative coordinates $(T, \mathbf{X})$, for which the expansion at future null infinity, $R \to +\infty$ with $U \equiv T - R/c = \text{const}$, takes the form

$$h_{\text{rad}(n)}^{\alpha\beta}(\mathbf{X}, T) = \sum_{k \ge 1} \frac{\hat{N}_L}{R^k} K_{L,k,n}^{\alpha\beta}(U) + \mathcal{O}\left(\frac{1}{R^N}\right). \qquad (90)$$

The functions $K_{L,k,n}^{\alpha\beta}$ are computable functionals of the source multipole moments. In radiative coordinates the retarded time $U$ is a null coordinate in the asymptotic far zone limit. Furthermore the metric

$$h_{\text{rad MPM}}^{\alpha\beta} = \sum_{n=1}^{+\infty} G^n h_{\text{rad}(n)}^{\alpha\beta}, \qquad (91)$$

is asymptotically simple in the sense of Penrose (1963,1965), Geroch and Horowitz (1978),<fnref n="27" /> perturbatively to any post-Minkowskian order.

The proof consists of performing the multipolar-post-Minkowskian construction directly in radiative coordinates, thus solving the Einstein field equations in vacuum, in a non-harmonic coordinate system:

$$\Box h^{\alpha\beta} - \partial\hat{H}^{\alpha\beta} = \Lambda^{\alpha\beta}. \qquad (92)$$

The source term is defined by Eqs. (55)–(56). In radiative coordinates the retarded time coordinate $u = t - r/c$ is null, i.e. satisfies<fnref n="28" />

$$(\eta^{\alpha\beta} + h_{\text{radMPM}}^{\alpha\beta}) \partial_{\alpha} u \partial_{\beta} u = 0, \qquad (93)$$

or at least, becomes a null coordinate in the asymptotic limit $r \to +\infty$ with $u = \text{const}$. Even at the linearized level, it is necessary to correct the harmonic coordinate metric in order to satisfy the requirement of asymptotically null retarded time. Consequently, the radiative MPM algorithm starts by performing a linear gauge transformation of the harmonic-coordinate metric defined by Eq. (83). A crucial point is that the multipole moments that parametrize the radiative algorithm will differ

<footnote n="27">²⁷ If the assumption of stationarity in the past is relaxed, there may be a loss of smoothness of the metric at future null infinity (Damour 1986).</footnote>

<footnote n="28">²⁸ Since we are iteratively constructing a coordinate system order by order, it is convenient to consider the coordinates as dummy variables and denote them by $(\mathbf{x}, t)$. In the end, when we obtain the full radiative metric, we shall denote the radiative coordinates by $(\mathbf{X}, T)$ as in Theorem 5.</footnote>
from their counterparts in the harmonic algorithm (for a given physical source). Hence we construct the radiative metric using multipole moments {$\overline{M}_L, \overline{S}_L$} which differ from the canonical moments {$M_L, S_L$} used in harmonic coordinates.<fnref n="29" /> At linear order we pose

$$h_{\text{rad}(1)}^{\alpha\beta} = h_{\text{can}(1)}^{\alpha\beta} [\overline{M}_L, \overline{S}_L] + \partial_{\xi_{(1)}}^{\alpha\beta}, \quad (94)$$

where $h_{\text{can}(1)}^{\alpha\beta}$ takes the same functional form as in (1), but is now computed with the moments {$\overline{M}_L, \overline{S}_L$}. The linear gauge transformation $\partial_{\xi_{(1)}}^{\alpha\beta}$ is defined by the gauge vector

$$\xi_{(1)}^{\alpha} = \frac{2M}{c^2} \eta^{0\alpha} \ln\left(\frac{r}{b_0}\right), \quad (95)$$

where M is the total ADM mass, $b_0$ denotes an arbitrary length scale, and we have $\eta^{0\alpha} = (-1, \mathbf{0})$ with our signature.

The effect of this linear gauge transformation is to correct for the well-known logarithmic deviation of the retarded time in harmonic coordinates, with respect to the true light cone. After the change of gauge, the coordinate $u = t - r/c$ coincides (asymptotically when $r \to \infty$) with a null coordinate at the linearized level. The linear radiative metric is no longer harmonic, since

$$H_{\text{rad}(1)}^{\alpha} \equiv \partial_{\mu} h_{\text{rad}(1)}^{\alpha\mu} = \Box \xi_{(1)}^{\alpha} = \frac{2M}{c^2 r^2} \eta^{0\alpha}. \quad (96)$$

Given any $n \ge 2$, let us recursively assume that: (i) we have obtained all the previous radiative post-Minkowskian coefficients $h_{\text{rad}(m)}^{\alpha\beta}$ for any $m \le n-1$; (ii) all of them admit an expansion as $r \to +\infty$ with $u = \text{const}$ in simple powers of $1/r$ (no logarithms); and (iii) all the previous coefficients satisfy the fall-off conditions ($\forall m \le n-1$)

$$k_{\alpha} k_{\beta} h_{\text{rad}(m)}^{\alpha\beta} = \mathcal{O}\left(\frac{1}{r^2}\right), \quad (97a)$$

$$H_{\text{rad}(m)}^{\alpha} = \mathcal{O}\left(\frac{1}{r^2}\right), \quad (97b)$$

where $k^\alpha = \eta^{\alpha\mu} k_\mu = (1, \mathbf{n})$ denotes the outgoing Minkowskian null vector. Denoting $z_{(m)}^{\alpha\beta}$ the leading $1/r$ coefficient of the metric when $r \to +\infty$, i.e.

<footnote n="29">²⁹ Although the two sets of moments {$M_L, S_L$} and {$\overline{M}_L, \overline{S}_L$} differ in general, the conserved mass monopole (as well as the mass and current dipoles) are in fact identical in the harmonic and radiative constructions. Hence, $\overline{M} = M$ for the ADM mass. See Eq. (5.4) of Trestini et al. (2023) for the relationship between the quadrupole moment $\overline{M}_{ij}$ in the radiative algorithm and the moments {$M_L, S_L$} in harmonic coordinates.</footnote>
$$h_{\text{rad}(m)}^{\alpha\beta} = \frac{1}{r} z_{(m)}^{\alpha\beta}(u, \mathbf{n}) + \mathcal{O}\left(\frac{1}{r^2}\right), \qquad (98)$$

we see that the two conditions (97) say that $k_\alpha k_\beta z_{(m)}^{\alpha\beta} = 0$ while $k_\mu z_{(m)}^{\alpha\mu} = 0$ is constant. For a stationary metric (independent of time) we know that the m-th post-Minkowskian coefficient decreases as $1/r^m$ when $r \to +\infty$. By our assumption of stationarity in the past, Eq. (60), this means that $h_{\text{rad}(m)}^{\alpha\beta} = \mathcal{O}(r^{-m})$ when $t \le -T$, hence the latter constant is zero when $m \ge 2$. Hence the recursive assumptions (97) are equivalent to

$$k_{\mu} z_{(1)}^{\alpha\mu} = \frac{2M}{c^2} k^{\alpha}, \qquad (99a)$$

$$k_{\mu} z_{(m)}^{\alpha\mu} = 0 \quad (\text{for } 2 \le m \le n-1). \qquad (99b)$$

One can easily check that the condition is satisfied by the linear metric (94).

The dominant piece when $r \to +\infty$ of the non-linear source term at the n-th order will be of order $1/r^2$ and will only be made of quadratic products of $h_{\text{rad}(m)}^{\alpha\beta}$, since each of the $h_{\text{rad}(m)}^{\alpha\beta}$'s behaves like $1/r$. Under our recursive assumptions, in particular (99), and from the structure of the source term at quadratic order, see Eq. (58a), one can prove that the n-th post-Minkowskian source term takes the form at leading order when $r \to +\infty$:

$$\Lambda_{\text{rad}(n)}^{\alpha\beta} = \frac{k^\alpha k^\beta}{r^2} \sigma_{(n)}(u, \mathbf{n}) + \mathcal{O}\left(\frac{1}{r^3}\right). \qquad (100)$$

This is the form of the stress-energy tensor of massless particles, i.e. gravitons in our case, with $\sigma_{(n)}$ being proportional to the n-th order contribution in the total power emitted by the massless radiation. The explicit expression in terms of the metric coefficient (98) is (the dot denoting the time derivative)

$$\sigma_{(n)} = \sum_{m=1}^{n-1} \left( \frac{1}{2} \eta_{\alpha\gamma} \eta_{\beta\delta} - \frac{1}{4} \eta_{\alpha\beta} \eta_{\gamma\delta} \right) \dot{z}_{(m)}^{\alpha\beta} \dot{z}_{(n-m)}^{\gamma\delta}. \qquad (101)$$

We know that logarithms in the asymptotic expansion when $r \to +\infty$ only arise due to the retarded integral of source terms that behave like $1/r^2$ (see Appendix A of Blanchet 1998c). Hence the dominant term written in Eq. (100) is the only piece of $\Lambda_{\text{rad}(n)}^{\alpha\beta}$ that can yield logarithms at order $n$; see indeed the integration formula (212), which behaves like $\ln r/r$ at infinity. But now, thanks to the particular structure of this term in (100), which follows from our recursive assumptions, we can gauge it away, thus constructing a coordinate system valid at the n-th post-Minkowskian order which avoids the appearance of logarithms. We find that an adequate gauge vector is

$$\xi_{(n)}^{\alpha} = \Box_{\text{ret}}^{-1} \left[ \frac{c k^{\alpha}}{2r^2} \int_{-\infty}^{u} dv \sigma_{(n)}(v, \mathbf{n}) \right]. \qquad (102)$$

The retarded integral is convergent at the bound $r \to 0$ and so there is no need to
include a finite part operation FP. With this choice of gauge vector, the logarithms that will be generated from the gauge transformation will cancel the logarithms coming from the retarded integral of the source term (100). Hence, similarly to the corresponding steps (84) in the harmonic algorithm, we construct the n-th post-Minkowskian metric by correcting for the new logarithms using the above gauge transformation:

$$ 
\left. 
\begin{aligned}
u_{\text{rad}(n)}^{\alpha\beta} &= \text{FP} \Box_{\text{ret}}^{-1} [\tilde{r}^B \Lambda_{\text{rad}(n)}^{\alpha\beta}] \\
w_{\text{rad}(n)}^{\alpha} &= \partial_{\mu} u_{\text{rad}(n)}^{\alpha\mu} \\
v_{\text{rad}(n)}^{\alpha\beta} &= \mathcal{V}^{\alpha\beta}[w_{\text{rad}(n)}]
\end{aligned}
\right\} \Rightarrow h_{\text{rad}(n)}^{\alpha\beta} = u_{\text{rad}(n)}^{\alpha\beta} + v_{\text{rad}(n)}^{\alpha\beta} + \partial_{\xi}(z_{(n)}^{\alpha\beta}). \quad (103)
$$

Note that the combination $u_{\text{rad}(n)}^{\alpha\beta} + v_{\text{rad}(n)}^{\alpha\beta}$ is divergenceless, so the radiative metric obeys by construction the non-harmonic gauge condition

$$ H_{\text{rad}(n)}^{\alpha} = \Box_{\xi(n)}^{\alpha\xi} = \frac{c k^{\alpha}}{2r^2} \int_{-\infty}^{t-r/c} dv \, \sigma_{(n)}(v, \mathbf{n}), \quad (104) $$

and the Einstein field equations (92) are trivially satisfied to order $n$. The far-zone expansion of the full non-linear radiative metric constructed by virtue of this procedure is free of any logarithms, and the retarded time $u = t - r/c$ in these coordinates tends asymptotically toward a null coordinate at future null infinity [i.e., satisfying (93)]. The property of asymptotic simplicity in the form given by Geroch and Horowitz (1978), is proved by introducing the conformal factor $\Omega = 1/r$ in radiative coordinates (Blanchet 1987). Finally, as stated in Theorem 5, the metric, as a general functional of the moments $\{\overline{M}_L, \overline{S}_L\}$, represents physically the most general solution to the vacuum field equations outside the source, as legitimate as the general harmonic-coordinate metric of Theorem 2.

Finally, having constructed the metric step by step in a radiative (or Bondi-type) coordinate system, we formally sum up the whole post-Minkowskian series from $n = 1$ up to $+\infty$, thus defining

$$ h_{\text{rad}}^{\alpha\beta} \equiv \sum_{n=1}^{+\infty} G^n h_{\text{rad}(n)}^{\alpha\beta}, \quad (105) $$

Posing then $z^{\alpha\beta}$ to be the leading $1/R$ coefficient of the radiative metric [we restore the names $(T, X)$ appropriate for radiative coordinates],

$$ h_{\text{rad}}^{\alpha\beta} = \frac{1}{R} z^{\alpha\beta}(U, N) + \mathcal{O}\left(\frac{1}{R^2}\right), \quad (106) $$

where $z^{\alpha\beta} = \sum_{n \ge 1} G^n z_{(n)}^{\alpha\beta}$, we obtain after resummation the leading $1/R^2$ coefficient of the source term in radiative coordinates as
$$A_{\text{rad}}^{\alpha\beta} = \frac{K^{\alpha}K^{\beta}}{R^2} \sigma(U, N) + \mathcal{O}\left(\frac{1}{R^3}\right), \qquad (107)$$

where $K^{\alpha} = (1, N)$ and the resummed energy distribution reads

$$\sigma(U, N) = \frac{1}{2} \dot{z}^{\alpha\beta} \dot{z}_{\alpha\beta} - \frac{1}{4} (\dot{z}_{\alpha}^{\alpha})^2 = \frac{1}{2} \perp_{ijkl} \dot{z}_{ij} \dot{z}_{kl} = \frac{1}{2} \dot{z}_{ij}^{\text{TT}} \dot{z}_{ij}^{\text{TT}}. \qquad (108)$$

Reminding from (99) that $K_{\mu\nu}\dot{z}^{\alpha\mu} = 0$, an easy calculation gives the second and third equalities in (108), where $\perp_{ijkl}(N)$ is the TT projection operator already introduced in Eq. (3), and we pose $\dot{z}_{ij}^{\text{TT}} = \perp_{ijab} \dot{z}_{ab}$. Finally the energy distribution in the waves is directly given by $\sigma(U, N)$ as

$$\left(\frac{dE}{dU d\Omega}\right)^{\text{GW}} = \frac{c^5}{16\pi G} \sigma(U, N) = \lim_{R \to +\infty} \frac{c^3 R^2}{32\pi G} \dot{h}_{ij}^{\text{TT}} \dot{h}_{ij}^{\text{TT}}, \qquad (109) \\ U = \text{const}$$

where $h_{ij}^{\text{TT}} = \perp_{ijab} h_{\text{rad ab}}$ denotes the TT waveform. The energy (and angular momentum) in the waves can be defined using the Isaacson (1968) stress-energy tensor

$$T_{\alpha\beta}^{\text{GW}} = \frac{c^4}{32\pi G} \langle \partial_{\alpha} h_{ij}^{\text{TT}} \partial_{\beta} h_{ij}^{\text{TT}} \rangle, \qquad (110)$$

where the brackets denote an average over the wavelength. This definition is general, and reduces asymptotically to

$$T_{\alpha\beta}^{\text{GW}} = \frac{c^4}{16\pi G} \frac{K_{\alpha}K_{\beta}}{R^2} \sigma + \mathcal{O}\left(\frac{1}{R^3}\right). \qquad (111)$$

We systematically investigate the properties of the TT waveform in Sect. 2.4.

## 2.2 Matching to a post-Newtonian source

By Theorem 2 we control the most general class of solutions of the vacuum equations outside the source, in the form of non-linear functionals of the source multipole moments. For instance, these solutions can describe the Schwarzschild and Kerr solutions for black holes, as well as all their perturbations, and neutron stars with arbitrary internal structure. By Theorem 5 we learned how to control the far zone expansion of the solutions; this will yield the notion of radiative moments at infinity, which constitute the observables of the radiation field far from the source (see Sect. 2.4). We now want to understand how a specific choice of matter stress-energy tensor $T^{\alpha\beta}$, i.e., a specific choice of some physical model describing the material source, selects a particular physical exterior solution among our general class, and therefore a given set of multipole moments for the source.
### 2.2.1 The matching equation

We shall provide the answer to that problem in the case of a post-Newtonian source for which the post-Newtonian parameter $\epsilon \sim 1/c$ defined by Eq. (1) is small. The fundamental fact that permits the connection of the exterior field to the inner field of the source is the existence of a “matching” region, in which both the multipole expansion and the post-Newtonian expansion are valid. This region is nothing but the exterior part of the near zone, such that $r > a$ (exterior) and $r \ll \lambda$ (near zone); it always exists around post-Newtonian sources whose radius is much less than the emitted wavelength, $a/\lambda \sim \epsilon \ll 1$, hence $a \ll \lambda$. In the present formalism the multipole expansion is defined by the multipolar-post-Minkowskian (MPM) solution; see Sect. 2.1. Matching together the post-Newtonian and MPM solutions in this overlapping region is an application of the method of matched asymptotic expansions, which has intensively been applied in the present context, both for radiation-reaction (Burke and Thorne 1970; Burke 1971; Anderson et al. 1982; Blanchet and Damour 1988; Blanchet 1993) and wave-generation (Blanchet and Damour 1989; Damour and Iyer 1991b; Blanchet 1995, 1998b) problems.

Let us denote by $M(h)$ the multipole expansion of the true metric $h$, valid everywhere inside and outside the source.<fnref n="30" /> By $M(h)$ we really mean the MPM exterior metric that we have constructed in Sect. 2.1. Specifically we adopt the general metric of Theorem 2, see Eq. (81), which is parametrized by six types of source multipole moments $I_L, J_L, W_L, X_L, Y_L, Z_L$:

$$M(h) = h_{\text{MPM}} = \sum_{n=1}^{+\infty} G^n h_{(n)} [I_L, \dots, Z_L]. \qquad (112)$$

An important point about the multipole expansion (not always well appreciated) is that it is a formal solution of the vacuum equations which is defined for any radius $r > 0$. Of course, the true solution $h$ agrees with its own multipole expansion in the exterior of the source, i.e.

$$r > a \implies M(h) = h. \qquad (113)$$

By contrast, inside the source, $h$ and $M(h)$ disagree with each other because $h$ is a fully-fledged solution of the field equations within the matter source, while $M(h)$ is a vacuum solution becoming singular at $r = 0$. Note also that the multipole expansion, physically valid when $r > a$, is to be distinguished from the far zone expansion of the field, valid when $r \gg \lambda$.

Now let us denote by $\bar{h}$ the post-Newtonian expansion of $h$. We have already anticipated the general structure of this expansion which is given in Eq. (88). In the matching region, where both the multipolar and post-Newtonian expansions are valid, we can write the numerical equality

$$a<r\ll\lambda \implies M(h)=\bar{h}. \qquad (114)$$

This “numerical” equality is viewed here in a sense of formal expansions, as we do

<footnote n="30">³⁰ For simplicity in this section, we do not write the space-time indices $\alpha\beta$ on the metric.</footnote>
not control the convergence of the series. In fact, we should be aware that such an equality, though quite natural and even physically obvious, has not really been justified mathematically within the approximation scheme, and we simply take it here as part of our fundamental assumptions.

We now transform Eq. (114) into a matching equation, by replacing in the left-hand side $\mathcal{M}(h)$ by its near-zone re-expansion $\overline{\mathcal{M}(h)}$, and in the right-hand side $\bar{h}$ by its multipole expansion $\mathcal{M}(\bar{h})$. The structure of the near-zone expansion ($r \to 0$) of the exterior multipolar field has been found in Theorem 4, see Eq. (87). We denote the corresponding infinite series $\overline{\mathcal{M}(h)}$ with the same overbar as for the post-Newtonian expansion because it is really an expansion when $r/c \to 0$, equivalent to an expansion when $c \to \infty$. Concerning the multipole expansion of the post-Newtonian metric, $\mathcal{M}(\bar{h})$, we simply postulate for the moment its existence, but we shall show later how to construct it explicitly. Therefore, the matching equation is the statement that

$$ \overline{\mathcal{M}(h)} = \mathcal{M}(\bar{h}), \qquad (115) $$

by which we really mean an infinite set of functional identities, valid $\forall (\mathbf{x}, t) \in \mathbb{R}_*^3 \times \mathbb{R}$, between the coefficients of the series in both sides of the equation. Note that such a meaning is different from that of a numerical equality like Eq. (114), which is valid only when $\mathbf{x}$ belongs to some limited spatial domain. The matching equation (115) tells us that the formal near-zone expansion of the multipole decomposition is identical, term by term, to the multipole expansion of the post-Newtonian solution. However, the latter expansion is nothing but the formal far-zone expansion, when $r \to \infty$, of each of the post-Newtonian coefficients. Most importantly, it is possible to write down, within the present formalism, the general structure of these identical expansions as a consequence of Eq. (87):

$$ \overline{\mathcal{M}(\bar{h})} = \sum \hat{n}_L r^m (\ln r)^p F_{L,m,p}(t) = \mathcal{M}(\bar{h}), \qquad (116) $$

where the functions $F_{L,m,p} = \sum_{n \ge 1} G^n F_{L,m,p,n}$. The latter expansion can be interpreted either as the singular re-expansion of the multipole decomposition when $r \to 0$—i.e., the first equality in Eq. (116), or the singular re-expansion of the post-Newtonian series when $r \to +\infty$—the second equality.

We recognize the beauty of singular perturbation theory, where two asymptotic expansions, taken formally outside their respective domains of validity, are matched together. Of course, the method works because there exists, physically, an overlapping region in which the two approximation series are expected to be numerically close to the exact solution. As we shall detail in Sects. 2.2.2 and 2.3.2, the matching equation (115), supplemented by the condition of no-incoming radiation [say in the form of Eq. (60)], permits determining all the unknowns of the problem: On the one hand, the external multipolar decomposition $\mathcal{M}(\bar{h})$, i.e., the explicit expressions of the multipole moments therein (see Sects. 2.2.2 and 2.2.3); on the other hand, the terms in the inner post-Newtonian expansion $\bar{h}$ that are associated with radiation-reaction effects, i.e., those terms which depend on the boundary conditions of the radiative field at infinity, and which correspond in the present case
to a post-Newtonian source which is isolated from other sources in the Universe; see Sect. 2.3.2.

### 2.2.2 General expression of the multipole expansion

**Theorem 6** (Blanchet 1995, 1998b) Under the hypothesis of matching, Eq. (115), the multipole expansion of the solution of the Einstein field equation outside a post-Newtonian source reads

$$ \mathcal{M}(h^{\alpha\beta}) = \underset{B=0}{\mathrm{FP}} \Box_{\mathrm{ret}}^{-1} [\tilde{r}^B \mathcal{M}(\Lambda^{\alpha\beta})] - \frac{4G}{c^4} \sum_{\ell=0}^{+\infty} \frac{(-1)^\ell}{\ell!} \partial_L \left\{ \frac{1}{r} \mathcal{H}_L^{\alpha\beta} (t-r/c) \right\}, \quad (117) $$

where the “multipole moments” are given by

$$ \mathcal{H}_L^{\alpha\beta}(u) = \underset{B=0}{\mathrm{FP}} \int d^3\mathbf{x} \, \tilde{r}^B x_L \bar{\tau}^{\alpha\beta}(\mathbf{x}, u). \quad (118) $$

Here, $\bar{\tau}^{\alpha\beta}$ denotes the post-Newtonian expansion of the stress-energy pseudo-tensor in harmonic coordinates as defined by Eq. (53).

First notice where the physical restriction of considering a post-Newtonian source enters this theorem: The multipole moments (118) depend on the post-Newtonian expansion $\bar{\tau}^{\alpha\beta}$ of the pseudo-tensor, rather than on $\tau^{\alpha\beta}$ itself. Consider $\Delta^{\alpha\beta}$, namely the difference between $h^{\alpha\beta}$, which is a solution of the field equations everywhere inside and outside the source, and the first term in Eq. (117), namely the finite part of the retarded integral of the multipole expansion $\mathcal{M}(\Lambda^{\alpha\beta})$:

$$ \Delta^{\alpha\beta} \equiv h^{\alpha\beta} - \underset{B=0}{\mathrm{FP}} \Box_{\mathrm{ret}}^{-1} [\tilde{r}^B \mathcal{M}(\Lambda^{\alpha\beta})]. \quad (119) $$

Now the true solution $h^{\alpha\beta}$, subject to the past-stationarity condition (60), is given by the retarded integral of the pseudo-tensor $\tau^{\alpha\beta}$,

$$ h^{\alpha\beta} = \frac{16\pi G}{c^4} \Box_{\mathrm{ret}}^{-1} \tau^{\alpha\beta}, \quad (120) $$

where we recall that the retarded inverse d'Alembertian operator is defined by Eq. (71). So we have,

$$ \Delta^{\alpha\beta} = \frac{16\pi G}{c^4} \Box_{\mathrm{ret}}^{-1} \tau^{\alpha\beta} - \underset{B=0}{\mathrm{FP}} \Box_{\mathrm{ret}}^{-1} [\tilde{r}^B \mathcal{M}(\Lambda^{\alpha\beta})]. \quad (121) $$

In the second term the finite part plays a crucial role because the multipole expansion $\mathcal{M}(\Lambda^{\alpha\beta})$ is singular at $r = 0$. By contrast, the first term in Eq. (121), as it stands, is well-defined because we are considering only some smooth matter distribution: $\tau^{\alpha\beta} \in C^\infty(\mathbb{R}^4)$. There is no need to include the finite part regularization in the first term, but a *contrario* there is no harm to add one in front of it, because for convergent integrals the finite part simply gives back the value of the integral. The
advantage of adding artificially the finite part in the first term is that we can re-write Eq. (121) into the more interesting form

$$ \Delta^{\alpha\beta} = \frac{16\pi G}{c^4} \mathrm{FP}_{B=0} \Box_{\mathrm{ret}}^{-1} [\tilde{r}^B (\tau^{\alpha\beta} - \mathcal{M}(\tau^{\alpha\beta}))], \quad (122) $$

in which we have also used the fact that $\mathcal{M}(\Delta^{\alpha\beta}) = \frac{16\pi G}{c^4} \mathcal{M}(\tau^{\alpha\beta})$ because $T^{\alpha\beta}$ has a compact support. The interesting point about Eq. (122) is that $\Delta^{\alpha\beta}$ appears now to be the (finite part of a) retarded integral of a source with spatially compact support. This follows from the fact that the pseudo-tensor agrees numerically with its own multipole expansion when $r > a$ [by the same equation as Eq. (114)]. Therefore, $\mathcal{M}(\Delta^{\alpha\beta})$ can be obtained from the known formula for the multipole expansion of the retarded solution of a wave equation with compact-support source. This formula, given in Appendix B of Blanchet and Damour (1989), yields the second term in Eq. (117),

$$ \mathcal{M}(\Delta^{\alpha\beta}) = -\frac{4G}{c^4} \sum_{\ell=0}^{+\infty} \frac{(-1)^{\ell}}{\ell!} \partial_L \left\{ \frac{1}{r} \mathcal{H}_L^{\alpha\beta}(u) \right\}, \quad (123) $$

but in which the moments do not yet match the result (118); instead,

$$ \mathcal{H}_L^{\alpha\beta} = \mathrm{FP}_{B=0} \int d^3\mathbf{x} \tilde{r}^B x_L [\tau^{\alpha\beta} - \mathcal{M}(\tau^{\alpha\beta})]. \quad (124) $$

The reason is that we have not yet applied the assumption of a post-Newtonian source. Such sources are entirely covered by their own near zone (i.e., $a \ll \lambda$), and, in addition, for them the integral (124) has a compact support limited to the domain of the source. In consequence, we can replace the integrand in Eq. (124) by its post-Newtonian expansion, valid over all the near zone:

$$ \mathcal{H}_L^{\alpha\beta} = \mathrm{FP}_{B=0} \int d^3\mathbf{x} \tilde{r}^B x_L \left[ \bar{\tau}^{\alpha\beta} - \overline{\mathcal{M}(\tau^{\alpha\beta})} \right]. \quad (125) $$

Strangely enough, we do not get the expected result because of the presence of the second term in Eq. (125). Actually, this term is a bit curious, because the object $\overline{\mathcal{M}(\tau^{\alpha\beta})}$ it contains is only known in the form of the formal series whose structure is given by the first equality in Eq. (116) (indeed $\tau^{\alpha\beta}$ and $h^{\alpha\beta}$ have the same type of structure). Happily—because we would not know what to do with this term in practice—we are going to prove that the second term in Eq. (125) is identically zero. The proof is based on the properties of the analytic continuation as applied to the formal structure (116) of $\overline{\mathcal{M}(\tau^{\alpha\beta})}$. Each term of this series yields a contribution to Eq. (125) that takes the form, after performing the angular integration, of the integral $\mathrm{FP}_{B=0} \int_0^{+\infty} dr r^{B+b} (\ln r)^p$, and multiplied by some function of time. We want to prove that the radial integral $\int_0^{+\infty} dr r^{B+b} (\ln r)^p$ is zero by analytic continuation ($\forall B \in \mathbb{C}$). First we can get rid of the logarithms by considering some repeated differentiations with respect to $B$; thus we need only to consider the simpler integral $\int_0^{+\infty} dr r^{B+b}$. We split the integral into a “near-zone” integral $\int_0^R dr r^{B+b}$ and a “far-
zone” one $\int_{\mathbb{R}}^{+\infty} dr \, r^{B+b}$, where $\mathbb{R}$ is some constant radius. When $\Re(B)$ is a large enough positive number, the value of the near-zone integral is $\mathbb{R}^{B+b+1}/(B+b+1)$, while when $\Re(B)$ is a large negative number, the far-zone integral reads the opposite, $-\mathbb{R}^{B+b+1}/(B+b+1)$. Both obtained values represent the unique analytic continuations of the near-zone and far-zone integrals for any $B \in \mathbb{C}$ except $-b-1$. The complete integral $\int_{0}^{+\infty} dr \, r^{B+b}$ is equal to the sum of these analytic continuations, and is therefore identically zero ($\forall B \in \mathbb{C}$, including the value $-b-1$). At last we have completed the proof of Theorem 6:

$$ \mathcal{H}_L^{\alpha\beta} = \underset{B=0}{\mathrm{FP}} \int d^3\mathbf{x} \, \tilde{r}^B x_L \bar{\tau}^{\alpha\beta}. \qquad (126) $$

The latter proof makes it clear how crucial the analytic-continuation finite part regularization is, which we recall is the same as in our iteration of the exterior post-Minkowskian field [see Eq. (75)]. Without a finite part, the multipole moment (126) would be strongly divergent, because the pseudo-tensor $\bar{\tau}^{\alpha\beta}$ has a non-compact support owing to the contribution of the gravitational field, and the multipolar factor $x_L$ behaves like $r^\ell$ when $r \to +\infty$. The latter divergence has plagued the field of post-Newtonian expansions of gravitational radiation for many years. In applications such as in Sect. 3 of this article, we must carefully follow the rules for handling the FP regularization.

The two terms in the right-hand side of Eq. (117) depend separately on the length scale $r_0$ that we have introduced into the definition of the finite part, through the analytic-continuation factor $\tilde{r}^B = (r/r_0)^B$ introduced in Eq. (72). However, the sum of these two terms, i.e., the exterior multipolar field $M(h^{\alpha\beta})$ itself, is independent of $r_0$. To see this, the simplest way is to differentiate formally $M(h^{\alpha\beta})$ with respect to $r_0$; the differentiations of the two terms of Eq. (117) cancel each other. The independence of the field upon $r_0$ is quite useful in applications, since in general many intermediate calculations do depend on $r_0$, and only in the final stage does the cancellation of $r_0$ occurs.

### 2.2.3 The source multipole moments

In principle, the bridge between the exterior gravitational field generated by the post-Newtonian source and its inner field is provided by Theorem 6; however, we still have to make the connection with the explicit construction of the general multipolar-post-Minkowskian metric in Sect. 2.1. Namely, we must find the expressions of the six STF source multipole moments $I_L, J_L, \dots, Z_L$ parametrizing the linearized metric (64)–(66) at the basis of that construction. Recall that in actual applications we need mostly the mass-type moment $I_L$ and current-type one $J_L$, because the other moments parametrize a linearized gauge transformation and give contributions to high PN order.

To do this we first find the equivalent of the multipole expansion given in Theorem 6, which was parametrized by non-trace-free multipole functions $\mathcal{H}_L^{\alpha\beta}$, in
terms of new multipole functions $\mathcal{F}_L^{\alpha\beta}$ that are STF in all their indices $L$. The result is (introducing the intermediate notation $u^{\alpha\beta}$ for convenience)

$$u^{\alpha\beta} \equiv \underset{B=0}{\mathrm{FP}} \Box_{\mathrm{ret}}^{-1} [\tilde{r}^B \mathcal{M}(\Lambda^{\alpha\beta})], \qquad (127\mathrm{a})$$

$$\mathcal{M}(h^{\alpha\beta}) = u^{\alpha\beta} - \frac{4G}{c^4} \sum_{\ell=0}^{+\infty} \frac{(-1)^\ell}{\ell!} \partial_L \left\{ \frac{1}{r} \mathcal{F}_L^{\alpha\beta}(t-r/c) \right\}, \qquad (127\mathrm{b})$$

where the STF multipole functions (witness the multipolar factor $\hat{x}_L \equiv \mathrm{STF}[x_L]$ therein) read

$$\mathcal{F}_L^{\alpha\beta}(u) = \underset{B=0}{\mathrm{FP}} \int d^3\mathbf{x} \tilde{r}^B \hat{x}_L \int_{-1}^{1} dz \, \delta_{\ell}(z) \bar{\tau}^{\alpha\beta}(\mathbf{x}, u + zr/c). \qquad (128)$$

Notice the presence of an extra integration variable $z$, ranging from $-1$ to $1$. The $z$-integration involves the weighting function

$$\delta_{\ell}(z) = \frac{(2\ell + 1)!!}{2^{\ell+1}\ell!} (1 - z^2)^{\ell}, \qquad (129)$$

which approaches the Dirac delta-function (hence its name) in the limit of large multiplicities, i.e. $\lim_{\ell \to +\infty} \delta_{\ell}(z) = \delta(z)$, and is normalized such that

$$\int_{-1}^{1} dz \, \delta_{\ell}(z) = 1. \qquad (130)$$

The next step is to impose the harmonic-gauge conditions (50) onto the multipole decomposition (127b), hence $\mathcal{M}(\partial_{\mu}h^{\alpha\mu}) = 0$. This implies

$$\partial_{\mu} u^{\alpha\mu} = \frac{4G}{c^4} \sum_{\ell=0}^{+\infty} \frac{(-1)^{\ell}}{\ell!} \partial_L \left\{ \frac{1}{r} \mathcal{G}_L^{\alpha}(t-r/c) \right\}, \qquad (131)$$

with STF multipole functions therein given by

$$\mathcal{G}_L^{\alpha}(u) = \underset{B=0}{\mathrm{FP}} B \int d^3\mathbf{x} \, \tilde{r}^B \frac{x_i}{r^2} \hat{x}_L \int_{-1}^{1} dz \, \delta_{\ell}(z) \bar{\tau}^{\alpha i}(\mathbf{x}, u + zr/c). \qquad (132)$$

Note the explicit factor $B$ in front of the integral, so that $\mathcal{G}_L^{\alpha}$ is actually zero unless the integral develops a pole when $B \to 0$.

The final step is to decompose the multipole functions $\mathcal{F}_L^{\alpha\beta}(u)$ and $\mathcal{G}_L^{\alpha}(u)$ into STF irreducible pieces with respect to both $L$ and their spatial indices contained into $\alpha\beta = 00, 0i, ij$ or $\alpha = 0, i$. This technical part of the calculation is identical to the one of the STF irreducible multipole moments of linearized gravity (Damour and Iyer 1991a). Thus the formulas needed in this decomposition read

$$\mathcal{F}_L^{00} = R_L, \qquad (133\mathrm{a})$$
$$ \mathcal{F}_L^{0i} = {}^{(+)}T_{iL} + \epsilon_{ai<i_l} {}^{(0)}T_{L-1>a} + \delta_{i<i_l} {}^{(-)}T_{L-1>}, \quad (133b) $$

$$ \begin{aligned} \mathcal{F}_L^{ij} = {}^{(+2)}T_{ijL} + \text{STF} \sum_{ij} \text{STF} \left[ \epsilon_{aii_l} {}^{(+1)}U_{ajL-1} + \delta_{ii_l} {}^{(0)}U_{jL-1} \right. \\ \left. + \delta_{ii_l} \epsilon_{aji_{l-1}} {}^{(-1)}U_{aL-2} + \delta_{ii_l} \delta_{ji_{l-1}} {}^{(-2)}U_{L-2} \right] + \delta_{ij} V_L, \end{aligned} \quad (133c) $$

where the ten tensors $R_L, {}^{(+)}T_{L+1}, {}^{(0)}T_L, \dots, {}^{(-2)}T_{L-2}, V_L$ are STF, and can be uniquely given in terms of the $\mathcal{F}_L^{\alpha\beta}$'s by some inverse formulas. Similarly, with four supplementary STF tensors $P_L, \dots, {}^{(-1)}Q_{L-1}$, we have

$$ \mathcal{G}_L^0 = P_L, \quad (134a) $$

$$ \mathcal{G}_L^i = {}^{(+)}Q_{iL} + \epsilon_{ai<i_l} {}^{(0)}Q_{L-1>a} + \delta_{i<i_l} {}^{(-)}Q_{L-1>}. \quad (134b) $$

Finally, the latter decompositions yield the following.

**Theorem 7** (Blanchet 1998b) The STF multipole moments $I_L$ and $J_L$ of a post-Newtonian source are given, formally up to any post-Newtonian order, by $(\ell \ge 2)^{31}$

$$ I_L(u) = \text{FP}_{B=0} \int d^3\mathbf{x} \, \tilde{r}^B \int_{-1}^{1} dz \left\{ \delta_{\ell}(z) \hat{x}_L \bar{\Sigma} - \frac{4(2\ell+1)}{c^2(\ell+1)(2\ell+3)} \delta_{\ell+1}(z) \hat{x}_{iL} \bar{\Sigma}_i^{(1)} \right. \\ \left. + \frac{2(2\ell+1)}{c^4(\ell+1)(\ell+2)(2\ell+5)} \delta_{\ell+2}(z) \hat{x}_{ijL} \bar{\Sigma}_{ij}^{(2)} \right\} (\mathbf{x}, u + zr/c), \quad (135a) $$

$$ J_L(u) = \text{FP}_{B=0} \int d^3\mathbf{x} \, \tilde{r}^B \int_{-1}^{1} dz \, \epsilon_{ab\langle i_l} \left\{ \delta_{\ell}(z) \hat{x}_{L-1}{}_a \bar{\Sigma}_b \right. \\ \left. - \frac{2\ell+1}{c^2(\ell+2)(2\ell+3)} \delta_{\ell+1}(z) \hat{x}_{L-1}{}_{ac} \bar{\Sigma}_{bc}^{(1)} \right\} (\mathbf{x}, u + zr/c). \quad (135b) $$

These moments are the ones that are to be inserted into the linearized metric $h_{(1)}^{\alpha\beta}$ which represents the lowest approximation to the post-Minkowskian field $h_{\text{MPM}}^{\alpha\beta} = \sum_{n \ge 1} G^n h_{(n)}^{\alpha\beta}$ defined in Eq.(81).

In Eqs. (7) some convenient source densities are defined from the post-Newtonian expansion (denoted by an overbar) of the pseudo-tensor $\tau^{\alpha\beta}$ by

$$ \bar{\Sigma} = \frac{\bar{\tau}^{00} + \bar{\tau}^{ii}}{c^2}, \quad \bar{\Sigma}_i = \frac{\bar{\tau}^{0i}}{c}, \quad \bar{\Sigma}_{ij} = \bar{\tau}^{ij}, \quad (136) $$

(where $\bar{\tau}^{ii} \equiv \delta_{ij}\bar{\tau}^{ij}$). As indicated in Eqs. (7) all these quantities are to be evaluated at the spatial point $\mathbf{x}$ and at time $u + zr/c$. Since the source densities (136) have a non-compact spatial support (because of the gravitational field contribution in the pseudo-tensor) the FP regularization in (7) plays a crucial role in curing IR type divergences.

<footnote n="31">³¹ The monopole and dipole cases $\ell = 0, 1$ are discussed in Blanchet (1998b).</footnote>
The source moments $I_L(u)$ and $J_L(u)$ are physically valid for post-Newtonian sources and make sense only in the form of a post-Newtonian expansion, so in practice we need to know how to expand the z-integrals as series when $c \to +\infty$. Here is the appropriate formula:

$$ \int_{-1}^{1} dz \delta_{\ell}(z) \bar{\Sigma}(\mathbf{x}, u + zr/c) = \sum_{k=0}^{+\infty} \frac{(2\ell + 1)!!}{2^k k! (2\ell + 2k + 1)!!} \left(\frac{r}{c} \frac{\partial}{\partial u}\right)^{2k} \bar{\Sigma}(\mathbf{x}, u). \quad (137) $$

Since the right-hand side involves only even powers of $1/c$, the same result holds equally well for the advanced variable $u + zr/c$ or the retarded one $u - zr/c$. Of course, in the Newtonian limit, the moments $I_L$ and $J_L$ reduce to the standard Newtonian expressions. For instance, $I_{ij}(u) = Q_{ij}(u) + \mathcal{O}(c^{-2})$ recovers the Newtonian quadrupole moment (4); similarly for the current quadrupole $J_{ij} = C_{ij} + \mathcal{O}(c^{-2})$ and mass octupole $I_{ijk} = O_{ijk} + \mathcal{O}(c^{-2})$ appearing in Eq. (35).

An interesting fact is that at the 1PN order the mass type source moments $I_L$ extend up only over the *compact support* of the source, hence we do not need the FP regularization in this case. We have (Blanchet and Damour 1989)

$$ I_L(u) = \int d^3\mathbf{x} \left\{ \hat{x}_L \sigma + \frac{r^2 \hat{x}_L}{2c^2(2\ell+3)} \frac{\partial^2 \sigma}{\partial t^2} - \frac{4(2\ell+1)\hat{x}_{iL}}{c^2(\ell+1)(2\ell+3)} \frac{\partial \sigma_i}{\partial t} \right\} (\mathbf{x}, u) + \mathcal{O}\left(\frac{1}{c^4}\right), \quad (138) $$

where $\sigma$ and $\sigma_i$ (and later $\sigma_{ij}$) denote the compact-support parts of the source densities (136), i.e. defined from the matter stress-energy tensor $T^{\alpha\beta}$ alone:

$$ \sigma \equiv \frac{T^{00} + T^{ii}}{c^2}, \quad \sigma_i \equiv \frac{T^{0i}}{c}, \quad \sigma_{ij} \equiv T^{ij}, \quad (139) $$

(with $T^{ii} \equiv \delta_{ij}T^{ij}$). Only from the 2PN order will the mass moments acquire some non-compact supported contributions—i.e., with some integrals extending up to infinity (Blanchet 1995). Concerning the current type moments $J_L$, they are given at Newtonian order by

$$ J_L = \int d^3\mathbf{x} \epsilon_{ab<i_i} \hat{x}_{L-1>a} \sigma_b + \mathcal{O}\left(\frac{1}{c^2}\right), \quad (140) $$

and the non compactness appearing already at the next 1PN order (Damour and Iyer 1991b; Blanchet 1995).

For completeness, we give also the formulas for the four auxiliary source moments $W_L, \dots, Z_L$, which parametrize the gauge vector $\varphi^{\alpha}_{(1)}$ as defined in Eqs. (66):
$$ W_L(u) = \underset{B=0}{\mathrm{FP}} \int d^3\mathbf{x} \, \tilde{\mathbf{r}}^B \int_{-1}^1 \mathrm{d}z \left\{ \frac{2\ell + 1}{(\ell + 1)(2\ell + 3)} \delta_{\ell+1} \hat{x}_{iL} \bar{\Sigma}_i - \frac{2\ell + 1}{2c^2(\ell + 1)(\ell + 2)(2\ell + 5)} \delta_{\ell+2} \hat{x}_{ijL} \bar{\Sigma}_{ij}^{(1)} \right\}, \quad (141a) $$

$$ X_L(u) = \underset{B=0}{\mathrm{FP}} \int d^3\mathbf{x} \, \tilde{\mathbf{r}}^B \int_{-1}^1 \mathrm{d}z \left\{ \frac{2\ell + 1}{2(\ell + 1)(\ell + 2)(2\ell + 5)} \delta_{\ell+2} \hat{x}_{ijL} \bar{\Sigma}_{ij} \right\}, \quad (141b) $$

$$ Y_L(u) = \underset{B=0}{\mathrm{FP}} \int d^3\mathbf{x} \, \tilde{\mathbf{r}}^B \int_{-1}^1 \mathrm{d}z \left\{ -\delta_{\ell} \hat{x}_L \bar{\Sigma}_{ii} + \frac{3(2\ell + 1)}{(\ell + 1)(2\ell + 3)} \delta_{\ell+1} \hat{x}_{iL} \bar{\Sigma}_i^{(1)} - \frac{2(2\ell + 1)}{c^2(\ell + 1)(\ell + 2)(2\ell + 5)} \delta_{\ell+2} \hat{x}_{ijL} \bar{\Sigma}_{ij}^{(2)} \right\}, \quad (141c) $$

$$ Z_L(u) = \underset{B=0}{\mathrm{FP}} \int d^3\mathbf{x} \, \tilde{\mathbf{r}}^B \int_{-1}^1 \mathrm{d}z \, \epsilon_{ab\langle i\ell} \left\{ -\frac{2\ell + 1}{(\ell + 2)(2\ell + 3)} \delta_{\ell+1} \hat{x}_{L-1} \rangle_{bc} \bar{\Sigma}_{ac} \right\}. \quad (141d) $$

As discussed in Sect. 2.4.6, one can always find two intermediate “packages” of multipole moments, namely the canonical moments $M_L$ and $S_L$, which are some non-linear functionals of the source moments (7) and (141), and such that the exterior field depends only on them, modulo a change of coordinates. However, the canonical moments $M_L$, $S_L$ do not admit general closed-form expressions like (7)–(141).<fnref n="32" />

Needless to say, the formalism becomes prohibitively difficult to apply at very high post-Newtonian approximations. Some PN order being given, we must first compute the relevant relativistic corrections to the pseudo stress-energy-tensor $\bar{\tau}^{\alpha\beta}$; this necessitates solving the field equations inside the matter source, which we shall investigate in the next Sect. 2.3. Then $\bar{\tau}^{\alpha\beta}$ is to be inserted into the source moments (7) and (141), where the formula (137) permits expressing all the terms up to that post-Newtonian order by means of more tractable integrals extending over $\mathbb{R}^3$. Given a specific model for the matter source we then have to find a way to compute all these spatial integrals; this is done in Sect. 3.4.1 for the case of point-mass binaries. Next, we must substitute the source multipole moments into the linearized metric (64)–(66), and iterate them until all the necessary multipole interactions taking place in the radiative moments $U_L$ and $V_L$ are under control. We shall work out these multipole interactions for general sources in Sect. 2.4.5 up to the 4.5PN order for the flux of circular orbits, and 4PN order for the dominant gravitational-wave (2, 2) mode. Only at this point does one have the physical radiation field at infinity, from which we can build the templates for the detection and analysis of gravitational waves. We advocate here that the complexity of the formalism simply reflects the complexity of the Einstein field equations. It is probably impossible to devise a different formalism, valid for general sources a priori devoid of symmetries, that would be substantially simpler.

<footnote n="32">³² The work Blanchet et al. (2005b) provided some alternative expressions for all the multipole moments (7) and (141), useful for some applications, in the form of surface integrals extending on the outer part of the source’s near zone.</footnote>
## 2.3 Interior field of a post-Newtonian source

**Theorem 7** solves the question of the generation of gravitational waves by extended post-Newtonian matter sources. However, notice that this result has still to be completed by the precise procedure, i.e., an explicit “algorithm”, for the post-Newtonian iteration of the near-zone field, analogous to the multipolar-post-Minkowskian algorithm we defined in Sect. 2.1. Such procedure will permit the systematic computation of the source multipole moments, which contain the full post-Newtonian expansion of the pseudo-tensor $\bar{\tau}^{\alpha\beta}$, and of the radiation reaction effects occurring within the matter source.

Before proceeding, let us recall that the “standard” post-Newtonian approximation, as it was used until, say, the early 1980s,<fnref n="33" /> was plagued with some apparently inherent difficulties, which cropped up at some high post-Newtonian order. Historically, these difficulties, even appearing at higher approximations, have cast a doubt on the actual soundness, from a theoretical point of view, of the post-Newtonian expansion. Practically speaking, they posed the question of the reliability of the approximation, when comparing the theory’s predictions with very precise experimental results. This was one of the reasons for the “radiation-reaction controversy” raging at the time of the binary pulsar data (Ehlers et al. 1976; Walker and Will 1980).

One can distinguish two main issues:

1. The first issue is that in higher approximations some divergent Poisson-type integrals appear. Indeed the post-Newtonian expansion replaces the resolution of a hyperbolic-like d'Alembertian equation by a perturbatively equivalent hierarchy of elliptic-like Poisson equations. Rapidly it is found during the post-Newtonian iteration that the right-hand side of the Poisson equations acquires a non-compact support (it is distributed all over space $\mathbb{R}^3$), and that as a result the standard Poisson integral diverges at the bound of the integral at spatial infinity, i.e., when $r \equiv |x| \to +\infty$, with $t = \text{const}$—an IR problem;

2. The second issue is related to the limitation of the post-Newtonian approximation to the near zone—the region surrounding the source of small extent with respect to the wavelength of the emitted radiation: $r \ll \lambda$. As we have seen, the PN expansion assumes from the start that all retardations $r/c$ are small, so it can rightly be viewed as a formal *near-zone* expansion, when $r \to 0$. Note that the fact which makes the Poisson integrals to become typically divergent, namely that the coefficients of the post-Newtonian series blow up at spatial infinity, when $r \to +\infty$, has nothing to do with the actual behaviour of the field at infinity. However, the serious consequence is that it is a *priori* impossible to implement within the post-Newtonian scheme alone the physical information that the matter system is isolated from the rest of the Universe. Most importantly, the no-incoming radiation condition, imposed at past null infinity, cannot be taken directly into account, *a priori*, into the post-Newtonian scheme. In this sense the

<footnote n="33">³³ See notably Anderson and DeCanio (1975), Ehlers (1980), Kerlick (1980a, b), Papapetrou and Linet (1981), and also the earlier works by Plebański and Bazański (1959), Chandrasekhar (1965), Chandrasekhar and Nutku (1969), Chandrasekhar and Esposito (1970).</footnote>
post-Newtonian approximation is not “self-supporting”, because it necessitates
some information taken from outside its own domain of validity.

The divergencies are linked to the fact that the PN expansion is actually a singular perturbation, in the sense that the coefficients of the successive powers of $1/c$ are not uniformly valid in space, since they typically blow up at spatial infinity like some powers of $r$. For instance the post-Newtonian metric cannot be “asymptotically flat” starting at the 2PN or 3PN level, depending on the adopted coordinate system (Rendall 1992). The consequence is that the standard Poisson integrals are in general badly-behaving at infinity. Trying to solve the post-Newtonian equations by means of the Poisson integral does not make sense. However, this does not mean that there are no solutions to the problem, but simply that the Poisson integral does not constitute the appropriate solution of the Poisson equation in the PN context.

Here we review a solution of both issues (Poujade and Blanchet 2002; Blanchet et al. 2005c), in the form of a general expression for the near-zone gravitational field, developed to any post-Newtonian order, which is determined as a solution of the matching equation (115). This solution is free of the divergences of Poisson-type integrals (issue 1), and yields, in particular, some general expression, valid up to any order, for the terms associated with the gravitational radiation reaction force inside the PN source (issue 2).

Other methods are possible in order to resolve the problems of the PN approximation. Notably, an alternative solution to the problem of divergencies has been proposed by Futamase and Schutz (1983), Futamase (1983), based on the initial-value formulation. In this method the problem of the appearance of divergencies is avoided because of the finiteness of the causal region of integration, between the initial Cauchy hypersurface and the considered field point. On the other hand, a different approach to the problem of radiation reaction, which does not use a matching procedure, is to work only within a post-Minkowskian iteration scheme without expanding the retardations, and thus to naturally incorporate the correct boundary conditions at infinity (Chicone et al. 2001).

### 2.3.1 Post-Newtonian iteration in the near zone

We perform the post-Newtonian iteration of the field equations in harmonic coordinates in the near zone of an isolated matter distribution. We deal with a general hydrodynamical fluid, whose stress-energy tensor is smooth, i.e., $T^{\alpha\beta} \in C^{\infty}(\mathbb{R}^4)$. Thus the scheme *a priori* excludes the presence of singularities and black holes; these will be dealt with in Sect. 3 of this article.

We shall now prove that the post-Newtonian expansion can be *indefinitely* iterated without divergences. Like in Eq. (118) we denote by means of an overline the formal (infinite) post-Newtonian expansion of the field inside the source's near-zone. The general structure of the post-Newtonian expansion is denoted as