# The Quadrupole Moment of Compact Binaries
to the Fourth post-Newtonian Order
I. Non-Locality in Time and Infra-Red Divergencies

François Larrouturou<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  Deutsches Elektronen-Synchrotron DESY, Notkestr.
85, 22607 Hamburg, Germany  Quentin Henry<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  Max-Planck-Institute for Gravitational Physics
(Albert-Einstein-Institute), Am Mühlenberg 1, 14476 Potsdam, Germany  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  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

###### Abstract

With the aim of providing high accuracy post-Newtonian (PN) templates for the
analysis of gravitational waves generated by compact binary systems, we
complete the analytical derivation of the source type mass quadrupole moment
of compact binaries (without spins) at the fourth PN order of general
relativity. Similarly to the case of the conservative 4PN equations of motion,
we show that the quadrupole moment at that order contains a non-local (in
time) contribution, arising from the tail-transported interaction entering the
conservative part of the dynamics. Furthermore, we investigate the infra-red
(IR) divergences of the quadrupole moment. In a previous work, this moment has
been computed using a Hadamard partie finie procedure for the IR divergences,
but the knowledge of the conservative equations of motion indicates that those
divergences have to be dealt with by means of dimensional regularization. This
work thus derives the difference between the two regularization schemes, which
has to be added on top of the previous result. We show that unphysical IR
poles start to appear at the 3PN order, and we determine all of these up to
the 4PN order. In particular, the non-local tail term comes in along with a
specific pole at the 4PN order. It will be proven in a companion paper that
the poles in the source-type quadrupole are cancelled in the physical
radiative type quadrupole moment measured at future null infinity.

###### pacs:

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

##  I Introduction

The theory of gravitational waves (GW) generated by binary systems of compact
objects has been developed using perturbative methods in classical general
relativity (see [1, 2, 3, 4] for reviews). In particular the post-Newtonian
(PN) approximation is a major and widely developed technique for computing
analytically the dynamics and GW emission of compact binaries. The state of
the art on the conservative dynamics are the equations of motion of point-
particle binaries at the 4PN (fourth-post-Newtonian) order [5, 6, 7, 8, 9, 10,
11, 12, 13, 14]. Recently progresses have been made so that the equations of
motion are now determined up to two unknown parameters at 5PN order, and up to
six unknowns at 6PN order [15, 16]. Currently the field is also evolving
thanks to new methods coming from effective field theory and scattering
amplitudes, naturally combined with the classical post-Minkowskian
approximation [17].

This paper is concerned with the GW emission aspect, which is directly related
to the data analysis of GW detectors. Here, the state of the art is the 3PN
approximation in the waveform, beyond the Einstein quadrupole formula [18, 19,
20, 21, 22]. Actually, the flux and orbital phase evolution due to
gravitational radiation are known to the 3.5PN order [23, 24], as well as the
dominant mass-quadrupole mode $(\ell,m)=(2,2)$ [25, 26], the current-
quadrupole mode $(2,1)$ [27] and the mass-octupole ones $(3,3)$ and $(3,1)$
[28]. The gravitational flux and mass quadrupole were confirmed at the 2PN
order by means of effective field theory techniques [29].

Extending the GW emission up to the 4PN (and even 4.5PN [30]) order is the
target of the present program. A central part of this program is of course the
control of the mass-type quadrupole moment of the system with the 4PN
precision. This computation faces subtle issues regarding the choice and
proper use of regularization schemes, both for the ultra-violet (UV) and
infra-red (IR) divergences. Recently, a preliminary calculation of the source
mass-type quadrupole moment of compact binaries (of spinless bodies) at the
4PN order has been tackled [31]. In this calculation, the UV divergences,
appearing because of the point-like structure of the source (modelling compact
objects with negligible internal structure), were properly treated by means of
the powerful dimensional regularization. However, the IR ones were regularized
with the Hadamard partie finie (PF) regularization procedure [32].
Furthermore, in the calculation of Ref. [31], the non-local-in-time
contributions, due to retarded correlations over arbitrarily large time spans
in the dynamics of the source [33], were neglected.

In the present paper and the next one [34], we complete the derivation of the
4PN source mass-quadrupole moment of compact binaries. More precisely, our
goals are two-fold:

  1. 1.

To derive the non-local (in time) effect in the source mass quadrupole moment
which is due to the radiation modes associated with propagating tails at
infinity. This effect is the analogue of the one occurring in the conservative
equations of motion and the Lagrangian/Hamiltonian at the 4PN order [5, 10,
7];

  2. 2.

To compute all the contributions to the quadrupole moment due to the IR
divergences, implementing a dimensional regularization scheme rather than the
Hadamard PF scheme adopted in [31]. Such procedure leads to the appearance of
specific IR poles which start to arise at the 3PN order and play a crucial
role at the 4PN order.

In the follow-up paper [34], we investigate the fate of these IR poles in the
radiative-type quadrupole moment, which represents the actual observable
moment at future null infinity, up to the 4PN order. We will find that they
are exactly canceled by radiation contributions due to non-linear propagation
effects (the most important ones being coined as “tails-of-tails” and “tails-
of-memory”) and that we can therefore safely define a three-dimensional
“renormalized” mass quadrupole moment at the 4PN order, which will constitute
a basic ingredient in the construction of 4PN GW templates.

The central object investigated in this work is thus the source-type mass-
quadrupole moment, defined in $d$ spatial dimensions by the expression [35,
31, 27]111The overbar denotes the formal PN expansion; superscript parenthesis
$(n)$ denote time derivatives; the hat refers to the symmetric-trace-free
(STF) product, e.g. $\hat{x}_{ija}\equiv\text{STF}[x_{i}x_{j}x_{a}]$;
$\text{PF}_{B=0}$ denotes the Hadamard partie finie with regulator
$(r/r_{0})^{B}$ and associated length scale $r_{0}$ ($B\in\mathbb{C}$); the
characteristic dimensional regularization length scale is $\ell_{0}$.

$\displaystyle I_{ij}$
$\displaystyle=\frac{d-1}{2(d-2)}\mathop{\mathrm{PF}}_{B=0}\int\frac{\mathrm{d}^{d}\mathbf{x}}{\ell_{0}^{d-3}}\left(\frac{r}{r_{0}}\right)^{B}\biggr{\\{}\hat{x}_{ij}\,\overline{\Sigma}_{[2]}-\frac{4(d+2)}{c^{2}\,d(d+4)}\,\hat{x}_{ija}\,\overline{\Sigma}^{(1)}_{a[3]}$
$\displaystyle\qquad\qquad+\frac{2(d+2)}{c^{4}\,d(d+1)(d+6)}\,\hat{x}_{ijab}\,\overline{\Sigma}^{(2)}_{ab[4]}-\frac{4B(d-3)(d+2)}{c^{2}\,(d-1)d(d+4)}\,\hat{x}_{ija}\,\frac{x_{b}}{r^{2}}\,\overline{\Sigma}_{ab[3]}\biggr{\\}}\,.$
(1)

We let the reader refer to [31] for a comprehensive review of the definitions
and properties of the quantities entering the source quadrupole moment, as
well as its computation using the IR Hadamard and UV dimensional
regularizations. Let us just emphasize a few points.

The main quantity over which the source quadrupole integrates is the pseudo
stress-energy tensor $\tau^{\mu\nu}$ in harmonic coordinates composed of a
matter and a gravitational part,

$\tau^{\mu\nu}\equiv|g|T^{\mu\nu}+\frac{c^{4}}{16\pi G}\,\Lambda^{\mu\nu}\,,$
(2)

where $T^{\mu\nu}$ is the matter stress-energy tensor and $\Lambda^{\mu\nu}$
the non-linear gravitational source term of the Einstein field equations. It
enters Eq. (I) through the following PN-expanded source densities:

$\overline{\Sigma}=\frac{2}{d-1}\frac{(d-2)\overline{\tau}^{00}+\overline{\tau}^{ii}}{c^{2}}\,,\qquad\overline{\Sigma}_{i}=\frac{\overline{\tau}^{0i}}{c}\,,\qquad\overline{\Sigma}_{ij}=\overline{\tau}^{ij}\,.$
(3)

The PN expansions implicit in (I) read as (with $\Gamma$ the Euler function)

$\overline{\Sigma}_{[\ell]}(\mathbf{x},t)\equiv\sum_{k=0}^{+\infty}\frac{1}{2^{2k}k!}\frac{\Gamma(\ell+\frac{d}{2})}{\Gamma(\ell+\frac{d}{2}+k)}\left(\frac{r}{c}\frac{\partial}{\partial
t}\right)^{2k}\overline{\Sigma}(\mathbf{x},t)\,.$ (4)

The partie finie procedure $\text{PF}_{B=0}$ comes from the matching between
near zone and exterior zone. It is crucial for the proper definition of the
multipole moments in 3 dimensions [36]. The derivation of the equations of
motion [8, 9] indicated that the PF operator should be kept even in $d$
dimensions. It was also shown there that, in $d$ dimensions, the limit $B\to
0$ is finite in physical quantities (without pole $1/B$); so, in the end the
whole procedure is equivalent to the usual dimensional regularization. We had
however to keep the PF operator as explicit in the expression of the
quadrupole moment (I). Here, similarly to the equations of motion [8, 9], we
apply first the PF process when $B\to 0$ on the $d$ dimensional expression
and, second, the usual dimensional regularization when $\varepsilon\equiv
d-3\to 0$. We call this mixed regularization the “$B\varepsilon$”
regularization. Note that the last term in (I), proportional to both $B$ and
$d-3$, will be shown to play no role with this particular regularization.

This paper is organized as follows. In the next Section II, we derive the non-
local in time part of the source quadrupole, relying on results derived in
[8]. We then perform the proper IR dimensional regularization for all the
various categories of terms composing the quadrupole in Section III. The final
result is presented in Section IV in the form of a pole followed by a finite
part contribution, which will be the starting point for the final
renormalization of the quadrupole in the next paper [34]. Appendix A present
technical formulas generalizing the Riesz formula in $d$ dimensions. Appendix
B gives the expression of the local part of the IR shift coming from the 4PN
equations of motion.

##  II Non-locality in time of the source quadrupole moment

Crucial to the completion of the ambiguity-free equations of motion at the 4PN
order was the proper inclusion of the tail effect, i.e. the non-local in time
back-scattering of emitted gravitational waves, modifying the conservative
dynamics of the system at the current time [8]. This effect enters at the 4PN
order in the near-zone metric and, as such, plays a key role in the
computation of the 4PN source mass quadrupole by introducing a non-local term,
together with a pole.

###  II.1 The tail effect in the conservative 4PN equations of motion

Let us first recall how the tail effect is included in the near-zone metric
and Fokker action for the conservative dynamics. In Ref. [8], it was found
that the PN-expanded gravitational field in harmonic coordinates
($h^{\mu\nu}\equiv\sqrt{-g}g^{\mu\nu}-\eta^{\mu\nu}$, which is such that
$\partial_{\nu}h^{\mu\nu}=0$) contains the following pieces responsible for
tails at the 4PN order:

$\displaystyle\overline{h}^{00ii}_{\text{tail}}$
$\displaystyle=\frac{8G^{2}M}{15c^{10}}\,x^{ij}\int_{0}^{+\infty}\mathrm{d}\tau\left[L_{\varepsilon}(\tau)+\frac{61}{60}\right]I^{(7)}_{ij}(t-\tau)+\mathcal{O}\left(\frac{1}{c^{12}}\right)\,,$
(5a) $\displaystyle\overline{h}^{0i}_{\text{tail}}$
$\displaystyle=-\frac{8G^{2}M}{3c^{9}}\,x^{j}\int_{0}^{+\infty}\mathrm{d}\tau\left[L_{\varepsilon}(\tau)+\frac{107}{120}\right]I^{(6)}_{ij}(t-\tau)+\mathcal{O}\left(\frac{1}{c^{11}}\right)\,,$
(5b) $\displaystyle\overline{h}^{ij}_{\text{tail}}$
$\displaystyle=\frac{8G^{2}M}{c^{8}}\int_{0}^{+\infty}\mathrm{d}\tau\left[L_{\varepsilon}(\tau)+\frac{4}{5}\right]I^{(5)}_{ij}(t-\tau)+\mathcal{O}\left(\frac{1}{c^{10}}\right)\,,$
(5c)

