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

HAL
open science

# Gravitational Radiation from Post-Newtonian Sources and Inspiralling Compact Binaries

Luc Blanchet

► To cite this version:

Luc Blanchet. Gravitational Radiation from Post-Newtonian Sources and Inspiralling Compact Binaries. Living Reviews in Relativity, 2014, 17, 10.12942/lrr-2014-2. insu-03645235

HAL Id: insu-03645235

https://insu.hal.science/insu-03645235

Submitted on 25 Apr 2022

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Gravitational Radiation from Post-Newtonian Sources
and Inspiralling Compact Binaries

Luc Blanchet

Institut d'Astrophysique de Paris
98<sup>bis</sup> Boulevard Arago

75014 Paris, France

email: blanchet@iap.fr

http://www2.iap.fr/users/blanchet

Accepted: 27 January 2014

Published: 13 February 2014

Revised: 16 February 2016

(Update of lrr-2006-4)

Abstract

To be observed and analyzed by the network of gravitational wave detectors on ground (LIGO, VIRGO, etc.) and by the future detectors in space (eLISA, etc.), inspiralling compact binaries – binary star systems composed of neutron stars and/or black holes in their late stage of evolution – require high-accuracy templates predicted by general relativity theory. The gravitational waves emitted by these very relativistic systems can be accurately modelled using a high-order post-Newtonian gravitational wave generation formalism. In this article, we present the current state of the art on post-Newtonian methods as applied to the dynamics and gravitational radiation of general matter sources (including the radiation reaction back onto the source) and inspiralling compact binaries. We describe the post-Newtonian equations of motion of compact binaries and the associated Lagrangian and Hamiltonian formalisms, paying attention to the self-field regularizations at work in the calculations. Several notions of innermost circular orbits are discussed. We estimate the accuracy of the post-Newtonian approximation and make a comparison with numerical computations of the gravitational self-force for compact binaries in the small mass ratio limit. The gravitational waveform and energy flux are obtained to high post-Newtonian order and the binary’s orbital phase evolution is deduced from an energy balance argument. Some landmark results are given in the case of eccentric compact binaries – moving on quasi-elliptical orbits with non-negligible eccentricity. The spins of the two black holes play an important role in the definition of the gravitational wave templates. We investigate their imprint on the equations of motion and gravitational wave phasing up to high post-Newtonian order (restricting to spin-orbit effects which are linear in spins), and analyze the post-Newtonian spin precession equations as well as the induced precession of the orbital plane.

**Keywords:** Gravitational radiation, Post-Newtonian approximations, Multipolar expansion, Inspiralling compact binary

<figure><img src="image_3.png" /></figure>
Imprint / Terms of Use

*Living Reviews in Relativity* is a peer-reviewed open access journal published by the Springer International Publishing AG, Gewerbestrasse 11, 6330 Cham, Switzerland. ISSN 1433-8351.

This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Figures that have been previously published elsewhere may not be reproduced without consent of the original copyright holders.

Luc Blanchet,

“Gravitational Radiation from Post-Newtonian Sources and Inspiralling Compact Binaries”,
*Living Rev. Relativity*, **17**, (2014), 2.

DOI 10.12942/lrr-2014-2.
Article Revisions

*Living Reviews* supports two ways of keeping its articles up-to-date:

**Fast-track revision.** A fast-track revision provides the author with the opportunity to add short notices of current research results, trends and developments, or important publications to the article. A fast-track revision is refereed by the responsible subject editor. If an article has undergone a fast-track revision, a summary of changes will be listed here.

**Major update.** A major update will include substantial changes and additions and is subject to full external refereeing. It is published with a new publication number.

For detailed documentation of an article's evolution, please refer to the history document of the article's online version at http://dx.doi.org/10.12942/lrr-2014-2.

**13 February 2014:** Major revision, updated and expanded. About 180 new references. Improved introduction with added Section 1.2 on quadrupole formalism; More extensive description of asymptotic waveform and relation between radiative and source moments (Section 3.3); More complete discussion on radiation reaction effects (5.2) and explicit radiation reaction potentials (5.4); New Section 7.5 on the PN metric in the near zone; New sections on conservative dynamics: stability of orbits (8.2), first law of binary dynamics (8.3), comparison with self-force computations (8.4); New Section 9.5 on spherical-harmonic modes of gravitational waves; New sections on eccentric binaries: structure of the motion (10.1), quasi-Keplerian representation (10.2), averaged fluxes (10.3); New sections on spin effects: Lagrangian formalism (11.1), equations of motion (11.2), spin-orbit terms in the phase (11.3).

**15 January 2016:** Minor update. Added 13 new references. More accurate quotations on spins. Minor improvements in wording. Added new PN subdominant gravitational-wave modes and SF comparisons.

**16 February 2016:** Minor update. Added 12 new references. Added a discussion on 4PN equations of motion.
# Contents

<table><tr><td>1</td><td>Introduction</td><td>7</td></tr><tr><td>1.1</td><td>Analytic approximations and wave generation formalism</td><td>8</td></tr><tr><td>1.2</td><td>The quadrupole moment formalism</td><td>12</td></tr><tr><td>1.3</td><td>Problem posed by compact binary systems</td><td>14</td></tr><tr><td>1.4</td><td>Post-Newtonian equations of motion</td><td>17</td></tr><tr><td>1.5</td><td>Post-Newtonian gravitational radiation</td><td>20</td></tr></table>

## Part A: Post-Newtonian Sources 23

<table><tr><td>2</td><td>Non-linear Iteration of the Vacuum Field Equations</td><td>23</td></tr><tr><td>2.1</td><td>Einstein's field equations</td><td>23</td></tr><tr><td>2.2</td><td>Linearized vacuum equations</td><td>25</td></tr><tr><td>2.3</td><td>The multipolar post-Minkowskian solution</td><td>28</td></tr><tr><td>2.4</td><td>Generality of the MPM solution</td><td>31</td></tr><tr><td>2.5</td><td>Near-zone and far-zone structures</td><td>33</td></tr></table>

<table><tr><td>3</td><td>Asymptotic Gravitational Waveform</td><td>35</td></tr><tr><td>3.1</td><td>The radiative multipole moments</td><td>35</td></tr><tr><td>3.2</td><td>Gravitational-wave tails and tails-of-tails</td><td>37</td></tr><tr><td>3.3</td><td>Radiative versus source moments</td><td>42</td></tr></table>

<table><tr><td>4</td><td>Matching to a Post-Newtonian Source</td><td>46</td></tr><tr><td>4.1</td><td>The matching equation</td><td>46</td></tr><tr><td>4.2</td><td>General expression of the multipole expansion</td><td>48</td></tr><tr><td>4.3</td><td>Equivalence with the Will-Wiseman formalism</td><td>50</td></tr><tr><td>4.4</td><td>The source multipole moments</td><td>51</td></tr></table>

<table><tr><td>5</td><td>Interior Field of a Post-Newtonian Source</td><td>53</td></tr><tr><td>5.1</td><td>Post-Newtonian iteration in the near zone</td><td>55</td></tr><tr><td>5.2</td><td>Post-Newtonian metric and radiation reaction effects</td><td>57</td></tr><tr><td>5.3</td><td>The 3.5PN metric for general matter systems</td><td>59</td></tr><tr><td>5.4</td><td>Radiation reaction potentials to 4PN order</td><td>62</td></tr></table>

## Part B: Compact Binary Systems 65

<table><tr><td>6</td><td>Regularization of the Field of Point Particles</td><td>65</td></tr><tr><td>6.1</td><td>Hadamard self-field regularization</td><td>66</td></tr><tr><td>6.2</td><td>Hadamard regularization ambiguities</td><td>68</td></tr><tr><td>6.3</td><td>Dimensional regularization of the equations of motion</td><td>72</td></tr><tr><td>6.4</td><td>Dimensional regularization of the radiation field</td><td>76</td></tr></table>

<table><tr><td>7</td><td>Newtonian-like Equations of Motion</td><td>78</td></tr><tr><td>7.1</td><td>The 3PN acceleration and energy for particles</td><td>78</td></tr><tr><td>7.2</td><td>Lagrangian and Hamiltonian formulations</td><td>84</td></tr><tr><td>7.3</td><td>Equations of motion in the center-of-mass frame</td><td>90</td></tr><tr><td>7.4</td><td>Equations of motion and energy for quasi-circular orbits</td><td>94</td></tr><tr><td>7.5</td><td>The 2.5PN metric in the near zone</td><td>97</td></tr></table>
<table>
  <tbody>
    <tr>
      <td><b>8 Conservative Dynamics of Compact Binaries</b></td>
      <td><b>101</b></td>
    </tr>
    <tr>
      <td>8.1 Concept of innermost circular orbit</td>
      <td>101</td>
    </tr>
    <tr>
      <td>8.2 Dynamical stability of circular orbits</td>
      <td>103</td>
    </tr>
    <tr>
      <td>8.3 The first law of binary point-particle mechanics</td>
      <td>108</td>
    </tr>
    <tr>
      <td>8.4 Post-Newtonian approximation versus gravitational self-force</td>
      <td>112</td>
    </tr>
    <tr>
      <td><b>9 Gravitational Waves from Compact Binaries</b></td>
      <td><b>116</b></td>
    </tr>
    <tr>
      <td>9.1 The binary's multipole moments</td>
      <td>117</td>
    </tr>
    <tr>
      <td>9.2 Gravitational wave energy flux</td>
      <td>120</td>
    </tr>
    <tr>
      <td>9.3 Orbital phase evolution</td>
      <td>123</td>
    </tr>
    <tr>
      <td>9.4 Polarization waveforms for data analysis</td>
      <td>125</td>
    </tr>
    <tr>
      <td>9.5 Spherical harmonic modes for numerical relativity</td>
      <td>131</td>
    </tr>
    <tr>
      <td><b>10 Eccentric Compact Binaries</b></td>
      <td><b>135</b></td>
    </tr>
    <tr>
      <td>10.1 Doubly periodic structure of the motion of eccentric binaries</td>
      <td>135</td>
    </tr>
    <tr>
      <td>10.2 Quasi-Keplerian representation of the motion</td>
      <td>137</td>
    </tr>
    <tr>
      <td>10.3 Averaged energy and angular momentum fluxes</td>
      <td>140</td>
    </tr>
    <tr>
      <td><b>11 Spinning Compact Binaries</b></td>
      <td><b>146</b></td>
    </tr>
    <tr>
      <td>11.1 Lagrangian formalism for spinning point particles</td>
      <td>148</td>
    </tr>
    <tr>
      <td>11.2 Equations of motion and precession for spin-orbit effects</td>
      <td>152</td>
    </tr>
    <tr>
      <td>11.3 Spin-orbit effects in the gravitational wave flux and orbital phase</td>
      <td>157</td>
    </tr>
    <tr>
      <td><b>References</b></td>
      <td><b>161</b></td>
    </tr>
  </tbody>
</table>

List of Tables

<table>
  <tbody>
    <tr>
      <td>1</td>
      <td>Numerically determined value of the 3PN coefficient for the SF part of the redshift observable defined by Eq. (292).</td>
      <td>115</td>
    </tr>
    <tr>
      <td>2</td>
      <td>Numerically determined values of higher-order PN coefficients for the SF part of the redshift observable defined by Eq. (292).</td>
      <td>115</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Post-Newtonian contributions to the accumulated number of gravitational-wave cycles N<sub>cycle</sub> for compact binaries detectable in the bandwidth of LIGO-VIRGO detectors.</td>
      <td>125</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Spin-orbit contributions to the number of gravitational-wave cycles N<sub>cycle</sub> for binaries detectable by ground-based detectors LIGO-VIRGO.</td>
      <td>159</td>
    </tr>
  </tbody>
</table>
<figure><img src="image_4.png" /></figure>
# 1 Introduction

The theory of gravitational radiation from isolated sources, in the context of general relativity, is a fascinating science that can be explored by means of what was referred to in the XVIIIth century France as *l'analyse sublime*: The analytical (i.e., mathematical) method, and more specifically the resolution of partial differential equations. Indeed, the field equations of general relativity, when use is made of the harmonic-coordinate conditions, take the form of a quasi-linear hyperbolic differential system of equations, involving the famous wave operator or d'Alembertian [140]. The resolution of that system of equations constitutes a *problème bien posé* in the sense of Hadamard [236, 104], 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 con-
temporary astrophysical observations, of binary pulsars like the historical Hulse–Taylor pulsar
PSR 1913+16 [250], and, in a forthcoming future, of gravitational waves produced by massive and
rapidly evolving systems such as inspiralling compact binaries. These should be routinely detected
on Earth by the network of large-scale laser interferometers, including the advanced versions of
the ground-based interferometers LIGO and VIRGO, with also GEO and the future cryogenic
detector KAGRA. The first direct detection of a coalescence of two black holes has been achieved
on September 14, 2015 by the advanced LIGO detector [1]. Further ahead, the space-based laser
interferometer LISA (actually, the evolved version eLISA) should be able to detect supermassive
black-hole binaries at cosmological distances.

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.

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.

The flagship of approximation methods is the *post-Newtonian* approximation, which has been
developed from the early days of general relativity [303]. 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 *N* point-like objects at the first post-Newtonian approximation level [184], is 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 [148, 147, 143, 142], has been experimentally verified by the observation of the secular acceleration of the orbital motion of the Hulse–Taylor binary pulsar PSR 1913+16 [399, 400, 398];

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.

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 [359, 116, 21]. We shall
refer to other review articles like Refs. [233, 187] for discussions of 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.

Part A 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 post-Minkowskian and multipolar 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 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.

Part B 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 third post-Newtonian (3PN) 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 3.5PN
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 (in the simple case of circular orbits).

All over the article we try to state the main results in a form that is simple enough to be
understood without the full details; however, we also outline some of the proofs when they present
some interest on their own. To emphasize the importance of some key results, we present them in
the form of mathematical theorems. In applications we generally show the most up-to-date results
up to the highest known post-Newtonian order.<fnref n="1" />

1.1 Analytic approximations and wave generation formalism

