# The Quadrupole Moment of Compact Binaries
to the Fourth post-Newtonian Order:
From Source to Canonical Moment

Luc Blanchet<EMAIL_ADDRESS>$\mathcal{G}\mathbb{R}\varepsilon{\mathbb{C}}\mathcal{O}$, Institut
d’Astrophysique de Paris,
UMR 7095, CNRS, Sorbonne Université,
98bis boulevard Arago, 75014 Paris, France  Institut de Physique Théorique,
Université Paris-Saclay,
CEA, CNRS, 91191 Gif-sur-Yvette, France  Guillaume Faye<EMAIL_ADDRESS>$\mathcal{G}\mathbb{R}\varepsilon{\mathbb{C}}\mathcal{O}$, Institut
d’Astrophysique de Paris,
UMR 7095, CNRS, Sorbonne Université,
98bis boulevard Arago, 75014 Paris, France  François Larrouturou
<EMAIL_ADDRESS>Deutsches Elektronen-Synchrotron DESY, Notkestr.
85, 22607 Hamburg, Germany
$\mathcal{G}\mathbb{R}\varepsilon{\mathbb{C}}\mathcal{O}$, Institut
d’Astrophysique de Paris,
UMR 7095, CNRS, Sorbonne Université,
98bis boulevard Arago, 75014 Paris, France

###### Abstract

As a crucial step towards the completion of the fourth post-Newtonian (4PN)
gravitational-wave generation from compact binary systems, we obtain the
expressions of the so-called “canonical” multipole moments of the source in
terms of the “source” and “gauge” moments. The canonical moments describe the
propagation of gravitational waves outside the source’s near zone, while the
source and gauge moments encode explicit information about the matter source.
Those two descriptions, in terms of two sets of canonical moments or in terms
of six sets of source and gauge moments, are isometric. We thus construct the
non-linear diffeomorphism between them up to the third post-Minkowskian order,
and we exhibit the concrete expression of the canonical mass-type quadrupole
moment at the 4PN order. This computation is one of the last missing pieces
for the determination of the gravitational-wave phasing of compact binary
systems at 4PN order.

###### pacs:

04.25.Nx, 04.30.-w, 97.60.Jd, 97.60.Lf

##  I Overview and result

Establishing accurate gravitational wave (GW) templates is crucial for modern
astronomy, as those constitute critical material for the signal analysis of
ground-based detectors and, in the future, airborne ones [1, 2, 3]. A major
technique is the post-Newtonian (PN) approximation, which allows depicting the
inspiralling phase of compact binaries, and constitutes the basis for
effective phenomenological methods such as EOB (effective-one-body) or IMR
(inspiral-merger-ringdown), able to describe the late inspiral and merger
phases (see [4, 5, 6, 7]).

Crucial to building PN waveforms is the knowledge of the mass-type quadrupole
moment at a high accuracy level. Such quantity has been computed, in the case
of non-spinning compact binaries, at the increasingly high 1PN [8, 9], 2PN
[10, 11, 12, 13], 3PN [14, 15, 16, 17] and finally 4PN [18, 19, 20] orders.
Similarly, the mass octupole and current quadrupole have been computed up to
3PN order [21, 22].

At the 4PN order, the mass quadrupole moment has been regularized by means of
dimensional regularization, and both UV and IR divergences have been properly
renormalized [19, 20]. An interesting feature of this 4PN accuracy is the non-
locality in time which appears in the near-zone quantities, due to the
conservative GW tail effect, _i.e._ the backscattering of the radiation onto
the static space-time curvature generated by the source.

Importantly, the previous computations concern the source-type moments, either
mass moments $\mathrm{I}_{L}$ or current ones $\mathrm{J}_{L}$ (where
$L=i_{1}\cdots i_{\ell}$ involves $\ell$ spatial indices). The source moments
are directly connected to the matter distribution in the source, being known
as explicit closed form integrals over the matter plus gravitation pseudo
stress-energy tensor. However, when considering gravitational waves emitted by
the source, it is more convenient to use a different set of multipole moments
called canonical and denoted $\mathrm{M}_{L}$ and $\mathrm{S}_{L}$. The
canonical moments $\mathrm{M}_{L}$ and $\mathrm{S}_{L}$ are directly
associated with the two usual polarization states of GR, which are the two
physical degrees of freedom in the GW propagation.

The goal of the present paper is to connect the canonical moments
$\\{\mathrm{M}_{L},\mathrm{S}_{L}\\}$ to the source moments
$\\{\mathrm{I}_{L},\mathrm{J}_{L}\\}$ and to four additional gauge moments
$\\{\mathrm{W}_{L},\mathrm{X}_{L},\mathrm{Y}_{L},\mathrm{Z}_{L}\\}$, which
parametrize a linear gauge transformation performed in the external part of
the source’s near zone. We refer to [23, 24, 25] for discussions and details
about the matching procedure we employ to link the near zone to the external
zone.

It is important to emphasize that the two descriptions in terms of canonical
moments $\\{\mathrm{M}_{L},\mathrm{S}_{L}\\}$ or in terms of source and gauge
moments
$\\{\mathrm{I}_{L},\mathrm{J}_{L},\mathrm{W}_{L},\mathrm{X}_{L},\mathrm{Y}_{L},\mathrm{Z}_{L}\\}$
are physically equivalent, i.e. describe the same physical matter source if
and only if the canonical moments are related in a precise way to the source
and gauge moments up to arbitrary high orders [23]. Hence there is a
coordinate transformation linking the two descriptions of the same source,
which is a non-linear deformation of the linear gauge transformation
parametrized by the gauge moments
$\\{\mathrm{W}_{L},\mathrm{X}_{L},\mathrm{Y}_{L},\mathrm{Z}_{L}\\}$.

In this paper, we complete a missing step towards the knowledge of the 4PN
orbital phasing for compact binary inspiral by obtaining the relation between
the canonical mass quadrupole moment $\mathrm{M}_{ij}$ and the corresponding
source one $\mathrm{I}_{ij}$ at the 4PN order. Such relation was known
previously at the leading 2.5PN order [23] and at the next-to-leading 3.5PN
order [25]. Up to the 3.5PN order, the correction terms in the canonical
quadrupole are quadratic in the multipole moments; the cubic corrections start
at 4PN, and are thus the aim of the present computation. Let us recapitulate
the complete result for the canonical quadrupole moment up to the 4PN
order111Hereafter, we denote with a capital letter $L$ a multi-index with
$\ell$ indices, _e.g._ $\mathrm{I}_{L}=\mathrm{I}_{i_{1}i_{2}\ldots
i_{\ell}}$, with angular brackets the symmetric-trace-free (STF) projection,
_e.g._ $\mathrm{W}_{a\langle i}\mathrm{I}_{j\rangle
a}\equiv\text{STF}_{ij}[\mathrm{W}_{ai}\mathrm{I}_{ja}]$; we systematically
use STF harmonics [26, 27, 28, 29]; in the spirit of [22], we use the
shorthands $\mathrm{J}_{i|L}\equiv\varepsilon_{ii_{\ell}a}\mathrm{J}_{aL-1}$
and $\mathrm{Z}_{i|L}\equiv\varepsilon_{ii_{\ell}a}\mathrm{Z}_{aL-1}$ for
current type moments; the superscript $(n)$ denotes $n$ time derivatives;
$(-)^{n}$ stands for $(-1)^{n}$ and $c$ is the speed of light and $G$ the
gravitational constant.

$\displaystyle\mathrm{M}_{ij}$
$\displaystyle=\mathrm{I}_{ij}+\frac{4G}{c^{5}}\bigg{[}\mathrm{W}^{(2)}\mathrm{I}_{ij}-\mathrm{W}^{(1)}\mathrm{I}^{(1)}_{ij}\bigg{]}$
$\displaystyle\
+\frac{4G}{c^{7}}\Bigg{\\{}\frac{4}{7}\mathrm{W}^{(1)}_{a\langle
i}\mathrm{I}^{(3)}_{j\rangle a}+\frac{6}{7}\mathrm{W}_{a\langle
i}\mathrm{I}^{(4)}_{j\rangle a}-\frac{1}{7}\mathrm{Y}^{(3)}_{a\langle
i}\mathrm{I}_{j\rangle a}-\mathrm{Y}_{a\langle i}\mathrm{I}^{(3)}_{j\rangle
a}-2\mathrm{X}\,\mathrm{I}^{(3)}_{ij}$
$\displaystyle\qquad\qquad-\frac{5}{21}\mathrm{W}^{(4)}_{a}\mathrm{I}_{ija}+\frac{1}{63}\mathrm{W}^{(3)}_{a}\mathrm{I}^{(1)}_{ija}-\frac{25}{21}\mathrm{Y}^{(3)}_{a}\mathrm{I}_{ija}-\frac{22}{63}\mathrm{Y}^{(2)}_{a}\mathrm{I}^{(1)}_{ija}+\frac{5}{63}\mathrm{Y}^{(1)}_{a}\mathrm{I}^{(2)}_{ija}$
$\displaystyle\qquad\qquad+2\mathrm{W}^{(3)}\mathrm{W}_{ij}+2\mathrm{W}^{(2)}\mathrm{W}^{(1)}_{ij}-\frac{4}{3}\mathrm{W}_{\langle
i}\mathrm{W}^{(3)}_{j\rangle}+2\mathrm{W}^{(2)}\mathrm{Y}_{ij}$
$\displaystyle\qquad\qquad-4\mathrm{W}_{\langle
i}\mathrm{Y}^{(2)}_{j\rangle}-\frac{1}{3}\mathrm{Z}^{(3)}_{a|\langle
i}\mathrm{I}_{j\rangle a}+\mathrm{Z}_{a|\langle i}\mathrm{I}^{(3)}_{j\rangle
a}-\frac{4}{9}\mathrm{W}^{(3)}_{a}\mathrm{J}_{\langle i|j\rangle
a}+\frac{4}{9}\mathrm{Y}^{(2)}_{a}\mathrm{J}_{\langle i|j\rangle
a}-\frac{8}{9}\mathrm{Y}^{(1)}_{a}\mathrm{J}^{(1)}_{\langle i|j\rangle
a}\Bigg{\\}}$
$\displaystyle+\frac{G^{2}\mathrm{M}}{c^{8}}\Bigg{\\{}\frac{106}{21}\mathrm{M}\mathrm{W}_{ij}^{(3)}+\frac{1084}{105}\mathrm{M}\mathrm{Y}_{ij}^{(2)}$
$\displaystyle\qquad\qquad+\frac{17642}{945}\mathrm{W}^{(3)}\mathrm{I}_{ij}-\frac{7018}{315}\mathrm{W}^{(2)}\mathrm{I}_{ij}^{(1)}+\frac{806}{45}\mathrm{W}^{(1)}\mathrm{I}_{ij}^{(2)}-\frac{2098}{135}\mathrm{W}\,\mathrm{I}_{ij}^{(3)}$
$\displaystyle\qquad\qquad+16\,\mathrm{I}_{ij}\int_{0}^{+\infty}\\!\\!\mathrm{d}\tau\,\ln\left(\frac{c\tau}{2r_{0}}\right)\mathrm{W}^{(4)}(t-\tau)-16\,\mathrm{I}_{ij}^{(1)}\int_{0}^{+\infty}\\!\\!\mathrm{d}\tau\,\ln\left(\frac{c\tau}{2r_{0}}\right)\mathrm{W}^{(3)}(t-\tau)\Bigg{\\}}$
$\displaystyle+{\mathcal{O}}\left(\frac{1}{c^{9}}\right)\,.$ (1)

Here $\mathrm{M}$ is the constant (ADM) total mass, $\mathrm{W}_{ij}$ and
$\mathrm{Y}_{ij}$ for instance denote the quadrupoles associated with the
series of gauge moments $\mathrm{W}_{L}$ and $\mathrm{Y}_{L}$, $\mathrm{W}$
and $\mathrm{W}_{i}$ are the monopole and dipole of $\mathrm{W}_{L}$, and so
on. All those quantities are evaluated at time $t$ except when specified
otherwise. In addition, the source is considered to be stationary in the
remote past so that all time derivatives of the moments vanish before some
instant $-\mathcal{T}$. This means in particular that the hereditary integrals
are well defined. The expression (I) is valid in the center-of-mass frame, for
which the mass dipole moment $\mathrm{I}_{i}$ vanishes. The 4PN cubic terms
are new with this paper, which extends the previously computed 3.5PN quadratic
relation [23, 25, 30]. Note that the analogous relations for the mass octupole
and current quadrupole moments (currently known only at the leading 2.5PN
order) do not receive 3PN corrections and can be found in [25].

Interestingly, the relation between the canonical and source/gauge moments is
not local, as clear from the 4PN tail integrals appearing in the last line of
(I). The associated scale $r_{0}$ is unphysical and should naturally disappear
from any physical results, such as the 4PN flux. This will be a stringent test
for our computation, which will however have to wait for the complete
calculation of the 4PN radiative-type moment, directly observable at future
null infinity. Notice however that the tail terms in (I) will be zero in the
case of circular orbits; see Eqs. (5.7) in [30].

Another worthy remark is that the subtleties arising from the use of an IR
dimensional regularization scheme for the mass quadrupole have been completely
tackled and solved in [19, 20]. We can thus safely perform the computation in
three dimensions, using the standard Hadamard regularization scheme.

Finally, after the result (I), only one last step remains before getting the
4PN mass quadrupole: the three-dimensional computation of cubic non-linear
terms called “tails-of-memory”, entering the relation between the radiative
quadrupole moment and the canonical one at 4PN order. This is left to future
work.

The plan of this paper is as follows. After reminders about the multipolar-
post-Minkowskian (MPM) formalism in Sec. II, we describe the general method
for relating the canonical moments to the source and gauge moments up to any
post-Minkowskian (PM) order in Sec. III (extending earlier works in [23, 25]).
Finally Sec. IV is devoted to the practical implementation that led to the
result (I), together with required formulas for retarded integrals of non-
linear source terms. The essential, but technical, near-zone expansion of tail
integrals is presented in details in App. A. The verification of our main
result (I) via an alternative procedure is relegated in App. B.

##  II The Multipolar-post-Minkowskian expansion

The MPM formalism [29, 31, 32, 33] obtains the general solution of the
Einstein field equations outside a matter source in the form of a post-
Minkowskian expansion, with each PM coefficient expanded as a formal
multipolar series. The gothic metric deviation from the Minkowski metric,
$h^{\mu\nu}\equiv\sqrt{-g}\,g^{\mu\nu}-\eta^{\mu\nu}$ where $g^{\mu\nu}$ is
the inverse of the usual covariant metric, $\eta^{\mu\nu}$, that of the
Minkowski metric, and $g\equiv\text{det}(g_{\mu\nu})$, obeys the Einstein’s
vacuum field equations in harmonic coordinates,

$\displaystyle\Box h^{\mu\nu}$ $\displaystyle=\Lambda^{\mu\nu}[h,\partial
h,\partial^{2}h]\,,$ (2a) $\displaystyle\partial_{\nu}h^{\mu\nu}$
$\displaystyle=0\,.$ (2b)

Here $\Box$ is the flat d’Alembertian operator and the gravitational source
term $\Lambda^{\mu\nu}$ is at least quadratic in $h$ and its first and second
partial derivatives. In this paper, we will reason to any PM order for the
general method, but in practical computations, as we are interested in the 3PM
interaction, we will only use the leading quadratic $\sim
h\partial^{2}h+\partial h\partial h$ and sub-leading cubic $\sim h\partial
h\partial h$ pieces in the source term.

###  II.1 The generic MPM algorithm