where we have posed $h^{00ii}\equiv\frac{2}{d-1}[(d-2)h^{00}+h^{ii}]$ and
introduced the short-hand notation for the logarithmic divergence with the
associated pole $\propto\varepsilon^{-1}\equiv(d-3)^{-1}$:

$L_{\varepsilon}(\tau)\equiv\ln\left(\frac{c\sqrt{\bar{q}}\,\tau}{2\ell_{0}}\right)-\frac{1}{2\varepsilon}\,,$
(6)

where $\bar{q}\equiv 4\pi\,\mathrm{e}^{\gamma_{\text{E}}}$ with
$\gamma_{\text{E}}$ the Euler constant, vanishingly small terms
$\mathcal{O}(\varepsilon)$ being neglected.

In our set-up and notation, the two and only two series of multipole moments
describing the radiation field generated by an isolated source are the so-
called canonical moments, $M_{L}$ and $S_{L}$. Those differ from the source
moments $I_{L}$ and $J_{L}$, e.g. (I) and, thus, in principle, the quadrupole
moment in (5) should rather be viewed as the canonical moment and denoted
$M_{ij}$. However, we take here advantage that $M_{ij}$ and the source moment
$I_{ij}$ are equivalent at Newtonian order (and even up to 2PN order).

Still in Ref. [8], we next applied a gauge transformation, at quadratic order,
for the particular interaction $M\times I_{ij}$, so designed as to transfer
all relevant tail terms in the “$00ii$” component of the metric. Namely, we
posed
${h^{\prime}}^{\mu\nu}=h^{\mu\nu}+\partial^{\mu}\epsilon_{\text{tail}}^{\nu}+\partial^{\nu}\epsilon_{\text{tail}}^{\mu}-\eta^{\mu\nu}\partial_{\rho}\epsilon_{\text{tail}}^{\rho}$,
with PN-expanded gauge vector given by the following tail pieces (used in [8]
but published here for the first time):

$\displaystyle\overline{\epsilon}^{0}_{\text{tail}}$
$\displaystyle=-\frac{2G^{2}M}{3c^{9}}\,x^{ij}\int_{0}^{+\infty}\mathrm{d}\tau\left[L_{\varepsilon}(\tau)+\frac{37}{60}\right]I^{(6)}_{ij}(t-\tau)+\mathcal{O}\left(\frac{1}{c^{11}}\right)\,,$
(7a) $\displaystyle\overline{\epsilon}^{i}_{\text{tail}}$
$\displaystyle=-\frac{4G^{2}M}{c^{8}}\,x^{j}\int_{0}^{+\infty}\mathrm{d}\tau\left[L_{\varepsilon}(\tau)+\frac{4}{5}\right]I^{(5)}_{ij}(t-\tau)+\mathcal{O}\left(\frac{1}{c^{10}}\right)\,.$
(7b)

In the new gauge, the 4PN tail effect is thus entirely described by the single
scalar potential ${h^{\prime}}^{00ii}_{\text{tail}}$ (or, equivalently, by the
$00$ component of the covariant metric ${g^{\prime}}_{00}^{\text{tail}}$),
which becomes

$\displaystyle{\overline{h}^{\prime}}^{00ii}_{\text{tail}}$
$\displaystyle=\frac{16G^{2}M}{5c^{10}}\,x^{ij}\int_{0}^{+\infty}\mathrm{d}\tau\left[L_{\varepsilon}(\tau)+\frac{41}{60}\right]I^{(7)}_{ij}(t-\tau)+\mathcal{O}\left(\frac{1}{c^{12}}\right)\,.$
(8)

This tail piece in the metric yields the tail term in the conservative Fokker
action, which is found to be manifestly symmetric under time reversal,

$\displaystyle S_{\text{tail}}$
$\displaystyle=\frac{G^{2}M}{5c^{8}}\int_{-\infty}^{+\infty}\mathrm{d}t\,I_{ij}^{(3)}(t)\int_{0}^{+\infty}\mathrm{d}\tau\biggl{[}L_{\varepsilon}(\tau)+\frac{41}{60}\biggr{]}\left(I_{ij}^{(4)}(t-\tau)-I_{ij}^{(4)}(t+\tau)\right)\,.$
(9)

An elegant alternative form is provided by the Hadamard partie finie (Pf)
integral

$\displaystyle
S_{\text{tail}}=\frac{G^{2}M}{5c^{8}}\,\mathop{\text{Pf}}_{\tau_{0}}\int\\!\\!\\!\int\frac{\mathrm{d}t\mathrm{d}t^{\prime}}{|t-t^{\prime}|}I_{ij}^{(3)}(t)\,I_{ij}^{(3)}(t^{\prime})\,,$
(10)

with the Hadamard regularization scale
$\tau_{0}=\frac{2\ell_{0}}{c\sqrt{\bar{q}}}\,\text{exp}[\frac{1}{2\varepsilon}-\frac{41}{60}]$.

###  II.2 Direct tail term in the 4PN mass quadrupole moment

Inserting Eqs. (5) into the quadrupole moment (I), it is straightforward to
see that the only 4PN effect comes from the first term, namely

$\displaystyle
I_{ij}^{\text{tail}}=\frac{d-1}{2(d-2)}\mathop{\mathrm{PF}}_{B=0}\int\mathrm{d}^{d}\mathbf{x}\,\left(\frac{r}{r_{0}}\right)^{B}\,\hat{x}_{ij}\,\overline{\Sigma}_{\text{tail}}+\cdots\,,$
(11)

where only the Newtonian term in the PN series (4) needs to be considered, and
where the ellipsis denote other terms that do not participate to the effect at
4PN order. Considering the order of appearance of the tail integrals in the
metric, the only terms in the non-linear source $\Lambda^{\mu\nu}$ which may
contribute to the effect at the 4PN order are
$\Lambda^{00}=-h^{ab}\partial_{ab}h^{00}+\frac{1}{4}\partial_{a}h^{00}\partial_{a}h^{bb}$
and $\Lambda^{ii}=-\frac{1}{4}(d-2)\partial_{a}h^{00}\partial_{a}h^{bb}$. Now,
only the components of the metric that are proportional to $1/c^{8}$, namely
$h_{\text{tail}}^{ab}$, can subsist at this accuracy level. They are directly
given by (5c) and are obviously trace-free. The tail sector of the effective
non-linear source in (11) thus reduces to

$\overline{\Sigma}_{\text{tail}}=-\frac{c^{2}}{16\pi
G}\frac{2(d-2)}{d-1}\,h_{\text{tail}}^{ab}\partial_{ab}h^{00}+\cdots\,.$ (12)

On the other hand, to the lowest order, we have
$h^{00}=-\frac{2(d-1)}{d-2}\frac{V}{c^{2}}+\cdots$, with $V$ being, at leading
order, the Newtonian potential in $d$ dimensions. Using also the fact that
$h_{\text{tail}}^{ab}$ is a function of time only, we find that the 4PN tail
contribution in the quadrupole moment reads

$\displaystyle I_{ij}^{\text{tail}}$ $\displaystyle=\frac{d-1}{8\pi
G(d-2)}\,h^{ab}_{\text{tail}}\mathop{\mathrm{PF}}_{B=0}\int\mathrm{d}^{d}\mathbf{x}\left(\frac{r}{r_{0}}\right)^{B}\hat{x}^{ij}\partial_{ab}V\,.$
(13)

Injecting the Newtonian potential
$V=\tilde{k}\,G\,m_{1}\frac{2(d-2)}{d-1}\,r_{1}^{2-d}+1\leftrightarrow 2$,
with $\tilde{k}=\pi^{1-\frac{d}{2}}\Gamma(\frac{d}{2}-1)$, the term (13) can
be computed by means of the “generalized Riesz integrals” presented in
Appendix A, more precisely Eq. (64a) with $\bm{y}_{2}=0$. Finally, applying
the operator PF at $B=0$ which reduces here to a simple limit when $B\to
0$,222The calculation boils down to the single simple elementary integral
$\frac{\tilde{k}}{2\pi}\mathop{\mathrm{PF}}_{B=0}\int\mathrm{d}^{d}\mathbf{x}\left(\frac{r}{r_{0}}\right)^{B}\hat{x}^{ij}\,r_{1}^{2-d}=-\frac{1}{d+4}\,\bm{y}_{1}^{2}\,y_{1}^{\langle
ij\rangle}\,.$ we find at 4PN order

$\displaystyle
I_{ij}^{\text{tail}}=-\frac{8G^{2}M}{c^{8}}\,\int_{0}^{+\infty}\mathrm{d}\tau\left[L_{\varepsilon}(\tau)+\frac{4}{5}\right]\left(\frac{4}{d+4}I_{k\langle
i}(t)I^{(5)}_{j\rangle
k}(t-\tau)+\frac{I(t)}{d}I^{(5)}_{ij}(t-\tau)\right)\,.$ (14)

Here, the STF quadrupole moment reads $I_{ij}=m_{1}\,y_{1}^{\langle
ij\rangle}+1\leftrightarrow 2+\mathcal{O}(c^{-2})$; we have also introduced
the Newtonian moment of inertia $I=m_{1}\,\bm{y}_{1}^{2}+1\leftrightarrow
2+\mathcal{O}(c^{-2})$.

###  II.3 Indirect tail contribution due to a 4PN shift and total tail effect

Now, an important point for our purpose is to remark that the gauge
transformation that has been applied to the 4PN equations of motion, with
gauge vector (7), will induce a spatial shift of the particles’ world-lines,
which has to be taken into account when evaluating the quadrupole moment.
Applying such a shift is crucial to ensure the coherence between the
coordinate systems used in the derivation of the equations of motion and the
multipole moments. Denoting the value on particle 1 of the gauge vector as
$\epsilon_{\text{tail\,1}}^{\mu}\equiv\epsilon_{\text{tail}}^{\mu}[t,\bm{y}_{1}(t)]$,
the spatial shift
$\bm{\zeta}_{1}(t)\equiv\bm{y}_{1}(t)-\bm{y}^{\prime}_{1}(t)$ will be given in
general by
$\zeta_{\text{1}}^{i}=-\epsilon^{i}_{\text{tail\,1}}+v_{1}^{i}\,\epsilon^{0}_{\text{tail\,1}}/c+\mathcal{O}(\epsilon^{2})$.
It is clear from the $1/c$ factors of Eq. (7) that only the first term
contributes to 4PN order; hence the shift reads

$\displaystyle\overline{\zeta}_{1}^{i}=\frac{4G^{2}M}{c^{8}}\,y_{1}^{j}\int_{0}^{+\infty}\mathrm{d}\tau\left[L_{\varepsilon}(\tau)+\frac{4}{5}\right]I^{(5)}_{ij}(t-\tau)+\mathcal{O}\left(\frac{1}{c^{10}}\right)\,,$
(15)

and the contribution brought about by this shift in the quadrupole moment at
4PN is

$\displaystyle\delta_{\zeta}I_{ij}\equiv 2m_{1}y_{1}^{\langle
i}\overline{\zeta}^{\,j\rangle}_{1}+1\leftrightarrow
2+\mathcal{O}\left(\frac{1}{c^{10}}\right)\,.$ (16)

Therefore, again denoting by $I=m_{1}\,\bm{y}_{1}^{2}+1\leftrightarrow
2+\mathcal{O}(c^{-2})$ the Newtonian moment of inertia, we obtain at 4PN order

$\displaystyle\delta_{\zeta}I_{ij}=\frac{8G^{2}M}{c^{8}}\,\int_{0}^{+\infty}\mathrm{d}\tau\left[L_{\varepsilon}(\tau)+\frac{4}{5}\right]\left(I_{k\langle
i}(t)I^{(5)}_{j\rangle
k}(t-\tau)+\frac{I(t)}{d}I^{(5)}_{ij}(t-\tau)\right)\,.$ (17)