The basic problem we face is to relate the asymptotic gravitational-wave form $h_{ij}$ generated by
some isolated source, 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="2" /> 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
in the distance to the source), and the perturbation in the small mass limit (ν-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

<footnote n="1">¹ A few errata have been published in this intricate field; all formulas take into account the latest changes.</footnote>

<footnote n="2">² In this article Greek indices αβ...μν... take space-time values 0, 1, 2, 3 and Latin indices ab...ij... spatial values 1, 2, 3. Cartesian coordinates are assumed throughout and boldface notation is often used for ordinary Euclidean vectors. In Section 11 upper Latin letters AB... will refer to tetrad indices 0, 1, 2, 3 with ab... the corresponding spatial values 1, 2, 3. Our signature is +2; hence the Minkowski metric reads ηαβ = diag(-1, +1, +1, +1) = ηAB. As usual G and c are Newton's constant and the speed of light.</footnote>
see Refs. [404, 142, 144, 145, 405, 421, 46, 52, 378]. For reviews on black-hole pertubations and the self-force approach see Refs. [348, 373, 177, 23].

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="3" /> 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\}, \qquad (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, following Chandrasekhar et al. [122, 124, 123], 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, still following Refs. [122, 124, 123], 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 $\frac{n}{2}\text{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)^4$.<fnref n="4" /> 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 Section 1.2 the utterly important “Newtonian” quadrupole moment formalism [183, 285]. Taking into account higher post-Newtonian corrections in a wave generation

<footnote n="3">³ Establishing the post-Newtonian expansion rigorously has been the subject of numerous mathematically oriented works, see e.g., [361, 362, 363].</footnote>

<footnote n="4">⁴ Note that for very eccentric binaries (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$ at the apastron of the orbit.</footnote>
formalism will mean including into the waveform the contributions of higher multipole moments, beyond the quadrupole. Post-Newtonian corrections of order $O(1/c^n)$ beyond the quadrupole formalism will still be denoted as $\frac{n}{2}$PN. 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 [219, 238, 384], 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. (4)] below to include the current-quadrupole and mass-octupole moments [332, 333], and obtained the corresponding formulas for linear momentum [332, 333, 30, 358] and angular momentum [339, 134]. The general structure of the infinite multipole series in the linearized theory was investigated by several works [369, 367, 343, 403], from which it emerged that the expansion is 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 [403]. The expressions of the multipole moments as integrals over the source, valid in the linearized theory but irrespective of a slow motion hypothesis, have been worked out in [309, 119, 118, 154]. In particular, Damour & Iyer [154] obtained the complete closed-form expressions for the time-dependent mass and spin multipole moments (in STF guise) of 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 $T - R/c = \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. [93] and Sachs [368]. The precise way under which such radiative space-times fall off asymptotically has been formulated geometrically by Penrose [337, 338] in the concept of an asymptotically simple space-time (see also Ref. [220]). The resulting 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 defined operationally by means of explicit integrals extending over the matter and gravitational fields.

An alternative way of defining the multipole expansion within the complete non-linear theory is that of Blanchet & Damour [57, 41], following pioneering works by Bonnor and collaborators [94, 95, 96, 251] and Thorne [403]. In this approach the basic multipole moments are the *source* moments, rather than the radiative ones. 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.<fnref n="5" /> 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 (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, we shall see that, in the current

<footnote n="5">⁵ Whereas, the direct attack of the post-Minkowskian expansion, valid at once inside and outside the source, faces some calculational difficulties [408, 136].</footnote>
situation, the closed-form expressions of the source multipole moments can be established only in the case where the source is post-Newtonian [44, 49]. 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 a standard application of the method of matched asymptotic expansions in general relativity [114, 113, 7, 357].

To be more precise, we shall show how a systematic multipolar and post-Minkowskian iteration scheme for the vacuum Einstein field equations yields the most general physically admissible solution of these equations [57]. The solution is specified once we give two and only two sets of time-varying (source) multipole moments. Some general theorems about the near-zone and far-zone expansions of that general solution will be stated. Notably, we show [41] that the asymptotic behaviour of the solution at future null infinity is in agreement with the findings of the Bondi-Sachs-Penrose [93, 368, 337, 338, 220] approach to gravitational radiation. However, checking that the asymptotic structure of the radiative field is correct is not sufficient by itself, because the ultimate aim, as we said, is to relate the far field to the properties of the source, and we are now obliged to ask: What are the multipole moments corresponding to a given stress-energy tensor $T^{\alpha\beta}$ describing the source? The general expression of the moments was obtained at the level of the second post-Newtonian (2PN) order in Ref. [44], and was subsequently proved to be in fact valid up to any post-Newtonian order in Ref. [49]. The source moments are given by some integrals extending over the post-Newtonian expansion of the total (pseudo) stress-energy tensor $\tau^{\alpha\beta}$, which is made of a matter part described by $T^{\alpha\beta}$ and a crucial non-linear gravitational source term $\Lambda^{\alpha\beta}$. These moments carry in front a particular operation of taking the finite part ($FP$ as we call it below), which makes them mathematically well-defined despite the fact that the gravitational part $\Lambda^{\alpha\beta}$ has a spatially infinite support, which would have made the bound of the integral at spatial infinity singular (of course the finite part is not added *a posteriori* to restore the well-definiteness of the integral, but is *proved* to be actually present in this formalism). The expressions of the moments had been obtained earlier at the 1PN level, albeit in different forms, in Ref. [59] for the mass-type moments [see Eq. (157a) below], and in Ref. [155] for the current-type ones.

The wave-generation formalism resulting from matching the exterior multipolar and post-Minkowskian field [57, 41] to the post-Newtonian source [44, 49] is able to take into account, in principle, any post-Newtonian correction to both the source and radiative multipole moments (for any multipolarity of the moments). The relationships between the radiative and source moments include many non-linear multipole interactions, because the source moments mix with each other as they “propagate” from the source to the detector. Such multipole interactions include the famous effects of wave tails, corresponding to the coupling between the non-static moments with the total mass M of the source. The non-linear multipole interactions have been computed within the present wave-generation formalism up to the 3.5PN order in Refs. [60, 50, 48, 74, 197]. Furthermore, the back-reaction of the gravitational-wave emission onto the source, up to the 1.5PN order relative to the leading order of radiation reaction, has also been studied within this formalism [58, 43, 47]. Now, recall that the leading radiation reaction force, which is quadrupolar, occurs already at the 2.5PN order in the source’s equations of motion. Therefore the 1.5PN “relative” order in the radiation reaction corresponds in fact to the 4PN order in the equations of motion, beyond the Newtonian acceleration. It has been shown that the gravitational-wave tails enter the radiation reaction at precisely the 1.5PN *relative* order, i.e., 4PN absolute order [58]. A systematic post-Newtonian iteration scheme for the near-zone field, formally taking into account all radiation reaction effects, has been obtained, fully consistent with the present formalism [357, 75].

A different wave-generation formalism has been devised by Will & Wiseman [424] (see also Refs. [422, 335, 336]), after earlier attempts by Epstein & Wagoner [185] and Thorne [403]. This formalism has exactly the same scope as the one of Refs. [57, 41, 44, 49], i.e., it applies to any isolated post-Newtonian sources, but it differs in the definition of the source multipole moments
and in many technical details when properly implemented [424]. In both formalisms, the moments are generated by the post-Newtonian expansion of the pseudo-tensor $\tau^{\alpha\beta}$, but in the Will-Wiseman formalism they are defined by some compact-support integrals terminating at some finite radius enclosing the source, e.g., the radius $\mathcal{R}$ of the near zone. By contrast, in Refs. [44, 49], the moments are given by some integrals covering the whole space ($\mathbb{R}^3$) and regularized by means of the finite part $\mathcal{F}\mathcal{P}$. Nevertheless, in both formalisms the source multipole moments, which involve a whole series of relativistic 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. We shall prove in Section 4.3 the complete equivalence, at the most general level, between the two formalisms.

## 1.2 The quadrupole moment formalism

The lowest-order wave generation formalism is the famous quadrupole formalism of Einstein [183] and Landau & Lifshitz [285]. This 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}^{\text{TT}}$ 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}^{\text{TT}} = \frac{2G}{c^4 R} \mathcal{P}_{ijab}(\mathbf{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 = |\mathbf{X}|$ is the distance to the source, $T - R/c$ is the retarded time, $\mathbf{N} = \mathbf{X}/R$ is the unit direction from the source to the far away observer, and $\mathcal{P}_{ijab} = \mathcal{P}_{ia}\mathcal{P}_{jb} - \frac{1}{2}\mathcal{P}_{ij}\mathcal{P}_{ab}$ is the TT projection operator, with $\mathcal{P}_{ij} = \delta_{ij} - N_i N_j$ being the projector onto the plane orthogonal to $\mathbf{N}$. The source’s quadrupole moment takes the familiar Newtonian form

$$Q_{ij}(t) = \int_{\text{source}} d^3\mathbf{x} \rho(\mathbf{x}, t) \left( x_i x_j - \frac{1}{3} \delta_{ij} \mathbf{x}^2 \right), \quad (3)$$

where $\rho$ is the Newtonian mass density. 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{dE}{dT} \right)^{\text{GW}} = \frac{G}{c^5} \left\{ \frac{1}{5} \frac{d^3 Q_{ab}}{dT^3} \frac{d^3 Q_{ab}}{dT^3} + \mathcal{O}\left(\frac{1}{c^2}\right) \right\}. \quad (4)$$

Our notation $\mathcal{F}$ stands for the total gravitational energy flux or gravitational “luminosity” of the source. Similarly, the total angular momentum flux is given by

$$\mathcal{G}_i \equiv \left( \frac{d\mathbf{J}_i}{dT} \right)^{\text{GW}} = \frac{G}{c^5} \left\{ \frac{2}{5} \epsilon_{iab} \frac{d^2 Q_{ac}}{dT^2} \frac{d^3 Q_{bc}}{dT^3} + \mathcal{O}\left(\frac{1}{c^2}\right) \right\}, \quad (5)$$

where $\epsilon_{abc}$ denotes the standard Levi-Civita symbol with $\epsilon_{123} = 1$.

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

<footnote n="6">⁶ The TT coordinate system can be extended to the near zone of the source as well; see for instance Ref. [282].</footnote>
pulsar and inspiralling compact binaries. At the position $(\mathbf{x}, t)$ in a particular coordinate system covering the source, the reaction force density can be written as [114, 113, 319]

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

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

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

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

The specific internal energy of the fluid is denoted $\Pi$, and obeys the usual thermodynamic 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 usual Eulerian equation of motion $\rho dv^i/dt = -\partial_i 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 (8a)$$

$$\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 (8b)$$

where one recognizes the fluxes at infinity given by Eqs. (4) and (5), 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 (see e.g., [103] for a discussion). Hence we obtain

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

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

where the brackets denote the time averaging over an orbit. These balance equations encode the secular decreases of energy and angular momentum by gravitational radiation emission.

The cardinal virtues of the Einstein-Landau-Lifshitz 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 [399, 400, 398]. Indeed let us apply the balance equations (9) to a system of two point masses moving on an eccentric orbit modelling the binary pulsar PSR 1913+16 – the classic references are [340, 339]; see also [186, 415]. We use the binary’s Newtonian energy and angular momentum,

$$E = -\frac{Gm_1 m_2}{2a}, \qquad (10a)$$
$$J = m_1 m_2 \sqrt{\frac{Ga(1 - e^2)}{m_1 + m_2}}, \qquad (10b)$$

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 (9a) we obtain first the secular evolution of a; next changing from a to the orbital period P using Kepler's third law,<fnref n="7" /> we get the secular evolution of the orbital period P as

$$\langle \frac{dP}{dt} \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 (11)$$

The last factor, depending on the eccentricity, comes out from the orbital average and is known as the Peters & Mathews [340] “enhancement” factor, 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 [399, 400, 398]. On the other hand the secular evolution of the eccentricity $e$ is deduced from the angular momentum balance equation (9b) [together with the previous result (11)], as

$$\langle \frac{de}{dt} \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 (12)$$

Interestingly, the system of equations (11)–(12) can be thoroughly integrated in closed analytic form. This yields the evolution of the eccentricity [339]:

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

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}$.

For a long while, it was thought that the various quadrupole formulas would be sufficient for sources of gravitational radiation to be observed directly on Earth – as they had proved to be amply sufficient in the case of the binary pulsar. However, further works [139]<fnref n="8" /> and [87, 138] showed that this is not true, as one has to include post-Newtonian corrections to the quadrupole formalism in order to prepare for the data analysis of future detectors, in the case of inspiralling compact binaries. From the beautiful test of the orbital decay (11) of the binary pulsar, we can say that the post-Newtonian corrections to the “Newtonian” quadrupole formalism – which we shall compute in this article – have already received a strong observational support.

## 1.3 Problem posed by compact binary systems

Inspiring compact binaries, containing neutron stars and/or black holes, are likely to become the bread-and-butter sources of gravitational waves for the detectors LIGO, VIRGO, GEO and KAGRA on ground, and also eLISA in space. The two compact objects steadily lose their orbital binding energy by emission of gravitational radiation; as a result, the orbital separation between them decreases, and the orbital frequency increases. Thus, the frequency of the gravitational-wave

<footnote n="7">⁷ 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 [319].</footnote>

<footnote n="8">⁸ This work entitled: “The last three minutes: Issues in gravitational-wave measurements of coalescing compact binaries” is sometimes coined the “3mn Caltech paper”.</footnote>
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 circularize the motion rapidly. 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. (12). 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 (13). 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 Kozai mechanism [283, 300] is one important scenario that produces eccentric binaries and involves the interaction between a pair of binaries in the dense cores of globular clusters [315]. If the mutual inclination angle of the inner binary is strongly tilted with respect to the outer compact star, 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 Section 10).

The main point about modelling the inspiralling compact binary is that a model made of two structureless point particles, characterized solely by two mass parameters $m_a$ and possibly two spins $S_a$ (with $a = 1, 2$ 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. The main justification for a model of point particles is that the effects due to the finite size of the compact bodies are small. Consider for instance the influence of the Newtonian quadrupole moments $Q_a$ induced by tidal interaction between two neutron stars. Let $a_a$ be the radius of the stars, and $r_{12}$ be the distance between the two centers of mass. We have, for tidal moments,

$$Q_1 = k_1 \frac{m_2 a_1^5}{r_{12}^3} \quad \text{and} \quad Q_2 = k_2 \frac{m_1 a_2^5}{r_{12}^3}, \qquad (14)$$

where $k_a$ are the star’s dimensionless (second) Love numbers [321], which depend on their internal structure, and are, typically, of the order unity. On the other hand, for compact objects, we can introduce their “compactness” parameters, defined by the dimensionless ratios

$$K_a \equiv \frac{Gm_a}{a_a c^2}, \qquad (15)$$

and equal $\sim 0.2$ for neutron stars (depending on their equation of state). The quadrupoles $Q_a$ will affect the Newtonian binding energy $E$ of the two bodies, and also the emitted total gravitational flux $\mathcal{F}$ as computed using the Newtonian quadrupole formula (4). 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 [276] 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 face each other at any instant during the inspiral), but for which the fluid motion has zero circulation in the inertial frame. In the Newtonian approximation, using the energy balance equation (9a), we find that within such a model (in the case of two identical neutron stars with same mass $m$, compactness $K$ and Love number $k$), the orbital phase reads

$$\phi^{\text{finite size}} - \phi_0 = -\frac{1}{8x^{5/2}} \left\{ 1 + \text{const} \, k \left( \frac{x}{K} \right)^5 \right\}, \qquad (16)$$
where “const” denotes a numerical coefficient of order one, $\phi_0$ is some initial phase, and $x \equiv (Gm\Omega/c^3)^{2/3}$ is a standard dimensionless post-Newtonian parameter of the order of $\sim 1/c^2$ (with $\Omega = 2\pi/P$ the orbital frequency). The first term in the right-hand side of Eq. (16) corresponds to the gravitational-wave damping of two point masses without internal structure; the second term is the finite-size effect, which appears as a relative correction, proportional to $(x/K)^5$, to the latter radiation damping effect. Because the finite-size effect is purely Newtonian, its relative correction $\sim (x/K)^5$ should not depend on the speed of light $c$; and indeed the factors $1/c^2$ cancel out in the ratio $x/K$. However, the compactness $K$ of neutron stars is of the order of 0.2 say, and by definition of compact objects we can consider that $K$ is formally of the order of unity or one half; therefore the factor $1/c^2$ it contains in (15) should not be taken into account when estimating numerically the effect. So the real order of magnitude of the relative contribution of the finite-size effect in Eq. (16) is given by the factor $x^5$ alone. This means that for compact objects the finite-size effect should roughly be comparable, numerically, to a post-Newtonian correction of magnitude $x^5 \sim 1/c^{10}$ namely 5PN order. This is a much higher post-Newtonian order than the one at which we shall investigate the gravitational effects on the phasing formula. Using $k \sim 1$, $K \sim 0.2$ and the bandwidth of detectors between 10 Hz and 1000 Hz, we find that the cumulative phase error due to the finite-size effect amounts to less that one orbital rotation over a total of $\sim 16\,000$ produced by the gravitational-wave damping of two neutron stars. The conclusion is that the finite-size effects can in general be neglected in comparison with purely gravitational-wave damping effects. 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 spins); hence this property has been called the “effacement” principle of general relativity [142]. But note that for non-compact or moderately compact objects (such as white dwarfs for instance) the Newtonian tidal interaction dominates over the radiation damping. The constraints on the nuclear equation of state and the tidal deformability of neutron stars which can be inferred from gravitational wave observations of neutron star binary inspirals have been investigated in Refs. [320, 200, 414]. For numerical computations of the merger of two neutron stars see Refs. [187, 249].

Inspiring compact binaries are ideally suited for application of a high-order post-Newtonian wave generation formalism. These systems are very relativistic, with orbital velocities as high as $0.5c$ in the last rotations (as compared to $\sim 10^{-3}c$ for the binary pulsar), so it is not surprising that the quadrupole-moment formalism (2)–(6) constitutes a poor description of the emitted gravitational waves, since many post-Newtonian corrections are expected to play a substantial role. This expectation has been confirmed by measurement-analyses [139, 137, 198, 138, 393, 346, 350, 284, 157], which have demonstrated that the post-Newtonian precision needed to implement successfully the optimal filtering technique for the LIGO/VIRGO detectors corresponds grossly, in the case of neutron-star binaries, to the 3PN approximation, or $1/c^6$ beyond the quadrupole moment approximation. Such a high precision is necessary because of the large number of orbital rotations that will be monitored in the detector’s frequency bandwidth, giving the possibility of measuring very accurately the orbital phase of the binary. Thus, the 3PN order is required mostly to compute the time evolution of the orbital phase, which depends, via Eq. (9a), on the center-of-mass binding energy $E$ and the total gravitational-wave energy flux $F$.

In summary, the theoretical problem is two-fold: On the one hand $E$, and on the other hand $F$, are to be computed with 3PN precision or better. To obtain $E$ we must control the 3PN equations of motion of the binary in the case of general, not necessarily circular, orbits; as for $F$ it necessitates the application of a 3PN wave generation formalism. It is remarkable that such high PN approximation is needed in preparation for the LIGO and VIRGO data analyses. As we shall see, the signal from compact binaries contains the signature of several non-linear effects which are specific to general relativity. We thus have the possibility of probing, experimentally, some aspects of the non-linear structure of Einstein’s theory [84, 85, 15, 14].
1.4 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 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 Lorentz & Droste [303].
Subsequently, Einstein, Infeld & Hoffmann [184] obtained the 1PN corrections for N particles 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 [201] (see also
Refs. [341, 330]).

The 2PN approximation was tackled by Ohta et al. [324, 327, 326, 325], who considered the post-Newtonian iteration of the Hamiltonian of *N* 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 & Deruelle [148, 147, 175, 141, 142], building on a non-linear (post-Minkowskian) iteration of the metric of two particles initiated in Ref. [31]. The corresponding result for the ADM-Hamiltonian of two particles at the 2PN order was given in Ref. [169] (see also Refs. [375, 376]). The 2.5PN equations of motion have also been derived in the case of two *extended* compact objects [280, 234]. The 2.5PN equations of two point masses as well as the near zone gravitational field in harmonic-coordinate were computed in Ref. [76].<fnref n="9" />

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 [148, 147, 175, 141, 142], have been used for the study of the radiation damping of the binary pulsar – its orbital $\dot{P}$ [142, 143, 173]. The result was in agreement with the prediction of the quadrupole formalism given by (11). 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 [142]. This is via the *effacing principle* for the internal structure of the compact bodies. As a result, the equations depend only on the “Schwarzschild” masses, $m_1$ and $m_2$, of the neutron stars. Notably their compactness parameters $K_1$ and $K_2$, defined by Eq. (15), do not enter the equations of motion. This has also been explicitly verified up to the 2.5PN order by Kopeikin et al. [280, 234], 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, Futamase & Asada [256, 257] in harmonic coordinates, using a variant (but, much simpler and more developed) of the surface-integral approach of Einstein et al. [184], 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 independently by several
groups, by means of different methods, and with equivalent results:

1. Jaranowski & Schäfer [261, 262, 263], and then with Damour [162, 164], employ the ADM-Hamiltonian canonical formalism of general relativity, following the line of research initiated in Refs. [324, 327, 326, 325, 169];

2. Blanchet & Faye [69, 71, 70, 72], and with de Andrade [174] and Iyer [79], founding their approach on the post-Newtonian iteration initiated in Ref. [76], compute directly the equations of motion (instead of a Hamiltonian) in harmonic coordinates;

3. Itoh & Futamase [255, 253] (see [213] for a review), continuing the surface-integral method

<footnote n="9">⁹ All the works reviewed in this section concern general relativity. However, let us mention here that the equations of motion of compact binaries in scalar-tensor theories are known up to 2.5PN order [318].</footnote>
of Refs. [256, 257], obtain the complete 3PN equations of motion in harmonic coordinates, without need of a self-field regularization;

4. Foffa & Sturani [203] derive the 3PN Lagrangian in harmonic coordinates within the effective field theory approach pioneered by Goldberger & Rothstein [223].

It has been shown [164, 174] that there exists a unique “contact” transformation of the dynamical variables that changes the harmonic-coordinates Lagrangian of Ref. [174] (identical to the ones issued from Refs. [255, 253] and [203]) into a new Lagrangian, whose Legendre transform coincides with the ADM-Hamiltonian of Ref. [162]. The equations of motion are therefore physically equivalent. For a while, however, they depended on one unspecified numerical coefficient, which is due to some incompleteness of the Hadamard self-field regularization method. This coefficient has been fixed by means of a better regularization, *dimensional regularization*, both within the ADM-Hamiltonian formalism [163], and the harmonic-coordinates equations of motion [61]. These works have demonstrated the power of dimensional regularization and its adequateness to the classical problem of interacting point masses in general relativity. By contrast, notice that, interestingly, the surface-integral method [256, 257, 255, 253] by-passes the need of a regularization. We devote Section 6 to questions related to the use of self-field regularizations.

The effective field theory (EFT) approach to the problems of motion and radiation of compact binaries, has been extensively developed since the initial proposal [223] (see [206] for a review). It borrows techniques from quantum field theory and consists of treating 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.<fnref n="10" /> 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. Like traditional approaches [163, 61] the EFT uses the dimensional regularization.