The “generic” non-linear MPM solution of the field equations is searched in
the form of a non-linear expansion in the field perturbation, labeled by the
gravitational constant $G$. Formally, it may be represented to arbitrary high
orders by the asymptotic series:

$h_{\text{gen}}^{\mu\nu}=\sum_{n=1}^{+\infty}G^{n}h^{\mu\nu}_{\text{gen}\,n}\,.$
(3)

The starting point is the most general solution of the linearized Einstein’s
vacuum equations in harmonic coordinates, $\Box
h_{\text{gen}\,1}^{\mu\nu}=\partial_{\nu}h_{\text{gen}\,1}^{\mu\nu}=0$, which
can be written in terms of six sets of symmetric-trace-free (STF) multipole
moments
$\\{\mathrm{I}_{L},\mathrm{J}_{L},\mathrm{W}_{L},\mathrm{X}_{L},\mathrm{Y}_{L},\mathrm{Z}_{L}\\}$,
dispatched between a simpler linear solution called “canonical”, and a linear
gauge transformation, as

$h_{\text{gen}\,1}^{\mu\nu}=h_{\text{can}\,1}^{\mu\nu}+\partial\varphi_{1}^{\mu\nu}\,.$
(4)

We employ the shorthand notation
$\partial\varphi_{1}^{\mu\nu}\equiv\partial^{\mu}\varphi_{1}^{\nu}+\partial^{\nu}\varphi_{1}^{\mu}-\eta^{\mu\nu}\,\partial_{\rho}\varphi_{1}^{\rho}$
for the linear gauge transformation. Denoting functionals of multipole moments
by means of capital calligraphic letters, the functional dependence of the two
terms in (4) are

$\displaystyle h_{\text{can}\,1}^{\mu\nu}$
$\displaystyle=\mathcal{H}_{\text{can}\,1}^{\mu\nu}\left[\mathrm{I}_{L},\mathrm{J}_{L}\right]\,,$
(5a) $\displaystyle\partial\varphi_{1}^{\mu\nu}$
$\displaystyle=\partial\Phi_{1}^{\mu\nu}\left[\mathrm{W}_{L},\mathrm{X}_{L},\mathrm{Y}_{L},\mathrm{Z}_{L}\right]\,.$
(5b)

For obvious reasons, the moments $\\{\mathrm{I}_{L},\mathrm{J}_{L}\\}$ (mass-
type $I_{L}$ and current-type $J_{L}$) are called the source moments while the
moments $\\{\mathrm{W}_{L},\mathrm{X}_{L},\mathrm{Y}_{L},\mathrm{Z}_{L}\\}$
are the gauge moments. The linear canonical solution, evaluated at field point
$\mathbf{x}$ and at time $t$, reads explicitly [34, 27, 28, 29]

$\displaystyle h_{\text{can}\,1}^{00}$
$\displaystyle=-\frac{4}{c^{2}}\sum_{\ell\geqslant
0}\frac{(-)^{\ell}}{\ell!}\partial_{L}\left[\frac{1}{r}\,\mathrm{I}_{L}\left(t-\frac{r}{c}\right)\right]\,,$
(6a) $\displaystyle h_{\text{can}\,1}^{0i}$
$\displaystyle=\frac{4}{c^{3}}\sum_{\ell\geqslant
1}\frac{(-)^{\ell}}{\ell!}\left\\{\partial_{L-1}\left[\frac{1}{r}\,\mathrm{I}_{iL-1}^{(1)}\left(t-\frac{r}{c}\right)\right]+\frac{\ell}{\ell+1}\partial_{L}\left[\frac{1}{r}\,\mathrm{J}_{i|L}\left(t-\frac{r}{c}\right)\right]\right\\}\,,$
(6b) $\displaystyle h_{\text{can}\,1}^{ij}$
$\displaystyle=-\frac{4}{c^{4}}\sum_{\ell\geqslant
2}\frac{(-)^{\ell}}{\ell!}\left\\{\partial_{L-2}\left[\frac{1}{r}\,\mathrm{I}_{ijL-2}^{(2)}\left(t-\frac{r}{c}\right)\right]+\frac{2\ell}{\ell+1}\partial_{L-1}\left[\frac{1}{r}\,\mathrm{J}_{(i|j)L-1}^{(1)}\left(t-\frac{r}{c}\right)\right]\right\\}\,,$
(6c)

where $r=|\mathbf{x}|$ represents the radial distance to the origin located in
the source. From the harmonic coordinate condition
$\partial_{\nu}h_{\text{can}\,1}^{\mu\nu}=0$, the mass monopole $I$ and
current dipole $J_{i}$ must be constant, while the mass dipole $I_{i}$ is
varying linearly with time. In applications, we choose a center-of-mass frame
for which $I_{i}=0$.222Note that the mass dipole $I_{i}$ is defined here in an
ADM sense and, as such, includes both matter and gravitational contributions.
Thus the displacement of the center-of-mass due to gravitational radiation,
which concerns the matter system and begins at the 3.5PN order (see e.g.
[35]), is already included in $I_{i}$ and is compensated by the contribution
of radiation. With these exceptions, the moments are arbitrary functions of
time encoding the properties of the source. The linearized gauge vector is
defined by

$\displaystyle\varphi_{1}^{0}=$ $\displaystyle\
\frac{4}{c^{3}}\sum_{\ell\geqslant
0}\frac{(-)^{\ell}}{\ell!}\partial_{L}\left[\frac{1}{r}\,\mathrm{W}_{L}\left(t-\frac{r}{c}\right)\right]\,,$
(7a) $\displaystyle\varphi_{1}^{i}=$
$\displaystyle-\frac{4}{c^{4}}\sum_{\ell\geqslant
0}\frac{(-)^{\ell}}{\ell!}\partial_{iL}\left[\frac{1}{r}\,\mathrm{X}_{L}\left(t-\frac{r}{c}\right)\right]$
$\displaystyle-\frac{4}{c^{4}}\sum_{\ell\geqslant
1}\frac{(-)^{\ell}}{\ell!}\left\\{\partial_{L-1}\left[\frac{1}{r}\,\mathrm{Y}_{iL-1}\left(t-\frac{r}{c}\right)\right]+\frac{\ell}{\ell+1}\partial_{L}\left[\frac{1}{r}\,\mathrm{Z}_{i|L}\left(t-\frac{r}{c}\right)\right]\right\\}\,.$
(7b)

The gauge moments
$\\{\mathrm{W}_{L},\mathrm{X}_{L},\mathrm{Y}_{L},\mathrm{Z}_{L}\\}$ are
arbitrary functions of time without restriction.

The MPM construction is defined by induction on the PM order $n$ [29]. Suppose
that for some given $n\geqslant 2$, one has obtained the first $n-1$ PM
coefficients $h_{\text{gen}\,m}^{\mu\nu}$, $\forall m\leqslant n-1$. Then, the
next order coefficient $h_{\text{gen}\,n}^{\mu\nu}$ is constructed as follows.
It satisfies

$\displaystyle\Box h_{\text{gen}\,n}^{\mu\nu}$
$\displaystyle=\Lambda_{\text{gen}\,n}^{\mu\nu}\equiv\Lambda^{\mu\nu}_{n}\left[h_{\text{gen}\,1},\cdots,h_{\text{gen}\,n-1}\right]\,,$
(8a) $\displaystyle\partial_{\nu}h_{\text{gen}\,n}^{\mu\nu}$
$\displaystyle=0\,,$ (8b)

where the source term $\Lambda_{\text{gen}\,n}^{\mu\nu}$, being at least
quadratic in $h$, depends only on the previous iterations as indicated. We
first construct a particular retarded solution of the wave equation $\Box
h_{\text{gen}\,n}^{\mu\nu}=\Lambda_{\text{gen}\,n}^{\mu\nu}$ as

$u_{\text{gen}\,n}^{\mu\nu}\equiv\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\left[\left(\frac{r}{r_{0}}\right)^{B}\Lambda_{\text{gen}\,n}^{\mu\nu}\right]\,,$
(9)

where $\Box^{-1}_{\text{ret}}$ denotes the usual retarded inverse
d’Alembertian operator, and the symbol $\mathop{\mathrm{FP}}_{B=0}$ refers to
a specific operation of taking the finite part (FP) in the Laurent expansion
when the complex parameter $B$ tends to zero. This finite part involves the
multiplication of the source term by the regularization factor
$(r/r_{0})^{B}$, where we introduce an arbitrary constant length scale
$r_{0}$. Such FP operation is required for dealing with source terms made with
multipolar expansions like in (6) that are singular at the origin $r=0$. More
generally, the regularized retarded integral operator
$\mathrm{FP}\,\Box^{-1}_{\mathrm{ret}}\equiv\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}(r/r_{0})^{B}$
is well defined when acting on a source term admitting an expansion when $r\to
0$ of the form, for any $N\in\mathbb{N}$,

$\displaystyle
f(\mathbf{x},t)=\sum_{\ell=0}^{+\infty}\,\sum_{a=a_{\text{min}}}^{N}\,\sum_{p=0}^{p_{\text{max}}}f^{L}_{a}(t)\,\hat{n}_{L}\,r^{a}\ln^{p}r+o\left(r^{N}\right)\,,$
(10)

where $\hat{n}_{L}=\text{STF}[n_{i_{1}}n_{i_{2}}\cdots n_{i_{\ell}}]$, with
$n_{i}=x^{i}/r$, is the STF spherical harmonics of order $\ell$, while the sum
boundaries $a_{\text{min}}$ and $p_{\text{max}}$ are integers (depending on
the PM order $n$). For any function in the class (10), we have
$\Box[\mathrm{FP}\,\Box^{-1}_{\mathrm{ret}}f]=f$. In the end of the recursive
process, the structure (10) turns out to be proved by induction.

Because of the regularization scheme, the object $u_{\text{gen}\,n}^{\mu\nu}$
does not satisfy the harmonic gauge condition, but a simple calculation using
the fact that the source term is divergenceless,
$\partial_{\nu}\Lambda_{\text{gen}\,n}^{\mu\nu}=0$, gives

$w_{\text{gen}\,n}^{\mu}\equiv\partial_{\nu}u_{\text{gen}\,n}^{\mu\nu}=\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\left[B\left(\frac{r}{r_{0}}\right)^{B}\frac{n_{i}}{r}\,\Lambda_{\text{gen}\,n}^{\mu
i}\right]\,.$ (11)

Due to the explicit factor $B$, this term is non zero only when the integral
generates a pole $\propto 1/B$ in the Laurent expansion as $B\to 0$. In turn,
the pole arises only from the singular behaviour of the source term as $r\to
0$, which is of the type (10). Furthermore, the coefficient of the pole is
necessarily a homogeneous retarded solution of the wave equation, since the
source term has no pole, hence $\Box w_{\text{gen}\,n}^{\mu}=0$. At this
stage, we apply the MPM “harmonicity” algorithm to construct from
$w_{\text{gen}\,n}^{\mu}$ another homogeneous retarded solution, say

$v_{\text{gen}\,n}^{\mu\nu}\equiv\mathcal{V}^{\mu\nu}\bigl{[}w_{\text{gen}\,n}\bigr{]}\,,$
(12)

satisfying
$\partial_{\mu}v_{\text{gen}\,n}^{\mu\nu}=-w_{\text{gen}\,n}^{\mu}$, together
with $\Box v_{\text{gen}\,n}^{\mu\nu}=0$. The formulas specifying the above
harmonicity algorithm $w^{\mu}\longrightarrow\mathcal{V}^{\mu\nu}[w]$ are
given by e.g. Eqs. (2.11)-(2.12) in [33]. Finally, the metric at the PM order
$n$, now satisfying the full Einstein vacuum equations in harmonic coordinates
at the PM order $n$, is naturally defined as

$h_{\text{gen}\,n}^{\mu\nu}=u_{\text{gen}\,n}^{\mu\nu}+v_{\text{gen}\,n}^{\mu\nu}\,.$
(13)

This construction yields the most general solution $h_{\text{gen}}^{\mu\nu}$
of the Einstein field equations in the vacuum region external to any isolated
matter system [29]. It is obtained in the form of a functional of six sets of
multipole moments,

$h^{\mu\nu}_{\text{gen}}=\sum_{n=1}^{+\infty}G^{n}h^{\mu\nu}_{\text{gen}\,n}=\mathcal{H}^{\mu\nu}_{\text{gen}}\left[\mathrm{I}_{L},\mathrm{J}_{L},\mathrm{W}_{L},\mathrm{X}_{L},\mathrm{Y}_{L},\mathrm{Z}_{L}\right]\,.$
(14)

###  II.2 The canonical MPM algorithm

However, we also know that the general field of an isolated system in GR can
be described by two and only two sets of STF multipole moments, which we call
the “canonical” moments and denote $\\{\mathrm{M}_{L},\mathrm{S}_{L}\\}$.
Notably, these moments describe the GW propagation far away from the source,
and they parametrize the two GW tensorial modes of GR [29]. The description of
the external field of the source in terms of
$\\{\mathrm{M}_{L},\mathrm{S}_{L}\\}$ starts at linear order just by the
canonical metric (6) instead of the generic metric (4), but, in this different
set up, the linear approximation is parametrized by the canonical moments
$\mathrm{M}_{L}$ and $\mathrm{S}_{L}$:

$h_{\text{can}\,1}^{\mu\nu}=\mathcal{H}_{\text{can}\,1}^{\mu\nu}\left[\mathrm{M}_{L},\mathrm{S}_{L}\right]\,.$
(15)

The canonical MPM algorithm proceeds then by induction over the PM order
$n\geqslant 2$ exactly in the same way as before, i.e. following the synthetic
steps

$\displaystyle\left.\begin{array}[]{l}\displaystyle
u_{\text{can}\,n}^{\mu\nu}=\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}\left(\frac{r}{r_{0}}\right)^{B}\Lambda_{\text{can}\,n}^{\mu\nu}\biggr{]}\\\\[5.69046pt]
\displaystyle
w_{\text{can}\,n}^{\mu}=\partial_{\nu}u_{\text{can}\,n}^{\mu\nu}=\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}B\left(\frac{r}{r_{0}}\right)^{B}\frac{n_{i}}{r}\,\Lambda_{\text{can}\,n}^{\mu
i}\biggr{]}\\\\[11.38092pt]
v_{\text{can}\,n}^{\mu\nu}=\mathcal{V}^{\mu\nu}\bigl{[}w_{\text{can}\,n}\bigr{]}\end{array}\right\\}\Longrightarrow~{}h_{\text{can}\,n}^{\mu\nu}=u_{\text{can}\,n}^{\mu\nu}+v_{\text{can}\,n}^{\mu\nu}\,.$
(19)

This results in a full non-linear metric, which is a functional of the
canonical moments only, and represents as well the most general solution of
the Einstein field equations outside the source:

$\displaystyle h^{\mu\nu}_{\text{can}}$
$\displaystyle=\sum_{n=1}^{+\infty}G^{n}h^{\mu\nu}_{\text{can}\,n}=\mathcal{H}^{\mu\nu}_{\text{can}}\left[\mathrm{M}_{L},\mathrm{S}_{L}\right]\,.$
(20)

The next section presents the general method to obtain the relationships
linking the canonical moments $\\{\mathrm{M}_{L},\mathrm{S}_{L}\\}$ to the set
of source and gauge moments
$\\{\mathrm{I}_{L},\mathrm{J}_{L},\cdots,\mathrm{Z}_{L}\\}$, valid in
principle to any PM order. We thus apply this general method to the practical
case of cubic interaction to derive the 4PN relation (I).

##  III Relation between canonical and source/gauge moments

###  III.1 General method