As the computation carried out in [31] was purely local, none of the non-local
effects (neither the direct effect of the tails nor the contribution of the
non-local shift) were included in this preliminary result. Combining the two
new contributions of this work (14) and (17), we thus have to add a non-local
in time term to the mass quadrupole moment previously computed in [31] which
reads explicitly at 4PN [neglecting terms $\mathcal{O}(\varepsilon)$]

$\displaystyle I_{ij}^{\text{non-loc}}\equiv
I_{ij}^{\text{tail}}+\delta_{\zeta}I_{ij}=\frac{24G^{2}M}{7c^{8}}\,I_{k\langle
i}(t)\int_{0}^{+\infty}\mathrm{d}\tau\left[L_{\varepsilon}(\tau)+\frac{74}{105}\right]I^{(5)}_{j\rangle
k}(t-\tau)\,.$ (18)

An interesting point is that the contribution of the moment of inertia $I(t)$
cancels out, but there remains a pole. The next section and the following
paper [34] will show how this pole finally combines with other poles coming
from IR dimensional regularization to disappear from the observable radiative
moment.

Note finally that $I_{ij}^{\text{non-loc}}$ contains both conservative and
dissipative effects, together with a purely instantaneous (non-tail) piece.
Those will be specified and computed on quasi-circular orbits in the companion
paper [34].

##  III IR regularization of the source quadrupole moment

The quadrupole moment in $d$ dimensions given by (I), when PN-expanded using
Eqs. (4), and after some suitable integrations by part, can be decomposed into
four different types of contributions (see App. C of [31] for the exhaustive
list of all those terms):

  1. 1.

Volume terms, where the integrands are made of products of derivatives of
elementary potentials. The complete list of potentials is provided in Appendix
A of [31];333The derivatives of those products of potentials are treated in
the context of the theory of distributions. They can thus include a
distributional sector, made with Dirac distributions. Those are naturally
regarded as compact terms.

  2. 2.

Compact-support terms, where the integrands are proportional to Dirac
distributions modeling the compact objects, multiplying products of
derivatives of elementary potentials;

  3. 3.

Surface terms, where the integrands are total spatial derivatives and can be
replaced by their expansions at spatial infinity;

  4. 4.

The “extra” term, which is the last term of (I), proportional to $B(d-3)$, and
formally zero in Hadamard’s sense, in 3 dimensions within a UV dimensional
regularization scheme.

Since the calculations concern only the IR bound, we only need the appropriate
accurate values for the potentials when expanded at spatial infinity.
Interestingly, the integral giving the value of a potential at the location of
particles extends up to infinity, thus this applies also to the case of
compact terms, which may be affected by the change of regularization [see Eq.
(40) below].

In Ref. [31], the applied IR regularization scheme was a pure Partie Finie
(PF) one as $B\rightarrow 0$ (with poles $1/B$ discarded) on the three-
dimensional expression of (I). Instead, we resort here to the mixed
“$B\varepsilon$” regularization which consists of computing first the limit
$B\to 0$ on the $d$-dimensional expression of the quadrupole (I), and only
then apply the dimensional regularization when $\varepsilon\to 0$, of course
keeping track of all the poles $1/\varepsilon$. By contrast with the pure PF
regularization, we expect from the work [8, 9] on equations of motion that the
first limit $B\to 0$, i.e., performed on the top of dimensional
regularization, will be finite in that case (no pole $1/B$), which is
confirmed by our calculations below.

Each different type of terms is to be dimensionally regularized (following the
$B\varepsilon$ scheme) with specific techniques, which are exposed in the rest
of the section. In a nutshell, we aim at computing

$\mathcal{D}I_{ij}\equiv
I_{ij}^{B\varepsilon}-I_{ij}^{\text{Had}}=\mathcal{D}I_{ij}^{\text{Vol}}+\mathcal{D}I_{ij}^{\text{Comp}}+\mathcal{D}I_{ij}^{\text{Surf}}+\mathcal{D}I_{ij}^{\text{extra}}\,,$
(19)

where $I_{ij}^{B\varepsilon}$ represents the source mass quadrupole properly
regularized following the $B\varepsilon$ prescription, $I_{ij}^{\text{Had}}$
is an abuse of notation for the result of [31], i.e. the source quadrupole
computed with dimensional regularization for the UV divergences and the
Hadamard regularization for the IR ones. The terms in the right-hand-side
(RHS) represent the differences for the four types of terms 1 to 4.

###  III.1 Volume terms

The most numerous terms to regularize are the volume terms. Their
regularizations consist in two parts: (i) one has to compute the general
expression of the difference between the two regularization schemes
$B\varepsilon$ and Had for each term, and then (ii) inject the accurate
$d$-dimensional values for the potentials expanded at spatial infinity.

####  III.1.1 Dimensional regularization of volume terms

Let us consider a generic volume term

$\mathcal{V}=\int\\!\frac{\mathrm{d}^{d}\mathbf{x}}{\ell_{0}^{\varepsilon}}\left(\frac{r}{r_{0}}\right)^{B}F(\mathbf{x},t)\,.$
(20)

The function $F(\mathbf{x},t)$ is some product of (derivatives of) potentials
and some $\hat{x}_{L}$; below, we will drop the time dependence as it plays no
role in the regularization procedure. As we investigate the difference between
IR regularization schemes, we can restrict the integral to $r>\mathcal{R}$,
where $\mathcal{R}$ is an arbitrary constant scale, significantly larger than
the distances of the particles to the origin, $|\bm{y}_{A}|$, so that we do
not have to consider the problem coming from the point-particle approximation,
already dealt with in [31].

In Ref. [31], we have computed the pure PF regularization

$\mathcal{V}^{\text{Had}}=\mathop{\mathrm{PF}}_{B=0}\int_{r>\mathcal{R}}\\!\\!\mathrm{d}^{3}\mathbf{x}\left(\frac{r}{r_{0}}\right)^{B}F^{(d=3)}(\mathbf{x})\,,$
(21)

where $F^{(d=3)}$ denotes naturally the three-dimensional limit of $F$
obtained by performing the PN iteration in 3 dimensions and given by (27). On
the other hand, in this work we consider the mixed regularization scheme
$B\varepsilon$,

$\mathcal{V}^{B\varepsilon}=\mathop{\mathrm{PF}}_{B=0}\int_{r>\mathcal{R}}\\!\\!\frac{\mathrm{d}^{d}\mathbf{x}}{\ell_{0}^{\varepsilon}}\,\left(\frac{r}{r_{0}}\right)^{B}F(\mathbf{x})\,.$
(22)

Thus, for each volume term, we are to compute the difference of regularization
schemes

$\mathcal{D}\mathcal{V}\equiv\mathcal{V}^{B\varepsilon}-\mathcal{V}^{\text{Had}}\,.$
(23)

Since, in the limit $\varepsilon\to 0$, the complementary integrals over
$r<\mathcal{R}$ agree with each other, this difference should be independent
of the cut-off scale $\mathcal{R}$. We conclude that $\mathcal{D}\mathcal{V}$
is the proper quantity that we have to add to the volume terms computed in
[31].

As it appears at the 4PN order, the functions $F$ we consider admit generic
far-zone (or multipolar) expansions in $d$ dimensions when
$r\equiv|\mathbf{x}|\rightarrow+\infty$, namely444By which we really mean
that, for any $N\in\mathbb{N}$, we can write
$F(\mathbf{x})=\sum_{p=-p_{0}}^{N}\sum_{q=q_{0}}^{q_{1}}\frac{\ell_{0}^{q\varepsilon}}{r^{p+q\varepsilon}}\,\varphi^{(\varepsilon)}_{p,q}(\mathbf{n})+o\left(\frac{1}{r^{N}}\right)\,,$
where $p_{0}\in\mathbb{N}$ indicates the maximal order of the IR divergence,
and $q_{0}$, $q_{1}\in\mathbb{Z}$ represent a finite range of values for $q$
depending implicitly on $N$.

$F(\mathbf{x})=\sum_{p,q}\frac{\ell_{0}^{q\varepsilon}}{r^{p+q\varepsilon}}\,\varphi^{(\varepsilon)}_{p,q}(\mathbf{n})\,,$
(24)

with coefficients that can contain poles $\propto 1/\varepsilon$, i.e. of the
type

$\varphi^{(\varepsilon)}_{p,q}(\mathbf{n})=\frac{1}{\varepsilon}\,\psi^{(-1)}_{p,q}(\mathbf{n})+\psi^{(\varepsilon)}_{p,q}(\mathbf{n})\,.$
(25)

We have verified that no double poles $\propto 1/\varepsilon^{2}$ appear at
the 4PN order. The coefficients $\psi^{(-1)}_{p,q}$ of the pole are naturally
defined with no dependence upon $\varepsilon$, while the
$\psi^{(\varepsilon)}_{p,q}$ are finite when $\varepsilon\to 0$. An important
point is that, despite the poles, the 3 dimensional limits of the functions
$F$ are finite, as clear from their expressions given in the App. C of [31]
and explicitly verified in our computation. This means that (for all $p$)

$\sum_{q}\psi^{(-1)}_{p,q}(\mathbf{n})=0\,.$ (26a) Furthermore, by posing
$\sum_{q}\psi^{(\varepsilon=0)}_{p,q}(\mathbf{n})\equiv\varphi_{p}(\mathbf{n})\,,\qquad\sum_{q}q\,\psi^{(-1)}_{p,q}(\mathbf{n})\equiv-\varphi^{\ln}_{p}(\mathbf{n})\,,$
(26b)

we get the corresponding logarithmic expansion in 3 dimensions

$F^{(d=3)}=\sum_{p}\frac{1}{r^{p}}\left[\varphi_{p}(\mathbf{n})+\varphi^{\ln}_{p}(\mathbf{n})\,\ln\left(\frac{r}{\ell_{0}}\right)\right]\,.$
(27)

Note that, as a confirmation of the absence of double poles in the
$d$-dimensional expressions, no squared logarithms appeared in the Hadamard
computation.

Using the relations (26) linking the $d$ and $3$ dimensional quantities, we
obtain

$\mathcal{D}\mathcal{V}=\sum_{q\neq
1}\frac{1}{\varepsilon}\left[\frac{1}{q-1}-\varepsilon\ln\left(\frac{r_{0}}{\ell_{0}}\right)\right]\\!\int\\!\mathrm{d}\Omega_{d-1}\,\varphi^{(\varepsilon)}_{3,q}(\mathbf{n})-\frac{1}{2}\ln^{2}\left(\frac{r_{0}}{\ell_{0}}\right)\int\\!\mathrm{d}\Omega_{2}\,\varphi^{\ln}_{3}(\mathbf{n})+\mathcal{O}\left(\varepsilon\right)\,,$
(28)

where $\mathrm{d}\Omega_{d-1}$ denotes the $(d-1)$-dimensional solid angle
element. As expected, we find that, modulo $\mathcal{O}(\varepsilon)$-terms,
this difference does not depend on the cut-off scale $\mathcal{R}$. The result
(28) generalizes that of Eq. (2.11) in [8], to include integrands containing
poles. Note that the value $q=1$ has to be excluded from the sum in (28). This
is _not_ an artificial requirement to ensure that the formula is well-defined:
it comes naturally out of the derivation. A last use of (26) yields the more
compact equivalent form

$\mathcal{D}\mathcal{V}=\frac{1}{\varepsilon}\sum_{q\neq
1}\frac{1}{q-1}\left(\frac{\ell_{0}}{r_{0}}\right)^{(q-1)\varepsilon}\\!\int\\!\mathrm{d}\Omega_{d-1}\,\varphi^{(\varepsilon)}_{3,q}(\mathbf{n})+\mathcal{O}\left(\varepsilon\right)\,.$
(29)

####  III.1.2 Far zone expansion of the potentials in $d$ dimensions

From (29), we see that the difference in regularization schemes depends on
specific orders in the far-zone expansion of the integrand of the volume
terms. As those are made of products of (derivatives of) potentials (defined
in App. A of [31]), we need the expansion of individual potentials in $d$
dimensions.