Computing the equations of motion and radiation field using Feynman diagrams in classical general relativity is not a new idea by itself: Bertotti & Plebanski [35] defined the diagrammatic tree-level perturbative expansion of the Green's functions in classical general relativity; Hari Dass & Soni [240]<fnref n="11" /> 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 [324, 327, 326, 325] and in scalar-tensor theories [151]. Nevertheless, the systematic EFT treatment has proved to be powerful and innovative for the field, e.g., with the introduction of a decomposition of the metric into “Kaluza–Klein type” potentials [277], the interesting link with the renormalization group equation [222], and the systematization of the computation of diagrams [203].

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 & Will [258, 259] for point-particle binaries in a general gauge, relying on energy and angular momentum balance equations and the known expressions of the 1PN fluxes. The latter works have been extended to 2PN order [226] and to include the leading spin-orbit effects [428]. The result has been then established from first principles (i.e., not relying on balance equations) in various works at 1PN order [260, 336, 278, 322, 254]. The 1PN radiation reaction force has also been obtained for general extended fluid systems in a particular gauge [43, 47]. Known also is the contribution of gravitational-wave tails in the

<footnote n="10">¹⁰ The effective action should be equivalent, in the tree-level approximation, to the Fokker action [207], for which the field degrees of freedom (i.e., the metric), that are solutions of the field equations derived from the original matter + field action with gauge-fixing term, have been inserted back into the action, thus defining the Fokker action for the sole matter fields.</footnote>

<footnote n="11">¹¹ This reference has an eloquent title: “Feynman graph derivation of the Einstein quadrupole formula”.</footnote>
equations of motion, which arises at the 4PN order and represents a 1.5PN modification of the radiation damping force [58]. This 1.5PN tail-induced correction to the radiation reaction force was also derived within the EFT approach [205, 215].

The state of the art on equations of motion is the 4PN approximation. Partial results on the equations of motion at the 4PN order have been obtained in [264, 265, 266] using the ADM Hamiltonian formalism, and in [204] using the EFT. The first derivation of the complete 4PN dynamics was accomplished in [166] by combining the local contributions [264, 265, 266] with a non-local contribution related to gravitational wave tails [58, 43], with the help of the result of an auxiliary analytical self-force calculation [36]. The non-local dynamics of [166] has been transformed in Ref. [167] into a local Hamiltonian containing an infinite series of even powers of the radial momentum. A second computation of the complete 4PN dynamics (including the same non-local interaction as in [166], but disagreeing on the local interaction) was accomplished in [33] using a Fokker Lagrangian in harmonic coordinates. Further works [39, 248] have given independent confirmations of the results of Refs. [166, 167]. More work is needed to understand the difference between the results of [166] and [33].

An important body of works concerns the effects of spins on the equations of motion of compact binaries. In this case we have in mind black holes rather than neutron stars, since astrophysical stellar-size black holes as well as super-massive galactic black holes have spins which can be close to maximal. The dominant effects are the spin-orbit (SO) coupling which is linear in spin, and the spin-spin (SS) coupling which is quadratic. For maximally spinning objects, and adopting a particular convention in which the spin is regarded as a 0.5PN quantity (see Section 11), the leading SO effect arises at the 1.5PN order while the leading SS effect appears at 2PN order. The leading SO and SS effects in the equations of motion have been determined by Barker & O'Connell [27, 28] and Kidder, Will & Wiseman [275, 271]. The next-to-leading SO effect, i.e., 1PN relative order corresponding to 2.5PN order, was obtained by Tagoshi, Ohashi & Owen [394], then confirmed and completed by Faye, Blanchet & Buonanno [194]. The results were also retrieved by two subsequent calculations, using the ADM Hamiltonian [165] and using EFT methods [292, 352]. The ADM calculation was later generalized to the N-body problem [241] and extended to the next-to-leading spin-spin effects (including both the coupling between different spins and spin square terms) in Refs. [387, 389, 388, 247, 243], and the next-to-next-to-leading SS interactions between different spins at the 4PN order [243]. In the meantime EFT methods progressed concurrently by computing the next-to-leading 3PN SS and spin-squared contributions [354, 356, 355, 293, 299], and the next-to-next-to-leading 4PN SS interactions for different spins [294] and for spin-squared [298]. Finally, 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 [242, 244] and in the harmonic-coordinates equations of motion [307, 90], with complete equivalence between the two approaches. Comparisons between the EFT and ADM Hamiltonian schemes for high-order SO and SS couplings can be found in Refs. [295, 299, 297]. We shall devote Section 11 to spin effects (focusing mainly on spin-orbit effects) in black hole binaries.

So far the status of post-Newtonian equations of motion is very satisfying. There is mutual agreement between all the results obtained by means of many different approaches and techniques, whenever they can be compared: point particles described by Dirac delta-functions or extended post-Newtonian fluids; surface-integrals methods; mixed post-Minkowskian and post-Newtonian expansions; direct post-Newtonian iteration and matching; EFT techniques versus traditional expansions; harmonic coordinates versus ADM-type coordinates; different processes or variants of the self-field regularization for point particles; different ways to including spins within the post-Newtonian approximation. In Part B of this article, we present complete results for the 3.5PN equations of motion (including the 1PN radiation reaction), and discuss the conservative part of the equations in the case of quasi-circular orbits. Notably, the conservative part of the dynamics is compared with numerical results for the gravitational self-force in Section 8.4.
## 1.5 Post-Newtonian gravitational radiation

The second problem, that of the computation of the gravitational waveform and the energy flux $F$, has to be solved by application of a wave generation formalism (see Section 1.1). The earliest computations at the 1PN level beyond the quadrupole moment formalism were done by Wagoner & Will [416], but based on some ill-defined expressions of the multipole moments [185, 403]. The computations were redone and confirmed by Blanchet & Schäfer [86] applying the rigorous wave generation formalism of Refs. [57, 60]. Remember that at that time the post-Newtonian corrections to the emission of gravitational waves had only a purely academic interest.

The energy flux of inspiralling compact binaries was then completed to the 2PN order by Blanchet, Damour & Iyer [64, 224], and, independently, by Will & Wiseman [424, 422], using their own formalism; see Refs. [66, 82] for joint reports of these calculations. The energy flux has been computed using the EFT approach in Ref. [221] with results agreeing with traditional methods.

At the 1.5PN order in the radiation field, appears the first contribution of “hereditary” terms, which are a priori sensitive to the entire past history of the source, i.e., which depend on all previous times up to $t \to -\infty$ in the past [60]. This 1.5PN hereditary term represents the dominant contribution of tails in the wave zone. It has been evaluated for compact binaries in Refs. [426, 87] by application of the formula for tail integrals given in Ref. [60]. Higher-order tail effects at the 2.5PN and 3.5PN orders, as well as a crucial contribution of tails generated by the tails themselves (the so-called “tails of tails”) at the 3PN order, were obtained in Refs [45, 48].

The 3PN approximation also involves, besides the tails of tails, many non-tail contributions coming from the relativistic corrections in the (source) multipole moments of the compact binary. Those have been almost completed in Refs. [81, 73, 80], in the sense that the result still involved one unknown numerical coefficient, due to the use of the Hadamard regularization. We shall review in Section 6 the computation of this parameter by means of dimensional regularization [62, 63], and shall present in Section 9 the most up-to-date results for the 3.5PN energy flux and orbital phase, deduced from the energy balance equation. In recent years all the results have been generalized to non-circular orbits, including both the fluxes of energy and angular momentum, and the associated balance equations [10, 9, 12]. The problem of eccentric orbits will be the subject of Section 10.

Besides the problem of the energy flux there is the problem of the gravitational waveform itself, which includes higher-order amplitude corrections and correlatively higher-order harmonics of the orbital frequency, consistent with the post-Newtonian order. Such full post-Newtonian waveform is to be contrasted with the so-called “restricted” post-Newtonian waveform which retains only the leading-order harmonic of the signal at twice the orbital frequency, and is often used in practical data analysis when searching the signal. However, for parameter estimation the full waveform is to be taken into account. For instance it has been shown that using the full waveform in the data analysis of future space-based detectors like eLISA will yield substantial improvements (with respect to the restricted waveform) of the angular resolution and the estimation of the luminosity distance of super-massive black hole binaries [16, 17, 410].

The full waveform has been obtained up to 2PN order in Ref. [82] by means of two independent wave generations (respectively those of Refs. [57, 44] and [424]), and it was subsequently extended up to the 3PN order in Refs. [11, 273, 272, 74]. At that order the signal contains the contributions of harmonics of the orbital frequency up to the eighth mode. The motivation is not only to build accurate templates for the data analysis of gravitational wave detectors, but also to facilitate the comparison and match of the high post-Newtonian prediction for the inspiral waveform with the numerically-generated waveforms for the merger and ringdown. For the latter application it is important to provide the post-Newtonian results in terms of a spin-weighted spherical harmonic decomposition suitable for a direct comparison with the results of numerical relativity. Recently the dominant quadrupole mode ($ℓ, m$) = (2, 2) in the spin-weighted spherical harmonic decomposition has been obtained at the 3.5PN order [197]. Available results will be provided in Sections 9.4
and 9.5.

At the 2.5PN order in the waveform appears the dominant contribution of another hereditary effect called the “non-linear memory” effect (or sometimes Christodoulou effect) [128, 427, 406, 60, 50]. This effect was actually discovered using approximation methods in Ref. [42] (see [60] for a discussion). It implies a permanent change in the wave amplitude from before to after a burst of gravitational waves, which can be interpreted as the contribution of gravitons in the known formulas for the linear memory for massless particles [99]. Note that the non-linear memory takes the form of a simple anti-derivative of an “instantaneous” term, and therefore becomes instantaneous (i.e., non-hereditary) in the energy flux which is composed of 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 the Newtonian level [427, 11]. It represents a part of the waveform whose amplitude steadily grows with time, but which is nearly constant over one orbital period. It is therefore essentially a zero-frequency effect (or DC effect), which has rather poor observational consequences in the case of the LIGO-VIRGO detectors, whose frequency bandwidth is always limited from below by some cut-off frequency $f_{\text{seismic}} > 0$. Non-linear memory contributions in the waveform of inspiralling compact binaries have been thoroughly computed by Favata [189, 192].