In Sec. II, we constructed two full non-linear MPM solutions, Eqs. (3) and
(20), which both represent the most general solution of the Einstein field
equations in the vacuum region outside the source. Requiring that these two
metrics describe the exterior field of the _same_ physical system, we now
impose that they are physically equivalent, i.e., just differ by a coordinate
transformation. This implies unique relations between the canonical moments
$\\{\mathrm{M}_{L},\mathrm{S}_{L}\\}$ and the source and gauge moments
$\\{\mathrm{I}_{L},\mathrm{J}_{L},\mathrm{W}_{L},\mathrm{X}_{L},\mathrm{Y}_{L},\mathrm{Z}_{L}\\}$,
which can be viewed as two physically equivalent sets of moments for
describing the source as seen from its exterior.

We thus look for a coordinate transformation $x^{\mu}\longrightarrow
x^{\prime\mu}$ such that

$h^{\mu\nu}_{\text{gen}}(x^{\prime})=\frac{1}{|J|}\,\frac{\partial
x^{\prime\mu}}{\partial x^{\rho}}\,\frac{\partial x^{\prime\nu}}{\partial
x^{\sigma}}\Bigl{(}h^{\rho\sigma}_{\text{can}}(x)+\eta^{\rho\sigma}\Bigr{)}-\eta^{\mu\nu}\,,$
(21)

where $J\equiv\text{det}(\partial x^{\prime}/\partial x)$ is the Jacobian of
the transformation. Eq. (21) immediately follows from the definition
$h^{\mu\nu}=\sqrt{-g}\,g^{\mu\nu}-\eta^{\mu\nu}$ and the law of transformation
of tensors [36]. Introducing a coordinate shift $\varphi^{\mu}$ such that
$x^{\prime\mu}=x^{\mu}+\varphi^{\mu}(x)$, we can rewrite the statement (21),
using the non-linear correction $\delta_{\varphi}h^{\mu\nu}_{\text{can}}(x)$
to the metric $h^{\mu\nu}_{\text{can}}(x)$ induced by the shift, as

$h^{\mu\nu}_{\text{gen}}(x)=h^{\mu\nu}_{\text{can}}(x)+\delta_{\varphi}h^{\mu\nu}_{\text{can}}(x)\,.$
(22)

It is implicit that we work perturbatively in both the vector
$\varphi^{\lambda}$ and the metric $h^{\lambda\rho}_{\text{can}}$, so that we
have for instance $h^{\mu\nu}_{\text{gen}}(x^{\prime})=\sum_{n\geqslant
0}\varphi^{\lambda_{1}}\cdots\varphi^{\lambda_{n}}\partial_{\lambda_{1}}\cdots\partial_{\lambda_{n}}h^{\mu\nu}_{\text{gen}}(x)/n!$.
To linear order, the correction reduces to
$\partial\varphi^{\mu\nu}=\partial^{\mu}\varphi^{\nu}+\partial^{\nu}\varphi^{\mu}-\eta^{\mu\nu}\,\partial_{\rho}\varphi^{\rho}$.
Let us then pose, to any order,

$\delta_{\varphi}h^{\mu\nu}_{\text{can}}=\partial\varphi^{\mu\nu}+\Omega^{\mu\nu}\left[\varphi,h_{\text{can}}\right]\,,$
(23)

where $\Omega^{\mu\nu}$ denotes a functional of $\varphi^{\lambda}$ and
$h^{\lambda\rho}_{\text{can}}$, as well as their derivatives, which is at
least quadratic and can be computed perturbatively up to any order using Eq.
(21). The harmonic gauge condition satisfied by both metrics
$h^{\mu\nu}_{\text{gen}}$ and $h^{\mu\nu}_{\text{can}}$ implies that (as a
consequence of the identity
$\partial_{\nu}\partial\varphi^{\mu\nu}=\Box\varphi^{\mu}$)

$\displaystyle\Box\varphi^{\mu}+\partial_{\nu}\Omega^{\mu\nu}=0\,.$ (24)

Let us now look for the coordinate shift in the form of a full PM expansion
series

$\varphi^{\mu}=\sum_{n=1}^{+\infty}G^{n}\varphi^{\mu}_{n}\,,$ (25)

and, conjointly, for the relations between the canonical moments and the
source/gauge moments in the same PM form

$\mathrm{M}_{L}=\sum_{n=1}^{+\infty}G^{n-1}\mathrm{M}_{n\,L}\,,\qquad\mathrm{S}_{L}=\sum_{n=1}^{+\infty}G^{n-1}\mathrm{S}_{n\,L}\,.$
(26)

Here $\mathrm{M}_{n\,L}$ and $\mathrm{S}_{n\,L}$ denote some $n$-th non-linear
functionals of the six types of source and gauge moments $\mathrm{I}_{K}$,
$\cdots$, $\mathrm{Z}_{K}$ (with, say, $K=i_{1}\cdots i_{k}$), which are to be
determined:

$\mathrm{M}_{n\,L}=\mathcal{M}_{n\,L}\left[\mathrm{I}_{K},\mathrm{J}_{K},\cdots,\mathrm{Z}_{K}\right]\,,\qquad\mathrm{S}_{n\,L}=\mathcal{S}_{n\,L}\left[\mathrm{I}_{K},\mathrm{J}_{K},\cdots,\mathrm{Z}_{K}\right]\,.$
(27)

At linear order, as we have seen with Eq. (4), by definition of the two linear
approximations for the two generic and canonical metrics,

$h_{\text{gen}\,1}^{\mu\nu}=\mathcal{H}_{\text{can}\,1}^{\mu\nu}\left[\mathrm{I}_{L},\mathrm{J}_{L}\right]+\partial\varphi_{1}^{\mu\nu}\,,$
(28)

where the functional $\mathcal{H}_{\text{can}\,1}^{\mu\nu}$ is explicitly
given by (6) and the gauge vector is parametrized by the gauge moments
$\\{\mathrm{W}_{L},\cdots,\mathrm{Z}_{L}\\}$ [see its expression (7)]. Eq.
(28) means that the relation (22) is satisfied at leading order provided that
(i)
$\varphi_{1}^{\mu}=\Phi_{1}^{\mu}[\mathrm{W}_{L},\mathrm{X}_{L},\mathrm{Y}_{L},\mathrm{Z}_{L}]$
and (ii) the moments $\\{\mathrm{M}_{L},\mathrm{S}_{L}\\}$ reduce to the
source moments $\\{\mathrm{I}_{L},\mathrm{J}_{L}\\}$ to leading order, i.e.

$\displaystyle\mathrm{M}_{1\,L}=\mathrm{I}_{L}\,,\qquad\mathrm{S}_{1\,L}=\mathrm{J}_{L}\,.$
(29)

Thus, we see that the looked-for coordinate transformation will be a non-
linear deformation of the linear gauge transformation associated with the
shift (7).

Our recurrence hypothesis will be that all the $\varphi^{\mu}_{m}$’s in (25)
are known up to a given PM order $n-1$, as are the functional relations (26)
up to the corresponding order. Hence we assume that we have already determined

$\displaystyle\varphi^{\mu}_{\leqslant n-1}$
$\displaystyle\equiv\sum_{m=1}^{n-1}G^{m}\varphi^{\mu}_{m}\,,$ (30a)
$\displaystyle\mathrm{M}_{\leqslant n-1\,L}$
$\displaystyle\equiv\sum_{m=1}^{n-1}G^{m-1}\mathcal{M}_{m\,L}\left[\mathrm{I}_{K},\mathrm{J}_{K},\cdots,\mathrm{Z}_{K}\right]\,,$
(30b) $\displaystyle\mathrm{S}_{\leqslant n-1\,L}$
$\displaystyle\equiv\sum_{m=1}^{n-1}G^{m-1}\mathcal{S}_{m\,L}\left[\mathrm{I}_{K},\mathrm{J}_{K},\cdots,\mathrm{Z}_{K}\right]\,,$
(30c)

in such a way that the equations (22)–(23) hold up to order $n-1$. This means
that, for all $m\leqslant n-1$,

$h^{\mu\nu}_{\text{gen}\,m}=h^{\mu\nu}_{\text{can}\,m}+\partial\varphi^{\mu\nu}_{m}+\Omega^{\mu\nu}_{m}\bigl{[}\varphi_{1}\cdots\varphi_{m-1};h_{\text{can}\,1}\cdots
h_{\text{can}\,m-1}\bigr{]}\,.$ (31)

Recall that $\Omega^{\mu\nu}_{m}$ is a non-linear, at least quadratic,
functional of the coordinate shift and metric, and is therefore known
following our induction hypothesis. Indeed, using (23) and (28), we get
$\Omega^{\mu\nu}_{1}\equiv 0$. Crucial in our hypothesis is the assumption
that the canonical metric depends on the moments so far determined to order
$n-1$:

$\displaystyle
h^{\mu\nu}_{\text{can}\,m}\equiv\mathcal{H}^{\mu\nu}_{\text{can}\,m}\left[\mathrm{M}_{\leqslant
n-1\,L},\mathrm{S}_{\leqslant n-1\,L}\right]\,.$ (32)

###  III.2 Implementation to $n$PM order

With our recurrence hypothesis, let us see how to the next PM order
$\varphi^{\mu}_{n}$ together with the functionals $\mathcal{M}_{n\,L}$ and
$\mathcal{S}_{n\,L}$ are uniquely determined. We thus want to find
$\varphi^{\mu}_{n}$ and $\mathcal{M}_{n\,L}$, $\mathcal{S}_{n\,L}$ such that

$h^{\mu\nu}_{\text{gen}\,n}=h^{\mu\nu}_{\text{can}\,n}+\partial\varphi^{\mu\nu}_{n}+\Omega^{\mu\nu}_{n}\,,$
(33)

where $h^{\mu\nu}_{\text{gen}\,n}$ and $h^{\mu\nu}_{\text{can}\,n}$ are
defined precisely by the two MPM constructions in Sec. II, and
$h^{\mu\nu}_{\text{can}\,n}$ is now a functional of $\mathrm{M}_{\leqslant
n\,L}$ and $\mathrm{S}_{\leqslant n\,L}$. First of all, note that
$\Omega^{\mu\nu}_{n}$ depends on $\varphi_{k}$ and $h_{\text{can}\,k}$ for
$k\leqslant n-1$ and is already known by our induction hypothesis. Second,
apply the harmonic coordinate conditions on Eq. (33): this shows that while
$\varphi^{\mu}_{n}$ is one of our unknowns, its d’Alembertian
$\Delta^{\mu}_{n}\equiv\Box\varphi^{\mu}_{n}$ is already determined as we have
[see Eq. (24)]

$\Delta^{\mu}_{n}+\partial_{\nu}\Omega^{\mu\nu}_{n}=0\,.$ (34)

The explicit expressions of $\Omega^{\mu\nu}_{n}$ and $\Delta^{\mu}_{n}$ at
quadratic and cubic orders are displayed in Eqs. (53), (58) and (59). At this
stage, we must use the specific definitions of the generic and canonical
metrics defined in Sec. II. Applying the d’Alembertian operator on (33), we
find that the two source terms $\Lambda^{\mu\nu}_{\text{gen}\,n}$ and
$\Lambda^{\mu\nu}_{\text{can}\,n}$ are related by

$\Lambda^{\mu\nu}_{\text{gen}\,n}=\Lambda^{\mu\nu}_{\text{can}\,n}+\partial\Delta^{\mu\nu}_{n}+\Box\Omega^{\mu\nu}_{n}\,.$
(35)

Hence, the particular retarded solutions of the two algorithms,
$u^{\mu\nu}_{\text{gen}\,n}$ and $u^{\mu\nu}_{\text{can}\,n}$, defined
respectively in (9) and (19), satisfy

$u^{\mu\nu}_{\text{gen}\,n}=u^{\mu\nu}_{\text{can}\,n}+\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}\left(\frac{r}{r_{0}}\right)^{B}\Bigr{(}\partial\Delta^{\mu\nu}_{n}+\Box\Omega^{\mu\nu}_{n}\Bigr{)}\biggr{]}\,.$
(36)

Next, we introduce a linear-looking gauge transformation with vector defined
by the retarded integral of $\Delta_{n}^{\mu}$, as

$\phi_{n}^{\mu}=\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}\left(\frac{r}{r_{0}}\right)^{B}\Delta_{n}^{\mu}\biggr{]}\,.$
(37)

This vector satisfies $\Box\phi_{n}^{\mu}=\Delta_{n}^{\mu}$ but is not yet our
looked-for vector $\varphi_{n}^{\mu}$. It a priori differs from it by an
homogeneous retarded solution of the wave equation. Thanks to (37) we can
advantageously rewrite (36) as

$u_{\mathrm{gen}\,n}^{\mu\nu}=u_{\mathrm{can}\,n}^{\mu\nu}+\partial\phi_{n}^{\mu\nu}+\Omega_{n}^{\mu\nu}+X_{n}^{\mu\nu}+Y_{n}^{\mu\nu}\,.$
(38)

The last two terms are the most interesting: they come from the non-
commutation of the finite part of the retarded integral
$\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}$ with the partial
derivative, due to the differentiation of the regularization factor
$(r/r_{0})^{B}$ therein. They are thus given by

$\displaystyle X_{n}^{\mu\nu}$
$\displaystyle\equiv\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\left[\left(\frac{r}{r_{0}}\right)^{B}\Box\Omega_{n}^{\mu\nu}\right]-\Omega_{n}^{\mu\nu}\,,$
(39a) $\displaystyle Y_{n}^{\mu\nu}$
$\displaystyle\equiv\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\left[\left(\frac{r}{r_{0}}\right)^{B}\partial\Delta_{n}^{\mu\nu}\right]-\partial\phi_{n}^{\mu\nu}\,.$
(39b)

Observing that, for the class of multipole-expanded functions $f$ we are
concerned with [see Eq. (10)], the statement
$\Box(\mathrm{FP}\,\Box^{-1}_{\mathrm{ret}}f)=f$ is always correct, we see
that $X_{n}^{\mu\nu}$ and $Y_{n}^{\mu\nu}$ represent the “commutators” of the
operators that appear inside the square brackets:

$\displaystyle X_{n}^{\mu\nu}$
$\displaystyle=\Bigl{[}\mathrm{FP}\,\Box^{-1}_{\mathrm{ret}},\,\Box\Bigr{]}\,\Omega_{n}^{\mu\nu}\,,$
(40a) $\displaystyle Y_{n}^{\mu\nu}$
$\displaystyle=\Bigl{[}\mathrm{FP}\,\Box^{-1}_{\mathrm{ret}},\,\partial\Bigr{]}\,\Delta_{n}^{\mu\nu}\,.$
(40b)

Since the differentiation of the regularization factor $(r/r_{0})^{B}$
produces an extra factor $B$, the quantities $X_{n}^{\mu\nu}$ and
$Y_{n}^{\mu\nu}$ will be non-zero only when the integral develops a pole $\sim
1/B$. In that case, they are necessarily homogeneous (retarded) solutions of
the wave equation: $\Box X_{n}^{\mu\nu}=\Box Y_{n}^{\mu\nu}=0$. It is easy to
figure out that if $X_{n}^{\mu\nu}$ and $Y_{n}^{\mu\nu}$ were actually zero,
the canonical moments $\\{\mathrm{M}_{L},\mathrm{S}_{L}\\}$ would simply agree
with their source counterparts $\\{\mathrm{I}_{L},\mathrm{J}_{L}\\}$. As a
result, the non trivial relations between those moments entirely follow from
the evaluation of the two quantities $X_{n}^{\mu\nu}$ and $Y_{n}^{\mu\nu}$. In
our practical calculations, we reshuffle the commutators (40) and use the
following expressions, exhibiting the explicit factor $B$ in front:

$\displaystyle X_{n}^{\mu\nu}$
$\displaystyle=\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}B\left(\frac{r}{r_{0}}\right)^{B}\biggl{(}-\frac{B+1}{r^{2}}\,\Omega_{n}^{\mu\nu}-\frac{2}{r}\,\partial_{r}\Omega_{n}^{\mu\nu}\biggr{)}\biggr{]}\,,$
(41a) $\displaystyle Y_{n}^{\mu\nu}$
$\displaystyle=\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}B\left(\frac{r}{r_{0}}\right)^{B}\frac{n_{i}}{r}\biggl{(}-\delta^{i\mu}\Delta_{n}^{\nu}-\delta^{i\nu}\Delta_{n}^{\mu}+\eta^{\mu\nu}\Delta_{n}^{i}\biggr{)}\biggr{]}\,.$
(41b)

Next, we carry on the MPM algorithm by computing the divergence of (38). With
Eq. (34), we obtain
$w_{\mathrm{gen}\,n}^{\mu}=w_{\mathrm{can}\,n}^{\mu}+W_{n}^{\mu}$, where we
have posed $U_{n}^{\mu\nu}\equiv X_{n}^{\mu\nu}+Y_{n}^{\mu\nu}$ and
$W_{n}^{\mu}\equiv\partial_{\nu}U_{n}^{\mu\nu}$. It is in fact necessary and
sufficient to apply the harmonicity algorithm $\mathcal{V}^{\mu\nu}$ only to
the divergence of the sum of the two commutators (41). We have

$\displaystyle\left.\begin{array}[]{l}\displaystyle
U_{n}^{\mu\nu}=X_{n}^{\mu\nu}+Y_{n}^{\mu\nu}\\\\[5.69046pt] \displaystyle
W_{n}^{\mu}=\partial_{\nu}U_{n}^{\mu\nu}\\\\[5.69046pt]
V_{n}^{\mu\nu}=\mathcal{V}^{\mu\nu}\bigl{[}W_{n}\bigr{]}\end{array}\right\\}\Longrightarrow~{}h_{\mathrm{gen}\,n}^{\mu\nu}=h_{\mathrm{can}\,n}^{\mu\nu}+\partial\phi_{n}^{\mu\nu}+\Omega_{n}^{\mu\nu}+U_{n}^{\mu\nu}+V_{n}^{\mu\nu}\,.$
(45)

The final step consists in remarking that $H_{n}^{\mu\nu}\equiv
U_{n}^{\mu\nu}+V_{n}^{\mu\nu}$ is not only divergenceless by definition of the
harmonicity algorithm $\mathcal{V}^{\mu\nu}$, but that it is also a retarded
homogeneous solution of the wave equation, since both $U_{n}^{\mu\nu}$ and
$V_{n}^{\mu\nu}$ are separately such. Hence $H_{n}^{\mu\nu}$ satisfies the
linearized vacuum Einstein field equations, i.e. $\Box
H_{n}^{\mu\nu}=\partial_{\nu}H_{n}^{\mu\nu}=0$, to which we know the general
solution. Namely, it can be decomposed in a unique way as

$H_{n}^{\mu\nu}=\mathcal{H}_{\mathrm{can}\,1}^{\mu\nu}\bigl{[}\mathrm{M}_{n\,L},\mathrm{S}_{n\,L}\bigr{]}+\partial\psi_{n}^{\mu\nu}\,,$
(46)

where $\mathcal{H}_{\mathrm{can}\,1}^{\mu\nu}$ denotes the linearized
functional (6) of the moments, but computed with certain moments
$\mathrm{M}_{n\,L}$ and $\mathrm{S}_{n\,L}$, and where
$\partial\psi_{n}^{\mu\nu}$ is some linear-looking gauge transformation. The
associated gauge “vector” $\psi_{n}^{\mu}$ is parametrized in a unique way by
some moments
$\\{\mathrm{W}_{n\,L},\mathrm{X}_{n\,L},\mathrm{Y}_{n\,L},\mathrm{Z}_{n\,L}\\}$.
The vector $\psi_{n}^{\mu}$ represents the homogeneous solution to be added to
$\phi_{n}^{\mu}$ as given by Eq. (37) in order to recover Eq. (33) with the
shift

$\varphi_{n}^{\mu}\equiv\phi_{n}^{\mu}+\psi_{n}^{\mu}\,.$ (47)

On the other hand, it is clear that
$\\{\mathrm{M}_{n\,L},\mathrm{S}_{n\,L}\\}$ represent the looked-for
corrections to the moments to the order $n$. Indeed, by the linearity of the
functional $\mathcal{H}_{\mathrm{can}\,1}^{\mu\nu}$, the first term in (46)
nicely combines with the linearized approximation in our induction hypothesis
(32) to give

$\mathcal{H}_{\mathrm{can}\,1}^{\mu\nu}\bigl{[}\mathrm{M}_{\leqslant
n-1\,L},\mathrm{S}_{\leqslant
n-1\,L}\bigr{]}+G^{n-1}\mathcal{H}_{\mathrm{can}\,1}^{\mu\nu}\bigl{[}\mathrm{M}_{n\,L},\mathrm{S}_{n\,L}\bigr{]}=\mathcal{H}_{\mathrm{can}\,1}^{\mu\nu}\bigl{[}\mathrm{M}_{\leqslant
n\,L},\mathrm{S}_{\leqslant n\,L}\bigr{]}\,.$ (48)

Finally, in any of the non-linear approximations in our induction hypothesis,
i.e. Eq. (32) with $m\geqslant 2$, we are entitled to consistently replace the
previous set of moments $\\{\mathrm{M}_{\leqslant n-1\,L},$
$\mathrm{S}_{\leqslant n-1\,L}\\}$ by the more accurate, newly determined set
$\\{\mathrm{M}_{\leqslant n\,L},\mathrm{S}_{\leqslant n\,L}\\}$, modulo
higher-order PM terms at least of order $\propto G^{n+1}$, which we can
discard to order $n$. Hence we have proved that the $n$PM contribution to the
canonical moments $\\{\mathrm{M}_{L},\mathrm{S}_{L}\\}$ is determined, as is
the $n$PM piece of the coordinate shift $\varphi^{\mu}$, and our recurrence
hypothesis is verified at the next order $n$. The practical implementation at
quadratic and cubic orders is described in Sec. IV. As a verification, we have
also followed an alternative approach, presented in App. B.

###  III.3 Extraction of physical multipole moments

A straightforward procedure permits reading off the expressions of the $n$-th
order corrections to the canonical moments
$\\{\mathrm{M}_{n\,L},\mathrm{S}_{n\,L}\\}$ from Eq. (46), where the left-hand
side $H^{\mu\nu}_{n}=U^{\mu\nu}_{n}+V^{\mu\nu}_{n}$ follows from the algorithm
(45). In fact, as we have seen that
$U_{n}^{\mu\nu}=X_{n}^{\mu\nu}+Y_{n}^{\mu\nu}$ is a retarded vacuum solution
of the wave equation, we can directly, and uniquely, give the desired moments
as functions of the ten sets of retarded STF moments composing
$U_{n}^{\mu\nu}$. We thus pose

$\displaystyle U_{n}^{00}$ $\displaystyle=\sum_{\ell\geqslant
0}\,\partial_{L}\left[\frac{1}{r}A_{L}\left(t-\frac{r}{c}\right)\right]\,,$
(49a) $\displaystyle U_{n}^{0i}$ $\displaystyle=\sum_{\ell\geqslant
0}\,\partial_{iL}\left[\frac{1}{r}B_{L}\left(t-\frac{r}{c}\right)\right]+\sum_{\ell\geqslant
1}\,\partial_{L-1}\left[\frac{1}{r}C_{iL-1}\left(t-\frac{r}{c}\right)\right]$
$\displaystyle+\sum_{\ell\geqslant
1}\,\partial_{L}\left[\frac{1}{r}D_{i|L}\left(t-\frac{r}{c}\right)\right]\,,$
(49b) $\displaystyle U_{n}^{ij}$ $\displaystyle=\sum_{\ell\geqslant
0}\,\partial_{ijL}\left[\frac{1}{r}E_{L}\left(t-\frac{r}{c}\right)\right]+\sum_{\ell\geqslant
0}\,\delta_{ij}\,\partial_{L}\left[\frac{1}{r}F_{L}\left(t-\frac{r}{c}\right)\right]+\sum_{\ell\geqslant
1}\,\partial_{L-1(i}\left[\frac{1}{r}G_{j)L-1}\left(t-\frac{r}{c}\right)\right]$
$\displaystyle+\sum_{\ell\geqslant
1}\,\partial_{L(i}\left[\frac{1}{r}H_{j)|L}\left(t-\frac{r}{c}\right)\right]+\sum_{\ell\geqslant
2}\,\partial_{L-2}\left[\frac{1}{r}K_{ijL-2}\left(t-\frac{r}{c}\right)\right]$
$\displaystyle+\sum_{\ell\geqslant
2}\,\partial_{L-1}\left[\frac{1}{r}L_{(i|j)L-1}\left(t-\frac{r}{c}\right)\right]\,,$
(49c)

and follow the steps (45), successively computing $W_{n}^{\mu}$,
$V_{n}^{\mu\nu}$ and $H_{n}^{\mu\nu}$, which is finally put into the form (46)
on which we read off the physical moments

$\displaystyle\mathrm{M}_{n\,L}$
$\displaystyle=-\frac{c^{2}(-)^{\ell}\ell!}{4}\biggl{[}A_{L}+4\,\frac{B^{(1)}_{L}}{c}+3\,\frac{E^{(2)}_{L}}{c^{2}}+3F_{L}+G_{L}\biggr{]}\,,$
$\displaystyle(\ell\geqslant 2)$ (50a) $\displaystyle\mathrm{S}_{n\,L}$
$\displaystyle=\frac{c^{3}(-)^{\ell}(\ell+1)!}{4\ell}\biggl{[}D_{L}+\frac{H^{(1)}_{L}}{2c}\biggr{]}\,.$
$\displaystyle(\ell\geqslant 2)$ (50b)

For completeness, we also give the corrections to the gauge moments composing
the gauge vector as $\psi_{n}^{\mu}$:

$\displaystyle\mathrm{W}_{n\,L}$
$\displaystyle=\frac{c^{3}(-)^{\ell}\,\ell!}{4}\biggl{[}B_{L}+\frac{E^{(1)}_{L}}{2c}\biggr{]}\,,$
$\displaystyle(\ell\geqslant 0)$ (51a) $\displaystyle\mathrm{X}_{n\,L}$
$\displaystyle=-\frac{c^{4}(-)^{\ell}\,\ell!}{8}\,E_{L}\,,$
$\displaystyle(\ell\geqslant 0)$ (51b) $\displaystyle\mathrm{Y}_{n\,L}$
$\displaystyle=\frac{3c^{4}(-)^{\ell}\,\ell!}{4}\biggl{[}\frac{B^{(1)}_{L}}{c}+\frac{E^{(2)}_{L}}{c^{2}}+F_{L}+\frac{G_{L}}{3}\biggr{]}\,,$
$\displaystyle(\ell\geqslant 1)$ (51c) $\displaystyle\mathrm{Z}_{n\,L}$
$\displaystyle=-\frac{c^{4}(-)^{\ell}\,(\ell+1)!}{8\ell}\,H_{L}\,.$
$\displaystyle(\ell\geqslant 1)$ (51d)

##  IV Practical implementation

Let us apply the previously described procedure to determine the quadratic and
cubic corrections $\mathcal{M}_{2\,ij}$ and $\mathcal{M}_{3\,ij}$ to the mass-
type quadrupole moment $\mathrm{M}_{ij}$ at the 4PN order [see Eqs. (26)].
Previous investigations [23, 25] focused on quadratic interactions and
determined the mass quadrupole at leading order 2.5PN and sub-leading order
3.5PN; such corrections are recalled in Eq. (I). However, in order to obtain
the 4PN correction, we need to derive the cubic interactions. A preliminary
dimensional analysis shows that at 4PN order and in the center-of-mass frame
(where the mass dipole $\mathrm{I}_{i}$ is vanishing), the only multipole
interactions between the source and gauge moments are cubic and necessarily of
the three types:

$\displaystyle\mathrm{M}\times\mathrm{M}\times\mathrm{W}_{ij}\,,\qquad\mathrm{M}\times\mathrm{M}\times\mathrm{Y}_{ij}\,,\qquad\mathrm{M}\times\mathrm{I}_{ij}\times\mathrm{W}\,,$
(52)

where $M$ is the constant ADM mass, $\mathrm{I}_{ij}$ is the source mass
quadrupole moment, and where the monopole $\mathrm{W}$ and the two quadrupoles
$\mathrm{W}_{ij}$ and $\mathrm{Y}_{ij}$ are gauge moments. Note that out of
the center-of-mass frame, a large number of additional interactions would
appear, such as $\mathrm{I}_{i}\times\mathrm{I}_{i}\times\mathrm{W}$,
$\mathrm{M}\times\mathrm{I}_{i}\times\mathrm{W}_{i}$, $\dots$, but those are
not needed in concrete applications.

###  IV.1 Controlling the cubic source terms

The first step towards the practical determination of the cubic couplings is
to successively construct the quadratic and cubic quantities
$\Omega_{n}^{\mu\nu}$ and $\Delta_{n}^{\mu}$ that enter Eqs. (41) for $n=2,3$.
To quadratic order we have

$\displaystyle\Omega_{2}^{\mu\nu}\equiv$
$\displaystyle-\partial_{\rho}\left[\varphi_{1}^{\rho}\left(h_{\text{can}\,1}^{\mu\nu}+\partial\varphi_{1}^{\mu\nu}\right)\right]+2\,\partial_{\rho}\varphi_{1}^{(\mu}\,h_{\text{can}\,1}^{\nu)\rho}$
(53a)
$\displaystyle+\partial^{\rho}\varphi_{1}^{(\mu}\,\partial_{\rho}\varphi_{1}^{\nu)}+\frac{1}{2}\eta^{\mu\nu}\left[\partial_{\rho}\varphi_{1}^{\sigma}\partial_{\sigma}\varphi_{1}^{\rho}-\partial_{\rho}\varphi_{1}^{\rho}\partial_{\sigma}\varphi_{1}^{\sigma}\right]\,,$
$\displaystyle\Delta_{2}^{\mu}\equiv$ $\displaystyle-
h_{\text{can}\,1}^{\rho\sigma}\,\partial_{\rho\sigma}\varphi_{1}^{\mu}\,.$
(53b)

This is valid “on-shell”, as we have used the facts that $\Box
h^{\mu\nu}_{\text{can}\,1}=\Box\varphi_{1}^{\mu}=0$, which hold at linear
order. One can directly verify that [see Eq. (34)]

$\partial_{\nu}\Omega^{\mu\nu}_{2}+\Delta^{\mu}_{2}=0\,.$ (54)

At the linear level, as is clear from Eqs. (6)–(7), $\mathrm{M}$ and
$\mathrm{I}_{ij}$ enter $h_{\text{can}\,1}^{\mu\nu}$ whereas $\mathrm{W}$,
$\mathrm{W}_{ij}$ and $\mathrm{Y}_{ij}$ enter $\varphi_{1}^{\mu}$. Thus,
before tackling the cubic interactions, we naturally require the knowledge of
the complete quadratic interactions $\mathrm{M}\times\mathrm{M}$,
$\mathrm{M}\times\mathrm{I}_{ij}$, $\mathrm{M}\times\mathrm{W}$, _etc_. The
canonical quadratic metrics $h_{\text{can}\,1}^{\mu\nu}$ corresponding to the
interactions $\mathrm{M}\times\mathrm{M}$ and
$\mathrm{M}\times\mathrm{I}_{ij}$ can be found in [32, 37]. The quadratic
couplings between source and gauge moments have been computed in [23, 25].
Notably, this yields [see (50a)]