The expressions of the simplest compact support potentials $V$, $V_{i}$, $K$
and the super-potentials (defined in Sec. III.B of [31]) in the whole
$d$-dimensional space are already available, so that we have just to expand
them when $r\rightarrow+\infty$. However, the non-compact support potential
$\hat{W}_{ij}$ is also required at 1PN order, and $\hat{X}$, $\hat{R}_{i}$ and
$\hat{Z}_{ij}$ at Newtonian order, none of those potentials being known in all
$d$-dimensional space. We thus have to compute their asymptotic behaviours by
iterating the propagator at infinity and, crucially, add an appropriate
homogeneous solution.

The far-zone expansion (or multipolar expansion, indicated by the operator
$\mathcal{M}$) of a potential $P$ with source $S$ (assumed to be PN expanded),
i.e., such that $\Box P=S=\overline{S}$ in $d$ dimensions, reads [8]

$\mathcal{M}(P)=\mathop{\mathrm{PF}}_{B=0}\,\Box^{-1}\left[\left(\frac{r}{r_{0}}\right)^{B}\mathcal{M}(S)\right]-\frac{1}{4\pi}\sum_{\ell\in\mathbb{N}}\frac{(-)^{\ell}}{\ell!}\partial_{L}S_{\star}^{L}\,.$
(30)

In principle, the first term is built from the PN-expanded retarded propagator
$\Box_{\text{R}}^{-1}$, but we trade it here for the symmetric one
$\Box^{-1}=\Box_{\text{S}}^{-1}$, since we are merely interested in even
orders. The second term of the RHS is a homogeneous solution constructed out
of the PN expansion of

$S_{\star}^{L}(t,r)\equiv\frac{\tilde{k}}{2r^{d-2}}\int_{1}^{+\infty}\\!\mathrm{d}z\,\gamma_{\frac{1-d}{2}}(z)\Bigl{[}\mathcal{S}_{L}\left(t-zr/c\right)+\mathcal{S}_{L}\left(t+zr/c\right)\Bigr{]}\,,$
(31)

where we recall that
$\tilde{k}\equiv\pi^{1-\frac{d}{2}}\Gamma(\frac{d}{2}-1)$, and that
$\gamma_{\frac{1-d}{2}}(z)$ is the kernel function entering the Green function
of the d’Alembertian equation in $d$ dimensions, defined by Eq. (3.3) in [8]
and recalled in Eq. (2.3) of the follow-up paper [34]. In other words,
$S_{\star}^{L}(t,r)$ is an elementary monopolar homogeneous solution of the
d’Alembertian equation, parametrized by the moments

$\mathcal{S}_{L}(u)\equiv\mathop{\mathrm{PF}}_{B=0}\int\\!\mathrm{d}^{d}\mathbf{x}\left(\frac{r}{r_{0}}\right)^{B}x_{L}\,S(\mathbf{x},u)\,.$
(32)

Note that these particular moments are chosen to be non-STF, with just
$x_{L}\equiv x_{i_{1}}\cdots x_{i_{\ell}}$. Performing explicitly the PN
expansion, Eq. (30) becomes

$\displaystyle\mathcal{M}(P)$
$\displaystyle=\sum_{k\in\mathbb{N}}\left(\frac{1}{c}\frac{\partial}{\partial
t}\right)^{2k}\mathop{\mathrm{PF}}_{B=0}\Delta^{-k-1}\left[\left(\frac{r}{r_{0}}\right)^{B}\mathcal{M}(S)\right]$
$\displaystyle\qquad-\frac{1}{4\pi^{\frac{d-1}{2}}}\sum_{k,\ell}\frac{(-)^{\ell}}{\ell!(2k)!}\frac{\Gamma\left(\frac{d}{2}-1-k\right)}{\Gamma\left(\frac{1}{2}-k\right)}\frac{1}{c^{2k}}\,\mathcal{S}_{L}^{(2k)}(t)\,\partial_{L}r^{2k+2-d}\,,$
(33)

plus odd powers of $1/c$ which can be ignored for the present purpose.

Computing each of those terms requires different techniques:

  * •

Particular solutions. There are no specific issues with the multipolar
expansion of the sources; so, computing the first term in (III.1.2) is
straightforwardly done by means of iterations of the $d$-dimensional
“Matthieu” formula [37]

$\Delta^{-1}\Bigl{[}r^{\alpha}\,\hat{n}_{L}\Bigr{]}=\frac{r^{\alpha+2}\,\hat{n}_{L}}{(\alpha+d+\ell)(\alpha+2-\ell)}\,.$
(34)

  * •

Homogeneous solutions. The delicate part in the computation of expanded
$d$-dimensional potentials is the control of the homogeneous solutions.
Indeed, the source integral $\mathcal{S}_{L}(u)$ of (32), of non-compact
support, is composed of terms like

$\int\\!\mathrm{d}^{d}\mathbf{x}\,x^{L}\,r^{B}\,r_{1}^{\alpha}\,r_{2}^{\beta}\,n_{1}^{P}n_{2}^{Q}\,.$
(35)

Those can be computed either by means of the generalized Riesz formulae
described in Appendix A below, or by implementing the nice method of App. A in
[38], relying on the use of prolate spheroidal coordinates. In addition to the
integrals (35), we had to deal with the cubic sector of the potential
$\hat{X}$, which is defined by Eq. (A.4f) in [31] and whose source contains
the delicate term

$W_{ij}\,\partial_{ij}V\,.$ (36)

This cubic part requires a priori the knowledge of the potential
$\hat{W}_{ij}$ all over the space in $d$ dimensions at Newtonian order, but
this has not been calculated yet. We could have done it relying on the
generalization of the Fock function $g=\ln(r_{1}+r_{2}+r_{12})$ in $d$
dimensions, which has been derived in [37]. However, its expression is only
given in an integral form, which is not very convenient in practice. Instead,
we employed the method of super-potentials [31] in order to replace
$\hat{W}_{ij}$ by the expression of its source, which is now, indeed, known
all over the space. Therefore, the asymptotic behaviours of required
potentials have been fully determined with the appropriate accuracy.

Nevertheless, these computations are heavy555For instance, it took more than
$160$ CPU hours to compute the 1PN $\hat{W}_{ij}$ up to $\mathcal{O}(r^{-4})$.
At this order, its expression contains more than $10\,000$ terms, of which
$\sim 8\,000$ come from the homogeneous part. and some consistency checks are
required. The first and most stringent one is that they were performed in a
“double-blind” fashion. In addition, we investigated the three-dimensional
limits of our potentials, confirming that they agree with the asymptotic
expressions that were used to compute the three-dimensional surface terms in
[31]. We have also verified that the harmonicity relations
$\partial_{\mu}h^{\mu\nu}=0$ hold, up to 1PN order and to the highest
achievable order in $1/r$. At the 1PN order, those conditions read [31]

$\displaystyle\partial_{\mu}\overline{h}^{\mu 0}=$
$\displaystyle\frac{d-1}{2(d-2)}\partial_{t}V+\partial_{i}V_{i}$ (37a)
$\displaystyle+\frac{1}{c^{2}}\left[\partial_{t}\left(-\frac{(d-1)(d-3)}{(d-2)^{2}}\,\hat{K}+\frac{\hat{W}}{2}+\frac{(d-1)^{2}}{2(d-2)^{2}}\,V^{2}\right)+\partial_{i}\left(2\hat{R}_{i}+\frac{d-1}{d-2}\,VV_{i}\right)\right]\,,$
$\displaystyle\partial_{\mu}\overline{h}^{\mu i}=$
$\displaystyle\partial_{t}V_{i}+\partial_{j}\left(\hat{W}_{ij}-\frac{\hat{W}}{2}\,\delta_{ij}\right)+\frac{1}{c^{2}}\left[\partial_{t}\left(2\hat{R}_{i}+\frac{d-1}{d-2}\,VV_{i}\right)+4\,\partial_{j}\left(\hat{Z}_{ij}-\frac{\hat{Z}}{2}\,\delta_{ij}\right)\right]\,.$
(37b)

The cancellation of $\partial_{\mu}h^{\mu 0}$ has been checked exactly at
Newtonian order, and up to $\mathcal{O}(r^{-6})$ at 1PN order. The fact that
$\partial_{\mu}h^{\mu i}$ vanishes has been checked up to
$\mathcal{O}(r^{-7})$ at Newtonian order, and to $\mathcal{O}(r^{-5})$ at 1PN
order.

###  III.2 Compact terms

An interesting and non-trivial feature of the IR regularization scheme is that
it affects the evaluation of the potentials at the location of the particles.
This is due to the fact that most of the potentials have a non-compact support
source: the non-linear source terms extend towards infinity, so that the value
of the potential at, say, $\bm{y}_{1}$ is sensitive to the IR regularization
process. (Naturally, this effect does not impact the compact-support
potentials, whose sources are proportional to Dirac distributions, and which
are thus sensitive to the UV regularization only.) Therefore, the IR
regularization scheme affects the compact terms that involve some non-compact
potentials evaluated at $\bm{y}_{1}$ or $\bm{y}_{2}$.

Let us consider a potential $P$ with $d$-dimensional source $S$. As in the
previous section, we are interested in the IR behaviour of $S$, hence we will
expand it in the far zone. In practice, none of the sources we are interested
in develop poles (or equivalently, logarithms in three dimensions); so, we can
safely consider that

$S(\mathbf{x},t)=\sum_{p,q}\frac{\ell_{0}^{q\varepsilon}}{r^{p+q\varepsilon}}\,\varphi^{(\varepsilon)}_{p,q}(\mathbf{n},t)\,,$
(38)

where $\varphi^{(\varepsilon)}_{p,q}$ has no pole. The source takes the three-
dimensional limit

$S^{(d=3)}(\mathbf{x},t)=\sum_{p}\frac{\varphi_{p}(\mathbf{n},t)}{r^{p}}\,,\quad\text{with}\quad\varphi_{p}\equiv\sum_{q}\varphi^{(\varepsilon=0)}_{p,q}\,.$
(39)

Let us first deal with the Newtonian case, i.e., let us compute the difference
induced by the change of IR regularization scheme in the Poisson integral
evaluated in $\bm{y}_{1}$,

$P(\bm{y}_{1})=-\frac{\tilde{k}}{4\pi}\mathop{\mathrm{PF}}_{B=0}\int_{r^{\prime}>\mathcal{R}}\\!\mathrm{d}^{d}\mathbf{x}^{\prime}\left(\frac{r^{\prime}}{r_{0}}\right)^{B}\frac{S\left(\mathbf{x}^{\prime}\right)}{|\bm{y}_{1}-\mathbf{x}^{\prime}|^{d-2}}\,,$
(40)

where we have safely replaced $\mathbf{x}$ by $\bm{y}_{1}$ in the kernel, as
we integrate over the domain
$r^{\prime}\equiv|\mathbf{x}^{\prime}|>\mathcal{R}\gg|\bm{y}_{1}|$, where we
are free from UV divergences. Expanding the kernel at spatial infinity
according to

$\frac{1}{|\bm{y}_{1}-\mathbf{x}^{\prime}|^{d-2}}=\sum_{\ell\in\mathbb{N}}\frac{2^{\ell}}{\ell!}\frac{\Gamma\left(\frac{d-2}{2}+\ell\right)}{\Gamma\left(\frac{d-2}{2}\right)}\,\frac{y_{1}^{L}\,\hat{n}_{L}}{{r^{\prime}}^{d-2+\ell}}\,,$
(41)

and using the machinery developed for the treatment of volume terms, we find
that the difference between the $B\varepsilon$ and Hadamard regularization
schemes turns out to be

$\mathcal{D}P^{\text{Newt}}(\bm{y}_{1})=-\sum_{\ell\in\mathbb{N}}\frac{(2\ell-1)!!}{4\pi\ell!}\,y_{1}^{L}\,\sum_{q\neq
0}\frac{1}{q}\left[\frac{1}{\varepsilon}-q\ln\left(\frac{r_{0}}{\ell_{0}}\right)+\sum_{k=0}^{\ell}\frac{1}{2k-1}\right]\int\\!\frac{\mathrm{d}\Omega_{d-1}}{\Omega_{d-1}}\,\hat{n}_{L}\,\varphi^{(\varepsilon)}_{2-\ell,q}(\mathbf{n})\,,$
(42)