The post-Newtonian results for the waveform and energy flux are in complete agreement (up to the 3.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. This is the test-mass limit $v \to 0$, in which we define the symmetric mass ratio to be the reduced mass divided by the total mass, $v \equiv \mu/m$ such that $v = 1/4$ for equal masses. Linear black-hole perturbations, triggered by the geodesic motion of a small particle around the black hole, have been applied to this problem by Poisson [345] at the 1.5PN order (following the pioneering work [216]), by Tagoshi & Nakamura [393], using a numerical code up to the 4PN order, and by Sasaki, Tagoshi & Tanaka [372, 395, 397] (see also Ref. [316]), analytically up to the 5.5PN order. More recently the method has been improved and extended up to extremely high post-Newtonian orders: 14PN [209] and even 22PN [210] orders – but still for linear black-hole perturbations.

To successfully detect the gravitational waves emitted by spinning black hole binaries and to estimate the binary parameters, it is crucial to include spins effects in the templates, most importantly the spin-orbit effect which is linear in spins. The spins will affect the gravitational waves through a modulation of their amplitude, phase and frequency. Notably the orbital plane will precess in the case where the spins are not aligned or anti-aligned with the orbital angular momentum, see e.g., Ref. [8]. The leading SO and SS contributions in the waveform and flux of compact binaries are known from Refs. [275, 271, 314]; the next-to-leading SO terms at order 2.5PN were obtained in Ref. [53] after a previous attempt in [328]; the 3PN SO contribution is due to tails and was computed in Ref. [54], after intermediate results at the same order (but including SS terms) given in [353]. Finally, 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 in Refs. [89, 306]. The next-to-leading 3PN SS and spin-squared contributions in the radiation field were derived in Ref. [88]. In Section 11 we shall give full results for the contributions of spins (at SO linear level) in the energy flux and phase evolution up to 4PN order.

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 [313]. 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 the host galaxy and effectively terminate the process. Recoils could eject coalescing black holes from dwarf galaxies or globular clusters. 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.

Post-Newtonian methods are not ideally suited to compute the recoil of binary black holes
because most of the recoil is generated in the strong field regime close to the coalescence [199].
Nevertheless, after earlier computations of the dominant Newtonian effect [30, 199]<fnref n="12" /> and the 1PN
relative corrections [425], the recoil velocity has been obtained up to 2PN order for point particle
binaries without spin [83], and is also known for the dominant spin effects [271]. Various estimations
of the magnitude of the kick include a PN calculation for the inspiraling phase together with a
treatment of the plunge phase [83], an application of the effective-one-body formalism [152], a
close-limit calculation with Bowen-York type initial conditions [385], and a close-limit calculation
with initial PN conditions for the ringdown phase [288, 290].

In parallel the problem of gravitational recoil of coalescing binaries has attracted considerable
attention from the numerical relativity community. These computations led to increasingly accu-
rate estimates of the kick velocity from the merger along quasicircular orbits of binary black holes
without spins [115, 20] and with spins [117]. In particular these numerical simulations revealed
the interesting result that very large kick velocities can be obtained in the case of spinning black
holes for particular spin configurations.

<footnote n="12">¹² In absence of a better terminology, we refer to the leading-order contribution to the recoil as “Newtonian”, although it really corresponds to a 3.5PN subdominant radiation-reaction effect in the binary’s equations of motion.</footnote>
# Part A: Post-Newtonian Sources

## 2 Non-linear Iteration of the Vacuum Field Equations

### 2.1 Einstein's field equations

The field equations of general relativity are obtained by varying the space-time metric $g_{\alpha\beta}$ in the famous Einstein-Hilbert action,

$$I_{\text{EH}} = \frac{c^3}{16\pi G} \int d^4 x \sqrt{-g} R + I_{\text{mat}}[\Psi, g_{\alpha\beta}]. \qquad (17)$$

They form a system of ten second-order partial differential equations obeyed by the metric,

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

where the Einstein curvature tensor $E^{\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 I_{\text{mat}}/\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,

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

The matter equations can also be obtained by varying the matter action in (17) with respect to the matter fields $\Psi$. 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.

In most of this paper we adopt the conditions of harmonic coordinates, sometimes also called de Donder coordinates. We define, as a basic variable, the gravitational-field amplitude

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

where $g^{\alpha\beta}$ denotes the contravariant metric (satisfying $g^{\alpha\mu}g_{\mu\beta} = \delta_{\beta}^{\alpha}$), 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 (19) corresponding to the conservation of the matter tensor, reads<fnref n="13" />

$$\partial_{\mu} h^{\alpha\mu} = 0. \qquad (21)$$

Equation (21) 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 (21) 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

$$\partial_{\mu} h^{\alpha\mu} = \sqrt{-g} \Box_g x^{\alpha},$$

<footnote n="13">¹³ 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 tells that the coordinates $x^{\alpha}$ themselves, considered as scalars, are harmonic, i.e., obey the vacuum (curved) d'Alembertian equation.</footnote>
to think of the problem using approximation methods. Indeed, the metric $\eta_{\alpha\beta}$, originally introduced in the coordinate condition (21), 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 Einstein field equations in harmonic coordinates can be written in the form of inhomogeneous flat d’Alembertian equations,

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

where $\Box \equiv \Box_{\eta} = \eta^{\mu\nu}\partial_{\mu}\partial_{\nu}$. 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 (23) $$

The exact expression of $\Lambda^{\alpha\beta}$ in harmonic coordinates, including all non-linearities, reads<fnref n="14" />

$$ \begin{aligned} \Lambda^{\alpha\beta} = & -h^{\mu\nu}\partial^2_{\mu\nu}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^{\alpha\tau}\partial_{\mu}h^{\nu\lambda} + g_{\mu\nu}g^{\lambda\tau}\partial_{\lambda}h^{\alpha\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 (24) $$

As is clear from this 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. 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^{\alpha\beta} = N^{\alpha\beta}[h, h] + M^{\alpha\beta}[h, h, h] + L^{\alpha\beta}[h, h, h, h] + \mathcal{O}(h^5). \qquad (25) $$

These various terms can be straightforwardly computed from expanding Eq. (24); for instance the leading quadratic piece is explicitly given by<fnref n="15" />

$$ \begin{aligned} N^{\alpha\beta} = & -h^{\mu\nu}\partial^2_{\mu\nu}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^{\beta\mu} + \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} \qquad (26) $$

As we said, the condition (21) is equivalent to the matter equations of motion, in the sense of the conservation of the total pseudo-tensor $\tau^{\alpha\beta}$,

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

In this article, we shall look for approximate solutions of the field equations (21)–(22) 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 = |\mathbf{x}|$ is the harmonic-coordinate radial distance. Outside the domain of the source, when $r > a$, the gravitational source term, according to Eq. (27), is divergence-free,

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

<footnote n="14">¹⁴ In $d+1$ space-time dimensions, only one coefficient in this expression is modified; see Eq. (175) below.</footnote>

<footnote n="15">¹⁵ See Eqs. (3.8) in Ref. [71] for the cubic and quartic terms. We denote e.g., $h_{\mu}^{\alpha} = \eta_{\mu\nu}h^{\alpha\nu}$, $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>
2. The matter distribution inside the source is smooth: $T^{\alpha\beta} \in C^{\infty}(\mathbb{R}^3)$.<fnref n="16" /> We have in mind a smooth hydrodynamical fluid system, without any singularities nor shocks (a priori), that is described by some Euler-type equations including high relativistic corrections. In particular, we exclude from the start the presence of any black holes; however, we shall return to this question in Part B 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 $-\mathcal{T}$ in the past, namely

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

The latter condition is a means to impose, by brute force, the famous no-incoming radiation condition, 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 (see Section 3.2) that our condition of stationarity in the past (29), although weaker than the ideal no-incoming radiation condition, does not entail any physical restriction on the general validity of the formulas we derive. Even more, the condition (29) 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 $-\mathcal{T}$ can be set right after the explosions of the two supernovæ yielding the formation of the compact binary system.

Subject to the past-stationarity condition (29), the differential equations (22) can be written equivalently into the form of the integro-differential equations

$$ h^{\alpha\beta} = \frac{16\pi G}{c^4} \Box_{\text{ret}}^{-1} \tau^{\alpha\beta}, \qquad (30) $$

containing the usual retarded inverse d'Alembertian integral operator, given by

$$ (\Box_{\text{ret}}^{-1} \tau)(\mathbf{x}, t) = -\frac{1}{4\pi} \iiint \frac{d^3\mathbf{x}'}{|\mathbf{x} - \mathbf{x}'|} \tau(\mathbf{x}', t - |\mathbf{x} - \mathbf{x}'|/c), \qquad (31) $$

extending over the whole three-dimensional space $\mathbb{R}^3$.

## 2.2 Linearized vacuum equations

In what follows we solve the field equations (21)–(22), 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. At the linearized level (or first-post-Minkowskian approximation), we write:

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

<footnote n="16">¹⁶ $\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>
where the subscript “ext” reminds us that the solution is valid 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. (32) 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, the harmonic-coordinate metric coefficient $h_{(1)}^{\alpha\beta}$ satisfies

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

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

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 and collaborators [94, 95, 96, 251], 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 [403], and looking for a definite algorithm for the approximation scheme [57]. The solution of the system of equations (33) takes the form of a series of retarded multipolar waves<fnref n="17" />

$$ h_{(1)}^{\alpha\beta} = \sum_{\ell=0}^{+\infty} \partial_L \left( \frac{K_L^{\alpha\beta}(t-r/c)}{r} \right), \qquad (34) $$

where $r = |\mathbf{x}|$, and where the functions $K_L^{\alpha\beta} \equiv K_{i_1...i_\ell}^{\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. (29). Since a monopolar wave satisfies $\Box(K_L(u)/r) = 0$ and the d’Alembertian commutes with the multi-derivative $\partial_L$, it is evident that Eq. (34) represents the most general solution of the wave equation (33a); but see Section 2 in Ref. [57] for a rigorous proof based on the Euler-Poisson-Darboux equation. The gauge condition (33b), 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 (33b)

<footnote n="17">¹⁷ Our notation is the following: $L = i_1 i_2 \dots i_\ell$ denotes a multi-index, made of $\ell$ (spatial) indices. Similarly, we write for instance $K = j_1 \dots j_k$ (in practice, we generally do not need to write explicitly the “carrier” letter $i$ or $j$), or $aL - 1 = a i_1 \dots i_{\ell-1}$. Always understood in expressions such as Eq. (34) 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_{i_1} \dots \partial_{i_\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_{...i_p...i_q...} = K_{...i_q...i_p...}$ and that $\delta_{i_p i_q} K_{...i_p...i_q...} = 0$ (see Ref. [403] and Appendices A and B in Ref. [57] 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_{\langle L \rangle}$. In particular, $\hat{n}_L = n_{\langle L \rangle}$ is the STF projection of the product of unit vectors $n_L = n_{i_1} \dots n_{i_\ell}$, for instance $\hat{n}_{ij} = n_{\langle ij \rangle} = n_{ij} - \frac{1}{3}\delta_{ij}$ and $\hat{n}_{ijk} = n_{\langle ijk \rangle} = n_{ijk} - \frac{1}{5}(\delta_{ij}n_k + \delta_{ik}n_j + \delta_{jk}n_i)$; an expansion into STF tensors $\hat{n}_L = \hat{n}_L(\theta, \phi)$ is equivalent to the usual expansion in spherical harmonics $Y_{lm} = Y_{lm}(\theta, \phi)$, see Eqs. (75) below. Similarly, we denote $x_L = x_{i_1} \dots x_{i_\ell} = r^l n_L$ where $r = |\mathbf{x}|$, and $\hat{x}_L = x_{\langle L \rangle} = \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>
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 [403, 57].

**Theorem 1.** The most general solution of the linearized field equations (33) outside some time-like world tube enclosing the source ($r > a$), and stationary in the past [see Eq. (29)], reads

$$h_{(1)}^{\alpha\beta} = k_{(1)}^{\alpha\beta} + \partial^{\alpha}\varphi_{(1)}^{\beta} + \partial^{\beta}\varphi_{(1)}^{\alpha} - \eta^{\alpha\beta}\partial_{\mu}\varphi_{(1)}^{\mu}. \qquad (35)$$

The first term 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

$$k_{(1)}^{00} = -\frac{4}{c^2} \sum_{\ell \ge 0} \frac{(-)^\ell}{\ell!} \partial_L \left( \frac{1}{r} I_L(u) \right), \qquad (36a)$$

$$k_{(1)}^{0i} = \frac{4}{c^3} \sum_{\ell \ge 1} \frac{(-)^\ell}{\ell!} \left\{ \partial_{L-1} \left( \frac{1}{r} I_{iL-1}^{(1)}(u) \right) + \frac{\ell}{\ell+1} \epsilon_{iab} \partial_{aL-1} \left( \frac{1}{r} J_{bL-1}(u) \right) \right\}, \qquad (36b)$$

$$k_{(1)}^{ij} = -\frac{4}{c^4} \sum_{\ell \ge 2} \frac{(-)^\ell}{\ell!} \left\{ \partial_{L-2} \left( \frac{1}{r} I_{ijL-2}^{(2)}(u) \right) + \frac{2\ell}{\ell+1} \partial_{aL-2} \left( \frac{1}{r} \epsilon_{ab(i} J_{j)bL-2}^{(1)}(u) \right) \right\}. \qquad (36c)$$

The other terms represent 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. (37)].

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$,<fnref n="18" /> 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 radiation.

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 components of the gauge-vector $\varphi_{(1)}^{\alpha}$ entering Eq. (35):

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

$$\varphi_{(1)}^{i} = -\frac{4}{c^4} \sum_{\ell \ge 0} \frac{(-)^\ell}{\ell!} \partial_{iL} \left( \frac{1}{r} X_L(u) \right) \qquad (37b)$$

$$-\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 (37c)$$

<footnote n="18">¹⁸ 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. (29). Hence, $P_i = 0$.</footnote>
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 Section 2.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 Section 2.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 Eqs. (97)–(98)]. The moments $M_L, S_L$ are called “canonical” moments; they play a useful role in intermediate calculations. All the multipole moments $I_L, J_L, W_L, X_L, Y_L, Z_L$ will be computed in Section 4.4.

## 2.3 The multipolar post-Minkowskian solution

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{ext}}^{\alpha\beta} = \sum_{n=1}^{+\infty} G^n h_{(n)}^{\alpha\beta}, \qquad (38)$$

where the first term is composed of the result given by Eqs. (35)–(37). In this article, we shall always understand the infinite sums such as the one in Eq. (38) 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 of that point (see, in the present context, Refs. [130, 171, 361, 362, 363]).

We substitute the post-Minkowski ansatz (38) into the vacuum Einstein field equations (21)–(22), i.e., with $\tau^{\alpha\beta}$ simply given by the gravitational source term $\Lambda^{\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 (39a)$$

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

The right-hand side of the wave equation (39a) 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 (39a) reads

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

$$\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 (40b)$$

$$\Box h_{(4)}^{\alpha\beta} = L^{\alpha\beta}[h_{(1)}, h_{(1)}, h_{(1)}, h_{(1)}] \\ + M^{\alpha\beta}[h_{(1)}, h_{(1)}, h_{(2)}] + M^{\alpha\beta}[h_{(1)}, h_{(2)}, h_{(1)}] + M^{\alpha\beta}[h_{(2)}, h_{(1)}, h_{(1)}] \\ + N^{\alpha\beta}[h_{(2)}, h_{(2)}] + N^{\alpha\beta}[h_{(1)}, h_{(3)}] + N^{\alpha\beta}[h_{(3)}, h_{(1)}], \qquad (40c)$$

...

...
The quadratic, cubic and quartic pieces of $\Lambda^{\alpha\beta}$ are defined by Eq. (25)–(26).

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. (39a), $\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. (39a) 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, following the method of Ref. [57], 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}$. 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 (35)–(37) in which the multipolar sums are actually finite.<fnref n="19" /> 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,

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

where the symbol $\Box_{\text{ret}}^{-1}$ stands for the retarded integral defined by Eq. (31). It is convenient to introduce inside the regularizing factor 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}. \qquad (42)$$

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 $I^{\alpha\beta}(B)$ on the complex plane, which was originally defined only when $\Re(B) > k_{\max} - 3$, admits a unique analytic continuation to all values of $B \in \mathbb{C}$ except at some integer values. Furthermore, the analytic continuation of $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 $I^{\alpha\beta}(B)$, when $B \to 0$, in the form

$$I^{\alpha\beta}(B) = \sum_{p=p_0}^{+\infty} \ell_p^{\alpha\beta} B^p, \qquad (43)$$

<footnote n="19">¹⁹ This assumption is justified because we are ultimately interested in the radiation field at some given finite post-Newtonian precision like 3PN, 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 (35)–(37), there is a maximal multipolarity $\ell_{\max}(n)$ at any post-Minkowskian order $n$, which grows linearly with $n$.</footnote>
where $p \in \mathbb{Z}$, and the various coefficients $\iota_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 $n$, refers to the maximal order of these poles. By applying the d'Alembertian operator onto both sides of Eq. (43), and equating the different powers of $B$, we arrive at

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

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

As we see, the case $p = 0$ shows that the finite-part coefficient in Eq. (43), namely $\iota_0^{\alpha\beta}$, is a particular solution of the requested equation: $\Box \iota_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 $\iota_0^{\alpha\beta}$, and denote, from now on, the latter solution by $u_{(n)}^{\alpha\beta} \equiv \iota_0^{\alpha\beta}$, or, in more explicit terms,

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

where the finite-part symbol $\mathcal{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 (39b); its divergence, say $w_{(n)}^{\alpha} = \partial_{\mu} u_{(n)}^{\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. (28)], we find instead

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

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. (46) 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 (from 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(u)$, $P_L(u)$, $Q_L(u)$ and $R_L(u)$, such that

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

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

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 simply 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)} + C_a^{(-2)} - 3P_a \right) \right], \qquad (48a)$$
$$v_{(n)}^{0i} = r^{-1}(-Q_{i}^{(-1)} + 3P_{i}^{(1)}) - \epsilon_{iab}\partial_{a}[r^{-1}R_{b}^{(-1)}] - \sum_{l=2}^{+\infty}\partial_{L-1}[r^{-1}N_{iL-1}], \quad (48b)$$

$$v_{(n)}^{ij} = -\delta_{ij}r^{-1}P + \sum_{l=2}^{+\infty} \left\{ \begin{aligned} & 2\delta_{ij}\partial_{L-1}[r^{-1}P_{L-1}] - 6\partial_{L-2}(i[r^{-1}P_{j}])_{L-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\}. \quad (48c)$$

Notice the presence of anti-derivatives, denoted e.g., by $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. (48) 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$. Finally, if we pose

$$h_{(n)}^{\alpha\beta} = u_{(n)}^{\alpha\beta} + v_{(n)}^{\alpha\beta}, \quad (49)$$

we see that we solve at once the d'Alembertian equation (39a) and the coordinate condition (39b). 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 (38) based on the linearized approximation $h_{(1)}^{\alpha\beta}$ given by Eqs. (35)-(37). The previous procedure constitutes an algorithm, which can be (and has recently been [74, 197]) implemented by an algebraic computer programme. Again, note that this algorithm permits solving the full Einstein field equations together with the gauge condition (i.e., not only the relaxed field equations).

## 2.4 Generality of the MPM solution

We have a solution, but is that a general solution? The answer, “yes”, is provided by the following result [57].

**Theorem 2.** 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{ext}}^{\alpha\beta} = \sum_{n=1}^{+\infty} G^n h_{(n)}^{\alpha\beta} [I_L, J_L, \dots, Z_L]. \quad (50)$$

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. (36), and on four supplementary functions $W_L(u), \dots, Z_L(u)$ parametrizing the gauge vector (37).

The proof is quite easy. With Eq. (49) we obtained a particular solution of the system of equations (39). To it we should add the most general solution of the corresponding homogeneous system of equations, which is obtained by setting $\Lambda_{(n)}^{\alpha\beta} = 0$ into Eqs. (39). But this homogeneous system of equations is nothing but the linearized vacuum field equations (33), to which we know the most general solution $h_{(1)}^{\alpha\beta}$ given by Eqs. (35)-(37). 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, \quad (51a)$$
$$Z_{L}^{\text{new}} = Z_{L} + G^{n-1} \delta Z_{L}. \qquad (51b)$$

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.

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. However, as we now discuss, it is always possible to find two, and only two, sets of 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 = \varphi_{(1)}^\alpha$, where $\varphi_{(1)}^\alpha$ is the gauge vector given by Eqs. (37). So, at the linearized level, we have only the two types of moments $M_L$ and $S_L$, parametrizing $k_{(1)}^{\alpha\beta}$ by the same formulas as in Eqs. (36). We must be careful to denote these moments with names different from $I_L$ and $J_L$ because they will ultimately correspond to a different physical source. Then we apply exactly the same post-Minkowskian algorithm, following the formulas (45)–(49) as we did above, but starting from the gauge-transformed linear metric $k_{(1)}^{\alpha\beta}$ instead of $h_{(1)}^{\alpha\beta}$. The result of the iteration is therefore some

$$k_{\text{ext}}^{\alpha\beta} = \sum_{n=1}^{+\infty} G^n k_{(n)}^{\alpha\beta} [M_L, S_L]. \qquad (52)$$

Obviously this post-Minkowskian algorithm yields some simpler calculations as we have only two multipole moments to iterate. The point is that one can show that the resulting metric (52) is *isometric* to the original one (50) if and only if the so-called canonical moments $M_L$ and $S_L$ are related to the source moments $I_L, J_L, \dots, Z_L$ by some (quite involved) non-linear equations. We shall give in Eqs. (97)–(98) the most up to date relations we have between these moments. Therefore, the most general solution of the field equations, modulo a coordinate transformation, can be obtained by starting from the linearized metric $k_{(1)}^{\alpha\beta}$ [$M_L, S_L$] instead of the more complicated $k_{(1)}^{\alpha\beta}$ [$I_L, J_L$] + $\partial^\alpha \varphi_{(1)}^\beta$ + $\partial^\beta \varphi_{(1)}^\alpha$ - $\eta^{\alpha\beta} \partial_\mu \varphi_{(1)}^\mu$, and continuing the post-Minkowskian calculation.

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 6 below 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$. 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 because they yield a simpler post-Minkowskian iteration. Then, one can generally come back to the more fundamental source-rooted moments by using the fact that $M_L$ and $S_L$ differ from the corresponding $I_L$ and $J_L$ only by high-order post-Newtonian terms like 2.5PN; see Eqs. (97)–(98) below. Indeed, this is to be expected because the physical difference between both types of moments stems only from non-linearities.
## 2.5 Near-zone and far-zone structures

In our presentation of the post-Minkowskian algorithm (45)–(49) 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 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)}$ has a certain structure (here we often omit the space-time indices $\alpha\beta$); this hypothesis is established as a theorem once the mathematical induction succeeds.

**Theorem 3.** 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},^{20}$

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

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 the source multipole moments $I_L, \dots, Z_L$.

For the proof see Ref. [57]. As we see, the near-zone expansion involves, besides the simple powers of $r$, some 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. (53) 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 necessarily of the type

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

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$ [57]. It is remarkable that there is no more complicated structure like for instance $\ln(\ln c)$.

Paralleling the structure of the near-zone expansion, we have a similar result concerning the structure of 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)}(\mathbf{x}, t) = \sum \frac{\hat{n}_L (\ln r)^p}{r^k} G_{L,k,p,n}(u) + o\left(\frac{1}{r^N}\right), \quad (55)$$

where $k, p \in \mathbb{N}$, with $1 \le k \le N$, and where, likewise in the near-zone expansion (53), some 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 [202]. One knows also that this is a coordinate effect, because the study of the “asymptotic” structure of space-time at future null infinity by Bondi et al. [93], Sachs [368], and Penrose [337, 338], has revealed the existence of other coordinate systems that avoid the appearance of any logarithms: the so-called *radiative* coordinates, 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 [252, 304, 41]. The following theorem, proved in Ref. [41], shows that our general construction of the metric in the exterior of the source, when developed at future null infinity, is consistent with the Bondi–Sachs–Penrose [93, 368, 337, 338] approach to gravitational radiation.

<footnote n="20">20 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>
**Theorem 4.** The most general multipolar-post-Minkowskian solution, stationary in the past [see Eq. (29)], 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_{(n)}(\mathbf{X}, T) = \sum \frac{\hat{N}_L}{R^k} K_{L,k,n}(U) + \mathcal{O}\left(\frac{1}{R^N}\right). \qquad (56)$$

The functions $K_{L,k,n}$ are computable functionals of the source multipole moments. In radiative coordinates the retarded time $U$ is a null coordinate in the asymptotic limit. The metric $H_{\text{ext}}^{\alpha\beta} = \sum_{n \ge 1} G^n H_{(n)}^{\alpha\beta}$ is asymptotically simple in the sense of Penrose [337, 338, 220], perturbatively to any post-Minkowskian order.

*Proof.* We introduce a linearized “radiative” metric by performing a gauge transformation of the harmonic-coordinates metric defined by Eqs. (35)–(37), namely

$$H_{(1)}^{\alpha\beta} = h_{(1)}^{\alpha\beta} + \partial^{\alpha} \xi_{(1)}^{\beta} + \partial^{\beta} \xi_{(1)}^{\alpha} - \eta^{\alpha\beta} \partial_{\mu} \xi_{(1)}^{\mu}, \qquad (57)$$

where the gauge vector $\xi_{(1)}^{\alpha}$ is

$$\xi_{(1)}^{\alpha} = \frac{2M}{c^2} \eta^{0\alpha} \ln \left( \frac{r}{r_0} \right). \qquad (58)$$

This gauge transformation is non-harmonic:

$$\partial_{\mu} H_{(1)}^{\alpha\mu} = \Box \xi_{(1)}^{\alpha} = \frac{2M}{c^2 r^2} \eta^{0\alpha}. \qquad (59)$$

Its effect is to correct for the well-known logarithmic deviation of the retarded time in harmonic coordinates, with respect to the true space-time characteristic or light cones. After the change of gauge, the coordinate $u = t - r/c$ coincides with a null coordinate at the linearized level.<fnref n="21" /> This is the requirement to be satisfied by a linearized metric so that it can constitute the linearized approximation to a full (post-Minkowskian) radiative field [304]. One can easily show that, at the dominant order when $r \to +\infty$,

$$k_{\mu} k_{\nu} H_{(1)}^{\mu\nu} = \mathcal{O}\left(\frac{1}{r^2}\right), \qquad (60)$$

where $k^{\mu} = \eta^{\mu\nu} k_{\nu} = (1, \mathbf{n})$ is the outgoing Minkowskian null vector. Given any $n \ge 2$, let us recursively assume that we have obtained all the previous radiative post-Minkowskian coefficients $H_{(m)}^{\alpha\beta}$, i.e. $\forall m \le n - 1$, and that all of them satisfy

$$k_{\mu} k_{\nu} H_{(m)}^{\mu\nu} = \mathcal{O}\left(\frac{1}{r^2}\right). \qquad (61)$$

From this induction hypothesis one can prove that the $n$-th post-Minkowskian source term $\Lambda_{(n)}^{\alpha\beta} = \Lambda_{(n)}^{\alpha\beta}(H_{(1)}, \dots, H_{(n-1)})$ is such that

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

To the leading order this term takes the classic form of the stress-energy tensor of massless particles, with $\sigma_{(n)}$ being proportional to the power in the massless waves. One can show that all the problems

<footnote n="21">²¹ In this proof the coordinates are considered as dummy variables denoted $(t, r)$. At the end, when we obtain the radiative metric, we shall denote the associated radiative coordinates by $(T, R)$.</footnote>
with the appearance of logarithms come from the retarded integral of the terms in Eq. (62) that behave like $1/r^2$: See indeed the integration formula (83), which behaves like $\ln r/r$ at infinity. But now, thanks to the particular index structure of the term (62), we can correct for the effect by adjusting the gauge at the n-th post-Minkowskian order. We pose, as a gauge vector,