$\mathrm{M}_{2\,ij}=\frac{4}{c^{5}}\left[\mathrm{W}^{(2)}\mathrm{I}_{ij}-\mathrm{W}^{(1)}\mathrm{I}_{ij}^{(1)}\right]+{\mathcal{O}}\left(\frac{1}{c^{7}}\right)\,,$
(55)

and we have already given the 3.5PN contribution in (I).

What remains to be computed, for insertion into the cubic quantities
$\Omega_{3}^{\mu\nu}$ and $\Delta_{3}^{\mu}$, are the quadratic couplings
contributing to the quadratic coordinate shift $\varphi_{2}^{\mu}$. As we have
seen in (47), the shift is composed of two parts: the first one,
$\phi_{2}^{\mu}$, has been defined in Eq. (37), and is computed by the usual
techniques (see for instance the Appendix A of [33]). The second part,
$\psi_{2}^{\mu}$, is extracted from (46) as quadratic corrections to the four
types of gauge moments; the general result has been provided in Eqs. (51).
Working out $\varphi_{2}^{\mu}=\phi_{2}^{\mu}+\psi_{2}^{\mu}$ for the needed
quadratic interactions $\mathrm{M}\times\mathrm{W}$,
$\mathrm{M}\times\mathrm{W}_{ij}$, $\mathrm{M}\times\mathrm{Y}_{ij}$ and
$\mathrm{I}_{ij}\times\mathrm{W}$, we find

$\displaystyle\varphi^{0}_{\mathrm{M}\times\mathrm{W}}$
$\displaystyle=-\frac{16\mathrm{M}}{c^{7}}\int_{1}^{+\infty}\\!\\!\mathrm{d}y\,Q_{0}(y)\mathrm{W}^{(2)}\left(t-\frac{yr}{c}\right)\,,$
(56a) $\displaystyle\varphi^{i}_{\mathrm{M}\times\mathrm{W}}$
$\displaystyle=0\,,$ (56b)
$\displaystyle\varphi^{0}_{\mathrm{M}\times\mathrm{W}_{ij}}$
$\displaystyle=-\frac{6\mathrm{M}}{c^{7}\,r^{2}}\hat{n}^{ij}\biggl{[}\mathrm{W}^{(2)}_{ij}\left(t-\frac{r}{c}\right)+\frac{r}{c}\,\mathrm{W}^{(3)}_{ij}\left(t-\frac{r}{c}\right)$
(56c)
$\displaystyle\qquad\qquad\qquad+\frac{4r^{2}}{3c^{2}}\int_{1}^{+\infty}\\!\\!\mathrm{d}y\,Q_{2}(y)\mathrm{W}^{(4)}_{ij}\left(t-\frac{yr}{c}\right)\biggr{]}\,,$
$\displaystyle\varphi^{i}_{\mathrm{M}\times\mathrm{W}_{ij}}$
$\displaystyle=0\,,$ (56d)
$\displaystyle\varphi^{0}_{\mathrm{M}\times\mathrm{Y}_{ij}}$
$\displaystyle=0\,,$ (56e)
$\displaystyle\varphi^{i}_{\mathrm{M}\times\mathrm{Y}_{ij}}$
$\displaystyle=-\frac{4\mathrm{M}}{c^{8}\,r}\hat{n}^{j}\left[\mathrm{Y}^{(2)}_{ij}\left(t-\frac{r}{c}\right)+\frac{2r}{c}\int_{1}^{+\infty}\\!\\!\mathrm{d}y\,Q_{1}(y)\mathrm{Y}^{(3)}_{ij}\left(t-\frac{yr}{c}\right)\right]\,,$
(56f) $\displaystyle\varphi^{0}_{\mathrm{I}_{ij}\times\mathrm{W}}$
$\displaystyle=-\frac{2}{c^{8}\,r}\hat{n}^{ij}\bigg{[}\mathrm{I}_{ij}\mathrm{W}^{(3)}+\mathrm{I}_{ij}^{(1)}\mathrm{W}^{(2)}+\mathrm{I}_{ij}^{(2)}\mathrm{W}^{(1)}+\mathrm{I}_{ij}^{(3)}\mathrm{W}\bigg{]}$
(56g)
$\displaystyle\quad-\frac{6}{c^{7}\,r^{2}}\hat{n}^{ij}\bigg{[}\mathrm{I}_{ij}\mathrm{W}^{(2)}-\frac{2}{3}\mathrm{I}_{ij}^{(1)}\mathrm{W}^{(1)}+\mathrm{I}_{ij}^{(2)}\mathrm{W}\bigg{]}\,,$
$\displaystyle\varphi^{i}_{\mathrm{I}_{ij}\times\mathrm{W}}$
$\displaystyle=0\,.$ (56h)

Note the presence of non-local tail integrals, involving $Q_{m}(y)$, the
Legendre functions of the second kind, here defined with a branch cut from
$-\infty$ to 1 and related to the Legendre polynomials $P_{m}(y)$ by

$Q_{m}(y)=\frac{1}{2}P_{m}(y)\,\mathrm{ln}\left(\frac{y+1}{y-1}\right)-\sum^{m}_{j=1}\frac{1}{j}P_{m-j}(y)P_{j-1}(y)\,.$
(57)

We inject the expressions (56) together with the expressions for the canonical
metric, notably $h^{\mu\nu}_{\text{can}\,2}$ corresponding to the interaction
$\mathrm{M}\times\mathrm{I}_{ij}$ which is also non-local and given by (B3) in
[32], into the cubic $\Omega_{3}^{\mu\nu}$ and $\Delta_{3}^{\mu}$ that define
(41). For convenience, we split $\Omega_{3}^{\mu\nu}$ into quadratic-type and
purely cubic interactions:
$\Omega_{3}^{\mu\nu}\equiv\Omega_{12}^{\mu\nu}+\Omega_{21}^{\mu\nu}+\Omega_{111}^{\mu\nu}$,
where

$\displaystyle\Omega_{12}^{\mu\nu}=$
$\displaystyle-\partial_{\rho}\left[\varphi_{1}^{\rho}\left(h_{\text{can}\,2}^{\mu\nu}+\partial\varphi_{2}^{\mu\nu}\right)\right]+2\,\partial_{\rho}\varphi_{1}^{(\mu}\,h_{\text{can}\,2}^{\nu)\rho}$
$\displaystyle+\partial^{\rho}\varphi_{1}^{(\mu}\,\partial_{\rho}\varphi_{2}^{\nu)}+\frac{1}{2}\eta^{\mu\nu}\left[\partial_{\rho}\varphi_{1}^{\sigma}\partial_{\sigma}\varphi_{2}^{\rho}-\partial_{\rho}\varphi_{1}^{\rho}\partial_{\sigma}\varphi_{2}^{\sigma}\right]\,,$
(58a) $\displaystyle\Omega_{21}^{\mu\nu}=$
$\displaystyle-\partial_{\rho}\left[\varphi_{2}^{\rho}\left(h_{\text{can}\,1}^{\mu\nu}+\partial\varphi_{1}^{\mu\nu}\right)\right]+2\,\partial_{\rho}\varphi_{2}^{(\mu}\,h_{\text{can}\,1}^{\nu)\rho}$
$\displaystyle+\partial^{\rho}\varphi_{2}^{(\mu}\,\partial_{\rho}\varphi_{1}^{\nu)}+\frac{1}{2}\eta^{\mu\nu}\left[\partial_{\rho}\varphi_{2}^{\sigma}\partial_{\sigma}\varphi_{1}^{\rho}-\partial_{\rho}\varphi_{2}^{\rho}\partial_{\sigma}\varphi_{1}^{\sigma}\right]\,,$
(58b) $\displaystyle\Omega_{111}^{\mu\nu}=$
$\displaystyle\frac{1}{2}\,\partial_{\rho\sigma}\left(\varphi_{1}^{\rho}\varphi_{1}^{\sigma}\,h_{\text{can}\,1}^{\mu\nu}\right)+h_{\text{can}\,1}^{\rho\sigma}\,\partial_{\rho}\varphi_{1}^{\mu}\,\partial_{\sigma}\varphi_{1}^{\nu}-2\partial_{\rho}\left(\varphi_{1}^{\rho}\,h_{\text{can}\,1}^{\sigma(\mu}\,\partial_{\sigma}\varphi_{1}^{\nu)}\right)$
$\displaystyle+\frac{1}{2}\,\partial_{\rho\sigma}\left(\varphi_{1}^{\rho}\varphi_{1}^{\sigma}\,\partial\varphi_{1}^{\mu\nu}\right)-\partial_{\rho}\left(\varphi_{1}^{\rho}\,\partial_{\sigma}\varphi_{1}^{\mu}\,\partial^{\sigma}\varphi_{1}^{\nu}\right)-\frac{1}{2}\eta^{\mu\nu}\,\varphi_{1}^{\lambda}\partial_{\lambda}\left(\partial_{\rho}\varphi_{1}^{\sigma}\partial_{\sigma}\varphi_{1}^{\rho}-\partial_{\rho}\varphi_{1}^{\rho}\partial_{\sigma}\varphi_{1}^{\sigma}\right)$
$\displaystyle+\frac{1}{3}\eta^{\mu\nu}\left(\partial_{\rho}\varphi_{1}^{\rho}\partial_{\sigma}\varphi_{1}^{\sigma}\partial_{\lambda}\varphi_{1}^{\lambda}-\partial_{\rho}\varphi_{1}^{\sigma}\partial_{\sigma}\varphi_{1}^{\lambda}\partial_{\lambda}\varphi_{1}^{\rho}\right)\,.$
(58c)

Recall that our computations are done on-shell, using the wave equations
satisfied at linear order, $\Box
h_{\text{can}\,1}^{\mu\nu}=\Box\varphi_{1}^{\mu}=0$, and quadratic order,
$\Box h_{\text{can}\,2}^{\mu\nu}=\Lambda_{\text{can}\,2}^{\mu\nu}$. An
important point is that $\Omega_{111}^{\mu\nu}$ only contains
$h_{\text{can}\,1}\times\varphi_{1}\times\varphi_{1}$ and
$\varphi_{1}\times\varphi_{1}\times\varphi_{1}$ terms, but there is no
$h_{\text{can}\,1}\times h_{\text{can}\,1}\times\varphi_{1}$ sector. This fact
allows us to discard it entirely in the practical implementation, since at the
4PN order we have only the three cubic interactions (52) which contain at most
one gauge moment. Similarly, we have for the coordinate shift
$\Delta_{3}^{\mu}\equiv\Delta_{12}^{\mu}+\Delta_{21}^{\mu}+\Delta_{111}^{\mu}$
with

$\displaystyle\Delta_{12}^{\mu}$
$\displaystyle=-h_{\text{can}\,2}^{\rho\sigma}\,\partial_{\rho\sigma}\varphi_{1}^{\mu}\,,$
(59a) $\displaystyle\Delta_{21}^{\mu}$
$\displaystyle=-h_{\text{can}\,1}^{\rho\sigma}\,\partial_{\rho\sigma}\varphi_{2}^{\mu}+\partial_{\rho}\left(\varphi_{1}^{\rho}\Delta_{2}^{\mu}\right)\,,$
(59b) $\displaystyle\Delta_{111}^{\mu}$
$\displaystyle=\partial_{\rho}\left(\varphi_{1}^{\rho}\,h_{\text{can}\,1}^{\rho\sigma}\partial_{\rho\sigma}\varphi_{1}^{\mu}\right)\,.$
(59c)

We have inserted $\Box\varphi_{2}^{\mu}=\Delta_{2}^{\mu}$, which comes from
(37), and used the fact that $\Box\psi_{2}^{\mu}=0$. The latter quantities
satisfy $\partial_{\nu}\Omega_{3}^{\mu\nu}+\Delta_{3}^{\mu}=0$, consistently
with the divergencelessness of the cubic source, and in fact, even separately,

$\partial_{\nu}\Omega_{12}^{\mu\nu}+\Delta_{12}^{\mu}=0\,,\qquad\partial_{\nu}\Omega_{21}^{\mu\nu}+\Delta_{21}^{\mu}=0\,,\qquad\partial_{\nu}\Omega_{111}^{\mu\nu}+\Delta_{111}^{\mu}=0\,.$
(60)

###  IV.2 Retarded integral of the source terms

At this stage, we control the cubic source terms that are required to evaluate
the “commutators” $X_{3}^{\mu\nu}$ and $Y_{3}^{\mu\nu}$ defined by the
formulas (41). Using canonical relations between the Legendre functions, we
find that the terms to be computed fall into two and only two classes:

$\displaystyle\mathcal{I}_{\text{inst}}\equiv\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}\left(\frac{r}{r_{0}}\right)^{B}\\!B^{b}\,\frac{\hat{n}_{L}}{r^{p}}\,F\left(t-\frac{r}{c}\right)\biggr{]}\,,$
(61a)
$\displaystyle\mathcal{I}_{\text{tail}}\equiv\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}\left(\frac{r}{r_{0}}\right)^{B}\\!B^{b}\,\frac{\hat{n}_{L}}{r^{p}}\,G\left(t-\frac{r}{c}\right)\int_{1}^{+\infty}\\!\\!\mathrm{d}y\,Q_{m}(y)\,F\left(t-\frac{yr}{c}\right)\biggr{]}\,.$
(61b)

Here, the functions $F$ and $G$ represent some products of (source or gauge)
multipole moments, with the function $G$ which can be constant when dealing
with the ADM mass. The second type of term $\mathcal{I}_{\text{tail}}$
integrates over a tail term which comes from the tail present in the quadratic
metric for the interaction $\mathrm{M}\times\mathrm{I}_{ij}$, see (B3) in
[32], and those which appeared in the coordinate shift [see Eqs. (56)]. The
crucial point here is the presence of the explicit factors $B$ and $B^{2}$ in
the definitions (41) hence we posed $b=1,2$ in Eqs. (61); $p$ in an integer
and the index of the Legendre function is $m\in[0,4]$; and $\hat{n}_{L}$ is
the usual STF harmonics.

Thanks to the factors $B$ or $B^{2}$, we know that the integrals (61) will
depend only on the behaviour of the source terms when $r\to 0$.333Recall that
in the MPM approach the $r\to+\infty$ limit is harmless as we assume that the
multipole moments, i.e. the functions $F$ and $G$, become stationary in the
remote past. Indeed the regularization process $\mathop{\mathrm{FP}}_{B=0}$
has been introduced to cope with the singular behaviour of the source when
$r\to 0$, and hence only that limit can generate poles $1/B$ or $1/B^{2}$
which can compensate the factors $B$, $B^{2}$ and lead to a finite part when
$B\to 0$. Therefore we are entitled to restrict ourselves to a ball of finite
and even infinitesimal size (say $r<\epsilon$) and replace the source terms in
(61) by their formal Taylor expansion when $r\to 0$ (i.e., the near zone or PN
expansion $r/c\to 0$). Once the near zone expansion of the source is known the
integration can be performed with standard techniques. The result for the
simpler case $\mathcal{I}_{\text{inst}}$ (with purely instantaneous source
term) was already given in Eq. (A.18) of [33]. It is non-zero only when
$p\geqslant\ell+3$ and $b=1$ (only a simple pole $1/B$ can appear in this
case), and reads