where we recall that the volume of the $(d-1)$-dimensional sphere is
$\Omega_{d-1}=2\pi^{d/2}/\Gamma(d/2)$. As expected, the scale $\mathcal{R}$
disappears from the difference in regularizations. The $q\neq 0$ criterion is
a natural consequence of the derivation. The main structural difference with
the formula for the volume terms (29) is the sum over $\ell$. Nevertheless,
this sum is finite by virtue of the form of the source term (38), since $\ell$
is bounded by $2-p_{0}$.

The formula (42), which is merely Newtonian, has to be generalized to higher
PN orders. For this purpose, instead of starting from the Poisson integral
(40), we have to use the $d$-dimensional propagator. Again, we can restrict
ourselves to an integration region $r^{\prime}>\mathcal{R}\gg|\bm{y}_{1}|$ in
which we are allowed to replace $\mathbf{x}$ by $\bm{y}_{1}$, so that

$P(\bm{y}_{1})=-\frac{\tilde{k}}{4\pi}\mathop{\mathrm{PF}}_{B=0}\int_{r^{\prime}>\mathcal{R}}\\!\frac{\mathrm{d}^{d}\mathbf{x}^{\prime}}{|\bm{y}_{1}-\mathbf{x}^{\prime}|^{d-2}}\left(\frac{r^{\prime}}{r_{0}}\right)^{B}\\!\int_{1}^{+\infty}\\!\\!\\!\mathrm{d}z\,\gamma_{\frac{1-d}{2}}(z)\,S\bigl{(}\mathbf{x}^{\prime},t-z|\bm{y}_{1}-\mathbf{x}^{\prime}|/c\bigr{)}\,,$
(43)

where the time-dependence of the source is now crucial. By PN-expanding this
propagator, it can be seen that the action of the full propagator, in a formal
PN sense, up to any PN order, is equivalent to the action of the mere Poisson
integral (40) but acting on the effective source

$S^{\text{eff}}\left(\mathbf{x}^{\prime},t,r_{1}^{\prime}\right)=\int_{1}^{+\infty}\\!\\!\\!\mathrm{d}z\,\gamma_{\frac{1-d}{2}}(z)\,S\left(\mathbf{x}^{\prime},t-zr_{1}^{\prime}/c\right)\,,$
(44)

with $\mathbf{x}^{\prime}$ just playing a spectator role. The PN expansion of
this effective source is given by

$\displaystyle\overline{S}^{\text{eff}}\equiv\overline{S}^{\text{eff}}_{\text{even}}+\overline{S}^{\text{eff}}_{\text{odd}}\,,$
(45a)
$\displaystyle\overline{S}^{\text{eff}}_{\text{even}}=\sum_{j\in\mathbb{N}}\frac{\sqrt{\pi}}{(2j)!}\frac{\Gamma\left(\frac{1+\varepsilon}{2}-j\right)}{\Gamma\left(\frac{1+\varepsilon}{2}\right)\Gamma\left(\frac{1}{2}-j\right)}\left(\frac{r_{1}}{c}\right)^{2j}S^{(2j)}(\mathbf{x},t)\,,$
(45b)
$\displaystyle\overline{S}^{\text{eff}}_{\text{odd}}=\sum_{j\in\mathbb{N}}\frac{(-)^{j}}{j!}\frac{2\sqrt{\pi}\,\Gamma\left(\varepsilon\right)}{\Gamma\left(\frac{1+\varepsilon}{2}\right)\Gamma\left(-j-\frac{\varepsilon}{2}\right)\Gamma\left(2j+2+\varepsilon\right)}\left(\frac{r_{1}}{c}\right)^{1+2j+\varepsilon}\int_{0}^{+\infty}\\!\mathrm{d}\tau\,\tau^{-\varepsilon}S^{(2j+2)}(\mathbf{x},t-\tau)\,.$
(45c)

These expansion series correspond to the PN-even and PN-odd parts of the
expansion of a homogeneous monopolar retarded solution of the wave equation in
$d$ dimensions (see App. A in [8] for more details). As expected, the $j=0$
term of $S^{\text{eff}}_{\text{even}}$ is exactly given by $S$. Note that the
PN-odd piece appears to be non-local in $d$ dimensions [8]. Nevertheless, as
already stated, we are merely interested in even terms666We will add the
already known 2.5 and 3.5PN terms [26] to our final result, together with the
4PN dissipative non-local tail term (see Eqs. (6.11) and (6.12) in the
companion paper [34]). and, thus, will not consider it. For a potential
entering the source at a given PN order, we can simply apply the Newtonian
result (42), but using the effective source $S^{\text{eff}}_{\text{even}}$
truncated at the appropriate PN order following Eqs. (45).

The required potentials that developp a non-vanishing difference at point
$\bm{y}_{1}$ are $\hat{W}_{1}^{2\text{PN}}$, $\hat{Z}_{1}^{1\text{PN}}$,
$\hat{X}_{1}^{1\text{PN}}$, $\hat{T}_{1}^{\text{N}}$ and
$\hat{M}_{1}^{\text{N}}$ (with the index 1 denoting the value at particle 1
and the superscript the PN order) while the ones of the super-potentials
vanish. The calculations reveal that all the non-vanishing corrections in our
potentials are connected. Indeed, we found the interesting, but probably not
very profound, relations

$\mathcal{D}\hat{W}_{1}^{2\text{PN}}=\frac{4}{c^{2}}\,\mathcal{D}\hat{Z}_{1}^{1\text{PN}}=-\frac{2\varepsilon}{c^{2}}\,\mathcal{D}\hat{X}_{1}^{1\text{PN}}=\frac{8}{3c^{4}}\,\mathcal{D}\hat{T}_{1}^{\text{N}}=-\frac{16\varepsilon}{3c^{4}}\,\mathcal{D}\hat{M}_{1}^{\text{N}}\,.$
(46)

None of the other potentials receives corrections at the required order. These
relations are valid up to the $\mathcal{O}(\varepsilon^{0})$ order only, as
the $\mathcal{O}(\varepsilon)$ remainders do not play any role in the IR
dimensional regularization of the compact-support terms. Note that only the
trace of the potential $\hat{M}_{ij}$ at Newtonian order (i.e.
$\hat{M}\equiv\hat{M}_{ii}$) as well as the potential $\hat{X}$ at 1PN order
develop poles. Moreover, only the “scalar” sector of the potentials is
affected: neither $\hat{R}_{i}$ nor $\hat{Y}_{i}$ are modified at the
particles’ positions; regarding the tensor potentials $\hat{Z}_{ij}$ and
$\hat{M}_{ij}$, only their traces are impacted. In addition to those
relations, the difference itself can be compactly written in terms of the
Newtonian moment of inertia $I\equiv m_{1}\bm{y}_{1}^{2}+m_{2}\bm{y}_{2}^{2}$
as

$\mathcal{D}\hat{M}_{1}^{\text{N}}=\left(\frac{3}{8}-\varepsilon\right)c^{2}\,\mathcal{D}\hat{X}_{1}^{1\text{PN}}=\frac{G^{2}(m_{1}+m_{2})}{4}\left[\ln\left(\frac{r_{0}\sqrt{\bar{q}}}{\ell_{0}}\right)-\frac{1}{2\varepsilon}+\frac{1}{2}\right]\,I^{(4)}\,,$
(47)

where we recall that $\bar{q}\equiv 4\pi\mathrm{e}^{\gamma_{E}}$. The
regularization induced differences of those potentials yield corrections in
the effective mass $\tilde{\mu}_{1}$ (see Eq. (2.17) in [31]), and in a few
compact-support terms that all enter the source mass quadrupole moment at 4PN
order.

###  III.3 Surface terms

Turning now to the surface terms, we will not calculate the difference between
their values obtained in the two regularization schemes, but rather directly
compute them in $d$ dimensions and show that they actually vanish. As
presented in [31], surface terms appearing in the mass quadrupole are of two
kinds: “Laplacian” and “divergence” terms.

  * •

A surface term of Laplacian type reads

$T_{L}=\mathop{\mathrm{PF}}_{B=0}\int_{r>\mathcal{R}}\\!\\!\frac{\mathrm{d}^{d}\mathbf{x}}{\ell_{0}^{\varepsilon}}\,\left(\frac{r}{r_{0}}\right)^{B}\,\hat{x}_{L}\,\Delta
G\,,$ (48)

where $G$ is a $d$-dimensional function (in practice a product of potentials).
As before, we have restricted ourselves to an integral over $r>\mathcal{R}$.
Like for the volume terms, the function $G$ we will consider can be expanded
near spatial infinity as:

$G(\mathbf{x})=\sum_{p,q}\frac{\ell_{0}^{q\varepsilon}}{r^{p+q\varepsilon}}\,\gamma^{(\varepsilon)}_{p,q}(\mathbf{n})\,,$
(49)

allowing the $\gamma^{(\varepsilon)}_{p,q}$ to contain poles. Inserting this
expansion into the integral (48), performing an integration by parts, using
$\Delta(r^{B}\hat{x}_{L})=B(B+d+2\ell-2)r^{B+\ell-2}\hat{n}_{L}$, and dropping
the integrated terms that are vanishing by analytic continuation in either $B$
or $\varepsilon$ near zero, we are led to

$T_{L}=-\left(2\ell+1+\varepsilon\right)\int\\!\mathrm{d}\Omega_{d-1}\hat{n}_{L}\,\gamma^{(\varepsilon)}_{\ell+1,1}(\mathbf{n})\,.$
(50)

  * •

A surface term of divergence type reads

$K=\mathop{\mathrm{PF}}_{B=0}\int_{r>\mathcal{R}}\\!\\!\frac{\mathrm{d}^{d}\mathbf{x}}{\ell_{0}^{\varepsilon}}\,\left(\frac{r}{r_{0}}\right)^{B}\,\partial_{i}H^{i}\,,$
(51)

where $H^{i}$ is a product of potentials or super-potentials. We again have
the far-zone expansion

$H^{i}(\mathbf{x})=\sum_{p,q}\frac{\ell_{0}^{q\varepsilon}}{r^{p+q\varepsilon}}\,\eta^{(\varepsilon)i}_{p,q}(\mathbf{n})\,,$
(52)

where the $\eta^{(\varepsilon)i}_{p,q}$ can contain poles. A similar procedure
yields

$K=\int\\!\mathrm{d}\Omega_{d-1}n_{i}\,\eta^{(\varepsilon)i}_{2,1}(\mathbf{n})\,.$
(53)

The noteworthy point is that the $d$-dimensional values of both types of
surface terms, in Eqs. (50) and (53), are non zero only in the special case
$q=1$. Now, as explicit in their expressions displayed in App. C of [31], the
integrands of those terms are made of non-linear products of potentials. As
such, the coefficients of their asymptotic expansions all bear $q\geqslant 2$,
which make them vanish.

More precisely, it is immediate to see, from the Green function
$\Delta^{-1}\delta^{(d)}(\mathbf{x}-\bm{y}_{1})=-\tilde{k}\,r_{1}^{-1-\varepsilon}/(4\pi)$,
that the expansions of compact support potentials have $q=1$. Now, the sources
of non-compact support potentials are made of products of compact support
ones, and the iteration of the Poisson integrals cannot reduce the number of
$q$’s. The latter fact can be understood from the Matthieu formula (34),
together with the fact that the expansion of the homogeneous solution bears
$q=1$, as explicitly shown by Eq. (III.1.2) where the second term is
$\propto\partial_{L}r^{2k+2-d}\propto r^{2k-\ell-1-\varepsilon}$. Thus, all
individual potentials $V$, $V_{i}$, $\hat{W}_{ij}$, etc. admit asymptotic
expansions with $q\geqslant 1$. A similar argument applies to the case of
super-potentials, which implies that all non-linear products of (derivatives
of) potentials or super-potentials have $q\geqslant 2$.

The conclusion is that the surface terms are vanishing in $d$ dimensions, so
that we have simply to subtract their Hadamard values from the final result:

$\mathcal{D}I_{ij}^{\text{Surf}}=-I_{ij}^{\text{Surf,Had}}\,.$ (54)

###  III.4 The “extra” term