$$ \xi_{(n)}^{\alpha} = \mathcal{F}\mathcal{P} \Box_{\text{ret}}^{-1} \left[ \frac{k^{\alpha}}{2r^2} \int_{-\infty}^{u} dv \, \sigma_{(n)}(v, \mathbf{n}) \right], \quad (63) $$

where $\mathcal{F}\mathcal{P}$ refers to the same finite part operation as in Eq. (45). This vector is such that the logarithms that will appear in the corresponding gauge terms cancel out the logarithms coming from the retarded integral of the source term (62); see Ref. [41] for the details. Hence, to the n-th post-Minkowskian order, we define the radiative metric as

$$ H_{(n)}^{\alpha\beta} = U_{(n)}^{\alpha\beta} + V_{(n)}^{\alpha\beta} + \partial^{\alpha} \xi_{(n)}^{\beta} + \partial^{\beta} \xi_{(n)}^{\alpha} - \eta^{\alpha\beta} \partial_{\mu} \xi_{(n)}^{\mu}. \quad (64) $$

Here $U_{(n)}^{\alpha\beta}$ and $V_{(n)}^{\alpha\beta}$ denote the quantities that are the analogues of $u_{(n)}^{\alpha\beta}$ and $v_{(n)}^{\alpha\beta}$, which were introduced into the harmonic-coordinates algorithm: See Eqs. (45)-(48). In particular, these quantities are constructed in such a way that the sum $U_{(n)}^{\alpha\beta} + V_{(n)}^{\alpha\beta}$ is divergence-free, so we see that the radiative metric does not obey the harmonic-gauge condition, but instead

$$ \partial_{\mu} H_{(n)}^{\alpha\mu} = \Box \xi_{(n)}^{\alpha} = \frac{k^{\alpha}}{2r^2} \int_{-\infty}^{u} dv \, \sigma_{(n)}(v, \mathbf{n}). \quad (65) $$

The far-zone expansion of the latter metric is of the type (56), i.e., is free of any logarithms, and the retarded time in these coordinates tends asymptotically toward a null coordinate at future null infinity. The property of asymptotic simplicity, in the form given by Geroch & Horowitz [220], is proved by introducing the usual conformal factor $\Omega = 1/R$ in radiative coordinates [41]. Finally, it can be checked that the metric so constructed, which is a functional of the source multipole moments $I_L, \dots, Z_L$ (from the definition of the algorithm), is as general as the general harmonic-coordinate metric of Theorem 2, since it merely differs from it by a coordinate transformation $(t, \mathbf{x}) \to (T, \mathbf{X})$, where $(t, \mathbf{x})$ are the harmonic coordinates and $(T, \mathbf{X})$ the radiative ones, together with a re-definition of the multipole moments.

# 3 Asymptotic Gravitational Waveform

## 3.1 The radiative multipole moments

The leading-order term $1/R$ of the metric in radiative coordinates $(T, \mathbf{X})$ as given in Theorem 4, neglecting $O(1/R^2)$, yields the operational definition of two sets of STF radiative multipole moments, mass-type $U_L(U)$ and current-type $V_L(U)$. As we have seen, radiative coordinates are such that the retarded time $U \equiv T - R/c$ becomes asymptotically a null coordinate at future null infinity. The radiative moments are defined from the spatial components $ij$ of the metric in a transverse-traceless (TT) radiative coordinate system. By definition, we have [403]

$$ H_{ij}^{\text{TT}}(U, \mathbf{X}) = \frac{4G}{c^2 R} \mathcal{P}_{ijab}(\mathbf{N}) \sum_{\ell=2}^{+\infty} \frac{1}{c^\ell \ell!} \left\{ N_{L-2} U_{abL-2}(U) - \frac{2\ell}{c(\ell+1)} N_{cL-2} \epsilon_{cd(a} V_{b)dL-2}(U) \right\} + \mathcal{O}\left(\frac{1}{R^2}\right). \quad (66) $$

We have formally re-summed the whole post-Minkowskian series in Eq. (56) from $n=1$ up to $+\infty$. As before we denote for instance $N_{L-2} = N_{i_1} \cdots N_{i_{\ell-2}}$ and so on, where $N_i = (\mathbf{N})_i$ and
$N = X/R$. The TT algebraic projection operator $P_{ijab}$ has already been defined at the occasion of the quadrupole-moment formalism in Eq. (2); and obviously the multipole decomposition (66) represents the generalization of the quadrupole formalism. Notice that the meaning of Eq. (66) is for the moment rather empty, because we do not yet know how to relate the radiative moments to the actual source parameters. Only at the Newtonian level do we know this relation, which is

$$U_{ij}(U) = Q_{ij}^{(2)}(U) + \mathcal{O}\left(\frac{1}{c^2}\right), \qquad (67)$$

where $Q_{ij}$ is the Newtonian quadrupole moment (3). Associated to the asymptotic waveform (66) we can compute by standard methods the total energy flux $\mathcal{F} = (dE/dU)^{\text{GW}}$ and angular momentum flux $\mathcal{G}_i = (dJ_i/dU)^{\text{GW}}$ in gravitational waves [403]:

$$\mathcal{F} = \sum_{\ell=2}^{+\infty} \frac{G}{c^{2\ell+1}} \left\{ \frac{(\ell+1)(\ell+2)}{(\ell-1)\ell\ell!(2\ell+1)!!} U_L^{(1)} U_L^{(1)} + \frac{4\ell(\ell+2)}{c^2(\ell-1)(\ell+1)!(2\ell+1)!!} V_L^{(1)} V_L^{(1)} \right\}. \qquad (68a)$$

$$\mathcal{G}_i = \epsilon_{iab} \sum_{\ell=2}^{+\infty} \frac{G}{c^{2\ell+1}} \left\{ \frac{(\ell+1)(\ell+2)}{(\ell-1)\ell!(2\ell+1)!!} U_{aL-1} U_{bL-1}^{(1)} + \frac{4\ell^2(\ell+2)}{c^2(\ell-1)(\ell+1)!(2\ell+1)!!} V_{aL-1} V_{bL-1}^{(1)} \right\}. \qquad (68b)$$

Next we introduce two unit polarization vectors $P$ and $Q$, orthogonal and transverse to the direction of propagation $N$ (hence $N_i N_j + P_i P_j + Q_i Q_j = \delta_{ij}$). Our convention for the choice of $P$ and $Q$ will be clarified in Section 9.4. Then the two “plus” and “cross” polarization states of the asymptotic waveform are defined by

$$h_+ = \frac{1}{2} (P_i P_j - Q_i Q_j) H_{ij}^{\text{TT}}, \qquad (69a)$$

$$h_{\times} = \frac{1}{2} (P_i Q_j + P_j Q_i) H_{ij}^{\text{TT}}. \qquad (69b)$$

Although the multipole decomposition (66) is completely general, it will also be important, having in view the comparison between the post-Newtonian and numerical results (see for instance Refs. [107, 34, 237, 97, 98]), to consider separately the various modes $(\ell, m)$ of the asymptotic waveform as defined with respect to a basis of spin-weighted spherical harmonics of weight $-2$. Those harmonics are function of the spherical angles $(\theta, \phi)$ defining the direction of propagation $N$, and given by

$$Y_{(-2)}^{\ell m} = \sqrt{\frac{2\ell + 1}{4\pi}} d^{\ell m}(\theta) e^{i m \phi}, \qquad (70a)$$

$$d^{\ell m} = \sum_{k=k_1}^{k_2} \frac{(-)^k}{k!} e_k^{\ell m} \left(\cos\frac{\theta}{2}\right)^{2\ell+m-2k-2} \left(\sin\frac{\theta}{2}\right)^{2k-m+2}, \qquad (70b)$$

$$e_k^{\ell m} = \frac{\sqrt{(\ell + m)!(\ell - m)!(\ell + 2)!(\ell - 2)!}}{(k - m + 2)!(\ell + m - k)!(\ell - k - 2)!}, \qquad (70c)$$

where $k_1 = \max(0, m - 2)$ and $k_2 = \min(\ell + m, \ell - 2)$. We thus decompose $h_+$ and $h_\times$ onto the basis of such spin-weighted spherical harmonics, which means (see e.g., [107, 272])

$$h_+ - i h_\times = \sum_{\ell=2}^{+\infty} \sum_{m=-\ell}^{\ell} h^{\ell m} Y_{(-2)}^{\ell m}(\theta, \phi). \qquad (71)$$
Using the orthonormality properties of these harmonics we can invert the latter decomposition and obtain the separate modes $h^{\ell m}$ from a surface integral,

$$h^{\ell m} = \int d\Omega [h_+ - i h_\times] \bar{Y}^{\ell m}_{(-2)}(\theta, \phi), \qquad (72)$$

where the overline refers to the complex conjugation. On the other hand, we can also relate $h^{\ell m}$ to the radiative multipole moments $U_L$ and $V_L$. The result is

$$h^{\ell m} = - \frac{G}{\sqrt{2} R c^{\ell+2}} \left[ U^{\ell m} - \frac{i}{c} V^{\ell m} \right], \qquad (73)$$

where $U^{\ell m}$ and $V^{\ell m}$ denote the radiative mass and current moments in standard (non-STF) guise. These are related to the STF moments by

$$U^{\ell m} = \frac{4}{\ell!} \sqrt{\frac{(\ell+1)(\ell+2)}{2\ell(\ell-1)}} \alpha_L^{\ell m} U_L, \qquad (74a)$$

$$V^{\ell m} = -\frac{8}{\ell!} \sqrt{\frac{\ell(\ell+2)}{2(\ell+1)(\ell-1)}} \alpha_L^{\ell m} V_L. \qquad (74b)$$

Here $\alpha_L^{\ell m}$ denotes the STF tensor connecting together the usual basis of spherical harmonics $Y^{\ell m}$ to the set of STF tensors $\hat{N}_L = N_{(i_1}...N_{i_\ell)}$ (where the brackets indicate the STF projection). Indeed both $Y^{\ell m}$ and $\hat{N}_L$ are basis of an irreducible representation of weight $\ell$ of the rotation group; the two basis are related by$^{22}$

$$\hat{N}_L(\theta, \phi) = \sum_{m=-\ell}^{\ell} \alpha_L^{\ell m} Y^{\ell m}(\theta, \phi), \qquad (75a)$$

$$Y^{\ell m}(\theta, \phi) = \frac{(2\ell + 1)!!}{4\pi l!} \bar{\alpha}_L^{\ell m} \hat{N}_L(\theta, \phi). \qquad (75b)$$

In Section 9.5 we shall present all the modes $(\ell, m)$ of gravitational waves from inspiralling compact binaries up to 3PN order, and even 3.5PN order for the dominant mode (2, 2).

## 3.2 Gravitational-wave tails and tails-of-tails

We learned from Theorem 4 the general method which permits the computation of the radiative multipole moments $U_L$, $V_L$ in terms of the source moments $I_L$, $J_L$, ..., $Z_L$, or in terms of the intermediate canonical moments $M_L$, $S_L$ discussed in Section 2.4. We shall now show that the relation between $U_L$, $V_L$ and $M_L$, $S_L$ (say) includes tail effects starting at the relative 1.5PN order.

Tails are due to the back-scattering of multipolar waves off the Schwarzschild curvature generated by the total mass monopole M of the source. They correspond to the non-linear interaction between M and the multipole moments $M_L$ and $S_L$, and are given by some non-local integrals, extending over the past history of the source. At the 1.5PN order we find [59, 44]

$$U_L(U) = M_L^{(\ell)}(U) + \frac{2GM}{c^3} \int_0^{+\infty} d\tau M_L^{(\ell+2)}(U-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \kappa_\ell \right] + \mathcal{O}\left(\frac{1}{c^5}\right), \qquad (76a)$$

<footnote>$^{22}$ The STF tensorial coefficient $\alpha_L^{\ell m}$ can be computed as $\alpha_L^{\ell m} = \int d\Omega \hat{N}_L \bar{Y}^{\ell m}$. Our notation is related to that used in Refs. [403, 272] by $Y_L^{\ell m} = \frac{(2\ell+1)!!}{4\pi l!} \bar{\alpha}_L^{\ell m}$.</footnote>
$$V_L(U) = S_L^{(\ell)}(U) + \frac{2GM}{c^3} \int_0^{+\infty} d\tau S_L^{(\ell+2)}(U-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \pi_\ell \right] + \mathcal{O}\left(\frac{1}{c^5}\right), \quad (76b)$$

where $r_0$ is the length scale introduced in Eq. (42), and the constants $\kappa_\ell$ and $\pi_\ell$ are given by

$$\kappa_{\ell} = \frac{2\ell^2 + 5\ell + 4}{\ell(\ell + 1)(\ell + 2)} + \sum_{k=1}^{\ell-2} \frac{1}{k}, \quad (77a)$$

$$\pi_{\ell} = \frac{\ell - 1}{\ell(\ell + 1)} + \sum_{k=1}^{\ell-1} \frac{1}{k}. \quad (77b)$$

Recall from the gauge vector $\xi_{(1)}^\alpha$ found in Eq. (58) that the retarded time $U = T - R/c$ in radiative coordinates is related to the retarded time $u = t - r/c$ in harmonic coordinates by

$$U = u - \frac{2GM}{c^3} \ln \left( \frac{r}{r_0} \right) + \mathcal{O}(G^2). \quad (78)$$

Inserting $U$ as given by Eq. (78) into Eqs. (76) we obtain the radiative moments expressed in terms of “source-rooted” harmonic coordinates $(t, r)$, e.g.,

$$U_L(U) = M_L^{(\ell)}(u) + \frac{2GM}{c^3} \int_0^{+\infty} d\tau M_L^{(\ell+2)}(u-\tau) \left[ \ln\left(\frac{c\tau}{2r}\right) + \kappa_\ell \right] + \mathcal{O}\left(\frac{1}{c^5}\right). \quad (79)$$

The remainder $\mathcal{O}(G^2)$ in Eq. (78) is negligible here. This expression no longer depends on the constant $r_0$, i.e., we find that $r_0$ gets replaced by $r$. If we now replace the harmonic coordinates $(t, r)$ to some new ones, such as, for instance, some “Schwarzschild-like” coordinates $(t', r')$ such that $t' = t$ and $r' = r + GM/c^2$ (and $u' = u - GM/c^3$), we get

$$U_L(U) = M_L^{(\ell)}(u') + \frac{2GM}{c^3} \int_0^{+\infty} d\tau M_L^{(\ell+2)}(u' - \tau) \left[ \ln\left(\frac{c\tau}{2r'}\right) + \kappa'_\ell \right] + \mathcal{O}\left(\frac{1}{c^5}\right), \quad (80)$$

where $\kappa'_\ell = \kappa_\ell + 1/2$. This shows that the constant $\kappa_\ell$ (and $\pi_\ell$ as well) depends on the choice of source-rooted coordinates $(t, r)$: For instance, we have $\kappa_2 = 11/12$ in harmonic coordinates from Eq. (77a), but $\kappa'_2 = 17/12$ in Schwarzschild coordinates [345].

The tail integrals in Eqs. (76) involve all the instants from $-\infty$ in the past up to the current retarded time $U$. However, strictly speaking, they do not extend up to infinite past, since we have assumed in Eq. (29) that the metric is stationary before the date $-\mathcal{T}$. The range of integration of the tails is therefore limited a priori to the time interval $[-\mathcal{T}, U]$. But now, once we have derived the tail integrals, thanks to the latter technical assumption of stationarity in the past, we can argue that the results are in fact valid in more general situations for which the field has never been stationary. We have in mind the case of two bodies moving initially on some unbound (hyperbolic-like) orbit, and which capture each other, because of the loss of energy by gravitational radiation, to form a gravitationally bound system around time $-\mathcal{T}$.

In this situation let us check, using a simple Newtonian model for the behaviour of the multipole moment $M_L(U-\tau)$ when $\tau \to +\infty$, that the tail integrals, when assumed to extend over the whole time interval $[-\infty, U]$, remain perfectly well-defined (i.e., convergent) at the integration bound $\tau = +\infty$. Indeed it can be shown [180] that the motion of initially free particles interacting gravitationally is given by $x^i(U-\tau) = V^i\tau + W^i \ln \tau + X^i + o(1)$, where $V^i, W^i$ and $X^i$ denote constant vectors, and $o(1) \to 0$ when $\tau \to +\infty$. From that physical assumption we find that the multipole moments behave when $\tau \to +\infty$ like

$$M_L(U - \tau) = A_L \tau^\ell + B_L \tau^{\ell-1} \ln \tau + C_L \tau^{\ell-1} + o(\tau^{\ell-1}), \quad (81)$$
where $A_L$, $B_L$ and $C_L$ are constant tensors. We used the fact that the moment $M_L$ will agree at the Newtonian level with the standard expression for the $\ell$-th mass multipole moment $Q_L$. The appropriate time derivatives of the moment appearing in Eq. (76a) are therefore dominantly like

$$M_L^{(\ell+2)}(U-\tau) = \frac{D_L}{\tau^3} + o(\tau^{-3}), \qquad (82)$$

which ensures that the tail integral is convergent. This fact can be regarded as an a posteriori justification of our a priori too restrictive assumption of stationarity in the past. Thus, this assumption does not seem to yield any physical restriction on the applicability of the final formulas. However, once again, we emphasize that the past-stationarity is appropriate for real astrophysical sources of gravitational waves which have been formed at a finite instant in the past.

To obtain the results (76), we must implement in details the post-Minkowskian algorithm presented in Section 2.3. Let us flash here some results obtained with such algorithm. Consider first the case of the interaction between the constant mass monopole moment M (or ADM mass) and the time-varying quadrupole moment $M_{ij}$. This coupling will represent the dominant non-static multipole interaction in the waveform. For these moments we can write the linearized metric using Eq. (35) in which by definition of the “canonical” construction we insert the canonical moments $M_{ij}$ in place of $I_{ij}$ (notice that $M=I$). We must plug this linearized metric into the quadratic-order part $N^{\alpha\beta}(h, h)$ of the gravitational source term (24)–(25) and explicitly given by Eq. (26). This yields many terms; to integrate these following the algorithm [cf. Eq. (45)], we need some explicit formulas for the retarded integral of an extended (non-compact-support) source having some definite multipolarity $\ell$. A thorough account of the technical formulas necessary for handling the quadratic and cubic interactions is given in the Appendices of Refs. [50] and [48]. For the present computation the most crucial formula, needed to control the tails, corresponds to a source term behaving like $1/r^2$:

$$\Box_{\text{ret}}^{-1} \left[ \frac{\hat{n}_L}{r^2} F(t-r) \right] = -\hat{n}_L \int_1^{+\infty} dx \, Q_\ell(x) F(t-rx), \qquad (83)$$

where F is any smooth function representing a time derivative of the quadrupole moment, and $Q_\ell$ denotes the Legendre function of the second kind.23 Note that there is no need to include a finite part operation $\mathcal{F}P$ in Eq. (83) as the integral is convergent. With the help of this and other formulas we obtain successively the objects defined in this algorithm by Eqs. (45)–(48) and finally obtain the quadratic metric (49) for that multipole interaction. The result is [60]$^{24}$

$$h_{(2)}^{00} = \frac{Mn_{ab}}{r^4} \left[ -21M_{ab} - 21rM_{ab}^{(1)} + 7r^2M_{ab}^{(2)} + 10r^3M_{ab}^{(3)} \right] \\ + 8Mn_{ab} \int_1^{+\infty} dx \, Q_2(x)M_{ab}^{(4)}(t-rx), \qquad (84a)$$

$$h_{(2)}^{0i} = \frac{Mn_{iab}}{r^3} \left[ -M_{ab}^{(1)} - rM_{ab}^{(2)} - \frac{1}{3}r^2M_{ab}^{(3)} \right] \\ + \frac{Mn_a}{r^3} \left[ -5M_{ai}^{(1)} - 5rM_{ai}^{(2)} + \frac{19}{3}r^2M_{ai}^{(3)} \right]$$

$$Q_{\ell}(x) = \frac{1}{2} \int_{-1}^{1} \frac{dz}{x-z} P_{\ell}(z) = \frac{1}{2} P_{\ell}(x) \ln \left( \frac{x+1}{x-1} \right) - \sum_{j=1}^{\ell} \frac{1}{j} P_{\ell-j}(x) P_{j-1}(x).$$

In the complex plane there is a branch cut from $-\infty$ to 1. The first equality is known as the Neumann formula for the Legendre function.

<footnote n="23">23 The function $Q_\ell$ is given in terms of the Legendre polynomial $P_\ell$ by</footnote>

<footnote n="24">24 We pose $c=1$ until the end of this section.</footnote>
$$ + 8Mn_a \int_{1}^{+\infty} dx Q_1(x) M_{ai}^{(4)}(t - rx), \qquad (84b) $$

$$ \begin{aligned} h_{(2)}^{ij} = & \frac{Mn_{ijab}}{r^4} \left[ -\frac{15}{2}M_{ab} - \frac{15}{2}rM_{ab}^{(1)} - 3r^2M_{ab}^{(2)} - \frac{1}{2}r^3M_{ab}^{(3)} \right] \\ & + \frac{M\delta_{ij}n_{ab}}{r^4} \left[ -\frac{1}{2}M_{ab} - \frac{1}{2}rM_{ab}^{(1)} - 2r^2M_{ab}^{(2)} - \frac{11}{6}r^3M_{ab}^{(3)} \right] \\ & + \frac{Mn_{a(i}}}{r^4} \left[ 6M_{j)a} + 6rM_{j)a}^{(1)} + 6r^2M_{j)a}^{(2)} + 4r^3M_{j)a}^{(3)} \right] \\ & + \frac{M}{r^4} \left[ -M_{ij} - rM_{ij}^{(1)} - 4r^2M_{ij}^{(2)} - \frac{11}{3}r^3M_{ij}^{(3)} \right] \\ & + 8M \int_{1}^{+\infty} dx Q_0(x) M_{ij}^{(4)}(t - rx). \end{aligned} \qquad (84c) $$