$\mathcal{I}_{\text{inst}}=\left\\{\begin{array}[]{ll}\displaystyle\frac{(-)^{p}2^{p-3}(p-3)!}{(p+\ell-2)!(p-\ell-3)!c^{p-\ell-3}}\,\hat{\partial}_{L}\left[\frac{F^{(p-\ell-3)}(t-r/c)}{r}\right]&\quad\text{($b=1$
and $p\geqslant\ell+3$)}\,,\\\\[17.07182pt] \displaystyle 0&\quad\text{($b=2$
or $p<\ell+3$)}\,.\end{array}\right.$ (62)

Next, we deal with the more difficult integral $\mathcal{I}_{\text{tail}}$. By
the previous argument we can replace the function $G(t-r/c)$ by its formal
Taylor expansion when $r\to 0$. The main problem is therefore the control of
the expansion series when $r\to 0$ of the tail integral
$\mathcal{F}_{m}\equiv\int_{1}^{+\infty}\\!\\!\mathrm{d}y\,Q_{m}(y)\,F(t-yr/c)$.
We devote the App. A to this not-so-easy question and state here the general
result for this expansion:

$\displaystyle\mathcal{F}_{m}$ $\displaystyle\stackrel{{\scriptstyle r/c\to
0}}{{~{}=~{}}}\,\sum_{i=0}^{+\infty}\beta_{i}^{m}\,\frac{(-)^{i}}{i!}\left(\frac{r}{c}\right)^{i}F^{(i)}(t)$
(63)
$\displaystyle\qquad+\sum_{j=0}^{+\infty}\frac{(-)^{m}\,c_{j}^{m}}{(m+2j)!}\left(\frac{r}{c}\right)^{m+2j}\int_{0}^{+\infty}\\!\\!\mathrm{d}\tau\biggl{[}\ln\left(\frac{c\tau}{2r}\right)-H_{m+j}+2H_{2m+2j+1}\biggr{]}\,F^{(m+2j+1)}\left(t-\tau\right)\,.$

Here, the symbol $\stackrel{{\scriptstyle r/c\to 0}}{{~{}=~{}}}$ indicates the
formal asymptotic expansion, the coefficients $c_{j}^{m}$ and $\beta_{i}^{m}$
are defined in (73) and (84), and $H_{n}$ denotes the usual harmonic number.
Thus, we have to multiply Eq. (63) by the Taylor expansion of $G(t-r/c)$ and
integrate term by term. Note that all those formal expansions are convergent,
due to the stationary of the source in the remote past. As we see from the
structure of the expansion (63), the final integration will boil down to the
control of just one type of term,

$\mathcal{J}\equiv\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}\left(\frac{r}{r_{0}}\right)^{B}\\!B^{b}\,\frac{\hat{n}_{L}}{r^{p}}\,(\ln
r)^{a}\,H(t)\biggr{]}\,.$ (64)

However, the function $H(t)$ here can be either an instantaneous function or a
non-local tail integral of the type given in (63). Furthermore, still because
of the tail term we must add the case where there is logarithm $\ln r$ in the
source, see again (63). Thus, we consider (64) with $a=0$ or $1$.

We give a few details on the calculation of $\mathcal{J}$. Writing (64) in
ordinary three-dimensional form we have, since as we said the integration is
limited to an infinitesimal ball $r<\epsilon$,

$\mathcal{J}=-\frac{1}{4\pi}\mathop{\mathrm{FP}}_{B=0}\int_{r<\epsilon}\frac{\mathrm{d}^{3}\mathbf{x}^{\prime}}{|\mathbf{x}-\mathbf{x}^{\prime}|}\left(\frac{r^{\prime}}{r_{0}}\right)^{B}B^{b}\,\frac{\hat{n}^{\prime}_{L}}{r^{\prime
p}}\,(\ln
r^{\prime})^{a}\,H\left(t-\frac{|\mathbf{x}-\mathbf{x}^{\prime}|}{c}\right)\,.$
(65)

Using the formal STF expansion of the multipolar factor when
$r^{\prime}\equiv|\mathbf{x}^{\prime}|\to 0$,

$\frac{H(t-|\mathbf{x}-\mathbf{x}^{\prime}|/c)}{|\mathbf{x}-\mathbf{x}^{\prime}|}=\sum_{q=0}^{+\infty}\frac{(-)^{q}}{q!}\sum_{j=0}^{+\infty}\,\frac{r^{\prime
2j}}{2^{j}j!}\,\frac{(2q+1)!!}{(2q+2j+1)!!}\,\hat{x}^{\prime}_{Q}\hat{\partial}_{Q}\left(\frac{H^{(2j)}(t-r/c)}{rc^{2j}}\right)\,,$
(66)

together with the angular integration performed using Eq. (A29a) in [29], we
end up with radial integrals of the type

$\mathop{\mathrm{FP}}_{B=0}\int_{0}^{\epsilon}\mathrm{d}r^{\prime}\,\left(\frac{r^{\prime}}{r_{0}}\right)^{B}B^{b}\,r^{\prime\ell+2j+2-p}\,(\ln
r^{\prime})^{a}=\mathop{\mathrm{FP}}_{B=0}\frac{B^{b}}{r_{0}^{B}}\left(\frac{\mathrm{d}}{\mathrm{d}B}\right)^{a}\left[\frac{\epsilon^{B+\ell+2j+3-p}}{B+\ell+2j+3-p}\right]\,.$
(67)

Here the $r^{\prime}=0$ boundary vanishes by analytic continuation on $B$.
With the factor $B^{b}$ in front the latter integral is zero unless there is a
pole, and the pole can come only when $p$ is of the form $p=\ell+3+2j$ (where
$j\in\mathbb{N}$). Hence the results are immediate.

When $a=0$ (source without $\ln r$) the integral is zero when $b=2$ as there
is only a simple pole $1/B$ in this case:

$\mathcal{J}_{a=0}=\left\\{\begin{array}[]{ll}\displaystyle\frac{(-)^{p}}{(p+\ell-2)!!(p-\ell-3)!!c^{p-\ell-3}}\,\hat{\partial}_{L}\left[\frac{H^{(p-\ell-3)}(t-r/c)}{r}\right]&~{}\text{($b=1$
and $p-\ell-3\in 2\mathbb{N}$)}\,,\\\\[17.07182pt] \displaystyle
0&~{}\text{($b=2$ or $p-\ell-3\notin 2\mathbb{N}$)}\,.\end{array}\right.$
(68a) One can check that this is perfectly consistent with the result for
$\mathcal{I}_{\text{inst}}$ in (62), in the sense that one can Taylor expand
the source term of $\mathcal{I}_{\text{inst}}$ when $r\to 0$ and recover the
same result by applying Eq. (68a) on each term of the Taylor series. By
contrast, when there is a $\ln r$ the integral is also non-zero when $b=2$
because of the presence of a double pole $1/B^{2}$:
$\mathcal{J}_{a=1}=\left\\{\begin{array}[]{ll}\displaystyle\frac{(-)^{p}\,\ln
r_{0}}{(p+\ell-2)!!(p-\ell-3)!!c^{p-\ell-3}}\,\hat{\partial}_{L}\left[\frac{H^{(p-\ell-3)}(t-r/c)}{r}\right]&~{}\text{($b=1$
and $p-\ell-3\in 2\mathbb{N}$)}\,,\\\\[11.38092pt]
\displaystyle\frac{(-)^{p+1}}{(p+\ell-2)!!(p-\ell-3)!!c^{p-\ell-3}}\,\hat{\partial}_{L}\left[\frac{H^{(p-\ell-3)}(t-r/c)}{r}\right]&~{}\text{($b=2$
and $p-\ell-3\in 2\mathbb{N}$)}\,,\\\\[17.07182pt] \displaystyle
0&~{}\text{($p-\ell-3\notin 2\mathbb{N}$)}\,.\end{array}\right.$ (68b)

Finally with those results in hand, we can implement all the terms up to cubic
non-linear order, using the _xAct_ library of the _Mathematica_ software [38].
Such computation ends up with the final result at the 4PN level already
recapitulated in Eq. (I). In the App. B, an alternative procedure is described
which permitted checking it independently.

###### Acknowledgements.

We acknowledge discussions with Laura Bernard, Quentin Henry and David
Trestini. F.L. received funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation program (grant
agreement No 817791).

##  Appendix A Near-zone expansion of the tail integral

This appendix is devoted to the determination of the near zone (or PN)
expansion when $r/c\to 0$ of the tail integral that enters the source of the
integral $\mathcal{I}_{\text{tail}}$ (61b). This will permit justifying the
claim of Eq. (63) concerning the structure of the near zone expansion, and
provide the explicit coefficients we need for the practical computation. Thus
we look for the expansion of

$\mathcal{F}_{m}(r,t)\equiv\int_{1}^{+\infty}\\!\\!\mathrm{d}y\,Q_{m}(y)\,F\left(t-\frac{yr}{c}\right)\,.$
(69)

If an explicit expression of the Legendre function of the second kind
$Q_{m}(y)$ is given in Eq. (57), we preferably use here the general expression
of the Legendre function for generic
$\mu\in\mathbb{C}\setminus\\{-1,-2,\cdots\\}$ in terms of the hypergeometric
function $F\equiv{}_{2}F_{1}$:

$Q_{\mu}(y)=\frac{\sqrt{\pi}}{2^{\mu+1}}\frac{\Gamma\left(\mu+1\right)}{\Gamma\left(\mu+\frac{3}{2}\right)}\,\frac{1}{y^{\mu+1}}\,F\left(\frac{\mu}{2}+1,\frac{\mu+1}{2};\mu+\frac{3}{2};\frac{1}{y^{2}}\right)\,,$
(70)

with $|y|>1$, $|\arg y|<\pi$. This representation of the Legendre function
reads explicitly, for any real argument such that $y>1$,

$\displaystyle Q_{\mu}(y)$
$\displaystyle=\sum_{j=0}^{+\infty}c_{j}^{\mu}\,y^{-1-\mu-2j}\,,$ (71a)
$\displaystyle\text{with}\quad c_{j}^{\mu}$
$\displaystyle\equiv\frac{2^{\mu}}{j!}\frac{\Gamma\left(\mu+2j+1\right)\Gamma\left(\mu+j+1\right)}{\Gamma\left(2\mu+2j+2\right)}\,.$
(71b)

It can naturally be regarded as an asymptotic expansion when $y\to+\infty$ on
the real axis. However, we stress that it is valid as soon as
$y\in]1,+\infty[$, as the series representation of the hypergeometric function
is well defined in that case.

Defining $\tau\equiv yr/c$ one can then rewrite (69) as the following series

$\displaystyle\mathcal{F}_{m}$
$\displaystyle=\sum_{j=0}^{+\infty}\,c_{j}^{m}\,\mathcal{F}_{m}^{j}\,,$ (72a)
$\displaystyle\text{where}\quad\mathcal{F}_{m}^{j}$
$\displaystyle\equiv\left(\frac{r}{c}\right)^{m+2j}\int_{r/c}^{+\infty}\\!\\!\mathrm{d}\tau\,\tau^{-m-2j-1}\,F(t-\tau)\,,$
(72b)

and the coefficients are just given by (71b) in the particular case where
$m\in\mathbb{N}$, i.e.,

$c_{j}^{m}=\frac{2^{m}}{j!}\frac{(m+2j)!(m+j)!}{(2m+2j+1)!}\,.$ (73)

Upon the expression of $\mathcal{F}_{m}^{j}$ we perform a series of
integrations by parts in order to increase the power of $\tau$ until we reach
a logarithm of $\tau$. The all-integrated terms contain the functions
$F^{(n)}(t-r/c)$ that we replace by their Taylor expansions when $r\to 0$. The
last integral containing $\ln\tau$ is split according to
$\int_{r/c}^{+\infty}=-\int_{0}^{r/c}+\int_{0}^{+\infty}$. In the first
integral, from $0$ to $r/c$, we are allowed to expand the integrand when
$\tau\to 0$, since by definition $r/c\to 0$ for the PN expansion. Then the
second integral from $0$ to $+\infty$ just gives a tail integral in the
ordinary sense. Therefore, we obtain the following asymptotic expansion when
$r/c\to 0$:

$\displaystyle\mathcal{F}_{m}^{j}$ $\displaystyle\stackrel{{\scriptstyle
r/c\to 0}}{{~{}=~{}}}\sum_{\begin{subarray}{c}i=0\\\ i\neq
m+2j\end{subarray}}^{+\infty}\frac{(-)^{i}}{i!\,(m+2j-i)}\left(\frac{r}{c}\right)^{i}F^{(i)}\left(t\right)$
$\displaystyle+\frac{(-)^{m}}{(m+2j)!}\left(\frac{r}{c}\right)^{m+2j}\int_{0}^{+\infty}\\!\\!\mathrm{d}\tau\biggl{[}\ln\left(\frac{c\tau}{r}\right)+H_{m+2j}\biggr{]}\,F^{(m+2j+1)}\left(t-\tau\right)\,.$
(74)

Very importantly, the value $i=m+2j$ is to be excluded from the first
summation. The last term is the tail integral, where $H_{n}$ denotes the usual
harmonic number. Resumming on $j$ yields

$\displaystyle\mathcal{F}_{m}\stackrel{{\scriptstyle r/c\to
0}}{{~{}=~{}}}\sum_{i=0}^{+\infty}\alpha_{i}^{m}\,\frac{(-)^{i}}{i!}\left(\frac{r}{c}\right)^{i}F^{(i)}(t)$
$\displaystyle\qquad+\sum_{j=0}^{+\infty}c_{j}^{m}\,\frac{(-)^{m}}{(m+2j)!}\left(\frac{r}{c}\right)^{m+2j}\int_{0}^{+\infty}\\!\\!\mathrm{d}\tau\biggl{[}\ln\left(\frac{c\tau}{r}\right)+H_{m+2j}\biggr{]}\,F^{(m+2j+1)}\left(t-\tau\right)\,.$
(75)

We already see the type of structure claimed in Eq. (63). The coefficients
$\alpha_{i}^{m}$ in the first sum (corresponding to instantaneous terms) are
still at this stage given by an infinite series:

$\alpha_{i}^{m}=\sum_{\begin{subarray}{c}j=0\\\ i\neq
m+2j\end{subarray}}^{+\infty}\frac{c_{j}^{m}}{m+2j-i}\,.$ (76)

Despite that, the result (A) can be dealt with as it is in practical
calculations. However, as it turns out the coefficient (76) can be resummed in
analytic closed form, and this will yield a more interesting and powerful
expression of the near zone expansion of $\mathcal{F}_{m}$.

To obtain such form of $\alpha_{i}^{m}$ we take advantage of the fact that the
coefficients $c_{j}^{m}$ can be generalized to any
$\mu\in\mathbb{C}\setminus\\{-1,-2,\cdots\\}$ through the general definition
of the Legendre function, see Eqs. (71). Therefore, one can extend the
definition of $\alpha_{i}^{m}$ to any generic value of $\mu$ by posing

$\alpha_{i}^{\mu}=\sum_{j=0}^{+\infty}\,\frac{c_{j}^{\mu}}{\mu+2j-i}\,,$ (77)

where the coefficients $c_{j}^{\mu}$ are now given by (71b). We assume that
$\mu$ is non-integral, so that the condition $i\neq m+2j$ is no longer
necessary and has been dropped. Now, from the very definition of the
coefficients $c_{j}^{\mu}$ in (71), the fact that the hypergeometric series is
absolutely convergent for $|y|>1$, and using the identity
$\int_{1}^{+\infty}\mathrm{d}y\,y^{i}\times y^{-1-\mu-2j}=(\mu+2j-i)^{-1}$
valid for $\Re(\mu)>i$, we prove that the coefficient $\alpha_{i}^{\mu}$ is
actually given by

$\alpha_{i}^{\mu}=\int_{1}^{+\infty}\\!\\!\mathrm{d}y\,y^{i}\,Q_{\mu}(y)\,.$
(78)

Furthermore, we know that [39]

$\int_{1}^{+\infty}\mathrm{d}y\,(y-1)^{\nu}\,Q_{\mu}(y)=2^{\nu}\frac{[\Gamma(\nu+1)]^{2}\,\Gamma(\mu-\nu)}{\Gamma(\mu+\nu+2)}\,.$
(79)

Hence we arrive at the closed-form expression

$\alpha_{i}^{\mu}=\sum_{k=0}^{i}{\genfrac{(}{)}{0.0pt}{}{i}{k}}\,2^{k}(k!)^{2}\frac{\Gamma(\mu-k)}{\Gamma(\mu+k+2)}\,,$
(80)

where ${\genfrac{(}{)}{0.0pt}{}{i}{k}}$ is the usual binomial coefficient. Its
validity may be extended to all non-integral values of $\mu$ by analytic
continuation. Still this expression is to be connected to the actual result
(76) we search for, as this result excludes the value $i=m+2j$ from the
summation. However, posing $\mu=m+\varepsilon$ one can also substract and then
re-add the contributions $i=m+2j$ in (77); in this way we rewrite
$\alpha_{i}^{m}$ as the following limit when $\varepsilon\to 0$:

$\alpha_{i}^{m}=\lim_{\varepsilon\to
0}\left\\{\begin{array}[]{ll}\displaystyle\alpha_{i}^{m+\varepsilon}&\qquad\text{when
$i\not=m+2j$}\,,\\\\[17.07182pt]
\displaystyle\alpha_{m+2j}^{m+\varepsilon}-\frac{1}{\varepsilon}\,c_{j}^{m+\varepsilon}&\qquad\text{when
$i=m+2j$}\,.\end{array}\right.$ (81)

In the case $i=m+2j$ an explicit pole $1/\varepsilon$ has to be added, and
which should cancel the pole present in $\alpha_{m+2j}^{m+\varepsilon}$ so
that the limit is finite. Furthermore since the coefficient
$c_{j}^{m+\varepsilon}$ is finite when $\varepsilon\to 0$, with limit
$c_{j}^{m}$ given by (73), we can expand it to order $\varepsilon$ and obtain

$\alpha_{i}^{m}=\lim_{\varepsilon\to
0}\left\\{\begin{array}[]{ll}\displaystyle\alpha_{i}^{m+\varepsilon}&\qquad\text{($i\not=m+2j$)}\,,\\\\[11.38092pt]
\displaystyle\alpha_{m+2j}^{m+\varepsilon}-c_{j}^{m}\Biggl{[}\frac{1}{\varepsilon}+\left(\frac{\mathrm{d}\ln
c_{j}^{\mu}}{\mathrm{d}\mu}\right)_{\mu=m}\Biggr{]}&\qquad\text{($i=m+2j$)}\,.\end{array}\right.$
(82)

Finally, the last step is to use the closed-form expression (80) for
$\alpha_{i}^{m+\varepsilon}$. The extra terms in (82) are easily computed with
the help of (71b) and (73). We verify that indeed the limit $\varepsilon\to 0$
is finite and obtain

$\alpha_{i}^{m}=\beta_{i}^{m}+\left\\{\begin{array}[]{ll}\displaystyle
0&\qquad\text{($i\not=m+2j$)}\,,\\\\[11.38092pt] \displaystyle-
c_{j}^{m}\Bigl{(}\ln
2+H_{m+2j}+H_{m+j}-2H_{2m+2j+1}\Bigr{)}&\qquad\text{($i=m+2j$)}\,,\end{array}\right.$
(83)

where the new coefficient $\beta_{i}^{m}$ reads explicitly

$\displaystyle\\!\\!\beta_{i}^{m}$
$\displaystyle=\sum_{k=0}^{m-1}{\genfrac{(}{)}{0.0pt}{}{i}{k}}\,2^{k}(k!)^{2}\frac{(m-k-1)!}{(m+k+1)!}+\sum_{k=m}^{i}{\genfrac{(}{)}{0.0pt}{}{i}{k}}\,2^{k}(k!)^{2}(-)^{m+k}\frac{H_{k-m}-H_{k+m+1}}{(k-m)!(k+m+1)!}\,,$
(84)

with the convention that $\genfrac{(}{)}{0.0pt}{}{i}{k}=0$ whenever $i<k$.
Hence our final result for the near zone expansion of $\mathcal{F}_{m}$ reads

$\displaystyle\mathcal{F}_{m}\stackrel{{\scriptstyle r/c\to
0}}{{~{}=~{}}}\,\sum_{i=0}^{+\infty}\beta_{i}^{m}\,\frac{(-)^{i}}{i!}\left(\frac{r}{c}\right)^{i}F^{(i)}(t)$
(85)
$\displaystyle\qquad+\sum_{j=0}^{+\infty}c_{j}^{m}\,\frac{(-)^{m}}{(m+2j)!}\left(\frac{r}{c}\right)^{m+2j}\int_{0}^{+\infty}\\!\\!\mathrm{d}\tau\biggl{[}\ln\left(\frac{c\tau}{2r}\right)-H_{m+j}+2H_{2m+2j+1}\biggr{]}\,F^{(m+2j+1)}\left(t-\tau\right)\,.$

##  Appendix B Verification of the result (I) via an alternative procedure

The general method exposed in Sec. III suggests a simplifying modification of
the MPM algorithm to compute the contributions to the metric of interactions
involving a least one gauge moment. In this alternative approach, the
construction of the generic MPM metric $h^{\mu\nu}_{\text{gen}}$, of the
coordinate shift $\varphi^{\mu}$ and of the functional relation (26), between
the canonical and the source/gauge moments, is achieved using a variant of the
recurrence procedure exposed in Sec. III.2.

The initial step remains unchanged. Let us now make the recurrence hypothesis
at order $n$. Namely, we assume that we have already determined
$h^{\mu\nu}_{\text{gen}\,m}$, $\varphi_{m}^{\mu}$, $\mathrm{M}_{m\,L}$, and
$\mathrm{S}_{m\,L}$, for all $m\leqslant n-1$. The equation (35) shows that
the generic MPM source of the field equations is a sum of three terms. The
first one is the source of the canonical metric
$\Lambda^{\mu\nu}_{\text{can}\,n}=\Lambda^{\mu\nu}_{n}[h_{\text{can}\,1},\cdots,h_{\text{can}\,n-1}]$.
While this term enters the derivation of
$\mathcal{H}^{\mu\nu}_{\text{can}\,n}[\mathrm{M}_{\leqslant
n-1\,L},\mathrm{S}_{\leqslant n-1\,L}]$, it does not play any role in the
calculation of $\mathrm{M}_{n\,L}$ and $\mathrm{S}_{n\,L}$, which are
determined by $X^{\mu\nu}_{n}$ and $Y^{\mu\nu}_{n}$ hence by
$\Delta^{\mu\nu}_{n}$ and $\Omega^{\mu\nu}_{n}$. The modification of the
recurrence will thus concern the treatment of the second and third terms,
which define the reduced source of order $n$:

$\displaystyle\delta\Lambda^{\mu\nu}_{\text{can}\,n}\equiv\Box\Omega^{\mu\nu}_{n}+\partial\Delta^{\mu\nu}_{n}\,.$
(86)

This source is at least quadratic in its arguments
$\varphi_{1}\cdots\varphi_{m-1},h_{\text{can}\,1}\cdots h_{\text{can}\,m-1}$,
but it is linear in the higher order pieces of the canonical metric or the
coordinate shift, $h^{\mu\nu}_{\text{can}\,n-1}$ and $\varphi^{\mu}_{n-1}$
respectively, which can only enter its quadratic part, through the combination

$\displaystyle\delta\Lambda^{\mu\nu}_{\text{quad}\,n}\equiv\Box\Omega^{\mu\nu}_{1,n-1}+\Box\Omega^{\mu\nu}_{n-1,1}+\partial\Delta^{\mu\nu}_{1,n-1}+\partial\Delta^{\mu\nu}_{n-1,1}\,,$
(87)

where the first two terms on the right-hand side belong to
$\Omega^{\mu\nu}_{n}$, and the next two ones belong to
$\partial\Delta^{\mu\nu}_{n}$. Explicitly, they read

$\displaystyle\Omega_{1,n-1}^{\mu\nu}=$
$\displaystyle-\partial_{\rho}\left[\varphi_{1}^{\rho}\left(h_{\text{can}\,n-1}^{\mu\nu}+\partial\varphi_{n-1}^{\mu\nu}\right)\right]+2\,\partial_{\rho}\varphi_{1}^{(\mu}\,h_{\text{can}\,n-1}^{\nu)\rho}$
$\displaystyle+\partial^{\rho}\varphi_{1}^{(\mu}\,\partial_{\rho}\varphi_{n-1}^{\nu)}+\frac{1}{2}\eta^{\mu\nu}\left[\partial_{\rho}\varphi_{1}^{\sigma}\partial_{\sigma}\varphi_{n-1}^{\rho}-\partial_{\rho}\varphi_{1}^{\rho}\partial_{\sigma}\varphi_{n-1}^{\sigma}\right]\,,$
(88a) $\displaystyle\Omega_{n-1,1}^{\mu\nu}=$
$\displaystyle-\partial_{\rho}\left[\varphi_{n-1}^{\rho}\left(h_{\text{can}\,1}^{\mu\nu}+\partial\varphi_{1}^{\mu\nu}\right)\right]+2\,\partial_{\rho}\varphi_{n-1}^{(\mu}\,h_{\text{can}\,1}^{\nu)\rho}$
$\displaystyle+\partial^{\rho}\varphi_{n-1}^{(\mu}\,\partial_{\rho}\varphi_{1}^{\nu)}+\frac{1}{2}\eta^{\mu\nu}\left[\partial_{\rho}\varphi_{n-1}^{\sigma}\partial_{\sigma}\varphi_{1}^{\rho}-\partial_{\rho}\varphi_{n-1}^{\rho}\partial_{\sigma}\varphi_{1}^{\sigma}\right]\,,$
(88b) $\displaystyle\Delta_{1,n-1}^{\mu}=$ $\displaystyle-
h_{\text{can}\,n-1}^{\rho\sigma}\,\partial_{\rho\sigma}\varphi_{1}^{\mu}\,,$
(88c) $\displaystyle\Delta_{n-1,1}^{\mu}=$ $\displaystyle-
h_{\text{can}\,1}^{\rho\sigma}\,\partial_{\rho\sigma}\varphi_{n-1}^{\mu}+\partial_{\rho}\left(\varphi_{1}^{\rho}\Delta_{n-1}^{\mu}\right)\,.$
(88d)

These relations are straightforward generalizations of Eqs. (58)–(59).

The action of the operator $\mathrm{FP}\,\Box^{-1}_{\mathrm{ret}}$ to
$\delta\Lambda^{\mu\nu}_{\text{quad}\,n}$ generates difficult integrals
characterized by the presence of $h^{\mu\nu}_{\text{can}\,n-1}$ or
$\varphi^{\mu}_{n-1}$ in the source. The aim of the method presented here is
to bypass their evaluation. As part of our recursive hypothesis, we assume
that the easier integration of source terms involving for instance
$\Omega^{\mu\nu}_{2,n-2}$, $\Omega^{\mu\nu}_{3,n-3}$ etc., has been solved in
previous recurrence steps.

Noticing that the partial source term
$\delta\Lambda^{\mu\nu}_{\text{quad}\,n}$ is divergenceless, we can apply to
it the same treatment as for the right-hand side of Eq. (35), described in
Sec. III.2, with the substitutions $h^{\mu\nu}_{\text{gen}\,n}\to\delta
h^{\mu\nu}_{\text{quad}\,n}$, $\Lambda^{\mu\nu}_{\text{can}\,n}\to 0$,
$\partial\Delta^{\mu\nu}_{n}+\Box\Omega^{\mu\nu}_{n}\to\delta\Lambda^{\mu\nu}_{\text{quad}\,n}$.
This yields

$\displaystyle\delta
h^{\mu\nu}_{\text{quad}\,n}=\partial\delta\phi^{\mu\nu}_{\text{quad}\,n}+\Omega^{\mu\nu}_{1,n-1}+\Omega^{\mu\nu}_{n-1,1}+\delta
U^{\mu\nu}_{\text{quad}\,n}+\delta V_{\text{quad}\,n}^{\mu\nu}\,,$ (89)

with the notations $\delta U^{\mu\nu}_{\text{quad}\,n}=\delta
X^{\mu\nu}_{\text{quad}\,n}+\delta Y^{\mu\nu}_{\text{quad}\,n}$ and

$\displaystyle\delta\phi_{\text{quad}\,n}^{\mu}$
$\displaystyle=\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}\left(\frac{r}{r_{0}}\right)^{B}(\Delta_{1,n-1}^{\mu}+\Delta_{n-1,1}^{\mu})\biggr{]}\,,$
(90a) $\displaystyle\delta X^{\mu\nu}_{\text{quad}\,n}$
$\displaystyle=\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}B\left(\frac{r}{r_{0}}\right)^{B}\biggl{(}-\frac{B+1}{r^{2}}\,\bigl{(}\Omega_{1,n-1}^{\mu\nu}+\Omega_{n-1,1}^{\mu\nu}\bigr{)}-\frac{2}{r}\,\partial_{r}\bigl{(}\Omega_{1,n-1}^{\mu\nu}+\Omega_{n-1,1}^{\mu\nu}\bigr{)}\biggr{)}\biggr{]}\,,$
(90b) $\displaystyle\delta Y^{\mu\nu}_{\text{quad}\,n}$
$\displaystyle=\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}B\left(\frac{r}{r_{0}}\right)^{B}\frac{n_{i}}{r}\biggl{(}-2\delta^{i(\mu}\bigl{(}\Delta_{1,n-1}^{\nu)}+\Delta_{n-1,1}^{\nu)}\bigr{)}+\eta^{\mu\nu}\bigl{(}\Delta_{1,n-1}^{i}+\Delta_{n-1,1}^{i}\bigr{)}\biggr{)}\biggr{]}\,;$
(90c)

the term $\delta V^{\mu\nu}_{\text{quad}\,n}$ is computed from the harmonicity
algorithm $\delta V_{\text{quad}\,n}^{\mu\nu}=\mathcal{V}^{\mu\nu}[\delta
W_{\text{quad}\,n}]$, with $\delta
W^{\mu}_{\text{quad}\,n}\equiv\partial_{\nu}\delta
U^{\mu\nu}_{\text{quad}\,n}$. The source of the coordinate shift
$\delta\phi^{\mu}_{\text{quad}\,n}$, which involves the $n-1$ order piece of
the metric $h^{\mu\nu}_{\text{can}\,n-1}$, is of the same undesirable type as
$\delta\Lambda^{\mu\nu}_{\text{quad}\,n}$. However, in the current procedure,
the canonical moments will be read off from the gravitational waveform, which
is not sensitive to linear-looking gauge transformations. On the other hand,
although the commutator contributions do contain terms proportional to
$h^{\mu\nu}_{\text{can}\,n-1}$ or $\varphi_{n-1}^{\mu}$, their brute force
integration is not required. Thus, it will be possible to determine
$\mathrm{M}_{L\,n}$ and $\mathrm{S}_{L\,n}$ without integrating any such term.
Once $\delta h^{\mu\nu}_{\text{quad}\,n}$ is obtained, we compute the rest of
the metric $h^{\mu\nu}_{\text{rest}\,n}\equiv
h^{\mu\nu}_{\text{gen}\,n}-\delta h^{\mu\nu}_{\text{quad}\,n}$ by solving the
equation