The argument developed in [31] to discard the contribution of the $B(d-3)$
piece of $I_{ij}$ [the last term in (I)] does not a priori hold with IR
dimensional regularization; further investigation is required. We can restrict
ourselves to an integration in the far zone $r>\mathcal{R}$, and recast

$I_{ij}^{\text{extra}}=C_{d}\,\mathop{\mathrm{PF}}_{B=0}\sum_{k\in\mathbb{N}}\frac{1}{2^{k}\,k!}\frac{\Gamma\left(\frac{d}{2}+3\right)}{\Gamma\left(\frac{d}{2}+3+k\right)}\left(\frac{1}{c}\frac{\mathrm{d}}{\mathrm{d}t}\right)^{2k}\int_{r>\mathcal{R}}\\!\\!\frac{\mathrm{d}^{d}\mathbf{x}}{\ell_{0}^{\varepsilon}}\,\left(\frac{r}{r_{0}}\right)^{B}\,B\,\frac{\hat{x}_{aij}\,x_{b}}{r^{2-2k}}\,\overline{\Sigma}_{ab}\,,$
(55)

with $C_{d}=-\frac{2}{c^{2}}\frac{(d-3)(d+2)}{d(d-2)(d+4)}$. The source
density $\overline{\Sigma}_{ab}$ (2)–(3) is composed of two pieces: The first
one, coming from the stress-energy tensor $T^{\mu\nu}$, involves Dirac
distributions and, thus, does not enter our computation, since we take
$r>\mathcal{R}$. The second piece, entailing the non-linear (NL) gravitational
interactions $\Lambda^{ij}$, is composed, as such, of products of potentials.
We expand it as

$\overline{\Sigma}_{ab}^{\text{NL}}(\mathbf{x},t)=\sum_{p,q}\frac{\ell_{0}^{q\varepsilon}}{r^{p+q\varepsilon}}\,\sigma_{p,q}^{(\varepsilon)ab}(\mathbf{n},t)\,,$
(56)

where the functions $\sigma_{p,q}^{(\varepsilon)ab}$ contain many PN orders
and can develop poles. Inserting (56) into (55), we can perform the radial
integration, which yields

$\displaystyle
I_{ij}^{\text{extra}}=-C_{d}\sum_{k\in\mathbb{N}}\frac{1}{2^{k}\,k!}\frac{\Gamma\left(\frac{d}{2}+3\right)}{\Gamma\left(\frac{d}{2}+3+k\right)}\left(\frac{1}{c}\frac{\mathrm{d}}{\mathrm{d}t}\right)^{2k}\left[\int\\!\mathrm{d}\Omega_{d-1}\,\hat{n}_{aij}n_{b}\,\sigma_{5+2k,1}^{(\varepsilon)ab}(\mathbf{n},t)\right]\,.$
(57)

Exactly as in the case of the surface terms, only the terms with $q=1$
contribute. However, recalling that $\sigma_{p,q}^{(\varepsilon)ab}$ is made
of products of potentials, it cannot involve $q=1$ terms. We see therefore
that, even with the present IR $B\varepsilon$ regularization scheme, the
“extra” piece does not contribute:

$\mathcal{D}I_{ij}^{\text{extra}}=0\,.$ (58)

##  IV Final expression of the source quadrupole moment

Summing up the pieces computed in the previous sections, and adding the local
IR shift $\chi^{i}_{1,2}$ coming from the cancellation of the remaining poles
in the conservative equations of motion [8] (as described in Appendix B), we
obtain the full dimensionally regularized source mass quadrupole moment at 4PN
order,

$I_{ij}=I_{ij}^{\text{Had}}+I_{ij}^{\text{non-
loc}}+\mathcal{D}I_{ij}+\delta_{\chi}I_{ij}\,.$ (59)

Here, $I_{ij}^{\text{Had}}$ stands for the end result of [31],777Note that,
after publishing our partial result in [31], we have spotted an error in the
$d$-dimensional computation of the value of the potential $\hat{R}_{i}$ at 1PN
order evaluated at $\bm{y}_{A}$. This error induces a small change in the
value of $I_{ij}^{\text{Had}}$ that has been taken into account in the final
result of the renormalized mass quadrupole, which will be displayed in the
companion paper [34]. then $I_{ij}^{\text{non-loc}}$ is the non-local part
computed in Sec. II and given by Eq. (18), $\mathcal{D}I_{ij}$ is the sum of
the four contributions (19) whose calculation has been detailed above, and
$\delta_{\chi}I_{ij}$ is the contribution of the IR shift of Appendix B.
Remember that the non-local sector $I_{ij}^{\text{non-local}}$ is composed of
both a direct tail contribution (14) and a shift contribution [see Eq. (17)].

As the complete expression for the regularized source moment is long and not
very enlightening nor interesting per se, we do not display it, but rather
discuss its most striking feature: the remaining poles. Indeed, the source
moment $I_{ij}$ contains poles, which will have to be compensated by the
proper treatment in dimensional regularization of the non-linear interactions
(such as tails-of-tails) relating the source and radiative moments. This will
be proven in the companion paper [34].

While the non-local tail $I_{ij}^{\text{non-loc}}$ and shifted pieces
$\delta_{\chi}I_{ij}$ are purely 4PN contributions, the effect of the
dimensional regularization $\mathcal{D}I_{ij}$ starts already at the 3PN
order. It can be written in a convenient form as

$\mathcal{D}I_{ij}^{\text{3PN}}=-\beta_{I}\,\frac{G^{2}M^{2}}{c^{6}}\left[\Pi_{\varepsilon}+\frac{246\,299}{44\,940}\right]\,I_{ij}^{(2)}+2\beta_{I}\,\frac{G^{2}M}{c^{6}}\left[\Pi_{\varepsilon}+\frac{252\,599}{44\,940}\right]\,P_{\langle
i}P_{j\rangle}\,,$ (60)

where we have used a notation similar to (6) for the “dressed” pole:

$\Pi_{\varepsilon}\equiv\ln\left(\frac{\sqrt{\bar{q}}\,r_{0}}{\ell_{0}}\right)-\frac{1}{2\varepsilon}\,,\qquad\bar{q}\equiv
4\pi\mathrm{e}^{\gamma_{E}}\,.$ (61)

In Eq. (60), $P_{i}$ is the constant linear momentum and
$\beta_{I}=-\frac{214}{105}$ is the coefficient associated with the
renormalization of the mass quadrupole moment [39]. The appearance of this
coefficient, which is known to be associated with the “tail-of-tail”
interaction [23], indicates the soundness of the removal of the poles by the
correct treatment in dimensional regularization of the non-linear interactions
in the radiative quadrupole.

As for the poles showing up at the 4PN order, once taken the 1PN correction of
Eq. (60) into account, they can be expressed in the center-of-mass (CoM) frame
as

$I_{ij}^{\text{pole,CoM,4PN}}=-\frac{G^{2}M}{2c^{8}\varepsilon}\left[\frac{12}{7}\,I_{a\langle
i}^{(2)}\,I_{j\rangle a}^{(2)}-\frac{24}{7}\,I_{a\langle i}^{(1)}\,I_{j\rangle
a}^{(3)}-\frac{4}{7}\,I_{a\langle i}\,I_{j\rangle
a}^{(4)}+\frac{4}{3}I_{a\langle i}^{(3)}J_{j\rangle|a}\right]\,,$ (62)

where $J_{i|j}$ is the $d$-dimensional constant angular momentum defined in
[27]. Once again, the fact that this pole can be recast as a combination of
non-linear interactions (in particular the expected coupling $M\times
I_{ij}\times I_{ij}$ corresponding to “tails-of-memory”) is a strong
indication that it will be compensated by the proper dimensional-
regularization treatment of the tail/memory effects. This will be the topic of
the next paper [34].

###### Acknowledgements.

We thank Gabriel Luz Almeida, Laura Bernard, Stefano Foffa, Sylvain Marsat,
Rafael Porto and Riccardo Sturani for interesting discussions. We are grateful
to Laura Bernard for providing us with the IR shift coming from the 4PN
equations of motion (see Appendix B). F.L. received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (grant agreement No 817791).

##  Appendix A Generalized Riesz formulae

The standard Riesz formula in arbitrary dimension $d$ (with
$a,b\in\mathbb{C}$) reads:

$\displaystyle\int\mathrm{d}^{d}x\,r_{1}^{a}r_{2}^{b}=\pi^{d/2}\frac{\Gamma(-\frac{a+b+d}{2})\Gamma(\frac{a+d}{2})\Gamma(\frac{b+d}{2})}{\Gamma(-\frac{a}{2})\Gamma(-\frac{b}{2})\Gamma(\frac{a+b}{2}+d)}r_{12}^{a+b+d}\,.$
(63)

It is a priori valid when the integral converges, i.e., for $a+d>0$, $b+d>0$
and $a+b+d<0$. However, the result can be extended by analytic continuation
everywhere but, possibly, a countable set of parameter values. On the other
hand, this formula may be generalized to similar integrals involving
additional STF angular factors $\hat{x}^{L}$ and various $r^{2}$. Among the
most useful ones, we have

$\displaystyle\int\mathrm{d}^{d}x\,\hat{x}^{L}r_{1}^{a}r_{2}^{b}=\pi^{d/2}\frac{\Gamma(-\frac{a+b+d}{2})}{\Gamma(-\frac{a}{2})\Gamma(-\frac{b}{2})}\,r_{12}^{a+b+d}\sum_{s=0}^{\ell}\genfrac{(}{)}{0.0pt}{}{\ell}{s}\frac{\Gamma(\frac{a+d}{2}+s)\Gamma(\frac{b+d}{2}+\ell-s)}{\Gamma(\frac{a+b}{2}+d+\ell)}y_{1}^{\langle
L-S}y_{2}^{S\rangle}\,,$ (64a)
$\displaystyle\int\mathrm{d}^{d}x\,r^{2}\hat{x}^{L}r_{1}^{a}r_{2}^{b}=-\pi^{d/2}\frac{\Gamma(-\frac{a+b+d}{2}-1)}{\Gamma(-\frac{a}{2})\Gamma(-\frac{b}{2})}\,r_{12}^{a+b+d}\sum_{s=0}^{\ell}\genfrac{(}{)}{0.0pt}{}{\ell}{s}\frac{\Gamma(\frac{a+d}{2}+s)\Gamma(\frac{b+d}{2}+\ell-s)}{\Gamma(\frac{a+b}{2}+d+\ell+1)}\times$
$\displaystyle\qquad\qquad~{}\times y_{1}^{\langle
L-S}y_{2}^{S\rangle}\Big{[}\Big{(}\frac{b}{2}-s+1\Big{)}\Big{(}\frac{b+d}{2}+\ell-s\Big{)}y_{1}^{2}+\Big{(}\frac{a}{2}-\ell+s+1\Big{)}\Big{(}\frac{a+d}{2}+s\Big{)}y_{2}^{2}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad+2\Big{(}\frac{a+d}{2}+s\Big{)}\Big{(}\frac{b+d}{2}+\ell-s\Big{)}(y_{1}y_{2})\Big{]}\,,$
(64b)

with $\genfrac{(}{)}{0.0pt}{}{\ell}{s}$ denoting the usual binomial
coefficient, and $y_{A}^{2}=\bm{y}_{A}^{2}$,
$(y_{1}y_{2})=\bm{y}_{1}\cdot\bm{y}_{2}$.

Let us sketch for instance the proof of the first one. The starting point
consists in rewriting $x^{L}$ as
$(r_{1}n_{1}+y_{1})^{L}\equiv(r_{1}n_{1}^{i_{1}}+y_{1}^{i_{1}})\cdots(r_{1}n_{1}^{i_{\ell}}+y_{1}^{i_{\ell}})$,
where $n^{i}_{1}$ is the unit vector $n^{i}_{1}=(x^{i}-y_{1}^{i})/r_{1}$, and
expanding the product. After permuting the integration and summation symbols,
we get a sum of elementary integrals of the form

$\int\mathrm{d}^{d}x\,r_{1}^{a+s}r_{2}^{b}\,\hat{n}_{1}^{L}\,,$ (65)