The metric is composed of two types of terms: “instantaneous” ones depending on the values of the quadrupole moment at the retarded time $u = t - r$, and “hereditary” tail integrals, depending on all previous instants $t - rx < u$.

Let us investigate now the cubic interaction between two mass monopoles M with the mass quadrupole $M_{ij}$. Obviously, the source term corresponding to this interaction will involve [see Eq. (40b)] cubic products of three linear metrics, say $h_M \times h_M \times h_{M_{ij}}$, and quadratic products between one linear metric and one quadratic, say $h_{M^2} \times h_{M_{ij}}$ and $h_M \times h_{MM_{ij}}$. The latter case is the most tricky because the tails present in $h_{MM_{ij}}$, which are given explicitly by Eqs. (84), will produce in turn some tails of tails in the cubic metric $h_{M^2 M_{ij}}$. The computation is rather involved [48] but can now be performed by an algebraic computer programme [74, 197]. Let us just mention the most difficult of the needed integration formulas for this calculation:<fnref n="25" />

$$ \mathcal{F}P_{\text{ret}}^{-1} \left[ \frac{\hat{n}_L}{r} \int_1^{+\infty} dx Q_m(x) F(t - rx) \right] = \hat{n}_L \int_1^{+\infty} dy F^{(-1)}(t - ry) \\ \times \left\{ Q_l(y) \int_1^y dx Q_m(x) \frac{dP_l}{dx}(x) + P_l(y) \int_y^{+\infty} dx Q_m(x) \frac{dQ_l}{dx}(x) \right\}, \qquad (85) $$

where $F^{(-1)}$ is the time anti-derivative of F. With this formula and others given in Ref. [48] we are able to obtain the closed algebraic form of the cubic metric for the multipole interaction $M \times M \times M_{ij}$, at the leading order when the distance to the source $r \to \infty$ with $u = \text{const}$. The result is<fnref n="26" />

$$ h_{(3)}^{00} = \frac{M^2 n_{ab}}{r} \int_{0}^{+\infty} d\tau M_{ab}^{(5)} \left[ -4 \ln^2 \left( \frac{\tau}{2r} \right) - 4 \ln \left( \frac{\tau}{2r} \right) + \frac{116}{21} \ln \left( \frac{\tau}{2r_0} \right) - \frac{7136}{2205} \right] \\ + o\left( \frac{1}{r} \right), \qquad (86a) $$

$$ h_{(3)}^{0i} = \frac{M^2 \hat{n}_{iab}}{r} \int_{0}^{+\infty} d\tau M_{ab}^{(5)} \left[ -\frac{2}{3} \ln \left( \frac{\tau}{2r} \right) - \frac{4}{105} \ln \left( \frac{\tau}{2r_0} \right) - \frac{716}{1225} \right] $$

$$ \frac{1}{2} \int_{-1}^{1} \frac{dz P_{\ell}(z)}{\sqrt{(xy-z)^2 - (x^2-1)(y^2-1)}} = Q_{\ell}(x)P_{\ell}(y), $$

<footnote n="25">²⁵ The equation (85) has been obtained using a not so well known mathematical relation between the Legendre functions and polynomials:</footnote>

<footnote>where $1 \le y < x$ is assumed. See Appendix A in Ref. [48] for the proof. This relation constitutes a generalization of the Neumann formula (see the footnote 23).</footnote>

<footnote n="26">²⁶ The neglected remainders are indicated by $o(1/r)$ rather than $O(1/r^2)$ because they contain powers of the logarithm of $r$; in fact they could be more accurately written as $o(r^{\epsilon-2})$ for some $\epsilon \ll 1$.</footnote>
$$ + \frac{M^2 n_a}{r} \int_0^{+\infty} d\tau M_{ai}^{(5)} \left[ -4 \ln^2 \left( \frac{\tau}{2r} \right) - \frac{18}{5} \ln \left( \frac{\tau}{2r} \right) + \frac{416}{75} \ln \left( \frac{\tau}{2r_0} \right) - \frac{22724}{7875} \right] \\ + o \left( \frac{1}{r} \right), \qquad (86b) $$

$$ h_{(3)}^{ij} = \frac{M^2 \hat{n}_{ijab}}{r} \int_0^{+\infty} d\tau M_{ab}^{(5)} \left[ -\ln\left(\frac{\tau}{2r}\right) - \frac{191}{210} \right] \\ + \frac{M^2 \delta_{ij} n_{ab}}{r} \int_0^{+\infty} d\tau M_{ab}^{(5)} \left[ -\frac{80}{21} \ln\left(\frac{\tau}{2r}\right) - \frac{32}{21} \ln\left(\frac{\tau}{2r_0}\right) - \frac{296}{35} \right] \\ + \frac{M^2 \hat{n}_{a(i}}}{r} \int_0^{+\infty} d\tau M_{j)a}^{(5)} \left[ \frac{52}{7} \ln\left(\frac{\tau}{2r}\right) + \frac{104}{35} \ln\left(\frac{\tau}{2r_0}\right) + \frac{8812}{525} \right] \\ + \frac{M^2}{r} \int_0^{+\infty} d\tau M_{ij}^{(5)} \left[ -4 \ln^2\left(\frac{\tau}{2r}\right) - \frac{24}{5} \ln\left(\frac{\tau}{2r}\right) + \frac{76}{15} \ln\left(\frac{\tau}{2r_0}\right) - \frac{198}{35} \right] \\ + o\left(\frac{1}{r}\right), \qquad (86c) $$

where all the moments $M_{ab}$ are evaluated at the instant $u-\tau = t-r-\tau$. Notice that the logarithms in Eqs. (86) contain either the ratio $\tau/r$ or $\tau/r_0$. We shall discuss in Eqs. (93)–(94) below the interesting fate of the arbitrary constant $r_0$.

From Theorem 4, the presence of logarithms of $r$ in Eqs. (86) is an artifact of the harmonic coordinates $x^\alpha$, and it is convenient to gauge them away by introducing radiative coordinates $X^\alpha$ at future null infinity. For controlling the leading $1/R$ term at infinity, it is sufficient to take into account the linearized logarithmic deviation of the light cones in harmonic coordinates: $X^\alpha = x^\alpha + G\xi_{(1)}^\alpha + \mathcal{O}(G^2)$, where $\xi_{(1)}^\alpha$ is the gauge vector defined by Eq. (58) [see also Eq. (78)]. With this coordinate change one removes the logarithms of $r$ in Eqs. (86) and we obtain the radiative (or Bondi-type [93]) logarithmic-free expansion

$$ H_{(3)}^{00} = \frac{M^2 N_{ab}}{R} \int_0^{+\infty} d\tau M_{ab}^{(5)} \left[ -4 \ln^2 \left( \frac{\tau}{2r_0} \right) + 3221 \ln \left( \frac{\tau}{2r_0} \right) - \frac{7136}{2205} \right] \\ + \mathcal{O} \left( \frac{1}{R^2} \right), \qquad (87a) $$

$$ H_{(3)}^{0i} = \frac{M^2 \hat{N}_{iab}}{R} \int_0^{+\infty} d\tau M_{ab}^{(5)} \left[ -\frac{74}{105} \ln \left( \frac{\tau}{2r_0} \right) - \frac{716}{1225} \right] \\ + \frac{M^2 N_a}{R} \int_0^{+\infty} d\tau M_{ai}^{(5)} \left[ -4 \ln^2 \left( \frac{\tau}{2r_0} \right) + \frac{146}{75} \ln \left( \frac{\tau}{2r_0} \right) - \frac{22724}{7875} \right] \\ + \mathcal{O} \left( \frac{1}{R^2} \right), \qquad (87b) $$

$$ H_{(3)}^{ij} = \frac{M^2 \hat{N}_{ijab}}{R} \int_0^{+\infty} d\tau M_{ab}^{(5)} \left[ -\ln\left(\frac{\tau}{2r_0}\right) - \frac{191}{210} \right] \\ + \frac{M^2 \delta_{ij} N_{ab}}{R} \int_0^{+\infty} d\tau M_{ab}^{(5)} \left[ -\frac{16}{3} \ln\left(\frac{\tau}{2r_0}\right) - \frac{296}{35} \right] \\ + \frac{M^2 \hat{N}_{a(i}}}{R} \int_0^{+\infty} d\tau M_{j)a}^{(5)} \left[ \frac{52}{5} \ln\left(\frac{\tau}{2r_0}\right) + \frac{8812}{525} \right] \\ + \frac{M^2}{R} \int_0^{+\infty} d\tau M_{ij}^{(5)} \left[ -4 \ln^2\left(\frac{\tau}{2r_0}\right) + \frac{4}{15} \ln\left(\frac{\tau}{2r_0}\right) - \frac{198}{35} \right] $$
$$ + \mathcal{O}\left(\frac{1}{R^2}\right), \qquad (87c) $$

where the moments are evaluated at time $U-\tau = T-R-\tau$. It is trivial to compute the contribution of the radiative moments corresponding to that metric. We find the “tail of tail” term which will be reported in Eq. (91) below.

## 3.3 Radiative versus source moments

We first give the result for the radiative quadrupole moment $U_{ij}$ expressed as a functional of the intermediate canonical moments $M_L, S_L$ up to 3.5PN order included. The long calculation follows from implementing the explicit MPM algorithm of Section 2.3 and yields various types of terms:

$$ U_{ij} = U_{ij}^{\text{inst}} + U_{ij}^{\text{tail}} + U_{ij}^{\text{tail-tail}} + U_{ij}^{\text{mem}} + \mathcal{O}\left(\frac{1}{c^8}\right). \qquad (88) $$

1. The instantaneous (i.e., non-hereditary) piece $U_{ij}^{\text{inst}}$ up to 3.5PN order reads

$$ \begin{aligned} U_{ij}^{\text{inst}} = & M_{ij}^{(2)} \\ & + \frac{G}{c^5} \left[ \frac{1}{7} M_{a\langle i}^{(5)} M_{j\rangle a} - \frac{5}{7} M_{a\langle i}^{(4)} M_{j\rangle a} - \frac{2}{7} M_{a\langle i}^{(3)} M_{j\rangle a}^{(2)} + \frac{1}{3} \epsilon_{ab\langle i} M_{j\rangle a}^{(4)} S_b \right] \\ & + \frac{G}{c^7} \left[ -\frac{64}{63} S_{a\langle i}^{(2)} S_{j\rangle a}^{(3)} + \frac{1957}{3024} M_{ijab}^{(3)} M_{ab}^{(4)} + \frac{5}{2268} M_{ab\langle i}^{(3)} M_{j\rangle ab}^{(4)} + \frac{19}{648} M_{ab}^{(3)} M_{ijab}^{(4)} \right. \\ & \left. \qquad + \frac{16}{63} S_{a\langle i}^{(1)} S_{j\rangle a}^{(4)} + \frac{1685}{1008} M_{ijab}^{(2)} M_{ab}^{(5)} + \frac{5}{126} M_{ab\langle i}^{(2)} M_{j\rangle ab}^{(5)} - \frac{5}{756} M_{ab}^{(2)} M_{ijab}^{(5)} \right. \\ & \left. \qquad + \frac{80}{63} S_{a\langle i} S_{j\rangle a}^{(5)} + \frac{5}{42} S_a S_{ija}^{(5)} + \frac{41}{28} M_{ijab}^{(1)} M_{ab}^{(6)} + \frac{5}{189} M_{ab\langle i}^{(1)} M_{j\rangle ab}^{(6)} \right. \\ & \left. \qquad + \frac{1}{432} M_{ab}^{(1)} M_{ijab}^{(6)} + \frac{91}{216} M_{ijab} M_{ab}^{(7)} - \frac{5}{252} M_{ab\langle i} M_{j\rangle ab}^{(7)} - \frac{1}{432} M_{ab} M_{ijab}^{(7)} \right. \\ & \left. \qquad + \epsilon_{ac\langle i} \left( \frac{32}{189} M_{j\rangle bc}^{(3)} S_{ab}^{(3)} - \frac{1}{6} M_{ab}^{(3)} S_{j\rangle bc}^{(3)} + \frac{3}{56} S_{j\rangle bc}^{(2)} M_{ab}^{(4)} + \frac{10}{189} S_{ab}^{(2)} M_{j\rangle bc}^{(4)} \right. \right. \\ & \left. \qquad + \frac{65}{189} M_{j\rangle bc}^{(2)} S_{ab}^{(4)} + \frac{1}{28} M_{ab}^{(2)} S_{j\rangle bc}^{(4)} + \frac{187}{168} S_{j\rangle bc}^{(1)} M_{ab}^{(5)} - \frac{1}{189} S_{ab}^{(1)} M_{j\rangle bc}^{(5)} \right. \\ & \left. \qquad - \frac{5}{189} M_{j\rangle bc}^{(1)} S_{ab}^{(5)} + \frac{1}{24} M_{ab}^{(1)} S_{j\rangle bc}^{(5)} + \frac{65}{84} S_{j\rangle bc} M_{ab}^{(6)} + \frac{1}{189} S_{ab} M_{j\rangle bc}^{(6)} \right. \\ & \left. \qquad - \frac{10}{63} M_{j\rangle bc} S_{ab}^{(6)} + \frac{1}{168} M_{ab} S_{j\rangle bc}^{(6)} \right] \end{aligned} \qquad (89) $$

The Newtonian term in this expression contains the Newtonian quadrupole moment $Q_{ij}$ and recovers the standard quadrupole formalism [see Eq. (67)];

2. The hereditary tail integral $U_{ij}^{\text{tail}}$ is made of the dominant tail term at 1.5PN order in agreement with Eq. (76a) above:

$$ U_{ij}^{\text{tail}} = \frac{2GM}{c^3} \int_{0}^{+\infty} d\tau \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{11}{12} \right] M_{ij}^{(4)}(U-\tau). \qquad (90) $$

The length scale $r_0$ is the one that enters our definition of the finite-part operation $\mathcal{FP}$ [see Eq. (42)] and it enters also the relation between the radiative and harmonic retarded times given by Eq. (78);
3. The hereditary tail-of-tail term appears dominantly at 3PN order [48] and is issued from the radiative metric computed in Eqs. (87):

$$ U_{ij}^{\text{tail-tail}} = 2 \left( \frac{GM}{c^3} \right)^2 \int_0^{+\infty} d\tau \left[ \ln^2 \left( \frac{c\tau}{2r_0} \right) + \frac{57}{70} \ln \left( \frac{c\tau}{2r_0} \right) + \frac{124627}{44100} \right] M_{ij}^{(5)}(U-\tau); \quad (91) $$

4. Finally the memory-type hereditary piece $U_{ij}^{\text{mem}}$ contributes at orders 2.5PN and 3.5PN and is given by

$$ U_{ij}^{\text{mem}} = \frac{G}{c^5} \left[ -\frac{2}{7} \int_0^{+\infty} d\tau M_{a\langle i}^{(3)} M_{j\rangle a}^{(3)}(U-\tau) \right] + \frac{G}{c^7} \left[ -\frac{32}{63} \int_0^{+\infty} d\tau S_{a\langle i}^{(3)} S_{j\rangle a}^{(3)}(U-\tau) - \frac{5}{756} \int_0^{+\infty} d\tau M_{ab}^{(4)} M_{ijab}^{(4)}(U-\tau) - \frac{20}{189} \epsilon_{ab\langle i} \int_0^{+\infty} d\tau S_{ac}^{(3)} M_{j\rangle bc}^{(4)}(U-\tau) + \frac{5}{42} \epsilon_{ab\langle i} \int_0^{+\infty} d\tau M_{ac}^{(3)} S_{j\rangle bc}^{(4)}(U-\tau) \right]. \quad (92) $$