$\displaystyle\Box
h^{\mu\nu}_{\text{rest}\,n}=\delta\Lambda^{\mu\nu}_{\text{rest}\,n}\equiv\Lambda^{\mu\nu}_{\text{can}\,n}-\delta\Lambda^{\mu\nu}_{\text{quad}\,n}\,,$
(91)

with the help of the standard MPM algorithm. By construction, the source term
$\delta\Lambda^{\mu\nu}_{\text{rest}\,n}$ does not depend on
$h^{\mu\nu}_{\text{can}\,n-1}$ nor $\varphi_{n-1}^{\mu}$ and is thus free of
the difficult contributions we wanted to avoid. At this stage, the $n$-th
order waveform may be built from $h^{\mu\nu}_{\text{gen}\,\leqslant
n}=Gh^{\mu\nu}_{\text{gen}\,1}+\cdots+G^{n}h^{\mu\nu}_{\text{gen}\,n}$, by
taking the limit $R\to+\infty$ for constant asymptotically null time
$U=T-R/c$, which corresponds to “radiative” coordinates
$(T=t-2G\mathrm{M}/c^{3}\ln(r/r_{0}),R=r)$ [5]. In those coordinates, the
leading order contribution to the metric, $H^{\mu\nu}_{\text{gen}\,\leqslant
n}$, admits the same expression as the one in harmonic gauge, but with the
logarithms $\ln r$ effectively replaced by $\ln r_{0}$, and with $(t,r)$
replaced by their radiative counterparts $(T,R)$; see, e.g., Ref. [21]. This
yields in particular [up to terms $\mathcal{O}(R^{-2})$]

$\displaystyle H^{\mu\nu}_{\text{gen}\,n}$
$\displaystyle=\Bigl{(}h^{\mu\nu}_{\text{rest}\,n}+\partial\delta\phi^{\mu\nu}_{\text{quad}\,n}\Bigr{)}\Big{|}_{\ln
r\to\ln r_{0}}+\Omega^{\mu\nu}_{1,n-1}+\Omega^{\mu\nu}_{n-1,1}+\delta
U^{\mu\nu}_{\text{quad}\,n}+\delta V_{\text{quad}\,n}^{\mu\nu}\,,$ (92)

where $\partial\delta\phi^{\mu\nu}_{\text{quad}\,n}$ is the linear gauge
transformation associated with $\delta\phi^{\mu}_{\text{quad}\,n}$. The $n$-th
order waveform $h^{\text{TT}}_{\leqslant n\,ij}$ is then the transverse trace-
free projection of the $1/R$ term in $H^{\mu\nu}_{\text{gen}\,\leqslant n}$.
As the TT projection of linear gauge transformations vanishes at order $1/R$,
the vector $\delta\phi^{\mu}_{\text{quad}\,n}$ is actually not required. The
result for $h^{\text{rad}}_{\leqslant n\,ij}$ is a certain functional of the
source/gauge moments:

$\displaystyle h^{\text{rad}}_{\leqslant
n\,ij}=\mathcal{H}^{\text{rad}}_{\leqslant
n\,ij}[\mathrm{I}_{L},\mathrm{J}_{L},\cdots,\mathrm{Z}_{L}]\,,$ (93a) and for
the canonical moments, we must also have, at the same time: $\displaystyle
h^{\text{rad}}_{\leqslant n\,ij}=\mathcal{H}^{\text{rad}}_{\leqslant
n\,ij}[\mathrm{M}_{L},\mathrm{S}_{L},0,0,0,0]\,.$ (93b)

The expressions of $\mathrm{M}_{n\,L}$ and $\mathrm{S}_{n\,L}$ are finally
found by guess work. We assume they are sums of terms involving the
source/gauge moments, with consistent index structures and physical
dimensions, but arbitrary coefficients. Those are fixed by identifying Eq.
(93a) with the outcome that ensues from inserting our ansatz into Eq. (93b).
Of course, if we wish to iterate the process to the next order $n+1$, we will
eventually need to tackle the difficult integrals of source terms containing
$h^{\mu\nu}_{\text{can}\,n-1}$, which arise in the calculation of
$h^{\mu\nu}_{\text{can}\,n}$ and $\delta\phi^{\mu}_{\text{quad}\,n}$.
Nonetheless, we have managed to push this step to the very end.

This strategy is particularly relevant to determine the canonical moments at
cubic order for two reasons: (i) The latter task does not demand computing
$h^{\mu\nu}_{\text{can}\,3}$ nor $\delta\phi^{\mu}_{3}$, which means that all
retarded integrals we have to consider are sourced by functions of
$h^{\mu\nu}_{\text{can}\,1}$, $\varphi^{\mu}_{1}$, or their derivatives; (ii)
The corresponding integrands have the form $f(t-r/c)\hat{n}_{L}/r^{k}$, with
$k\in\mathbb{N}\setminus\\{0,1\\}$, whose finite part retarded integral are
explicitly known [33]. In practice, we build the cubic source
$\Lambda^{\mu\nu}_{3}$, subtract the “difficult” part
$\delta\Lambda^{\mu\nu}_{3}$, and apply the MPM algorithm to the rest, which
leads to $h^{\mu\nu}_{\text{rest}}$. At last, we compute the commutators
$\delta U^{\mu\nu}_{\text{quad}\,3}$ with the method developed in Sec. IV.2,
from which we can infer $\delta V^{\mu\nu}_{\text{quad}\,3}$. The cubic
waveform follows from the effective metric

$\displaystyle h^{\mu\nu}_{\text{eff}\,\leqslant
3}=Gh^{\mu\nu}_{\text{gen}\,1}+G^{2}h^{\mu\nu}_{\text{gen}\,2}+G^{3}\left[h^{\mu\nu}_{\text{rest}\,3}+\Omega^{\mu\nu}_{12}+\Omega^{\mu\nu}_{21}+\delta
U^{\mu\nu}_{3}\right]\,.$ (94)

Our final result (I) for the canonical quadrupole moment was obtained
following the general method in Secs. III.1–III.2, and has then been entirely
checked using this approach.

## References

  * Aasi  _et al._ [2015] J. Aasi _et al._ (LIGO Scientific), Advanced LIGO, Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ [2015] F. Acernese _et al._ (VIRGO), Advanced Virgo: a second-generation interferometric gravitational wave detector, Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Akutsu  _et al._ [2019] T. Akutsu _et al._ (KAGRA), KAGRA: 2.5 Generation Interferometric Gravitational Wave Detector, Nature Astron. 3, 35 (2019), arXiv:1811.08079 [gr-qc] .
  * Maggiore [2008] M. Maggiore,  _Gravitational waves: Volume 1: Theory and experiments_ , Vol. 1 (Oxford university press, 2008).
  * Blanchet [2014] L. Blanchet, Gravitational radiation from post-Newtonian sources and inspiralling compact binaries, Living Rev. Relativ. 17, 2 (2014), arXiv:1310.1528 [gr-qc] .
  * Buonanno and Sathyaprakash [2015] A. Buonanno and B. Sathyaprakash, Sources of gravitational waves: Theory and observations, in  _General Relativity and Gravitation: A Centennial Perspective_ , edited by A. Ashtekar, B. Berger, J. Isenberg, and M. MacCallum (2015) p. 513, arXiv:1410.7832 [gr-qc] .
  * Porto [2016] R. A. Porto, The effective field theorist’s approach to gravitational dynamics, Phys. Rept. 633, 1 (2016), arXiv:1601.04914 [hep-th] .
  * Wagoner and Will [1976] R. Wagoner and C. Will, Post-Newtonian gravitational radiation from orbiting point masses, Astrophys. J. 210, 764 (1976).
  * Blanchet and Schäfer [1989] L. Blanchet and G. Schäfer, Higher order gravitational radiation losses in binary systems, Mon. Not. Roy. Astron. Soc. 239, 845 (1989).
  * Blanchet  _et al._ [1995a] L. Blanchet, T. Damour, and B. R. Iyer, Gravitational waves from inspiralling compact binaries: Energy loss and wave form to second post-Newtonian order, Phys. Rev. D 51, 5360 (1995a), gr-qc/9501029 .
  * Blanchet  _et al._ [1995b] L. Blanchet, T. Damour, B. R. Iyer, C. M. Will, and A. G. Wiseman, Gravitational radiation damping of compact binary systems to second post-Newtonian order, Phys. Rev. Lett. 74, 3515 (1995b), gr-qc/9501027 .
  * Will and Wiseman [1996] C. Will and A. Wiseman, Gravitational radiation from compact binary systems: Gravitational waveforms and energy loss to second post-Newtonian order, Phys. Rev. D 54, 4813 (1996), gr-qc/9608012 .
  * Leibovich  _et al._ [2020] A. K. Leibovich, N. T. Maia, I. Z. Rothstein, and Z. Yang, Second post-Newtonian order radiative dynamics of inspiralling compact binaries in the Effective Field Theory approach, Phys. Rev. D 101, 084058 (2020), arXiv:1912.12546 [gr-qc] .
  * Blanchet  _et al._ [2002] L. Blanchet, B. R. Iyer, and B. Joguet, Gravitational waves from inspiralling compact binaries: Energy flux to third post-Newtonian order, Phys. Rev. D 65, 064005 (2002), erratum Phys. Rev. D, 71:129903(E), 2005, gr-qc/0105098 .
  * Blanchet and Iyer [2004] L. Blanchet and B. R. Iyer, Hadamard regularization of the third post-Newtonian gravitational wave generation of two point masses, Phys. Rev. D 71, 024004 (2004), gr-qc/0409094 .
  * Blanchet  _et al._ [2004] L. Blanchet, T. Damour, G. Esposito-Farèse, and B. R. Iyer, Gravitational radiation from inspiralling compact binaries completed at the third post-Newtonian order, Phys. Rev. Lett. 93, 091101 (2004), gr-qc/0406012 .
  * Blanchet  _et al._ [2005] L. Blanchet, T. Damour, G. Esposito-Farèse, and B. R. Iyer, Dimensional regularization of the third post-Newtonian gravitational wave generation of two point masses, Phys. Rev. D 71, 124004 (2005), gr-qc/0503044 .
  * Marchand  _et al._ [2020] T. Marchand, Q. Henry, F. Larrouturou, S. Marsat, G. Faye, and L. Blanchet, The mass quadrupole moment of compact binary systems at the fourth post-Newtonian order, Class. Quant. Grav. 37, 215006 (2020), arXiv:2003.13672 [gr-qc] .
  * Larrouturou  _et al._ [2022a] F. Larrouturou, Q. Henry, L. Blanchet, and G. Faye, The quadrupole moment of compact binaries to the fourth post-Newtonian order: I. non-locality in time and infra-red divergencies, Class. Quant. Grav. 39, 115007 (2022a), arXiv:2110.02240 [gr-qc] .
  * Larrouturou  _et al._ [2022b] F. Larrouturou, L. Blanchet, Q. Henry, and G. Faye, The quadrupole moment of compact binaries to the fourth post-Newtonian order: II. dimensional regularization and renormalization, Class. Quant. Grav. 39, 115008 (2022b), arXiv:2110.02243 [gr-qc] .
  * Faye  _et al._ [2015] G. Faye, L. Blanchet, and B. R. Iyer, Non-linear multipole interactions and gravitational-wave octupole modes for inspiralling compact binaries to third-and-a-half post-Newtonian order, Class. Quant. Grav. 32, 045016 (2015), arXiv:1409.3546 [gr-qc] .
  * Henry  _et al._ [2021] Q. Henry, G. Faye, and L. Blanchet, The current-type quadrupole moment and gravitational-wave mode ($\ell$, m) = (2, 1) of compact binary systems at the third post-Newtonian order, Class. Quant. Grav. 38, 185004 (2021), arXiv:2105.10876 [gr-qc] .
  * Blanchet [1996] L. Blanchet, Energy losses by gravitational radiation in inspiralling compact binaries to five halves post-Newtonian order, Phys. Rev. D 54, 1417 (1996), erratum Phys. Rev. D 71, 129904(E) (2005), gr-qc/9603048 .
  * Blanchet [1998a] L. Blanchet, On the multipole expansion of the gravitational field, Class. Quant. Grav. 15, 1971 (1998a), gr-qc/9801101 .
  * Blanchet  _et al._ [2008] L. Blanchet, G. Faye, B. R. Iyer, and S. Sinha, The third post-Newtonian gravitational wave polarisations and associated spherical harmonic modes for inspiralling compact binaries in quasi-circular orbits, Class. Quant. Grav. 25, 165003 (2008), arXiv:0802.1249 [gr-qc] .
  * Sachs [1961] R. Sachs, Gravitational waves in general relativity. vi. outgoing radiation condition, Proc. R. Soc. London, Ser. A 264, 309+ (1961).
  * Pirani [1964] F. Pirani, Introduction to gravitational radiation theory, in  _Lectures on General Relativity_ , Brandeis Summer Institute in Theoretical Physics, Vol. 1, edited by A. Trautman, F. Pirani, and H. Bondi (Prentice-Hall, Englewood Cliffs, 1964) pp. 249–373.
  * Thorne [1980] K. Thorne, Multipole expansions of gravitational radiation, Rev. Mod. Phys. 52, 299 (1980).
  * Blanchet and Damour [1986] L. Blanchet and T. Damour, Radiative gravitational fields in general relativity. i. general structure of the field outside the source, Phil. Trans. Roy. Soc. Lond. A 320, 379 (1986).
  * Faye  _et al._ [2012] G. Faye, S. Marsat, L. Blanchet, and B. R. Iyer, The third and a half post-Newtonian gravitational wave quadrupole mode for quasi-circular inspiralling compact binaries, Class. Quant. Grav. 29, 175004 (2012), arXiv:1204.1043 [gr-qc] .
  * Blanchet [1987] L. Blanchet, Radiative gravitational fields in general relativity. ii. asymptotic behaviour at future null infinity, Proc. Roy. Soc. Lond. A 409, 383 (1987).
  * Blanchet and Damour [1992] L. Blanchet and T. Damour, Hereditary effects in gravitational radiation, Phys. Rev. D 46, 4304 (1992).
  * Blanchet [1998b] L. Blanchet, Quadrupole-quadrupole gravitational waves, Class. Quant. Grav. 15, 89 (1998b), gr-qc/9710037 .
  * Sachs and Bergmann [1958] R. Sachs and P. Bergmann, Structure of particles in linearized gravitational theory, Phys. Rev. 112, 674 (1958).
  * Blanchet and Faye [2019] L. Blanchet and G. Faye, Flux-balance equations for linear momentum and center-of-mass position of self-gravitating post-Newtonian systems, arXiv:1811.08966 [gr-qc] (2019).
  * Misner  _et al._ [1973] C. Misner, K. Thorne, and J. Wheeler,  _Gravitation_ (Freeman, San Francisco, 1973).
  * Blanchet [1998c] L. Blanchet, Gravitational-wave tails of tails, Class. Quant. Grav. 15, 113 (1998c), gr-qc/9710038 .
  * Martín-García  _et al._ [2012] J. M. Martín-García, A. García-Parrado, A. Stecchina, B. Wardell, C. Pitrou, D. Brizuela, D. Yllanes, G. Faye, L. Stein, R. Portugal, and T. Bäckdahl, xAct: Efficient tensor computer algebra for Mathematica (GPL 2002–2012), http://www.xact.es/.
  * Gradshteyn and Ryzhik [1980] I. Gradshteyn and I. Ryzhik,  _Table of Integrals, Series and Products_ (Academic Press, 1980).