for a summation index $0\leqslant s\leqslant\ell$, with $n_{1}^{L}\equiv
n_{1}^{i_{1}}\cdots n_{1}^{i_{\ell}}$. For each integral, the parameters $a$,
$b$ and $d$ are chosen in a domain of $\mathbb{C}^{3}$ where the convergence
is guaranteed. We can always manage to avoid the situation where one of them
is an integer. We then express the factors $r_{1}^{a+s}\hat{n}_{1}^{L}$ as a
multi-derivative using the relation

$\displaystyle\hat{\partial}_{L}r_{1}^{\alpha}=(-)^{\ell}\hat{\partial}_{1L}r_{1}^{\alpha}=2^{\ell}\frac{\Gamma(\alpha/2+1)}{\Gamma(\alpha/2-\ell+1)}r_{1}^{\alpha-\ell}\hat{n}_{1}^{L}\,,$
(66)

with $\partial_{1i}=\partial/\partial y_{1}^{i}$, valid as long as
$\hat{\partial}_{1L}r_{1}^{\alpha}$ does not involve any distributional
contribution, which is indeed the case when $\alpha$ is a non-integer. Next,
we commute the derivatives and the integral (taking for instance $2s<d$,
$2s<-a<d$ and $a+d<-b<d$ to bypass any convergence issue). We find

$\displaystyle\int\mathrm{d}^{d}x\,\hat{x}^{L}r_{1}^{a}r_{2}^{b}=\sum_{s=0}^{\ell}\genfrac{(}{)}{0.0pt}{}{\ell}{s}\frac{(-)^{s}}{2^{s}}\frac{\Gamma(a/2+1)}{\Gamma(a/2+s+1)}y_{1}^{\langle
L-S}\partial_{1S\rangle}\int\mathrm{d}^{d}x\,r_{1}^{a+2s}r_{2}^{b}\,.$ (67)

We now use the standard Riesz formula (63), which produces a dimensional
factor $r_{12}^{a+b+d+2s}$, and apply the derivative $\hat{\partial}_{1S}$ as
in Eq. (66). At this stage, we simplify some factors under the sum by noticing
that, for $s\in\mathbb{N}$,

$\displaystyle\Gamma(z+s+1)\Gamma(-z-s)=\frac{\pi}{\sin[-(z+s)\pi]}=(-)^{s}\Gamma(z+1)\Gamma(-z)\,,$
(68)

and expand the position-dependent factor $r_{12}^{a+b+d+s}n_{12}^{\langle
S}\,y_{1}^{L-S\rangle}=r_{12}^{a+b+d}(y_{1}-y_{2})^{\langle
S}y_{1}^{L-S\rangle}$. The result is obtained by resummation:888It is a mere
consequence of the relation
${}_{2}F_{1}(a,b;c;1)=\frac{\Gamma(c)\Gamma(c-a-b)}{\Gamma(c-a)\Gamma(c-b)}$
for $\Re(c-a-b)>0$, where ${}_{2}F_{1}(a,b;c;z)$ is the Gaussian
hypergeometric function, with standard notations.

$\displaystyle\sum_{s=0}^{+\infty}\frac{(-)^{s}}{(\ell-
k-s)!s!}\frac{\Gamma(\frac{a+d}{2}+s+k)}{\Gamma(\frac{a+b}{2}+d+s+k)}$
$\displaystyle\qquad\qquad\qquad=\frac{\Gamma(\frac{a+d}{2}+k)}{(\ell-k)!\Gamma(\frac{a+b}{2}+d+k)}\frac{\Gamma(\frac{a+b}{2}+d+k)\Gamma(\frac{b+d}{2}+\ell-k)}{\Gamma(\frac{a+b}{2}+d+\ell)\Gamma(\frac{b+d}{2})}\,,$
(69)

for $k$ non-negative integer lower than $\ell$.

A third particularly useful generalization of Eq. (63) reads

$\displaystyle\int\mathrm{d}^{d}x\,r_{1}^{a}\,r_{2}^{b}\,n_{1}^{L}n_{2\,K}=\pi^{d/2}\frac{\Gamma(-\frac{a+b+\ell+k+d}{2})\Gamma(\frac{a+b+\ell+k+d}{2}+1)}{\Gamma(-\frac{a-\ell}{2})\Gamma(-\frac{b-k}{2})\Gamma(\frac{a+b+\ell+k}{2}+d)}r_{12}^{a+b+d}\times$
$\displaystyle\qquad\qquad~{}\times\sum_{i=0}^{[\frac{\ell+k}{2}]}\sum_{j=\max(0,2i-k)}^{\min{(2i,\ell)}}\sum_{p=\max(0,j-i)}^{[\frac{j}{2}]}\\!\\!\\!\\!\\!\frac{\Gamma(\frac{a+\ell+d}{2}+i-j)\Gamma(\frac{b+k+d}{2}+j-i)}{\Gamma(\frac{a+b-\ell-k+d}{2}+i+1)}\times$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\times\frac{(-)^{k+i+j}\ell!k!\,\delta_{(2I2P-2J}n_{12\,KJ-2I}\delta_{J-2P)}^{(J-2P}n_{12}^{L-J}\delta^{2P)}}{2^{2i+2p-j}(\ell-j)!(k+j-2i)!p!(i+p-j)!(j-2p)!}\,,$
(70)

with the convention that the contravariant indices belong to
$\\{L\\}=\\{i_{1},\cdots,i_{\ell}\\}$ and the covariant ones to
$\\{K\\}=\\{j_{1},\cdots,j_{k}\\}$, while
$\delta^{J-2P}_{J-2P}=\delta^{a_{1}}_{b_{1}}\cdots\delta^{a_{j-2p}}_{b_{j-2p}}$,
$\delta^{2P}=\delta^{a_{j-2p+1}a_{j-2p+2}}\cdots\delta^{a_{j-1}a_{j}}$, and
$\delta_{2I2P-2J}=\delta_{b_{j-2p+1}b_{j-2p+2}}\cdots\delta_{b_{2i-j-1}b_{2i-j}}$.
The upper (lower) symmetrization only affects the contravariant (covariant)
indices, respectively.

The most delicate task in the derivation is the symmetrization over the multi-
indices $L$ and $K$ of the multi-derivative

$\displaystyle\delta^{(2S}\partial_{1}^{L-2S)}\delta_{(2R}\partial_{2\,K-2R)}r_{12}^{a+b+\ell+k-2(s+r)+d}\,.$
(71)

It is achieved with the help of the convenient identity

$\displaystyle\delta^{(2I}n_{12}^{L-2S}n_{12\,K-2R-2I)}$
$\displaystyle=i!\frac{(\ell+k-2s-2r-2i)!}{(\ell+k-2s-2r)!}\sum_{j=\max(0,2i+2r-k)}^{\min(2i,\ell-2s)}\sum_{p=\max(0,j-i)}^{[\frac{j}{2}]}\frac{1}{2^{2p-j}p!}\times$
$\displaystyle\qquad\times\frac{(\ell-2s)!}{(\ell-2s-j)!}\frac{(k-2r)!}{(k+j-2r-2i)!}\frac{1}{(i+p-j)!(j-2p)!}\times$
$\displaystyle\qquad\times\delta_{(2P2I2R-2J}n_{12\,KJ-2R-2I}\delta_{J-2P)}^{(J-2P}n_{12}^{L-2S-J}\delta^{2P)}\,.$
(72)

In this manner, we find that the integral on the left-hand side of Eq. (A),
after the appropriate change of indices, may be put in the form

$\displaystyle\sum_{i=0}^{[\frac{\ell+k}{2}]}\sum_{j=\max(0,2i-k)}^{\min(2i,\ell)}\sum_{p=\max(0,j-i)}^{[\frac{j}{2}]}\sum_{s=0}^{p}\sum_{r=0}^{p+i-j}r_{12}^{a+b+d}\,\lambda_{s,r}^{i,j,p}\,\delta_{(2I2P-2J}n_{12\,KJ-2I}\delta_{J-2P)}^{(J-2P}n_{12}^{L-J}\delta^{2P)}\,,$
(73)

with

$\displaystyle\lambda_{s,r}^{i,j,p}$
$\displaystyle=\frac{(-)^{s+r+\ell}2^{j-2p-2i}}{(\ell-j)!(k+j-2i)!(p-s)!(i-r+p-j)!(j-2p)!}\frac{\ell!k!}{s!r!}\frac{\Gamma(\frac{a-\ell}{2}+1)}{\Gamma(\frac{a+\ell}{2}+1-s)}\times$
$\displaystyle\times\frac{\Gamma(\frac{b-k}{2}+1)}{\Gamma(\frac{b+k}{2}+1-r)}\frac{\Gamma(\frac{a+b+\ell+k+d}{2}-s-r+1)}{\Gamma(\frac{a+b-\ell-k+d}{2}+i+1)}\bigg{[}\int\mathrm{d}^{d}x\,r_{1}^{a+\ell-2s}r_{2}^{b+k-2r}\bigg{]}_{r_{12}=1}\,.$
(74)

Resorting to techniques similar to those employed to compute the previous
integrals, we can perform explicitly the sums over $s$ and $r$,

$\displaystyle\sum_{s=0}^{+\infty}\sum_{r=0}^{+\infty}\lambda_{s,r}^{i,j,p}$
$\displaystyle=\frac{\ell!k!}{2^{2i+2p-j}(\ell-j)!(k+j-2i)!}\frac{\pi^{d/2}(-)^{k+i+j}}{p!(i+p-j)!(j-2p)!}\times$
$\displaystyle\times\frac{\Gamma(-\frac{a+b+\ell+k+d}{2})\Gamma(\frac{a+b+\ell+k+d}{2}+1)}{\Gamma(-\frac{a-\ell}{2})\Gamma(-\frac{b-k}{2})\Gamma(\frac{a+b+\ell+k}{2}+d)}\frac{\Gamma(\frac{a+\ell+d}{2}+i-j)\Gamma(\frac{b+k+d}{2}+j-i)}{\Gamma(\frac{a+b-\ell-k+d}{2}+i+1)}\,,$
(75)

which yields the desired result.

##  Appendix B The local piece of the IR shift

In addition to the non-local shift due to tails, as obtained in Eq. (15), the
source quadrupole has to be shifted by the exact same IR shift that was
applied to the conservative sector to compensate the remaining poles in the
equations of motion [8]. This IR shift $\chi^{i}_{1,2}$ is local and starts at
the 4PN order. Thus, its contribution to the source quadrupole moment is
simply given by $\delta_{\chi}I_{ij}=2m_{1}\,y_{1}^{\langle
i}\chi_{1}^{j\rangle}+1\leftrightarrow 2$. The local IR shift $\chi_{1}^{i}$
can be decomposed as

$\chi_{1}^{i}=\frac{1}{c^{8}}\left(\frac{1}{\varepsilon}\,\chi_{1}^{i\,(-1)}+\chi_{1}^{(0,y_{1})}\,y^{i}_{1}+\chi_{1}^{(0,n_{12})}\,n^{i}_{12}+\chi_{1}^{(0,v_{1})}\,v^{i}_{1}+\chi_{1}^{(0,v_{12})}\,v^{i}_{12}\right)\,.$
(76)

Recalling that $\bar{q}\equiv 4\pi e^{\gamma_{E}}$, $m\equiv m_{1}+m_{2}$ and
$\nu\equiv m_{1}m_{2}/m^{2}$, the various terms of Eq. (76) are given by