The 2.5PN non-linear memory integral – the first term inside the coefficient of $G/c^5$ – has been obtained using both post-Newtonian methods [42, 427, 406, 60, 50] and rigorous studies of the field at future null infinity [128]. The expression (92) is in agreement with the more recent computation of the non-linear memory up to any post-Newtonian order in Refs. [189, 192].

Be careful to note that the latter post-Newtonian orders correspond to “relative” orders when counted in the local radiation-reaction force, present in the equations of motion: For instance, the 1.5PN tail integral in Eq. (90) is due to a 4PN radiative effect in the equations of motion [58]; similarly, the 3PN tail-of-tail integral is expected to be associated with some radiation-reaction terms occurring at the 5.5PN order.

Note that $U_{ij}$, when expressed in terms of the intermediate moments $M_L$ and $S_L$, shows a dependence on the (arbitrary) length scale $r_0$; cf. the tail and tail-of-tail contributions (90)–(91). Most of this dependence comes from our definition of a radiative coordinate system as given by (78). Exactly as we have done for the 1.5PN tail term in Eq. (79), we can remove most of the $r_0$'s by inserting $U = u - \frac{2GM}{c^3} \ln(r/r_0)$ back into (89)–(92), and expanding the result when $c \to \infty$, keeping the necessary terms consistently. In doing so one finds that there remains a $r_0$-dependent term at the 3PN order, namely

$$ U_{ij} = M_{ij}^{(2)}(u) - \frac{214}{105} \ln\left(\frac{r}{r_0}\right) \left(\frac{GM}{c^3}\right)^2 M_{ij}^{(4)}(u) + \text{terms independent of } r_0. \quad (93) $$

However, the latter dependence on $r_0$ is fictitious and should *in fine* disappear. The reason is that when we compute explicitly the mass quadrupole moment $M_{ij}$ for a given matter source, we will find an extra contribution depending on $r_0$ occurring at the 3PN order which will cancel out the one in Eq. (93). Indeed we shall compute the source quadrupole moment $I_{ij}$ of compact binaries at the 3PN order, and we do observe on the result (300)–(301) below the requested terms depending on $r_0$, namely$^{27}$

$$ M_{ij} = Q_{ij} + \frac{214}{105} \ln \left( \frac{r_{12}}{r_0} \right) \left( \frac{Gm}{c^3} \right)^2 Q_{ij}^{(2)} + \text{terms independent of } r_0. \quad (94) $$

where $Q_{ij} = \mu \hat{x}_{ij}$ denotes the Newtonian quadrupole, $r_{12}$ is the separation between the particles, and $m$ is the total mass differing from the ADM mass $M$ by small post-Newtonian corrections.

<footnote>$^{27}$ The canonical moment $M_{ij}$ differs from the source moment $I_{ij}$ by small 2.5PN and 3.5PN terms; see Eq. (97).</footnote>
Combining Eqs. (93) and (94) we see that the $r_0$-dependent terms cancel as expected. The appearance of a logarithm and its associated constant $r_0$ at the 3PN order was pointed out in Ref. [7]; it was rederived within the present formalism in Refs. [58, 48]. Recently a result equivalent to Eq. (93) was obtained by means of the EFT approach using considerations related to the renormalization group equation [222].

The previous formulas for the 3.5PN radiative quadrupole moment permit to compute the dominant mode (2, 2) of the waveform up to order 3.5PN [197]; however, to control the full waveform one has also to take into account the contributions of higher-order radiative moments. Here we list the most accurate results we have for all the moments that permit the derivation of the waveform up to order 3PN [74]:<fnref n="28" />

$$ 
\begin{aligned}
U_{ijk}(U) = M_{ijk}^{(3)}(U) &+ \frac{2GM}{c^3} \int_0^{+\infty} d\tau \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{97}{60} \right] M_{ijk}^{(5)}(U-\tau) \\
&+ \frac{G}{c^5} \left\{ \int_0^{+\infty} d\tau \left[ -\frac{1}{3} M_{a\langle i}^{(3)} M_{jk\rangle a}^{(4)} - \frac{4}{5} \epsilon_{ab\langle i} M_{ja}^{(3)} S_{k\rangle b}^{(3)} \right] (U-\tau) \right. \\
&- \frac{4}{3} M_{a\langle i}^{(3)} M_{jk\rangle a}^{(3)} - \frac{9}{4} M_{a\langle i}^{(4)} M_{jk\rangle a}^{(2)} + \frac{1}{4} M_{a\langle i}^{(2)} M_{jk\rangle a}^{(4)} - \frac{3}{4} M_{a\langle i}^{(5)} M_{jk\rangle a}^{(1)} + \frac{1}{4} M_{a\langle i}^{(1)} M_{jk\rangle a}^{(5)} \\
&+ \frac{1}{12} M_{a\langle i}^{(6)} M_{jk\rangle a} + \frac{1}{4} M_{a\langle i} M_{jk\rangle a}^{(6)} + \frac{1}{5} \epsilon_{ab\langle i} \left[ -12 S_{ja}^{(2)} M_{k\rangle b}^{(3)} - 8 M_{ja}^{(2)} S_{k\rangle b}^{(3)} - 3 S_{ja}^{(1)} M_{k\rangle b}^{(4)} \right. \\
&\left. - 27 M_{ja}^{(1)} S_{k\rangle b}^{(4)} - S_{ja} M_{k\rangle b}^{(5)} - 9 M_{ja} S_{k\rangle b}^{(5)} - \frac{9}{4} S_a M_{jk\rangle b}^{(5)} \right] + \frac{12}{5} S_{\langle i} S_{jk\rangle}^{(4)} \right\} \\
&+ \mathcal{O}\left(\frac{1}{c^6}\right),
\end{aligned}
\quad (95a) $$

$$ 
\begin{aligned}
V_{ij}(U) = S_{ij}^{(2)}(U) &+ \frac{2GM}{c^3} \int_0^{+\infty} d\tau \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{7}{6} \right] S_{ij}^{(4)}(U-\tau) \\
&+ \frac{G}{7c^5} \left\{ 4S_{a\langle i}^{(2)} M_{j\rangle a}^{(3)} + 8M_{a\langle i}^{(2)} S_{j\rangle a}^{(3)} + 17S_{a\langle i}^{(1)} M_{j\rangle a}^{(4)} - 3M_{a\langle i}^{(1)} S_{j\rangle a}^{(4)} + 9S_{a\langle i} M_{j\rangle a}^{(5)} \right. \\
&- 3M_{a\langle i} S_{j\rangle a}^{(5)} - \frac{1}{4} S_a M_{ija}^{(5)} - 7\epsilon_{ab\langle i} S_a S_{j\rangle b}^{(4)} + \frac{1}{2}\epsilon_{ac\langle i} \left[ 3M_{ab}^{(3)} M_{j\rangle bc}^{(3)} + \frac{353}{24} M_{j\rangle bc}^{(2)} M_{ab}^{(4)} \right. \\
&\left. - \frac{5}{12} M_{ab}^{(2)} M_{j\rangle bc}^{(4)} + \frac{113}{8} M_{j\rangle bc}^{(1)} M_{ab}^{(5)} - \frac{3}{8} M_{ab}^{(1)} M_{j\rangle bc}^{(5)} + \frac{15}{4} M_{j\rangle bc} M_{ab}^{(6)} + \frac{3}{8} M_{ab} M_{j\rangle bc}^{(6)} \right] \\
&+ \mathcal{O}\left(\frac{1}{c^6}\right).
\end{aligned}
\quad (95b) $$

$$ 
\begin{aligned}
U_{ijkl}(U) = M_{ijkl}^{(4)}(U) &+ \frac{G}{c^3} \left\{ 2M \int_0^{+\infty} d\tau \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{59}{30} \right] M_{ijkl}^{(6)}(U-\tau) \right. \\
&\left. + \frac{2}{5} \int_0^{+\infty} d\tau M_{\langle ij}^{(3)} M_{kl\rangle}^{(3)}(U-\tau) - \frac{21}{5} M_{\langle ij}^{(5)} M_{kl\rangle} - \frac{63}{5} M_{\langle ij}^{(4)} M_{kl\rangle}^{(1)} - \frac{102}{5} M_{\langle ij}^{(3)} M_{kl\rangle}^{(2)} \right\} \\
&+ \mathcal{O}\left(\frac{1}{c^5}\right),
\end{aligned}
\quad (95c) $$

$$ 
\begin{aligned}
V_{ijk}(U) = S_{ijk}^{(3)}(U) &+ \frac{G}{c^3} \left\{ 2M \int_0^{+\infty} d\tau \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{5}{3} \right] S_{ijk}^{(5)}(U-\tau) \right. \\
&\left. + \frac{1}{10} \epsilon_{ab\langle i} M_{ja}^{(5)} M_{k\rangle b} - \frac{1}{2} \epsilon_{ab\langle i} M_{ja}^{(4)} M_{k\rangle b}^{(1)} - 2S_{\langle i} M_{jk\rangle}^{(4)} \right\}
\end{aligned}
$$

<footnote n="28">²⁸ In all formulas below the STF projection $\langle \rangle$ applies only to the “free” indices denoted $ijkl...$ carried by the moments themselves. Thus the dummy indices such as $abc...$ are excluded from the STF projection.</footnote>
$$ + \mathcal{O}\left(\frac{1}{c^5}\right). \qquad (95d) $$

$$ \begin{aligned} U_{ijklm}(U) = M_{ijklm}^{(5)}(U) &+ \frac{G}{c^3} \left\{ 2M \int_0^{+\infty} d\tau \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{232}{105} \right] M_{ijklm}^{(7)}(U-\tau) \right. \\ &+ \frac{20}{21} \int_0^{+\infty} d\tau M_{\langle ij}^{(3)} M_{klm\rangle}^{(4)}(U-\tau) - \frac{710}{21} M_{\langle ij}^{(3)} M_{klm\rangle}^{(3)} - \frac{265}{7} M_{\langle ijk}^{(2)} M_{lm\rangle}^{(4)} - \frac{120}{7} M_{\langle ij}^{(2)} M_{klm\rangle}^{(4)} \\ &\left. - \frac{155}{7} M_{\langle ijk}^{(1)} M_{lm\rangle}^{(5)} - \frac{41}{7} M_{\langle ij}^{(1)} M_{klm\rangle}^{(5)} - \frac{34}{7} M_{\langle ijk} M_{lm\rangle}^{(6)} - \frac{15}{7} M_{\langle ij} M_{klm\rangle}^{(6)} \right\} \\ &+ \mathcal{O}\left(\frac{1}{c^4}\right), \end{aligned} \qquad (95e) $$

$$ \begin{aligned} V_{ijkl}(U) = S_{ijkl}^{(4)}(U) &+ \frac{G}{c^3} \left\{ 2M \int_0^{+\infty} d\tau \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{119}{60} \right] S_{ijkl}^{(6)}(U-\tau) \right. \\ &- \frac{35}{3} S_{\langle ij}^{(2)} M_{kl\rangle}^{(3)} - \frac{25}{3} M_{\langle ij}^{(2)} S_{kl\rangle}^{(3)} - \frac{65}{6} S_{\langle ij}^{(1)} M_{kl\rangle}^{(4)} - \frac{25}{6} M_{\langle ij}^{(1)} S_{kl\rangle}^{(4)} - \frac{19}{6} S_{\langle ij} M_{kl\rangle}^{(5)} \\ &- \frac{11}{6} M_{\langle ij} S_{kl\rangle}^{(5)} - \frac{11}{12} S_{\langle i} M_{jkl\rangle}^{(5)} + \frac{1}{6} \epsilon_{ab\langle i} \left[ -5M_{ja}^{(3)} M_{kl\rangle b}^{(3)} - \frac{11}{2} M_{ja}^{(4)} M_{kl\rangle b}^{(2)} - \frac{5}{2} M_{ja}^{(2)} M_{kl\rangle b}^{(4)} \right. \\ &\left. - \frac{1}{2} M_{ja}^{(5)} M_{kl\rangle b}^{(1)} + \frac{37}{10} M_{ja}^{(1)} M_{kl\rangle b}^{(5)} + \frac{3}{10} M_{ja}^{(6)} M_{kl\rangle b}^{(6)} + \frac{1}{2} M_{ja} M_{kl\rangle b}^{(6)} \right] \\ &+ \mathcal{O}\left(\frac{1}{c^4}\right). \end{aligned} \qquad (95f) $$

For all the other multipole moments in the 3PN waveform, it is sufficient to assume the agreement between the radiative and canonical moments, namely

$$ U_L(U) = M_L^{(\ell)}(U) + \mathcal{O}\left(\frac{1}{c^3}\right), \qquad (96a) $$

$$ V_L(U) = S_L^{(\ell)}(U) + \mathcal{O}\left(\frac{1}{c^3}\right). \qquad (96b) $$

In a second stage of the general formalism, we must express the canonical moments {$M_L, S_L$} in terms of the six types of source moments {$I_L, J_L, W_L, X_L, Y_L, Z_L$}. For the control of the (2, 2) mode in the waveform up to 3.5PN order, we need to relate the canonical quadrupole moment $M_{ij}$ to the corresponding source quadrupole moment $I_{ij}$ up to that accuracy. We obtain [197]

$$ \begin{aligned} M_{ij} = I_{ij} &+ \frac{4G}{c^5} [W^{(2)}I_{ij} - W^{(1)}I_{ij}^{(1)}] \\ &+ \frac{4G}{c^7} \left[ \frac{4}{7} W_{a\langle i}^{(1)} I_{j\rangle a}^{(3)} + \frac{6}{7} W_{a\langle i} I_{j\rangle a}^{(4)} - \frac{1}{7} I_{a\langle i} Y_{j\rangle a}^{(3)} - Y_{a\langle i} I_{j\rangle a}^{(3)} - 2XI_{ij}^{(3)} - \frac{5}{21} I_{ija} W_a^{(4)} \right. \\ &+ \frac{1}{63} I_{ija}^{(1)} W_a^{(3)} - \frac{25}{21} I_{ija} Y_a^{(3)} - \frac{22}{63} I_{ija}^{(1)} Y_a^{(2)} + \frac{5}{63} Y_a^{(1)} I_{ija}^{(2)} + 2W_{ij} W^{(3)} \\ &+ 2W_{ij}^{(1)} W^{(2)} - \frac{4}{3} W_{\langle i} W_{j\rangle}^{(3)} + 2Y_{ij} W^{(2)} - 4W_{\langle i} Y_{j\rangle}^{(2)} \\ &\left. + \epsilon_{ab\langle i} \left( \frac{1}{3} I_{j\rangle a} Z_b^{(3)} + \frac{4}{9} J_{j\rangle a} W_b^{(3)} - \frac{4}{9} J_{j\rangle a} Y_b^{(2)} + \frac{8}{9} J_{j\rangle a}^{(1)} Y_b^{(1)} + Z_a I_{j\rangle b}^{(3)} \right) \right] + \mathcal{O}\left(\frac{1}{c^8}\right). \end{aligned} \qquad (97) $$

Here, for instance, W denotes the monopole moment associated with the moment $W_L$, and $Y_i$ is the dipole moment corresponding to $Y_L$. Notice that the difference between the canonical and source moments starts at the relatively high 2.5PN order. For the control of the full waveform
up to 3PN order we need also the moments $M_{ijk}$ and $S_{ij}$, which admit similarly some correction terms starting at the 2.5PN order:

$$M_{ijk} = I_{ijk} + \frac{4G}{c^5} [W^{(2)}I_{ijk} - W^{(1)}I_{ijk}^{(1)} + 3I_{\langle ij \rangle k} Y_{\langle k \rangle}^{(1)}] + \mathcal{O}\left(\frac{1}{c^6}\right), \qquad (98a)$$

$$S_{ij} = J_{ij} + \frac{2G}{c^5} \left[ \epsilon_{ab\langle i} \left( -I_{j\rangle b}^{(3)} W_a - 2I_{j\rangle b} Y_a^{(2)} + I_{j\rangle b}^{(1)} Y_a^{(1)} \right) + 3J_{\langle i} Y_{j\rangle}^{(1)} - 2J_{ij}^{(1)} W^{(1)} \right] \\ + \mathcal{O}\left(\frac{1}{c^6}\right). \qquad (98b)$$

The remainders in Eqs. (98) are consistent with the 3PN approximation for the full waveform. Besides the mass quadrupole moment (97), and mass octopole and current quadrupole moments (98), we can state that, with the required 3PN precision, all the other moments $M_L$, $S_L$ agree with their source counterparts $I_L$, $J_L$:

$$M_L = I_L + \mathcal{O}\left(\frac{1}{c^5}\right), \qquad (99a)$$

$$S_L = J_L + \mathcal{O}\left(\frac{1}{c^5}\right). \qquad (99b)$$

With those formulas we have related the radiative moments {$U_L, V_L$} parametrizing the asymptotic waveform (66) to the six types of source multipole moments {$I_L, J_L, W_L, X_L, Y_L, Z_L$}. What is missing is the explicit dependence of the source moments as functions of the actual parameters of some matter source. We come to grips with this important question in the next section.

# 4 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 include the Schwarzschild and Kerr solutions for black holes, as well as all their perturbations. By Theorem 4 we learned how to construct the radiative moments at infinity, which constitute the observables of the radiation field at large distances from the source, and we obtained in Section 3.3 explicit relationships between radiative and source moments. 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.

## 4.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$. In our formalism the multipole expansion is defined by the multipolar-post-Minkowskian (MPM) solution; see Section 2. Matching together the post-Newtonian and MPM solutions in this overlapping region is an application of the method of matched asymptotic expansions, which has frequently been applied in the present context, both for radiation-reaction [114, 113, 7, 58, 43] and wave-generation [59, 155, 44, 49] problems.
Let us denote by $\mathcal{M}(h)$ the multipole expansion of $h$ (for simplicity, we suppress the space-time indices). By $\mathcal{M}(h)$ we really mean the MPM exterior metric that we have constructed in Sections 2.2 and 2.3:

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

This solution is formally 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 \quad \Rightarrow \quad \mathcal{M}(h) = h. \qquad (101) $$

By contrast, inside the source, $h$ and $\mathcal{M}(h)$ disagree with each other because $h$ is a fully-fledged solution of the field equations within the matter source, while $\mathcal{M}(h)$ is a vacuum solution becoming singular at $r = 0$. 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. (54). In the matching region, where both the multipolar and post-Newtonian expansions are valid, we write the numerical equality

$$ a < r \ll \lambda \quad \Rightarrow \quad \mathcal{M}(h) = \bar{h}. \qquad (102) $$

This “numerical” equality is viewed here in a sense of formal expansions, as we do not control the convergence of the series. In fact, we should be aware that such an equality, though quite natural and even physically obvious, is probably not really justified within the approximation scheme (mathematically speaking), and we simply take it here as part of our fundamental assumptions.