$\displaystyle\chi_{1}^{i\,(-1)}$
$\displaystyle=\frac{8G^{4}m^{3}\nu}{5r_{12}{}^{4}}\bigg{[}m(n_{12}y_{1})n_{12}^{i}+r_{12}(2m_{1}-m_{2})n_{12}^{i}-\frac{m}{3}\,y_{1}^{i}\bigg{]}$
$\displaystyle-\frac{4G^{3}m^{2}\nu\,n_{12}^{i}}{5r_{12}^{3}}\bigg{[}2v_{12}^{2}\left(r_{12}(m_{1}-2m_{2})-3m\left(n_{12}y_{1}\right)\right)-8mr_{12}\left(v_{1}v_{12}\right)$
$\displaystyle\hskip
56.9055pt+\left(n_{12}v_{12}\right)^{2}\left(30m\left(n_{12}y_{1}\right)+9m_{2}r_{12}\right)+6m\left(n_{12}v_{12}\right)\left(2r_{12}\left(n_{12}v_{1}\right)-3\left(y_{1}v_{12}\right)\right)\bigg{]}$
$\displaystyle+\frac{4G^{3}m^{2}\nu\,v_{12}^{i}}{5r_{12}{}^{3}}\bigg{[}8m\left(r_{12}\left(n_{12}v_{1}\right)-(y_{1}v_{12})\right)+\left(n_{12}v_{12}\right)\left(18m\left(n_{12}y_{1}\right)+5m_{2}r_{12}\right)\bigg{]}$
$\displaystyle-\frac{8G^{3}m^{3}\nu\,y_{1}^{i}}{15r_{12}^{3}}\bigg{[}3\left(n_{12}v_{12}\right)^{2}-v_{12}^{2}\bigg{]}-\frac{16G^{3}m^{3}\nu\left(n_{12}v_{12}\right)}{15r_{12}^{2}}\,v_{1}^{i}\,,$
(77a) $\displaystyle\chi_{1}^{(0,y_{1})}=$
$\displaystyle\frac{G^{4}m^{4}\nu}{75r_{12}^{4}}\left[80\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+80\ln\left(\frac{r_{12}}{r_{0}}\right)+73\right]$
$\displaystyle+\frac{G^{3}m^{3}\nu}{75r_{12}^{3}}\bigg{[}\left(n_{12}v_{12}\right)^{2}\left(180\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+120\ln\left(\frac{r_{12}}{r_{0}}\right)+359\right)$
$\displaystyle\hskip 56.9055pt-
v_{12}^{2}\left(60\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+40\ln\left(\frac{r_{12}}{r_{0}}\right)+133\right)\bigg{]}\,,$
(77b) $\displaystyle\chi_{1}^{(0,n_{12})}=$
$\displaystyle-\frac{G^{4}m^{3}\nu}{2100r_{12}^{4}}\bigg{[}r_{12}\left(6720(2m_{1}-m_{2})\left(\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+\ln\left(\frac{r_{12}}{r_{0}}\right)\right)+8861m_{1}-5183m_{2}\right)$
$\displaystyle\hskip
56.9055pt+42m\left(n_{12}y_{1}\right)\left(160\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+160\ln\left(\frac{r_{12}}{r_{0}}\right)+61\right)\bigg{]}$
$\displaystyle+\frac{G^{3}m^{2}\nu}{4200r_{12}^{3}}\bigg{[}2r_{12}v_{12}^{2}\left((m_{1}-2m_{2})\left(5040\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+3360\ln\left(\frac{r_{12}}{r_{0}}\right)\right)-121916m_{1}-37995m_{2}\right)$
$\displaystyle\hskip
56.9055pt-84mv_{12}^{2}\left(n_{12}y_{1}\right)\left(360\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+240\ln\left(\frac{r_{12}}{r_{0}}\right)+223\right)$
$\displaystyle\hskip
56.9055pt+2100m\left(n_{12}v_{12}\right)^{2}\left(n_{12}y_{1}\right)\left(72\ln\left(\bar{q}r_{0}{}^{2}\right)+48\ln\left(\frac{r_{12}}{r_{0}}\right)+35\right)$
$\displaystyle\hskip
56.9055pt+3r_{12}\left(n_{12}v_{12}\right)^{2}\left(15120m_{2}\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+10080m_{2}\ln\left(\frac{r_{12}}{r_{0}}\right)+126959m_{1}+45908m_{2}\right)$
$\displaystyle\hskip
56.9055pt+21r_{12}\left(n_{12}v_{1}\right)\left(n_{12}v_{12}\right)(m_{1}+m_{2})\left(2880\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+1920\ln\left(\frac{r_{12}}{r_{0}}\right)-9661\right)$
$\displaystyle\hskip
56.9055pt-252m\left(n_{12}v_{12}\right)\left(y_{1}v_{12}\right)\left(360\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+240\ln\left(\frac{r_{12}}{r_{0}}\right)+223\right)$
$\displaystyle\hskip
56.9055pt-42mr_{12}\left(v_{1}v_{12}\right)\left(960\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+640\ln\left(\frac{r_{12}}{r_{0}}\right)-3007\right)\bigg{]}\,,$
(77c) $\displaystyle\chi_{1}^{(0,v_{1})}=$
$\displaystyle\frac{G^{3}m^{3}\nu\left(n_{12}v_{12}\right)}{600r_{12}^{2}}\left[960\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+640\ln\left(\frac{r_{12}}{r_{0}}\right)+2113\right]\,,$
(77d) $\displaystyle\chi_{1}^{(0,v_{12})}=$
$\displaystyle-\frac{G^{3}m^{2}\nu}{4200r_{12}^{3}}\bigg{[}5r_{12}\left(n_{12}v_{12}\right)\left(5040m_{2}\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+3360m_{2}\ln\left(\frac{r_{12}}{r_{0}}\right)+55132m_{1}+41681m_{2}\right)$
$\displaystyle\hskip
56.9055pt+42mr_{12}\left(n_{12}v_{1}\right)\left(960\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+640\ln\left(\frac{r_{12}}{r_{0}}\right)-3007\right)$
$\displaystyle\hskip
56.9055pt+252m\left(n_{12}v_{12}\right)\left(n_{12}y_{1}\right)\left(360\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+240\ln\left(\frac{r_{12}}{r_{0}}\right)+223\right)$
$\displaystyle\hskip
56.9055pt-366m\left(y_{1}v_{12}\right)\left(120\ln\left(\frac{r_{0}^{2}\bar{q}}{\ell_{0}^{2}}\right)+80\ln\left(\frac{r_{12}}{r_{0}}\right)+101\right)\bigg{]}\,.$
(77e)

## References

  * 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] .
  * Damour  _et al._ [2014] T. Damour, P. Jaranowski, and G. Schäfer, Non-local-in-time action for the fourth post-Newtonian conservative dynamics of two-body systems, Phys. Rev. D 89, 064058 (2014), arXiv:1401.4548 [gr-qc] .
  * Damour  _et al._ [2016] T. Damour, P. Jaranowski, and G. Schäfer, On the conservative dynamics of two-body systems at the fourth post-Newtonian approximation of general relativity, Phys. Rev. D 93, 084014 (2016), arXiv:1601.01283 [gr-qc] .
  * Bernard  _et al._ [2016] L. Bernard, L. Blanchet, A. Bohé, G. Faye, and S. Marsat, Fokker action of non-spinning compact binaries at the fourth post-Newtonian approximation, Phys. Rev. D 93, 084037 (2016), arXiv:1512.02876 [gr-qc] .
  * Bernard  _et al._ [2017] L. Bernard, L. Blanchet, A. Bohé, G. Faye, and S. Marsat, Dimensional regularization of the IR divergences in the Fokker action of point-particle binaries at the fourth post-Newtonian order, Phys. Rev. D 96, 104043 (2017), arXiv:1706.08480 [gr-qc] .
  * Marchand  _et al._ [2018] T. Marchand, L. Bernard, L. Blanchet, and G. Faye, Ambiguity-free completion of the equations of motion of compact binary systems at the fourth post-Newtonian order, Phys. Rev. D 97, 044023 (2018), arXiv:1707.09289 [gr-qc] .
  * Galley  _et al._ [2016] C. R. Galley, A. K. Leibovich, R. A. Porto, and A. Ross, Tail effect in gravitational radiation reaction: Time nonlocality and renormalization group evolution, Phys. Rev. D 93, 124010 (2016), arXiv:arXiv:1511.07379 [gr-qc] [gr-qc] .
  * Porto and Rothstein [2017] R. Porto and I. Rothstein, On the apparent ambiguities in the post-Newtonian expansion for binary systems, Phys. Rev. D 96, 024061 (2017), arXiv:1703.06433 [gr-qc] .
  * Foffa and Sturani [2019] S. Foffa and R. Sturani, Conservative dynamics of binary systems to fourth post-Newtonian order in the EFT approach I: Regularized Lagrangian, Phys. Rev. D 100, 024047 (2019), arXiv:1903.05113 [gr-qc] .
  * Foffa  _et al._ [2019] S. Foffa, R. Porto, I. Rothstein, and R. Sturani, Conservative dynamics of binary systems to fourth post-Newtonian order in the EFT approach II: Renormalized Lagrangian, Phys. Rev. D 100, 024048 (2019), arXiv:1903.05118 [gr-qc] .
  * Blümlein  _et al._ [2020] J. Blümlein, A. Maier, P. Marquard, and G. Schäfer, Fourth post-Newtonian hamiltonian dynamics oftwo-body systems from an effective field theory approach, 2003.01692 [hep-th] (2020), arXiv preprint.
  * Blümlein  _et al._ [2021] J. Blümlein, A. Maier, P. Marquard, and G. Schäfer, The fifth-order post-Newtonian Hamiltonian dynamics of two-body systems from an effective field theory approach: Potential contributions, Nuclear Physics B 965, 115352 (2021), arXiv:2021.13672 [gr-qc] .
  * Bini  _et al._ [2020] D. Bini, T. Damour, and A. Geralico, Sixth post-Newtonian local-in-time dynamics of binary systems, Physical Review D 102, 024061 (2020).
  * Bern  _et al._ [2019] Z. Bern, C. Cheung, R. Roiban, C.-H. Shen, M. P. Solon, and M. Zeng, Scattering Amplitudes and the Conservative Hamiltonian for Binary Systems at Third Post-Minkowskian Order, Phys. Rev. Lett. 122, 201603 (2019), arXiv:1901.04424 [hep-th] .
  * Blanchet  _et al._ [1995] 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 (1995), gr-qc/9501029 .
  * 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 .
  * 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 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._ [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 .
  * Blanchet [1998a] L. Blanchet, Gravitational-wave tails of tails, Class. Quant. Grav. 15, 113 (1998a), gr-qc/9710038 .
  * Blanchet  _et al._ [2002] L. Blanchet, G. Faye, B. R. Iyer, and B. Joguet, Gravitational-wave inspiral of compact binary systems to 7/2 post-Newtonian order, Phys. Rev. D 65, 061501(R) (2002), erratum Phys. Rev. D, 71:129902(E), 2005, gr-qc/0105099 .
  * 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] .
  * 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] .
  * 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] .
  * 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] .
  * 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] .
  * Marchand  _et al._ [2016] T. Marchand, L. Blanchet, and G. Faye, Gravitational-wave tail effects to quartic non-linear order, Class. Quant. Grav. 33, 244003 (2016), arXiv:1607.07601 [gr-qc] .
  * 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] .
  * Hadamard [1932] J. Hadamard,  _Le problème de Cauchy et les équations aux dérivées partielles linéaires hyperboliques_ (Hermann, Paris, 1932).
  * Blanchet and Damour [1988] L. Blanchet and T. Damour, Tail-transported temporal correlations in the dynamics of a gravitating system, Phys. Rev. D 37, 1410 (1988).
  * Larrouturou  _et al._ [2021] 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 (2021), companion paper.
  * Blanchet  _et al._ [2004a] 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 (2004a), gr-qc/0406012 .
  * Blanchet [1998b] L. Blanchet, On the multipole expansion of the gravitational field, Class. Quant. Grav. 15, 1971 (1998b), gr-qc/9801101 .
  * Blanchet  _et al._ [2004b] L. Blanchet, T. Damour, and G. Esposito-Farèse, Dimensional regularization of the third post-Newtonian dynamics of point particles in harmonic coordinates, Phys. Rev. D 69, 124007 (2004b), gr-qc/0311052 .
  * Hartung  _et al._ [2013] J. Hartung, J. Steinhoff, and G. Schäfer, Next-to-next-to-leading order post-Newtonian linear-in-spin binary Hamiltonians, Annalen der Physik 525, 359 (2013), arXiv:1302.6723 [gr-qc] .
  * Goldberger and Ross [2010] W. Goldberger and A. Ross, Gravitational radiative corrections from effective field theory, Phys. Rev. D 81, 124015 (2010), arXiv:0912.4254 [gr-qc] .