We now transform Eq. (102) 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 3, see Eq. (53). 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 (103) $$

by which we really mean an infinite set of functional identities, valid $\forall (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 somewhat different from that of a numerical equality like Eq. (102), which is valid only when $x$ belongs to some limited spatial domain. The matching equation (103) 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 former 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. (53):

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

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. (104) –, 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 Sections 4.2 and 5.2, the matching equation (103), supplemented by the condition of no-incoming radiation [say in the form of Eq. (29)], permits determining all the unknowns of the problem: On the one hand, the external multipolar decomposition $M(h)$, i.e., the explicit expressions of the multipole moments therein (see Sections 4.2 and 4.4); 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 Section 5.2.

## 4.2 General expression of the multipole expansion

**Theorem 5.** Under the hypothesis of matching, Eq. (103), the multipole expansion of the solution of the Einstein field equation outside a post-Newtonian source reads

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

where the “multipole moments” are given by

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

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

*Proof* (see Refs. [44, 49]): First notice where the physical restriction of considering a post-Newtonian source enters this theorem: The multipole moments (106) 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. (105), namely the finite part of the retarded integral of the multipole expansion $M(\Lambda^{\alpha\beta})$:

$$ \Delta^{\alpha\beta} \equiv h^{\alpha\beta} - \mathcal{FP} \Box_{\text{ret}}^{-1} [M(\Lambda^{\alpha\beta})]. \quad (107) $$

From now on we shall generally abbreviate the symbols concerning the finite-part operation at $B=0$ by a mere $\mathcal{FP}$. According to Eq. (30), $h^{\alpha\beta}$ is given by the retarded integral of the pseudo-tensor $\tau^{\alpha\beta}$. So,

$$ \Delta^{\alpha\beta} = \frac{16\pi G}{c^4} \Box_{\text{ret}}^{-1} \tau^{\alpha\beta} - \mathcal{FP} \Box_{\text{ret}}^{-1} [M(\Lambda^{\alpha\beta})]. \quad (108) $$

In the second term the finite part plays a crucial role because the multipole expansion $M(\Lambda^{\alpha\beta})$ is singular at $r=0$. By contrast, the first term in Eq. (108), as it stands, is well-defined because we are considering only some smooth field distributions: $\tau^{\alpha\beta} \in C^\infty(\mathbb{R}^4)$. There is no need to include a finite part $\mathcal{FP}$ 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 $\mathcal{FP}$ in the first term is that we can re-write Eq. (108) into the more interesting form

$$ \Delta^{\alpha\beta} = \frac{16\pi G}{c^4} \mathcal{FP} \Box_{\text{ret}}^{-1} [\tau^{\alpha\beta} - \mathcal{M}(\tau^{\alpha\beta})], \quad (109) $$

in which we have also used the fact that $M(\Lambda^{\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. (109) 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. (102)]. 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 Ref. [59], yields the second term in Eq. (105),

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

but in which the moments do not yet match the result (106); instead,<fnref n="29" />

$$ \mathcal{H}_L^{\alpha\beta} = \mathcal{F}\mathcal{P} \int d^3\mathbf{x} x_L [\tau^{\alpha\beta} - \mathcal{M}(\tau^{\alpha\beta})]. \quad (111) $$

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 (111) has a compact support limited to the domain of the source. In consequence, we can replace the integrand in Eq. (111) by its post-Newtonian expansion, valid over all the near zone:

$$ \mathcal{H}_L^{\alpha\beta} = \mathcal{F}\mathcal{P} \int d^3\mathbf{x} x_L [\bar{\tau}^{\alpha\beta} - \overline{\mathcal{M}(\tau^{\alpha\beta})}]. \quad (112) $$

Strangely enough, we do not get the expected result because of the presence of the second term in Eq. (112). 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. (104) (indeed $\tau$ and $h$ have the same type of structure). Happily – because we would not know what to do with this term in applications – we are now going to prove that the second term in Eq. (112) is in fact identically zero. The proof is based on the properties of the analytic continuation as applied to the formal structure (104) of $\overline{\mathcal{M}(\tau^{\alpha\beta})}$. Each term of this series yields a contribution to Eq. (112) that takes the form, after performing the angular integration, of the integral $\mathcal{F}\mathcal{P}_{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^\mathcal{R} dr r^{B+b}$ and a “far-zone” one $\int_\mathcal{R}^{+\infty} dr r^{B+b}$, where $\mathcal{R}$ is some constant radius. When $\Re(B)$ is a large enough positive number, the value of the near-zone integral is $\mathcal{R}^{B+b+1}/(B+b+1)$, while when $\Re(B)$ is a large negative number, the far-zone integral reads the opposite, $-\mathcal{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 5:

$$ \mathcal{H}_L^{\alpha\beta} = \mathcal{F}\mathcal{P} \int d^3\mathbf{x} x_L \bar{\tau}^{\alpha\beta}. \quad (113) $$

The latter proof makes it clear how crucial the analytic-continuation finite part $\mathcal{F}\mathcal{P}$ is, which we recall is the same as in our iteration of the exterior post-Minkowskian field [see Eq. (45)]. Without a finite part, the multipole moment (113) 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

<footnote n="29">²⁹ Recall that our abbreviated notation $\mathcal{F}\mathcal{P}$ includes the crucial regularization factor $\bar{\tau}^B$.</footnote>
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 Part B of this article, we must carefully follow the rules for handling the $FP$ operator.

The two terms in the right-hand side of Eq. (105) 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. (42). However, the sum of these two terms, i.e., the exterior multipolar field $M(h)$ itself, is independent of $r_0$. To see this, the simplest way is to differentiate formally $M(h)$ with respect to $r_0$; the differentiations of the two terms of Eq. (105) 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 the $r_0$'s occur. For instance, we have already seen in Eqs. (93)–(94) that the source quadrupole moment $I_{ij}$ depends on $r_0$ starting from the 3PN level, but that this $r_0$ is compensated by another $r_0$ coming from the non-linear “tails of tails” at the 3PN order.

## 4.3 Equivalence with the Will–Wiseman formalism

Will & Wiseman [424] (see also Refs. [422, 335]), extending previous work of Epstein & Wagoner [185] and Thorne [403], have obtained a different-looking multipole decomposition, with different definitions for the multipole moments of a post-Newtonian source. They find, instead of our multipole decomposition given by Eq. (105),

$$ \mathcal{M}(h^{\alpha\beta}) = \Box_{\text{ret}}^{-1} [\mathcal{M}(\Lambda^{\alpha\beta})]_{\mathcal{R}} - \frac{4G}{c^4} \sum_{\ell=0}^{+\infty} \frac{(-)^\ell}{\ell!} \partial_L \left\{ \frac{1}{r} W_L^{\alpha\beta}(t-r/c) \right\}. \quad (114) $$

There is no $FP$ operation in the first term, but instead the retarded integral is *truncated*, as indicated by the subscript $\mathcal{R}$, to extend only in the “far zone”: i.e., $|x'| > \mathcal{R}$ in the notation of Eq. (31), where $\mathcal{R}$ is a constant radius enclosing the source ($\mathcal{R} > a$). The near-zone part of the retarded integral is thereby removed, and there is no problem with the singularity of the multipole expansion $M(\Lambda^{\alpha\beta})$ at the origin. The multipole moments $W_L$ are then given, in contrast with our result (106), by an integral extending over the “near zone” only:

$$ W_L^{\alpha\beta}(u) = \int_{|x'|<\mathcal{R}} d^3\mathbf{x} \, x_L \bar{\tau}^{\alpha\beta}(\mathbf{x}, u). \quad (115) $$

Since the integrand is compact-supported there is no problem with the bound at infinity and the integral is well-defined (no need of a $FP$).

Let us show that the two different formalisms are equivalent. We compute the difference between our moment $H_L$ defined by Eq. (106), and the moment $W_L$ given by Eq. (115). For the comparison we split $H_L$ into far-zone and near-zone pieces corresponding to the radius $\mathcal{R}$. Since the finite part $FP$ present in $H_L$ deals only with the bound at infinity, it can be removed from the near-zone piece, which is then seen to reproduce $W_L$ exactly. So the difference between the two moments is simply given by the far-zone piece:

$$ H_L^{\alpha\beta}(u) - W_L^{\alpha\beta}(u) = FP \int_{|x'|>\mathcal{R}} d^3\mathbf{x} \, x_L \bar{\tau}^{\alpha\beta}(\mathbf{x}, u). \quad (116) $$

We transform next this expression. Successively we write $\bar{\tau}^{\alpha\beta} = M(\bar{\tau}^{\alpha\beta})$ because we are outside the source, and $M(\bar{\tau}^{\alpha\beta}) = \overline{M}(\tau^{\alpha\beta})$ thanks to the matching equation (103). At this stage, we recall from our reasoning right after Eq. (112) that the finite part of an integral over the whole space $\mathbb{R}^3$ of a quantity having the same structure as $\overline{M}(\tau^{\alpha\beta})$ is identically zero by analytic continuation. The main ingredient of the present proof is made possible by this fact, as it allows us to transform
the far-zone integration $|\mathbf{x}| > \mathcal{R}$ in Eq. (116) into a near-zone one $|\mathbf{x}| < \mathcal{R}$, at the price of changing the overall sign in front of the integral. So,

$$ \mathcal{H}_L^{\alpha\beta}(u) - \mathcal{W}_L^{\alpha\beta}(u) = -\mathcal{F}\mathcal{P} \int_{|\mathbf{x}|<\mathcal{R}} d^3\mathbf{x} x_L \overline{\mathcal{M}(\tau^{\alpha\beta})}(\mathbf{x}, u). \quad (117) $$

Finally, it is straightforward to check that the right-hand side of this equation, when summed up over all multipolarities $\ell$, accounts exactly for the near-zone part that was removed from the retarded integral of $\mathcal{M}(\Lambda^{\alpha\beta})$ in the first term in Eq. (114), so that the “complete” retarded integral as given by the first term in our own definition (105) is exactly reconstituted. In conclusion, the formalism of Ref. [424] is equivalent to the one of Refs. [44, 49].

## 4.4 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 5; however, we still have to make the connection with the explicit construction of the general multipolar and post-Minkowskian metric in Section 2. Namely, we must find the expressions of the six STF source multipole moments $I_L, J_L, \dots, Z_L$ parametrizing the linearized metric (35)–(37) at the basis of that construction.30

To do this we first find the equivalent of the multipole expansion given in Theorem 5, 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

$$ \mathcal{M}(h^{\alpha\beta}) = \mathcal{F}\mathcal{P} \Box_{\text{ret}}^{-1} [\mathcal{M}(\Lambda^{\alpha\beta})] - \frac{4G}{c^4} \sum_{\ell=0}^{+\infty} \frac{(-)^\ell}{\ell!} \partial_L \left\{ \frac{1}{r} \mathcal{F}_L^{\alpha\beta}(t-r/c) \right\}, \quad (118) $$

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

$$ \mathcal{F}_L^{\alpha\beta}(u) = \mathcal{F}\mathcal{P} \int d^3\mathbf{x} \hat{x}_L \int_{-1}^{1} dz \delta_{\ell}(z) \bar{\tau}^{\alpha\beta}(\mathbf{x}, u + zr/c). \quad (119) $$

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}, \quad (120) $$

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

$$ \int_{-1}^{1} dz \delta_{\ell}(z) = 1. \quad (121) $$

The next step is to impose the harmonic-gauge conditions (21) onto the multipole decomposition (118), and to decompose the multipole functions $\mathcal{F}_L^{\alpha\beta}(u)$ into STF irreducible pieces with respect to both $L$ and their spatial indices contained into $\alpha\beta = 00, 0i, ij$. This technical part of the calculation is identical to the one of the STF irreducible multipole moments of linearized gravity [154]. The formulas needed in this decomposition read

$$ \mathcal{F}_L^{00} = R_L, \quad (122a) $$

<footnote n="30">30 Recall that in actual applications we need mostly the mass-type moment $I_L$ and current-type one $J_L$, because the other moments simply parametrize a linearized gauge transformation.</footnote>
$$ \mathcal{F}_L^{0i} = {}^{(+)}T_{iL} + \epsilon_{ai<i\ell} {}^{(0)}T_{L-1>a} + \delta_{i<i\ell} {}^{(-)}T_{L-1>}, \quad (122b) $$

$$ \mathcal{F}_L^{ij} = {}^{(+2)}U_{ijL} + \text{STF} \text{STF} [\epsilon_{aii\ell} {}^{(+1)}U_{ajL-1} + \delta_{ii\ell} {}^{(0)}U_{jL-1} \quad (122c) $$

$$ + \delta_{ii\ell} \epsilon_{aji\ell-1} {}^{(-1)}U_{aL-2} + \delta_{ii\ell} \delta_{ji\ell-1} {}^{(-2)}U_{L-2}] + \delta_{ij} V_L, \quad (122d) $$

where the ten tensors $R_L$, ${}^{(+)}T_{L+1}$, ..., ${}^{(-2)}U_{L-2}$, $V_L$ are STF, and are uniquely given in terms of the $\mathcal{F}_L^{\alpha\beta}$'s by some inverse formulas. Finally, the latter decompositions yield the following.

**Theorem 6.** 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)$

$$ I_L(u) = \mathcal{F}\mathcal{P} \int d^3\mathbf{x} \int_{-1}^{1} dz \left\{ \delta_{\ell} \hat{x}_L \Sigma - \frac{4(2\ell+1)}{c^2(\ell+1)(2\ell+3)} \delta_{\ell+1} \hat{x}_{iL} \Sigma_i^{(1)} \right. \\ \left. \qquad + \frac{2(2\ell+1)}{c^4(\ell+1)(\ell+2)(2\ell+5)} \delta_{\ell+2} \hat{x}_{ijL} \Sigma_{ij}^{(2)} \right\} (\mathbf{x}, u+zr/c), \quad (123a) $$

$$ J_L(u) = \mathcal{F}\mathcal{P} \int d^3\mathbf{x} \int_{-1}^{1} dz \epsilon_{ab\langle i\ell} \left\{ \delta_{\ell} \hat{x}_{L-1}{}_{a} \Sigma_{b} - \frac{2\ell+1}{c^2(\ell+2)(2\ell+3)} \delta_{\ell+1} \hat{x}_{L-1}{}_{ac} \Sigma_{bc}^{(1)} \right\} (\mathbf{x}, u+zr/c). \quad (123b) $$

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

In these formulas the notation is as follows: Some convenient source densities are defined from the post-Newtonian expansion (denoted by an overbar) of the pseudo-tensor $\tau^{\alpha\beta}$ by

$$ \Sigma = \frac{\bar{\tau}^{00} + \bar{\tau}^{ii}}{c^2}, \quad (124a) $$

$$ \Sigma_i = \frac{\bar{\tau}^{0i}}{c}, \quad (124b) $$

$$ \Sigma_{ij} = \bar{\tau}^{ij}, \quad (124c) $$

(where $\bar{\tau}^{ii} \equiv \delta_{ij}\bar{\tau}^{ij}$). As indicated in Eqs. (123) all these quantities are to be evaluated at the spatial point $\mathbf{x}$ and at time $u + zr/c$.

For completeness, we give also the formulas for the four auxiliary source moments $W_L, \dots, Z_L$, which parametrize the gauge vector $\varphi_1^\alpha$ as defined in Eqs. (37):

$$ W_L(u) = \mathcal{F}\mathcal{P} \int d^3\mathbf{x} \int_{-1}^{1} dz \left\{ \frac{2\ell+1}{(\ell+1)(2\ell+3)} \delta_{\ell+1} \hat{x}_{iL} \Sigma_i - \frac{2\ell+1}{2c^2(\ell+1)(\ell+2)(2\ell+5)} \delta_{\ell+2} \hat{x}_{ijL} \Sigma_{ij}^{(1)} \right\}, \quad (125a) $$

$$ X_L(u) = \mathcal{F}\mathcal{P} \int d^3\mathbf{x} \int_{-1}^{1} dz \left\{ \frac{2\ell+1}{2(\ell+1)(\ell+2)(2\ell+5)} \delta_{\ell+2} \hat{x}_{ijL} \Sigma_{ij} \right\}, \quad (125b) $$

$$ Y_L(u) = \mathcal{F}\mathcal{P} \int d^3\mathbf{x} \int_{-1}^{1} dz \left\{ -\delta_{\ell} \hat{x}_L \Sigma_{ii} + \frac{3(2\ell+1)}{(\ell+1)(2\ell+3)} \delta_{\ell+1} \hat{x}_{iL} \Sigma_i^{(1)} - \frac{2(2\ell+1)}{c^2(\ell+1)(\ell+2)(2\ell+5)} \delta_{\ell+2} \hat{x}_{ijL} \Sigma_{ij}^{(2)} \right\}, \quad (125c) $$

$$ Z_L(u) = \mathcal{F}\mathcal{P} \int d^3\mathbf{x} \int_{-1}^{1} dz \epsilon_{ab\langle i\ell} \left\{ -\frac{2\ell+1}{(\ell+2)(2\ell+3)} \delta_{\ell+1} \hat{x}_{L-1}{}_{bc} \Sigma_{ac} \right\}. \quad (125d) $$
As discussed in Section 2, 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 (123) and (125), 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 (123) – (125).<fnref n="31" />

These source moments 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) \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} \Sigma(\mathbf{x}, u). \quad (126) $$

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$ (and also $M_L$ and $S_L$) reduce to the standard Newtonian expressions. For instance, $I_{ij}(u) = Q_{ij}(u) + O(1/c^2)$ recovers the Newtonian quadrupole moment (3).<fnref n="32" />

Needless to say, the formalism becomes prohibitively difficult to apply at very high post-Newtonian approximations. Some post-Newtonian 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, which we shall investigate in the next Section 5. Then $\bar{\tau}^{\alpha\beta}$ is to be inserted into the source moments (123) and (125), where the formula (126) 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 Section 9.1 for the case of point-mass binaries. Next, we must substitute the source multipole moments into the linearized metric (35) – (37), and iterate them until all the necessary multipole interactions taking place in the radiative moments $U_L$ and $V_L$ are under control. In fact, we have already worked out these multipole interactions for general sources in Section 3.3 up to the 3PN order in the full waveform, and 3.5PN order for the dominant (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 devoid of symmetries, that would be substantially simpler.

# 5 Interior Field of a Post-Newtonian Source

Theorem 6 solves in principle 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 Section 2. 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 1980’s (see for instance Refs. [6, 181, 269, 270, 334] and also the earlier works [344, 122, 124, 123]), was plagued with some apparently inherent difficulties, which cropped

<footnote n="31">³¹ The work [65] provided some alternative expressions for all the multipole moments (123) – (125), useful for some applications, in the form of surface integrals extending on the outer part of the source’s near zone.</footnote>

<footnote n="32">³² The moments $W_L, \dots, Z_L$ have also a Newtonian limit, but which is not particularly illuminating.</footnote>