# The Quadrupole Moment of Compact Binaries
to the Fourth post-Newtonian Order
II. Dimensional Regularization and Renormalization

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  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  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
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

The regularization and renormalization of the radiative mass-type quadrupole
moment of inspiralling compact binaries (without spins) is investigated at the
fourth post-Newtonian (4PN) approximation of general relativity. As clear from
the conservative 4PN equations of motion, a dimensional regularization has to
be implemented in order to properly treat the non-linear interactions
experienced by gravitational waves during their propagation toward future null
infinity. By implementing such procedure, we show that the poles coming from
the source moment (computed in a companion paper) are exactly cancelled in the
radiative moment, as expected for a physical quantity. We thus define and
obtain a “renormalized” source quadrupole, three-dimensional by nature, which
is an important step towards the computation of the gravitational-wave flux
with 4PN accuracy. Furthermore, we explicitly prove the equivalence between
the dimensional regularization and the previously used Hadamard partie finie
scheme up to the 3PN order.

###### pacs:

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

##  I Introduction

The mass quadrupole moment of compact binary systems is a crucial ingredient
toward the definition of accurate gravitational wave (GW) templates in the
post-Newtonian (PN) approximation. Such PN templates represent the main
technique for detection and analysis of binary neutron star signals in the
current network of detectors on ground, and they are at the basis of effective
phenomenological methods such as EOB (effective-one-body) and IMR (inspiral-
merger-ringdown), able to describe the late inspiral and merger of binary
black holes (see more details in the reviews [1, 2, 3, 4]).

Previous computations of the quadrupole moment in the case of non-spinning
compact binaries achieved the 2PN order [5, 6, 7, 8], then, the 3PN order [9,
10, 11, 12], and the state-of-the-art is currently the 4PN order. A
preliminary calculation was done at the 4PN order [13] where all the terms
were computed using a series of techniques (extensively documented in [13]);
however, as recognized in this paper, the computation was incomplete because
the infra-red (IR) divergences appearing at the 4PN order were regularized
thanks to the Hadamard partie finie regularization instead of the required
dimensional regularization. By contrast, the ultra-violet (UV) divergences,
due to the point-like nature of the source (model of compact objects by point
masses without internal structure) were correctly treated with dimensional
regularization [13].

Recently, we analyzed in Ref. [14] the IR divergences of the (source type)
mass quadrupole moment of compact binaries at 4PN order by means of
dimensional regularization. This study led to the presence of poles in the
dimension, i.e. $\propto\varepsilon^{-1}\equiv(d-3)^{-1}$, arising already
from the 3PN order, and, of course, also contributing at the 4PN order. The
thorough computation of IR poles and the finite terms following the poles was
done in [14].

Furthermore, we also obtained [14] the contribution of propagating GW tails in
the source quadrupole itself. This effect is due to retarded correlations over
arbitrarily large time spans in the dynamics of the source [15]. It is exactly
the analogue of the GW tail effect which adds a non-local piece to the
conservative dynamics (equations of motion and Lagrangian/Hamiltonian) at the
4PN order [16, 17, 18, 19, 20, 21, 22, 23, 24] and beyond [25, 26]. As such,
it induces a non-locality in time in the mass quadrupole moment at the 4PN
order.

The present work is the follow-up of the companion paper [14]. We prove that
the IR poles are cancelled by UV poles coming from the contributions of non-
linear effects, essentially the “tails-of-tails” and the “tails-of-memory”,
arising in the radiative type quadrupole moment, which constitutes the
physical observable at future null infinity.

More precisely, at the 3PN order, the non-linearities are made of the tail-of-
tails due to cubic interaction between two masses and the quadrupole [27], and
the interactions between one mass and two mass dipoles (the latter naturally
vanishes in the center-of-mass frame). Taking into account the pole term and
the finite term beyond the pole, we find that the contributions from
dimensional regularization in the radiative moment exactly cancel out at this
order. This proves that at 3PN order the quadrupole moment can be computed
using the Hadamard regularization for the IR divergences (but the dimensional
regularization for the UV) and, therefore, we totally confirm the previous
calculations of the 3PN quadrupole moment [9, 10, 11, 12].

Furthermore, we prove that the same cancellation holds for the 3PN current
quadrupole moment, thus supporting the recent computation of this moment at
3PN order using the Hadamard regularization for the IR [28], and for the 3PN
mass octupole moment, also confirming the previous calculation done in [29].

Finally, we demonstrate that all the poles cancel in the radiative type
quadrupole moment up to 4PN order. At the 4PN order, in addition to tails-of-
tails, the main non-linear effects to consider in the radiative moment are the
so-called tails-of-memory, which are cubic interations between one mass
monopole and two quadrupole moments. However, at the 4PN order, the net
contribution of the dimensional regularization is not vanishing and, thus, we
determine some finite contributions beyond the poles which are crucial to
include into the final expression of the mass quadrupole moment.

To control the occurrence of poles from the non-linear multipole interactions,
we have to obtain a general expression for the retarded integral of some
extended source in $d$ dimensions up to cubic order, and compute from it the
difference up to order $\mathcal{O}(\varepsilon)$ between the dimensional
regularization or, more precisely, the so-called $B\varepsilon$ regularization
which we systematically employ [19, 20, 14], and the previously used Hadamard
regularization.

The plan of this paper is as follows. We first present the derivations of the
main technical ingredients used for the analysis done in this work: the
generic solution of the $d$-dimensional wave equation in Sec. II, and the
difference induced by the change of regularization scheme in Sec. III. Sec. IV
then sketches the MPM (Multipolar-post-Minkowskian) method used to concretely
implement the computation of this difference. The results at 3PN order are
displayed in Sec.V, together with the proof that, at this order, both
regularization schemes are equivalent. Finally, Sec. VI displays the 4PN
results, introduces our notion of “renormalized mass quadrupole” and presents
its expression on quasi-circular orbits. Appendix A contains the proof that
the generic solution of the $d$-dimensional wave equation correctly reduces to
the known results in the three-dimensional limit.

##  II The retarded integral of a multipolar extended source

In this technical section we provide an explicit expression for the retarded
solution of the wave equation in any $d$ space dimensions, in the case where
the non-compact (extended) source term has a definite multipolarity. This
object is the core of the iteration presented in Sec. IV, that we use to
compute the difference between regularizations in non-linear interactions
entering the gravitational wave propagation. We thus seek for the generic
solution $h(\mathbf{x},t)$ solving the wave equation

$\Box h(\mathbf{x},t)=N(\mathbf{x},t)\,,$ (1)

where $\Box\equiv\Box_{\eta}$ is the flat d’Alembertian operator and
$N(\mathbf{x},t)$ is some (non-compact support) source term. While generally
the solution is provided in the Fourier domain, here we work in the physical
space. The retarded solution of Eq. (1) in the real domain reads (see e.g.
[19])111We pose $c=1$ in Secs. II, III and App. A.

$\displaystyle
h(\mathbf{x},t)=-\frac{\tilde{k}}{4\pi}\int_{1}^{+\infty}\\!\\!\mathrm{d}z\,\gamma_{\frac{1-d}{2}}(z)\int\mathrm{d}^{d}\mathbf{x}^{\prime}\,\frac{N(\mathbf{x}^{\prime},t-z|\mathbf{x}-\mathbf{x}^{\prime}|)}{|\mathbf{x}-\mathbf{x}^{\prime}|^{d-2}}\,,$
(2)

where $\tilde{k}=\Gamma(\frac{d}{2}-1)/\pi^{\frac{d}{2}-1}$ is the constant
entering the Green’s function of the Laplace operator, i.e.,
$\Delta\bigl{(}\tilde{k}\,r^{2-d}\bigr{)}=-4\pi\delta(\mathbf{x})$, $\Gamma$
is the Eulerian function, and we have posed

$\displaystyle\gamma_{\frac{1-d}{2}}(z)=\frac{2\sqrt{\pi}}{\Gamma(\frac{3-d}{2})\Gamma(\frac{d}{2}-1)}\,\big{(}z^{2}-1\bigr{)}^{\frac{1-d}{2}}\,,$
(3)

whith chosen normalization such that
$\int_{1}^{+\infty}\mathrm{d}z\,\gamma_{\frac{1-d}{2}}(z)=1$ and
$\gamma_{-1}(z)=\delta(z-1)$.

We split the spatial integration over the source point $\mathbf{x}^{\prime}$
in (2) into an inner domain corresponding to $r^{\prime}<r$, where we denote
$r^{\prime}\equiv|\mathbf{x}^{\prime}|$ and $r\equiv|\mathbf{x}|$ with
$\mathbf{x}$ being the field point, and the outer domain for which
$r^{\prime}>r$. Consider for instance the integration over the inner domain
$r^{\prime}<r$, say

$\displaystyle
h^{<}\equiv-\frac{\tilde{k}}{4\pi}\int_{1}^{+\infty}\\!\\!\mathrm{d}z\,\gamma_{\frac{1-d}{2}}(z)\int_{r^{\prime}<r}\mathrm{d}^{d}\mathbf{x}^{\prime}\,\frac{N(\mathbf{x}^{\prime},t-z|\mathbf{x}-\mathbf{x}^{\prime}|)}{|\mathbf{x}-\mathbf{x}^{\prime}|^{d-2}}\,.$
(4)

We obtain the multipole expansion in this domain by applying the formal Taylor
expansion when $r^{\prime}\to 0$. The symmetric-trace-free (STF) form of that
expansion reads

$\displaystyle\frac{N(\mathbf{x}^{\prime},t-z|\mathbf{x}-\mathbf{x}^{\prime}|)}{|\mathbf{x}-\mathbf{x}^{\prime}|^{d-2}}=\sum_{m=0}^{+\infty}\frac{(-)^{m}}{m!}\sum_{j=0}^{+\infty}\frac{\Gamma\left(\frac{d}{2}+m\right)}{\Gamma\left(\frac{d}{2}+m+j\right)}\frac{\hat{n}^{\prime}_{M}{r^{\prime}}^{2j+m}}{2^{2j}j!}\hat{\partial}_{M}\Delta^{j}\\!\left[\frac{N(\mathbf{x}^{\prime},t-zr)}{r^{d-2}}\right]\,,$
(5)

where $M$ denotes a running multi-index $j_{1}\cdots j_{m}$ with $m$ indices,
$\hat{n}^{\prime}_{M}\equiv\text{STF}(n^{\prime}_{j_{1}}\cdots
n^{\prime}_{j_{m}})$ is the STF product of $m$ unit vectors
$n^{\prime}_{j}\equiv x^{\prime}_{j}/r^{\prime}$, and we denote
$\hat{\partial}_{M}\equiv\text{STF}(\partial_{j_{1}}\cdots\partial_{j_{m}})$
where $\partial_{j}\equiv\partial/\partial x^{j}$. Note that
$\hat{\partial}_{M}$ and $\Delta^{j}$ in (5) act on $\mathbf{x}$, but
$\mathbf{x}^{\prime}$ is just spectator.

We consider the above general formula in the case where the source term
$N(\mathbf{x},t)$ itself has a definite multipolarity $\ell$, i.e. is of the
type

$N(\mathbf{x},t)\equiv\hat{n}_{L}\,N(r,t)\,,$ (6)

where $L$ denotes the multi-index $i_{1}\cdots i_{\ell}$ and
$\hat{n}_{L}\equiv\text{STF}(n_{i_{1}}\cdots n_{i_{\ell}})$ with $n_{i}\equiv
x_{i}/r$. Clearly, after summing up all the multipolar pieces, while assuming
the convergence of the multipolar series, there is no restriction on the
generality of the solution. For a source term of the type (6) the angular
integration in the volume element
$\mathrm{d}^{d}\mathbf{x}^{\prime}\equiv{r^{\prime}}^{d-1}\mathrm{d}{r^{\prime}}\mathrm{d}\Omega^{\prime}_{d-1}$
can be evaluated in closed form using

$\int\mathrm{d}\Omega^{\prime}_{d-1}\hat{n}^{\prime}_{L}\,\hat{n}^{\prime}_{M}=\frac{\ell!}{2^{\ell-1}}\,\frac{\pi^{\frac{d}{2}}}{\Gamma\left(\frac{d}{2}+\ell\right)}\,\delta_{LM}\,,$
(7)

where $\delta_{LM}\equiv\delta_{\ell
m}\delta_{(i_{1}\underline{j}{}_{1}}\cdots\delta_{i_{\ell})j_{\ell}}$
(symmetrization over the indices $L=i_{1}\cdots i_{\ell}$); we recall the
volume of the sphere $\Omega_{d-1}=2\pi^{\frac{d}{2}}/\Gamma(\frac{d}{2})$.
After angular integration, the part of the integral over the domain
$r^{\prime}<r$, let us now call it $h_{L}^{<}$ since it depends on the multi-
index $L$, becomes

$\displaystyle
h_{L}^{<}=\frac{1}{(-2)^{\ell+1}}\sum_{j=0}^{+\infty}\frac{\Gamma\left(\frac{d}{2}-1\right)}{\Gamma\left(\frac{d}{2}+\ell+j\right)}\frac{1}{2^{2j}j!}\int_{1}^{+\infty}\\!\\!\mathrm{d}z\,\gamma_{\frac{1-d}{2}}(z)\int_{0}^{r}\mathrm{d}r^{\prime}r^{\prime
B+2j+\ell+d-1}\hat{\partial}_{L}\\!\left[\frac{N^{(2j)}(r^{\prime},t-zr)}{r^{d-2}}\right]\,.$
(8)

Notice that in Eq. (8) the action of the iterated Laplacian $\Delta^{j}$ has
reduced to $2j$ time derivatives of the source term, as indicated by the
superscript in $N^{(2j)}$. Furthermore we dispose of the following useful
lemma (already employed in [19] but not given there)

$\displaystyle\int_{1}^{+\infty}\\!\\!\mathrm{d}z\,\gamma_{\frac{1-d}{2}}(z)\,\hat{\partial}_{L}\\!\left(\frac{F(t-zr)}{r^{d-2}}\right)=(-2)^{\ell}\frac{\Gamma\left(\frac{d}{2}+\ell-1\right)}{\Gamma\left(\frac{d}{2}-1\right)}\frac{\hat{n}_{L}}{r^{\ell+d-2}}\int_{1}^{+\infty}\\!\\!\mathrm{d}z\,\gamma_{\frac{1-d}{2}-\ell}(z)F(t-zr)\,,$
(9)

which permits to rewrite (8) into the interesting alternative form

$\displaystyle
h_{L}^{<}=-\frac{1}{2}\sum_{j=0}^{+\infty}\frac{\Gamma\left(\frac{d}{2}+\ell-1\right)}{\Gamma\left(\frac{d}{2}+\ell+j\right)}\frac{1}{2^{2j}j!}\frac{\hat{n}_{L}}{r^{\ell+d-2}}\int_{1}^{+\infty}\\!\\!\mathrm{d}z\,\gamma_{\frac{1-d}{2}-\ell}(z)\int_{0}^{r}\mathrm{d}{r^{\prime}}{r^{\prime}}^{B+\ell+2j+d-1}N^{(2j)}({r^{\prime}},t-zr)\,.$
(10)

Very important in our approach, we use the “$B\varepsilon$” regularization
scheme [19, 20, 14], which is a variant of the dimensional regularization
(with $\varepsilon=d-3$), in which a regularization factor $(r/r_{0})^{B}$ is
inserted first in order to protect against the usual divergence of the
multipole expansion when $r\to 0$, where $r_{0}$ is an arbitrary constant.
Such factor $(r/r_{0})^{B}$ is the same as the one used in many previous works
in 3 dimensions, and here it is also used on the “top” of dimensional
regularization. This procedure is not an ad hoc additional regularization
scheme, but comes from the matching between near zone and exterior zone and is
crucial for the proper definition of the MPM algorithm in three dimensions
[30]. The calculations of the tail sector in the conservative equations of
motion [19, 20] have indicated that the $(r/r_{0})^{B}$ factor should be kept
even in $d$ dimensions, at least in intermediate calculations. Therefore, in
Eqs. (8) and (10) we have replaced the source term by $N(r,t)\longrightarrow
r^{B}N(r,t)$. For convenience we pose $r_{0}=1$ but we will restore the
constant $r_{0}$ (as well as the dimensional regularization length scale
$\ell_{0}$) in our final results. Following the $B\varepsilon$ regularization
a finite part (or Partie Finie PF [31]) at $B=0$ is always understood before
considering the limit $\varepsilon\to 0$. We expect that the finite part at
$B=0$ just reduces to a finite limit when $B\to 0$.

The result (8)–(10) is still in the form of a formal series, and we now show
how this series can be re-summed so that we obtain the corresponding “exact”
result. To this end we use the expression of the near-zone (or PN) expansion
of an antisymmetric multipolar solution of the wave equation, i.e. of the type
retarded minus advanced. In the monopolar case ($\ell=0$) this expansion is
given by Eq. (A13) of [19], and we see that in $d$ dimensions it involves a
non-local integral. See also the discussion in Sec. IV and the PN expansion
given by (42)–(43). Let us rewrite Eq. (A13) of [19] in the form (with
$\varepsilon=d-3$)

$\displaystyle\frac{\tilde{k}}{r^{d-2}}\int_{1}^{+\infty}\\!\\!\mathrm{d}y\,\gamma_{\frac{1-d}{2}}(y)\Bigl{[}\widehat{F}_{\varepsilon}(t-yr)-\widehat{F}_{\varepsilon}(t+yr)\Bigr{]}=-\frac{\pi^{2-\frac{d}{2}}}{2^{d-3}}\,\sum_{j=0}^{+\infty}\frac{r^{2j}}{2^{2j}j!}\,\frac{F^{(2j+1)}(t)}{\Gamma\left(\frac{d}{2}+j\right)}\,.$
(11)

For any smooth function $F(t)$ with compact support on $\mathbb{R}$, we have
introduced its “$\varepsilon$-transform” $\widehat{F}_{\varepsilon}(t)$ (also
smooth and with compact support) to be the function defined by analytic
continuation in $\varepsilon$ as

$\displaystyle\widehat{F}_{\varepsilon}(t)=\frac{1}{2\cos(\frac{\pi\varepsilon}{2})\Gamma(\varepsilon)}\int_{-\infty}^{+\infty}\mathrm{d}\tau\,|\tau|^{\varepsilon-1}\,F(t+\tau)\,.$
(12)

This $\varepsilon$-transform satisfies two remarkable properties.

  * •

The main one is that its inverse transform is given by the same expression but
corresponding to the parameter $-\varepsilon$, thus

$\displaystyle
F(t)=\frac{1}{2\cos(\frac{\pi\varepsilon}{2})\Gamma(-\varepsilon)}\int_{-\infty}^{+\infty}\mathrm{d}\tau\,|\tau|^{-\varepsilon-1}\,\widehat{F}_{\varepsilon}(t+\tau)\,.$
(13)

This nice property can easily be proved by going to the Fourier domain, using
e.g. (3.17)–(3.18) in [19].222Hence we can check that Eq. (11) is indeed
equivalent to Eq. (A13) of [19].

  * •

The second important property of the $\varepsilon$-transform is that it
reduces to the identity in the limit when $\varepsilon\to 0$:

$\displaystyle\lim_{\varepsilon\to 0}\widehat{F}_{\varepsilon}(t)=F(t)\,.$
(14)

However, the transform (12)–(13) is valid by analytic continuation for any
$\varepsilon\in\mathbb{C}$.

Next we obtain the expansion for any multipolar antisymmetric wave by applying
on Eq. (11) the STF multi-spatial derivative operator $\hat{\partial}_{L}$,
and we obtain333Using the fact that
$\hat{\partial}_{L}r^{2j}=\frac{2^{\ell}j!}{(j-\ell)!}\,\hat{n}_{L}\,r^{2j-\ell}$
when $j\geqslant\ell$ and is zero when $0\leqslant j\leqslant\ell-1$.

$\displaystyle\tilde{k}\int_{1}^{+\infty}\\!\\!\mathrm{d}y\,\gamma_{\frac{1-d}{2}}(y)\,\hat{\partial}_{L}\Biggl{[}\frac{\widehat{F}_{\varepsilon}(t-yr)-\widehat{F}_{\varepsilon}(t+yr)}{r^{d-2}}\Biggr{]}=-\frac{\pi^{2-\frac{d}{2}}}{2^{\ell+d-3}}\,\sum_{j=0}^{+\infty}\frac{\hat{x}_{L}r^{2j}}{2^{2j}j!}\,\frac{F^{(2j+2\ell+1)}(t)}{\Gamma\left(\frac{d}{2}+\ell+j\right)}\,.$
(15)

Thanks to the relation (11) with the substitution $d\rightarrow d+2\ell$, one
can finally re-express the result for $h_{L}^{<}$ as given by (8) in the
following closed-analytic “exact” form

$\displaystyle h_{L}^{<}$
$\displaystyle=\tilde{k}\,\Gamma_{\ell}\int_{1}^{+\infty}\mathrm{d}y\,\gamma_{\frac{1-d}{2}}(y)\int_{1}^{+\infty}\mathrm{d}z\,\gamma_{\frac{1-d}{2}-\ell}(z)$
(16)
$\displaystyle\qquad\qquad\times\int_{0}^{r}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1}\hat{\partial}_{L}\\!\left[\frac{\widehat{N}_{\varepsilon}^{(-2\ell-1)}\bigl{(}{r^{\prime}},t-yr-z{r^{\prime}}\bigr{)}-\widehat{N}_{\varepsilon}^{(-2\ell-1)}\bigl{(}{r^{\prime}},t-yr+z{r^{\prime}}\bigr{)}}{r^{d-2}}\right]\,,$

where the superscript $(-2\ell-1)$ indicates $2\ell+1$ time anti-
derivatives444As usual in our formalism, we assume stationarity in the remote
past, and anti-derivatives are defined to be those which vanish when
$t\to-\infty$. and we have posed $\Gamma_{\ell}\equiv
2^{\ell+d-4}\pi^{\frac{d}{2}-2}\,\Gamma(\frac{d}{2}+\ell-1)$, such that
$\lim_{\varepsilon\to 0}\Gamma_{\ell}=\frac{(2\ell-1)!!}{2^{\ell-1}}$. We can
essentially repeat the same steps of the previous calculation with only some
simple adaptations, for the outer part of the solution $h_{L}^{>}$
corresponding to $r^{\prime}>r$, and we obtain

$\displaystyle h_{L}^{>}$
$\displaystyle=\tilde{k}\,\Gamma_{\ell}\int_{1}^{+\infty}\mathrm{d}y\,\gamma_{\frac{1-d}{2}}(y)\int_{1}^{+\infty}\mathrm{d}z\,\gamma_{\frac{1-d}{2}-\ell}(z)$
(17)
$\displaystyle\qquad\qquad\times\int_{r}^{+\infty}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1}\hat{\partial}_{L}\\!\left[\frac{\widehat{N}_{\varepsilon}^{(-2\ell-1)}\bigl{(}{r^{\prime}},t-yr-z{r^{\prime}}\bigr{)}-\widehat{N}_{\varepsilon}^{(-2\ell-1)}\bigl{(}{r^{\prime}},t+yr-z{r^{\prime}}\bigr{)}}{r^{d-2}}\right]\,.$

A finite part PF when $B\to 0$ is to be understood in front of both Eqs.
(16)–(17). Recalling the definition (12) of the $\varepsilon$-transform we see
that (16)–(17) involve three non-local integrals over $y$, $z$ and $\tau$,
that are specifically due to $d$ dimensions and reduce to trivial local
results in 3 dimensions. Notice that in the complete solution the first terms
in (16) and (17) merge to form a single integral from 0 to $+\infty$. Thus we
may write the complete solution $h_{L}=h_{L}^{<}+h_{L}^{>}$ as

$\displaystyle
h_{L}=\tilde{k}\,\Gamma^{\prime}_{\ell}\,\frac{\hat{n}_{L}}{r^{d+\ell-2}}\int_{1}^{+\infty}\mathrm{d}y\,\gamma_{\frac{1-d}{2}-\ell}(y)\int_{1}^{+\infty}\mathrm{d}z\,\gamma_{\frac{1-d}{2}-\ell}(z)$
(18)
$\displaystyle\qquad\times\Biggl{\\{}\int_{0}^{+\infty}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1}\widehat{N}_{\varepsilon}^{(-2\ell-1)}\bigl{(}{r^{\prime}},t-yr-z{r^{\prime}}\bigr{)}$
$\displaystyle\qquad-\int_{0}^{r}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1}\widehat{N}_{\varepsilon}^{(-2\ell-1)}\bigl{(}{r^{\prime}},t-yr+z{r^{\prime}}\bigr{)}-\int_{r}^{+\infty}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1}\widehat{N}_{\varepsilon}^{(-2\ell-1)}\bigl{(}{r^{\prime}},t+yr-z{r^{\prime}}\bigr{)}\Biggr{\\}}\,,$

where we pose
$\Gamma^{\prime}_{\ell}\equiv(-)^{\ell}2^{2\ell+d-4}\pi^{\frac{d}{2}-2}\,[\Gamma(\frac{d}{2}+\ell-1)]^{2}/\Gamma(\frac{d}{2}-1)$,
and we have chosen to apply the lemma (9) on the $y$ integration in order to
express the result. Note that the first term in (18) is an homogeneous
retarded solution of the wave equation in $d$ dimensions. In the Appendix A,
we explicitly check the limit of our result (18) in three dimensions, and we
gladly recover a formula derived previously in [32].

We end up this section by two remarks. For the computation of the conservative
Lagrangian and equations of motion at 4PN order, and handling the conservative
part of the tail effect at that order, it was shown in Ref. [19] that the
relevant quantity to be considered is $h_{L}^{0}$, defined by $h_{L}^{>}$ as
given by (17) but in which the radial integral from $r$ to $+\infty$ is
replaced by the integral from 0 to $+\infty$:555This solution is given in
expanded form by Eq. (3.20) in [19].

$\displaystyle h_{L}^{0}$
$\displaystyle=\tilde{k}\,\Gamma_{\ell}\int_{1}^{+\infty}\mathrm{d}y\,\gamma_{\frac{1-d}{2}}(y)\int_{1}^{+\infty}\mathrm{d}z\,\gamma_{\frac{1-d}{2}-\ell}(z)$
(19)
$\displaystyle\qquad\qquad\times\int_{0}^{+\infty}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1}\hat{\partial}_{L}\\!\left[\frac{\widehat{N}_{\varepsilon}^{(-2\ell-1)}\bigl{(}{r^{\prime}},t-yr-z{r^{\prime}}\bigr{)}-\widehat{N}_{\varepsilon}^{(-2\ell-1)}\bigl{(}{r^{\prime}},t+yr-z{r^{\prime}}\bigr{)}}{r^{d-2}}\right]\,.$

Indeed this $h_{L}^{0}$ is an homogeneous solution of the wave equation, $\Box
h_{L}^{0}=0$, and furthermore is regular when $r\to 0$, i.e. is of the type
retarded minus advanced. Therefore it was proved that this solution can be
extended inside the source by matching, and it permitted the computation of
the 4PN tail effect in the conservative equations of motion (with result in
perfect agreement with the one derived by the effective-field-theory approach
[22]).

Similarly we will find in the next section that in order to compute the
contribution from tails in the far zone, it is sufficient to define
$h_{L}^{\infty}$ by $h_{L}^{<}$ as given by (16), but in which the radial
integral from 0 to $r$ is replaced by the integral from 0 to $+\infty$:

$\displaystyle h_{L}^{\infty}$
$\displaystyle=\tilde{k}\,\Gamma_{\ell}\int_{1}^{+\infty}\mathrm{d}y\,\gamma_{\frac{1-d}{2}}(y)\int_{1}^{+\infty}\mathrm{d}z\,\gamma_{\frac{1-d}{2}-\ell}(z)$
(20)
$\displaystyle\qquad\qquad\times\int_{0}^{+\infty}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1}\hat{\partial}_{L}\\!\left[\frac{\widehat{N}_{\varepsilon}^{(-2\ell-1)}\bigl{(}{r^{\prime}},t-yr-z{r^{\prime}}\bigr{)}-\widehat{N}_{\varepsilon}^{(-2\ell-1)}\bigl{(}{r^{\prime}},t-yr+z{r^{\prime}}\bigr{)}}{r^{d-2}}\right]\,.$

Again this is an homogeneous solution of the wave equation, $\Box
h_{L}^{\infty}=0$, but this time a retarded solution, therefore regular at
infinity, when $r\to+\infty$. This is that solution $h_{L}^{\infty}$ which
will be responsible for the contribution of far zone tails in the difference
between the dimensional and Hadamard regularizations, as we will show in the
next section.

##  III The difference between two regularization schemes

Up to vanishingly small terms when $\varepsilon\to 0$, the difference between
the dimensional regularization and the Hadamard one comes only from the bound
${r^{\prime}}=0$ in the general formula for the retarded integral, i.e. it
comes only from the part of the integral $h_{L}^{<}$ given by (16). This is
clear because only the bound ${r^{\prime}}=0$ can give a pole in dimensional
regularization, and any integration over an interval excluding
${r^{\prime}}=0$ will have a finite limit when $\varepsilon\to 0$, and so
cancels out between the two regularization in the limit $\varepsilon\to 0$.
This further means that for the computation of the difference we can replace
the integral from 0 to $r$ in (16) by the integral from 0 up to $+\infty$ (or
any finite constant), i.e. we can just consider the homogeneous retarded
solution $h_{L}^{\infty}$ given by (20).

Thus our statement is that modulo terms of order $\varepsilon$ the difference
is

$\displaystyle\mathcal{D}h_{L}=\mathcal{D}h^{\infty}_{L}+\mathcal{O}(\varepsilon)\,,$
(21)

where the difference operator $\mathcal{D}$ is formally defined as the
commutator between the three-dimensional limit (viz. $\lim_{3d}$) and the PF
operator:

$\mathcal{D}A\equiv\left[\lim_{3d}\,,\mathop{\mathrm{PF}}_{B=0}\right]\,A=\lim_{3d}\left\\{\mathop{\mathrm{PF}}_{B=0}\,A\right\\}-\mathop{\mathrm{PF}}_{B=0}\left\\{\lim_{3d}\,A\right\\}=A^{B\varepsilon}-A^{\text{Had}}\,.$
(22)

Here, by “3d” limit, we really mean the expansion $\varepsilon\to 0$, while
crucially keeping all the divergent poles $\varepsilon^{-1}$, and neglecting
terms $\mathcal{O}(\varepsilon)$. By contrast, $\mathop{\mathrm{PF}}_{B=0}$
discards the poles $B^{-1}$, so the difference $\mathcal{D}A$ is actually a
function of the poles $\varepsilon^{-1}$ but does not depend on $B$.

Since $h^{\infty}_{L}$ is a retarded homogeneous solution of the wave equation
in $d$ dimensions, it can be written in multipolar form as

$\displaystyle
h_{L}^{\infty}=\hat{\partial}_{L}\,\biggl{[}\frac{\tilde{k}}{r^{d-2}}\,\int_{1}^{+\infty}\mathrm{d}y\,\gamma_{\frac{1-d}{2}}(y)\,H(t-yr)\biggr{]}\,,$
(23)

where we pose

$\displaystyle H(t)$
$\displaystyle=\Gamma_{\ell}\int_{1}^{+\infty}\mathrm{d}z\,\gamma_{\frac{1-d}{2}-\ell}(z)\int_{0}^{+\infty}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1}\,\left[\widehat{N}_{\varepsilon}^{(-2\ell-1)}\bigl{(}{r^{\prime}},t-z{r^{\prime}}\bigr{)}-\widehat{N}_{\varepsilon}^{(-2\ell-1)}\bigl{(}{r^{\prime}},t+z{r^{\prime}}\bigr{)}\right]\,.$
(24)

We recall that the “$\varepsilon$-transform” operation is defined by (12).

As we consider only the limit ${r^{\prime}}\to 0$ we insert into (24) the
near-zone expansion of the source term $N(r,t)$ in $d$ dimensions, which will
turn out to be always of the generic type (as shown in Sec. IV)

$\displaystyle N(r,t)=\sum_{p,q}r^{-p-q\varepsilon}f_{p,q}(t)\,,$ (25)

where the coefficients $f_{p,q}(t)$ may be complicated non-local integrals
involving poles, and the singularity is not essential, i.e. the maximal
divergence when $r\to 0$ is finite, say $p\leqslant p_{0}$. Furthermore, we
have explicitly checked that $q$ takes values in a finite range, say
$q_{0}\leqslant q\leqslant q_{1}$. The $d$-dimensional coefficients
$f_{p,q}(t)$ match the corresponding three-dimensional coefficients when
$\varepsilon\to 0$, in the sense that

$\displaystyle\lim_{\varepsilon\to
0}\,\sum_{q}r^{-q\varepsilon}\,f_{p,q}(t)=f^{(3d)}_{p}(t)+f^{(3d,\text{ln})}_{p}(t)\ln\left(\frac{r}{\ell_{0}}\right)\,.$
(26)

The three-dimensional logarithm comes from possible poles in $f_{p,q}$, see
also Eq. (3.9) of [14], but in our practical calculations the limit (26) is
finite. Furthermore, in our case, the source does not bear double poles so
that only simple logarithms show up in the three-dimensional limit. By
inserting (25) into (24) we obtain

$\displaystyle H(t)\\!$
$\displaystyle=\\!\sum_{p,q}\Gamma_{\ell}\\!\int_{1}^{+\infty}\\!\\!\mathrm{d}z\,\gamma_{-1-\ell-\frac{\varepsilon}{2}}(z)\\!\int_{0}^{+\infty}\\!\\!\mathrm{d}{r^{\prime}}{r^{\prime}}^{B-\ell+1-p-q\varepsilon}\,\biggl{[}{}_{\varepsilon}\widehat{f}_{p,q}^{(-2\ell-1)}\bigl{(}t-z{r^{\prime}}\bigr{)}\\!-\\!{}_{\varepsilon}\widehat{f}_{p,q}^{(-2\ell-1)}\bigl{(}t+z{r^{\prime}}\bigr{)}\biggr{]}\,,$
(27)

and by posing $\rho=zr^{\prime}$, we find that the integral over $z$
factorizes out and can nicely be computed in closed analytic form.666For that
we use
$\int_{1}^{+\infty}\mathrm{d}z\,z^{\alpha}\,\gamma_{s}(z)=\sqrt{\pi}\frac{\Gamma(-s-\frac{\alpha+1}{2})}{\Gamma(-s-\frac{1}{2})\Gamma(\frac{1-\alpha}{2})}.$
Finally, we are able to perform a series of integrations by parts (with all
integrated terms at infinity vanishing) to arrive at

$\displaystyle H(t)\\!$
$\displaystyle=\\!-2^{\ell-1+\varepsilon}\pi^{\frac{\varepsilon}{2}}\sum_{p,q}\int_{0}^{+\infty}\\!\\!\mathrm{d}\rho\,Q_{\ell}^{p,q}\bigl{(}B,\varepsilon,\rho\bigr{)}\,\biggl{[}(-)^{\ell+p}{}_{\varepsilon}\widehat{f}_{p,q}^{(-\ell-2+p)}\bigl{(}t-\rho\bigr{)}+{}_{\varepsilon}\widehat{f}_{p,q}^{(-\ell-2+p)}\bigl{(}t+\rho\bigr{)}\biggr{]}\,,$
(28)

in which we have inserted $\Gamma_{\ell}\equiv
2^{\ell+d-4}\pi^{\frac{d}{2}-2}\Gamma(\frac{d}{2}+\ell-1)$, and posed

$\displaystyle
Q_{\ell}^{p,q}\bigl{(}B,\varepsilon,\rho\bigr{)}\equiv\frac{\Gamma\left(\frac{B+\ell+3-p-(q-1)\varepsilon}{2}\right)\Gamma\left(\text{$-B+q\varepsilon$}\right)}{\Gamma\left(\frac{B-\ell+3-p-q\varepsilon}{2}\right)\Gamma\left(-B+\ell-1+p+q\varepsilon\right)}\,\rho^{B-q\varepsilon}\,.$
(29)

With this result we straightforwardly obtain the difference $\mathcal{D}H$
corresponding to the calculation with the $B\varepsilon$ method minus the
Hadamard regularization which is just the finite part at $B=0$ in three
dimensions. Of course we use the link between $d$ dimensions and 3 dimensions
saying that the functions $f_{p,q}(t)$ match their counterparts in 3
dimensions, see (26). From (28), it is clear that the difference reads as

$\displaystyle\mathcal{D}H(t)$
$\displaystyle=-2^{\ell-1+\varepsilon}\pi^{\frac{\varepsilon}{2}}\sum_{p,q}\int_{0}^{+\infty}\\!\mathrm{d}\rho\,\mathcal{D}Q_{\ell}^{p,q}\bigl{(}\varepsilon,\rho\bigr{)}\,\biggl{[}(-)^{\ell+p}{}_{\varepsilon}\widehat{f}_{p,q}^{(-\ell-2+p)}\bigl{(}t-\rho\bigr{)}+{}_{\varepsilon}\widehat{f}_{p,q}^{(-\ell-2+p)}\bigl{(}t+\rho\bigr{)}\biggr{]}\,,$
(30)

where we recall that
$\mathcal{D}Q_{\ell}^{p,q}\bigl{(}\varepsilon,\rho\bigr{)}$ denotes the
difference of the quantities (29) evaluated with the two regularizations, and
given explicitly by (22) as

$\displaystyle\mathcal{D}Q_{\ell}^{p,q}\bigl{(}\varepsilon,\rho\bigr{)}\equiv\lim_{3d}\left\\{\mathop{\mathrm{PF}}_{B=0}\bigl{[}Q_{\ell}^{p,q}\bigl{(}B,\varepsilon,\rho\bigr{)}\bigr{]}\right\\}-\mathop{\mathrm{PF}}_{B=0}\bigl{[}Q_{\ell}^{p,q}\bigl{(}B,0,\rho\bigr{)}\bigr{]}\,.$
(31)

Finally it is a straightforward (although fairly tedious) matter to evaluate
explicitly the difference (31) neglecting vanishing terms in the limit
$\varepsilon\to 0$. We find that there are two cases where the difference is
not zero:

  1. 1.

When $p=\ell+2i+3$ where $i\in\mathbb{N}$ and $q\neq 1$, in which case we have

$\displaystyle\mathcal{D}Q_{\ell}^{p,q}=\frac{(-)^{\ell}(\ell+i)!}{i!(2\ell+2i+1)!}\frac{1}{q-1}\biggl{[}\frac{1}{\varepsilon}-\sum_{k=0}^{\ell+i}\frac{1}{2k+1}-\frac{1}{2}\gamma_{\text{E}}-\ln\left(\frac{\rho\,r_{0}^{q-1}}{\ell_{0}^{q}}\right)\biggr{]}+\mathcal{O}\left(\varepsilon\right)\,.$
(32)

Note the important point that the exclusion of the $q=1$ is not an artificial
requirement to ensure that the formula is well-defined. Indeed, when $q=1$ the
difference $\mathcal{D}Q_{\ell}^{p,1}$ is exactly vanishing; this point will
be of uttermost importance in Sec. IV.

  2. 2.

When $p=-\ell+2i+2$ where $i\in\mathbb{N}$ and $q\neq 0$, where we have

$\displaystyle\mathcal{D}Q_{\ell}^{p,q}=(-)^{i}\frac{(2\ell-2i-1)!!}{2^{\ell+i}i!}\frac{1}{q}\biggl{[}\frac{1}{\varepsilon}+\sum_{k=0}^{|\ell-i|-1}\frac{1}{2k+1}-\frac{1}{2}\gamma_{\text{E}}+\ln\left(\frac{\ell_{0}^{q}}{2r_{0}^{q}}\right)\biggr{]}+\mathcal{O}\left(\varepsilon\right)\,.$
(33)

with the convention that $(-2n-1)!!(2n+1)!!\equiv(-)^{n}(2n+1)$ for any
integer $n\in\mathbb{N}$.

All other cases give zero (no difference). In case 2, since (33) does not
depend on $\rho$ the difference is local, and turns out to be vanishing.
Indeed, due to the even nature of $p+\ell$, the two integrals corresponding to
the two terms in the square brackets of (30) exactly cancel out. On the other
hand, the case 1 seems to be non-local as there is a $\ln\rho$ remaining and
therefore a non-local integral associated. However we find that this non-local
integral exactly cancels out the one present in the coefficients
${}_{\varepsilon}\widehat{f}_{p,q}$ up to order $\mathcal{O}(\varepsilon)$,
cf. Eq. (12) of the $\varepsilon$-transform, so that the final result does not
contain a remaining integral over $\rho$. Note that we have not considered the
$\mathcal{O}(\varepsilon)$ term, as we have explicitly verified that none of
the contributing $f_{p,q}$ (i.e. such that $p-\ell-3\in 2\mathbb{N}$) bears
poles. If it were the case, the difference would turn out to be non-local.
Fortunately this is not our case and thus we are able to write down our final
result as (with $\bar{q}=4\pi\mathrm{e}^{\gamma_{\text{E}}}$)

$\displaystyle\mathcal{D}H(t)$
$\displaystyle=(-2)^{\ell}\sum_{i=0}^{+\infty}\frac{(i+\ell)!}{i!(2\ell+2i+1)!}\sum_{p,q\not=1}\frac{1}{q-1}\biggl{[}\frac{1}{\varepsilon}-\sum_{k=0}^{\ell+i}\frac{1}{2k+1}+\ln\left(\sqrt{\bar{q}}\,\frac{\ell_{0}^{q-1}}{r_{0}^{q-1}}\right)\biggr{]}\frac{f_{p,q}^{(2i)}(t)}{c^{2i}}\,,$
(34)

plus the neglected remainder $\mathcal{O}(\varepsilon)$. It is understood that
the sum over $i$ is actually finite because the singularity in the source term
(25) is finite (i.e. $p\leqslant p_{0}$ for some $p_{0}\in\mathbb{N}$). In
(34) we have restored the relevant Hadamard regularization scale $r_{0}$ and
dimensional regularization scale $\ell_{0}$. Again we emphasize that this
result is local, as long as the $f_{p,q}(t)$ are themselves local. As
presented in the Sec. IV, this is always the case, and thus the complicated
difference between the two regularization schemes of non-linear interactions
is local. This appears to be a non-trivial result, considering that the tails
and their iterations are by essence non-local.

To summarize, we have computed the difference between regularization schemes
in the metric in harmonic coordinates. From (21) and (23) we have expressed
this difference, for each component of the metric and each multipolarity
$\ell$, see (6), in the form of a retarded solution of the d’Alembertian
equation in $d$ dimensions,

$\displaystyle\mathcal{D}h_{L}=\hat{\partial}_{L}\,\biggl{[}\frac{\tilde{k}}{r^{d-2}}\,\int_{1}^{+\infty}\mathrm{d}y\,\gamma_{\frac{1-d}{2}}(y)\,\mathcal{D}H(t-yr)\biggr{]}\,,$
(35)

where the function $\mathcal{D}H$ is given by (34) for any source term of the
form (25)

##  IV The MPM algorithm for non-linear interactions

The goal of this work is to apply the formulas derived in the previous section
to the case of the non-linear interactions entering the mass quadrupole moment
at 4PN. A quick dimensional analysis shows that such interactions are at most
cubic in $G$. Thus, we need to compute the cubic source term entering the
right-hand side of the wave equation (1). To this end, we rely on the so-
called Multipolar-Post-Minkowskian (MPM) algorithm [32, 30].

The MPM algorithm builds the sources (and solutions) of the wave equation (1)
in an iterative fashion, orders by orders in $G$, hence the name. Its starting
point is the so-called canonical metric, which is a generalization in $d$
dimensions of Thorne’s linearized metric [33]. The recent calculation in $d$
dimensions [28] writes it as

$\displaystyle h^{00}_{\text{can\,1}}=-\frac{4}{c^{2}}\sum_{\ell\geqslant
0}\frac{(-)^{\ell}}{\ell!}\,\hat{\partial}_{L}\,\widetilde{\mathcal{M}}_{L}\,,$
(36a) $\displaystyle h^{0i}_{\text{can\,1}}=\frac{4}{c^{3}}\sum_{\ell\geqslant
1}\frac{(-)^{\ell}}{\ell!}\,\left[\hat{\partial}_{L-1}\,\widetilde{\mathcal{M}}_{iL-1}^{(1)}+\frac{\ell}{\ell+1}\hat{\partial}_{L}\widetilde{\mathcal{S}}_{i|L}\right]\,,$
(36b) $\displaystyle
h^{ij}_{\text{can\,1}}=-\frac{4}{c^{4}}\sum_{\ell\geqslant
2}\frac{(-)^{\ell}}{\ell!}\,\left[\hat{\partial}_{L-2}\,\widetilde{\mathcal{M}}_{ijL-2}^{(2)}+\frac{2\ell}{\ell+1}\hat{\partial}_{L-1}\widetilde{\mathcal{S}}^{(1)}_{(i|j)L-1}+\frac{\ell-1}{\ell+1}\hat{\partial}_{L}\widetilde{\mathcal{K}}_{ij|L}\right]\,,$
(36c)

where we have defined the “tilded” versions of the moments as, e.g.,

$\widetilde{\mathcal{M}}_{L}(r,t)\equiv\frac{\tilde{k}}{r^{d-2}}\int_{1}^{+\infty}\\!\\!\mathrm{d}y\,\gamma_{\frac{1-d}{2}}(y)\,M_{L}\left(t-\frac{yr}{c}\right)\,,$
(37)

and where $M_{L}$ and $S_{i|L}$ denote the mass and current canonical moments,
together with the additional moment $K_{ij|L}$, which does not exist in 3
dimensions.

We let the reader refer to [2] for a comprehensive review of the MPM
formalism.777Note that the procedure is formally similar in three-dimensions
(as presented in [2]) and in $d$-dimensions (as should be implemented here);
however the three-dimensional moments have to be replaced by their
$d$-dimensional tilded counterparts, e.g.
$M_{L}(t-r/c)/r\longrightarrow\widetilde{\mathcal{M}}_{L}(r,t)$ as defined in
(37). We will only note that, due to its iterative nature, we need the
knowledge of the quadratic solution to the wave equation in order to build up
the cubic source. However, let us recall that the difference between
regularization schemes is a near zone effect [as only the ${r^{\prime}}\to 0$
limit of the integral (20) plays a role], and so we can focus on the near zone
behaviour of the source only. This notably implies that we do not require the
knowledge of the whole quadratic interactions to build up the cubic source,
but only their near zone behaviour, which we can iterate following the MPM
procedure.

We discuss the near zone expansion of a quadratic metric, taking first the
example of the quadratic interaction $M\times M_{ij}$. Later we will say a
word on the more difficult interaction $M_{ij}\times M_{ij}$, but the
generalization to other interactions is almost straightforward. The source of
a $d$-dimensional quadratic interaction like $M\times M_{ij}$ (for which the
monopole $\widetilde{\mathcal{M}}=\tilde{k}Mr^{2-d}$ is static) takes the
generic form

$N(\mathbf{x},t)=\hat{n}_{L}\,\frac{\ell_{0}^{q\varepsilon}}{r^{p+q\varepsilon}}\int_{1}^{+\infty}\\!\\!\mathrm{d}z\,\gamma_{\frac{1-d}{2}}(z)\,z^{k}\,F\left(t-\frac{zr}{c}\right)\,,$
(38)

and $F$ involves some time derivative of $M_{ij}$. The exact prescription to
be applied would naturally be the propagator (2), as we discussed in details
in Sec. II. Nevertheless, we do not need to use it for the quadratic
interactions, as we only seek for the near zone behaviour. Following [19], the
formal PN expansion of the quadratic solution is made of two terms,

$\displaystyle\overline{h}_{L}=\text{PF}\,\Box_{\text{R}}^{-1}\overline{N}-\mathop{\mathrm{PF}}_{B=0}\sum_{m\in\mathbb{N}}\frac{(-)^{m}}{m!}\sum_{j\in\mathbb{N}}\Delta^{-j}\hat{x}_{M}\int\\!\mathrm{d}^{d}\mathbf{x}^{\prime}\left(\frac{r^{\prime}}{r_{0}}\right)^{B}\hat{\partial}_{M}^{\prime}\left[\frac{\widetilde{\mathcal{N}}^{(2j)}(\mathbf{y},t-zr^{\prime}/c)}{c^{2j}\,r^{\prime}{}^{d-2}}\right]_{\mathbf{y}=\mathbf{x}^{\prime}}\,,$
(39)

where the “tilde” operation (with respect to variables $r^{\prime},t$) is
defined by (37) and we denoted

$\Delta^{-j}\hat{x}_{L}\equiv\frac{\Gamma\left(\frac{d}{2}+\ell\right)}{\Gamma\left(\frac{d}{2}+\ell+j\right)}\frac{r^{2j}\hat{x}_{L}}{2^{2j}j!}\,.$
(40)

The first term is a particular solution built up (in the sense of the PF) from
the formal PN expansion $\overline{N}$ of the source term (38). Introducing
$\tau=zr/c$ and using the asymptotic expansion when $z\to+\infty$ of the
function $\gamma_{\frac{1-d}{2}}$ given by (3), i.e.

$\gamma_{\frac{1-d}{2}}(z)\underset{z\to+\infty}{\sim}\
\sum_{j=0}^{+\infty}\frac{(-)^{j}}{j!}\frac{2\sqrt{\pi}\,z^{-2-2j-\varepsilon}}{\Gamma\left(-j-\frac{\varepsilon}{2}\right)\Gamma\left(\frac{1+\varepsilon}{2}\right)}\,,$
(41)

we can decompose $\overline{N}$ into “even” and “odd” parts given by

$\displaystyle\overline{N}^{\text{even}}=\sum_{j\in\mathbb{N}}A_{j}^{k}\,\frac{r^{2j+\kappa-
p-q\varepsilon}}{c^{2j+\kappa}}\,F^{(2j+\kappa)}(t)\,\hat{n}_{L}\,,$ (42a)
$\displaystyle\overline{N}^{\text{odd}}=\sum_{j\in\mathbb{N}}B_{j}^{k}\,\frac{r^{1+2j-k-p+(1-q)\varepsilon}}{c^{1+2j-k+\varepsilon}}\int_{0}^{\infty}\\!\\!\mathrm{d}\tau\,\tau^{-\varepsilon}\,F^{(2j+2-k)}(t-\tau)\,\hat{n}_{L}\,,$
(42b)

where $\kappa\equiv k-2[\frac{k}{2}]$, i.e. $\kappa=0$ if $k$ is even and $1$
if $k$ is odd, and the coefficients that enter those expressions are

$\displaystyle
A_{j}^{k}\equiv\frac{(-)^{j}}{2^{2j+\kappa}j!}\frac{\Gamma\left(\frac{1}{2}-j-\kappa\right)}{\Gamma\left(\frac{1+\varepsilon}{2}\right)}\frac{\Gamma\left(\frac{1-\kappa-k+\varepsilon}{2}-j\right)}{\Gamma\left(\frac{1-\kappa-k}{2}-j\right)}\,,$
(43a) $\displaystyle
B_{j}^{k}\equiv\frac{(-)^{k}}{2^{2j}j!}\frac{2\sqrt{\pi}\,\Gamma\left(2-k+\varepsilon\right)}{\Gamma\left(\frac{1+\varepsilon}{2}\right)\Gamma\left(1-\varepsilon\right)\Gamma\left(\frac{-\varepsilon}{2}\right)}\frac{\Gamma\left(1+j+\frac{\varepsilon}{2}\right)}{\Gamma\left(1+\frac{\varepsilon}{2}\right)}\frac{\Gamma\left(1-\frac{k}{2}+\frac{\varepsilon}{2}\right)}{\Gamma\left(1+j-\frac{k}{2}+\frac{\varepsilon}{2}\right)}\frac{\Gamma\left(\frac{3-k+\varepsilon}{2}\right)}{\Gamma\left(j+\frac{3-k+\varepsilon}{2}\right)}\,.$
(43b)

The expressions (42)–(43) generalize those of [19] in the case of a non-
vanishing $k$.

Once integrated with the PN-expanded propagator, using the “Matthieu” formula,
i.e. Eq. (3.16) in [14], we find that Eqs. (42) yield poles. Fortunately, the
three-dimensional limit is non-pathological, since those poles then disappear
by mutual cancellation. They simply generate logarithmic terms as emphasized
by (26). Moreover, such poles do not contribute to the difference between
regularization schemes in the iterated cubic solution, as already advertized
(they do not enter the terms bearing $p-\ell-3\in 2\mathbb{N}$).

As for the homogeneous solution, which is the second term of (39), after a
series of transformations following [19], it becomes

$\displaystyle\overline{h}_{L}^{\text{hom}}=\frac{(-)^{p+\ell}}{d+2\ell-2}\mathop{\mathrm{PF}}_{B=0}\,\frac{\Gamma\left(q\varepsilon-B\right)}{\Gamma\left(p+\ell-1+q\varepsilon-B\right)}\,C^{k,p,q}_{\ell}\,\sum_{j\in\mathbb{N}}\Delta^{-j}\hat{x}_{L}\int_{0}^{\infty}\\!\\!\mathrm{d}\tau\frac{\tau^{B-q\varepsilon}}{r_{0}^{B}}\frac{F^{(2j+\ell+p-1)}(t-\tau)}{c^{2j+\ell+p+q\varepsilon-B}}\,,$
(44)

where we have denoted

$C^{k,p,q}_{\ell}\equiv\int_{1}^{+\infty}\\!\\!\mathrm{d}y\,\gamma_{\frac{1-d}{2}-\ell}(y)\int_{1}^{+\infty}\\!\\!\mathrm{d}z\,\gamma_{\frac{1-d}{2}}(z)\,z^{k}(y+z)^{\ell-2+p+q\varepsilon-B}\,.$
(45)

The coefficients $C^{k,p,q}_{\ell}$ are generalizations for $q\in\mathbb{Z}$
of the ones introduced in [19]. They can be computed by methods similar to
those described in the App. D of [19].

The three types of terms we have met: integrated “even” and “odd” PN
expansions, as well as the homogeneous solutions, are to be inserted into the
cubic source which, again, can be PN-expanded and will take the same general
form as in (42), which confirms our claim on the structure of the cubic source
(25).

An important point to note is that only the “even” contributions (42a) will
contribute to the final result (34). Indeed we have explicitly checked that,
due to the particular powers of $r$, the “odd” terms given by (42b) will bear
$p-\ell-3=2i+1$, $i\in\mathbb{N}$, thus not contributing to the final
difference. Concerning the homogeneous solution, it is apparent from Eqs. (44)
and (40) that it bears $q=0$. As the linear metric bears $q=1$, the
homogeneous solution, once inserted into the cubic source, will produce only
terms with $q=1$, thus evading the condition for a term to contribute to the
difference. This crucial remark that only “even” terms contribute makes the
cubic source become local in the end, once expanded in the near zone.
Therefore, our final result (34) will itself be local, as claimed at the end
of Sec. III.

The only quadratic interaction that does not fall within the previous
investigation is the memory type one $M_{ij}\times M_{ij}$, which is needed to
compute the cubic source of the 4PN tail-of-memory (see Sec. VI.1). This
interaction involves two time-dependent moments, so that its source reads

$N(\mathbf{x},t)=\hat{n}_{L}\,\frac{\ell_{0}^{q\varepsilon}}{r^{p+q\varepsilon}}\int_{1}^{+\infty}\\!\\!\mathrm{d}y\,\gamma_{\frac{1-d}{2}}(y)\,y^{s}\int_{1}^{+\infty}\\!\\!\mathrm{d}z\,\gamma_{\frac{1-d}{2}}(z)\,z^{k}\,F\left(t-\frac{yr}{c}\right)\,G\left(t-\frac{zr}{c}\right)\,.$
(46)

In this case, the PN expansion of the solution involves five types of terms:
the squares “even-even” and “odd-odd”, the cross products “even-odd” and “odd-
even”, and the homogeneous solution, which keeps the same structure as Eq.
(44). Nevertheless, by arguments similar as for the tail interaction, only the
“even-even” terms contribute to the final result. Note that, although the
“odd-odd” terms could contribute as they contain terms with $p-\ell-3\in
2\mathbb{N}$, it turns out that they have $q=0$. Therefore, the same argument
as the one used to discard the contribution of the homogeneous solution
applies. From (42a), we see again that the result for the difference at cubic
order will be purely local.

Finally, applying the formula (35) on each terms of the (PN-expanded) cubic
source, we obtain the difference between regularizations in the form of a
piece $\mathcal{D}h^{\mu\nu}$ in the metric to the cubic order,

$\displaystyle\mathcal{D}h^{\mu\nu}=\sum_{\ell\in\mathbb{N}}\hat{\partial}_{L}\,\biggl{[}\frac{\tilde{k}}{r^{d-2}}\,\int_{1}^{+\infty}\mathrm{d}y\,\gamma_{\frac{1-d}{2}}(y)\,\mathcal{D}H^{\mu\nu}\left(t-\frac{yr}{c}\right)\biggr{]}\,.$
(47)

We are almost done but, while this piece clearly obeys
$\Box\mathcal{D}h^{\mu\nu}=0$, it does not a priori satisfy the harmonic
coordinate condition by itself, $\partial_{\nu}\mathcal{D}h^{\mu\nu}\not=0$.
Next, we compute from it another object, $\mathcal{D}h^{\mu\nu}_{\text{gen}}$,
which happens to be a general solution of the linearized Einstein field
equations in vacuum (and $d$ dimensions), i.e. which satisfies at once
$\Box\mathcal{D}h^{\mu\nu}_{\text{gen}}=0$ and
$\partial_{\nu}\mathcal{D}h^{\mu\nu}_{\text{gen}}=0$. In doing so, we are
exactly following the MPM procedure [32, 30], which can be summarized as

$\displaystyle
u^{\mu\nu}\equiv\mathcal{D}h^{\mu\nu}~{}\longrightarrow~{}w^{\mu}=\partial_{\nu}u^{\mu\nu}~{}\longrightarrow~{}v^{\mu\nu}=\mathcal{H}(w^{\mu})~{}\longrightarrow~{}\mathcal{D}h^{\mu\nu}_{\text{gen}}\equiv
u^{\mu\nu}+v^{\mu\nu}\,,$ (48)

$v^{\mu\nu}$ being here a homogeneous solution such that
$\partial_{\nu}v^{\mu\nu}=-w^{\mu}$, built from $w^{\mu}$ with the help of the
“harmonicity” algorithm $\mathcal{H}$ given by Eqs. (2.11)–(2.12) in [34].
Since $\mathcal{D}h^{\mu\nu}_{\text{gen}}$ solves the linearized vacuum field
equations, it can be written in a unique way as

$\displaystyle\mathcal{D}h^{\mu\nu}_{\text{gen}}=\mathcal{D}h^{\mu\nu}_{\text{can}}+\text{(linearized
gauge transformation)}\,,$ (49)

where the gauge transformation can be ignored as it plays no role when looking
for the invariants in the metric. The physically relevant quantity is another
object, $\mathcal{D}h^{\mu\nu}_{\text{can}}$, given as some simpler canonical
solution, which naturally takes the same form as the linearized metric (36),
i.e.

$\displaystyle\mathcal{D}h^{00}_{\text{can}}=-\frac{4}{c^{2}}\sum_{\ell\geqslant
0}\frac{(-)^{\ell}}{\ell!}\,\hat{\partial}_{L}\,\delta\widetilde{\mathcal{M}}_{L}\,,$
(50a)
$\displaystyle\mathcal{D}h^{0i}_{\text{can}}=\frac{4}{c^{3}}\sum_{\ell\geqslant
1}\frac{(-)^{\ell}}{\ell!}\,\left[\hat{\partial}_{L-1}\,\delta\widetilde{\mathcal{M}}_{iL-1}^{(1)}+\frac{\ell}{\ell+1}\hat{\partial}_{L}\delta\widetilde{\mathcal{S}}_{i|L}\right]\,,$
(50b)
$\displaystyle\mathcal{D}h^{ij}_{\text{can}}=-\frac{4}{c^{4}}\sum_{\ell\geqslant
2}\frac{(-)^{\ell}}{\ell!}\,\left[\hat{\partial}_{L-2}\,\delta\widetilde{\mathcal{M}}_{ijL-2}^{(2)}+\frac{2\ell}{\ell+1}\hat{\partial}_{L-1}\delta\widetilde{\mathcal{S}}^{(1)}_{(i|j)L-1}+\frac{\ell-1}{\ell+1}\hat{\partial}_{L}\delta\widetilde{\mathcal{K}}_{ij|L}\right]\,,$
(50c)

where we have straightforwardly defined

$\delta\widetilde{\mathcal{M}}_{L}(r,t)\equiv\frac{\tilde{k}}{r^{d-2}}\int_{1}^{+\infty}\\!\\!\mathrm{d}y\,\gamma_{\frac{1-d}{2}}(y)\,\mathcal{D}M_{L}\left(t-\frac{yr}{c}\right)\,,$
(51)

and where $\mathcal{D}M_{L}$ and $\mathcal{D}S_{i|L}$ denote the corresponding
corrections to the mass and current canonical moments $M_{L}$ and $S_{i|L}$,
together with those of the additional moment $K_{ij|L}$, which does not exist
in 3 dimensions. In conclusion, when written into the form (50), our result is
directly formatted to read off the non-linear corrections to the canonical
moments. From the detailed calculations of the cubic interactions (in
particular $M\times M_{ij}\times M_{ij}$) outlined below, and the comparison
with (50)–(51), we will find $\delta\widetilde{\mathcal{M}}_{ij}(r,t)$ and,
thus, the looked-for correction to the canonical quadrupole moment
$\mathcal{D}M_{ij}(t)$. This result will play a crucial role in the definition
of the “renormalized mass quadrupole” in Sec. VI.1.

Finally, we have used two methods to extract the correction to the canonical
moment $\delta\widetilde{\mathcal{M}}_{ij}$ from the metric. One method
consists of considering the leading order $1/r^{d-2}$ term in the metric when
$r\to+\infty$ and performing a suitable angular integration (in $d$
dimensions) to obtain $\delta\widetilde{\mathcal{M}}_{ij}^{(2)}$, then
$\delta\widetilde{\mathcal{M}}_{ij}$. The second method was to extract it
directly from the gauge invariant components of the linearized Riemann tensor

$\mathcal{R}_{0i0j}=\frac{G}{2}\left(\partial_{ij}\mathcal{D}h^{00}_{\text{can}}+\frac{2}{c}\,\partial_{t}\partial^{(i}\mathcal{D}h^{j)0}_{\text{can}}+\frac{1}{c^{2}}\,\partial_{t}^{2}\mathcal{D}h^{ij}_{\text{can}}+\frac{\partial_{ij}\mathcal{D}h_{\text{can}}}{d-1}-\frac{1}{c^{2}}\frac{\delta_{ij}}{d-1}\,\partial_{t}^{2}\mathcal{D}h_{\text{can}}\right)\,,$
(52)

where
$\mathcal{D}h_{\text{can}}\equiv-\mathcal{D}h^{00}_{\text{can}}+\mathcal{D}h^{ii}_{\text{can}}$.
Both methods perfectly agree with each other.

##  V Regularization of non-linear interactions at 3PN order

When computing the radiative mass quadrupole moment at 3PN order [12], a
Hadamard regularization scheme for the IR sector has been used for the source
moment and the non-linear interactions have been computed in ordinary 3
dimensions. However, the crucial need for a dimensional regularization scheme
in the 4PN equations of motion [19] implies that one should trade the Hadamard
regularization for dimensional regularization when computing the mass
quadrupole at the 4PN order. This has been achieved in the companion paper
[14]. It was found that, when using a dimensional regularization scheme, the
source quadrupole moment receives corrections including poles
$\propto\varepsilon^{-1}$ already starting at the 3PN order.

To describe the poles, it is convenient to “dress” them into the particular
combination

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

Here $r_{0}$ is the scale associated to the Partie Finie and $\ell_{0}$ the
one associated to the dimensional regularization. Now, the detailed
calculations in [14] conclude that, at the 3PN order, the dimensional
regularization corrections can be recast as

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

where $I_{ij}$ is the source quadrupole moment and $P_{i}$ is the conserved
linear momentum, i.e. $P_{i}=I^{(1)}_{i}$, with $I_{i}$ being the linearly
varying mass dipole. Note that this formula is valid in an arbitrary frame,
not necessarily the center-of-mass (CoM) frame for which $P_{i}=0$. Obviously,
$I_{ij}$ and $P_{i}$ are merely Newtonian with the 3PN accuracy of (54) but
will contain crucial 1PN corrections when this equation is applied to 4PN
order.

We recognize $\beta_{I}=-\frac{214}{105}$ in the coefficient of the first term
in (54). This coefficient is well known from “tail-of-tail” calculations [27].
It can be interpreted as the beta function coefficient associated with the
logarithmic renormalization of the mass quadrupole moment [35].

The 3PN correction to the source quadrupole (54) bears a pole, which looks
dreadful. Fortunately, the source quadrupole is not a physical quantity,
contrary to the radiative quadrupole defined at infinity. We thus have to
ensure that the correction (54) is exactly compensated by the regularization
of the 3PN non-linear interactions in the radiative moment. Otherwise, it
would show that the formalism used previously at 3PN order is inconsistent
with dimensional regularization.

A trivial dimensional analysis shows that only two interactions can arise at
3PN order: (i) the tail-of-tail $M\times M\times M_{ij}$, where $M$ is the
constant ADM mass and $M_{ij}$ the canonical moment (which agrees with the
source moment to 2PN order), and (ii) the interaction $M\times M_{i}\times
M_{i}$ where $M_{i}$ is the mass dipole moment, varying linearly with time and
vanishing in the CoM frame. As we have proved in Sec. III [see (34)], the
regularization of the non-linear interactions linking the canonical and
radiative quadrupoles yields purely local corrections. Moreover, since in
general relativity the propagation of GW is described entirely by the two
types of multipole moments $M_{L}$ and $S_{i|L}$ (called the canonical moments
in our terminology) we have to express the non-linear terms as corrections in
the canonical moments, as discussed in Sec. IV [see (50)].

Implementing the non-linear iteration of the requested multipole interactions,
as described in Sec. IV, we computed the 3PN correction to the canonical
quadrupole:

$\mathcal{D}M_{ij}^{\text{3PN}}=-\frac{214}{105}\,\frac{G^{2}M^{2}}{c^{6}}\left(\Pi_{\varepsilon}+\frac{246\,299}{44\,940}\right)\,M_{ij}^{(2)}+\frac{428}{105}\,\frac{G^{2}M}{c^{6}}\left(\Pi_{\varepsilon}+\frac{252\,599}{44\,940}\right)\,M^{(1)}_{\langle
i}M^{(1)}_{j\rangle}\,.$ (55)

As the canonical and source moments are equivalent at Newtonian order, we see
that the contribution of the non-linear interactions exactly cancels the
correction of the source quadrupole given by Eq. (54). Again, this result is
valid in a general frame, not necessarily in the CoM frame, and it is true not
only for the pole part, but also for the finite part following the pole.

We conclude that the corrections due to dimensional IR regularization exactly
vanish in the radiative quadrupole to 3PN order; so, using Hadamard or
dimensional regularization for the IR is equivalent at that order. A priori,
the same cancellation is not expected to happen at 4PN order, because we found
for the conservative equations of motion that the dimensional and Hadamard
regularization schemes are no more equivalent [19]. Correcting the quadrupole
moment to account for dimensional regularization at 4PN order is the main
motivation for this work.

Two other radiative moments are very important as they have to be computed to
3PN order in order to control the 4PN flux: these are the mass octupole
moment, obtained in [29], and the current quadrupole moment, computed in [28].
We have to verify the equivalence of both regularization schemes for these 3PN
moments too. Following the methods described in [14], we have computed the
correction due to the change in regularization schemes for both the source
mass octupole $I_{ijk}$ and the (dual of the) current quadrupole $J_{i|jk}$,
with results

$\displaystyle\mathcal{D}I_{ijk}^{\text{3PN}}$
$\displaystyle=\frac{26}{21}\frac{G^{2}}{c^{6}}\left(\Pi_{\varepsilon}+\frac{9281}{2730}\right)M^{2}I_{ijk}^{(2)}+\frac{92}{35}\frac{G^{2}}{c^{6}}\left(\Pi_{\varepsilon}-\frac{161\,597}{28\,980}\right)I_{\langle
i}P_{j}P_{k\rangle}$
$\displaystyle-\frac{52}{7}\frac{G^{2}}{c^{6}}\left(\Pi_{\varepsilon}+\frac{14\,387}{4\,095}\right)MI_{\langle
ij}^{(1)}P_{k\rangle}+\frac{12}{5}\frac{G^{2}}{c^{6}}\left(\Pi_{\varepsilon}+\frac{63\,421}{7560}\right)MI_{\langle
ij}^{(2)}I_{k\rangle}\,,$ (56a)
$\displaystyle\mathcal{D}J_{i|jk}^{\text{3PN}}$
$\displaystyle=\frac{214}{105}\frac{G^{2}}{c^{6}}\left(\Pi_{\varepsilon}+\frac{54\,989}{44\,940}\right)M^{2}J_{i|jk}^{(2)}$
$\displaystyle+\frac{856}{105}\frac{G^{2}}{c^{6}}\left(\Pi_{\varepsilon}+\frac{229\,289}{44\,940}\right)MI_{\\{i|j}^{(2)}P_{k\\}}-\frac{106}{7}\frac{G^{2}}{c^{6}}\,MI_{\\{i|j}^{(3)}I_{k\\}}\,.$
(56b)

Here
$\\{i|jk\\}\equiv\underset{ij}{\mathcal{A}}\,\underset{ijk}{\text{TF}}\,\underset{jk}{\text{STF}}$
denotes the appropriate symmetries of the $d$-dimensional current quadrupole
derived in [28] [see Eq. (2.36) there]. In Eqs. (56), we again recognize the
usual beta function coefficients $\beta_{O}=-26/21$ and $\beta_{J}=-214/105$
($=\beta_{I}$)888These coefficients are known for both mass and current tail-
of-tail interactions up to any multipolar order [15, 36]. in front of the
tail-of-tail contribution to the octupole and current quadrupole. On the other
hand, using the formulas of section III, we derived the corrections to those
canonical moments as

$\displaystyle\mathcal{D}M_{ijk}^{\text{3PN}}$
$\displaystyle=-\frac{26}{21}\frac{G^{2}}{c^{6}}\left(\Pi_{\varepsilon}+\frac{9281}{2730}\right)M^{2}M_{ijk}^{(2)}-\frac{92}{35}\frac{G^{2}}{c^{6}}\left(\Pi_{\varepsilon}-\frac{161\,597}{28\,980}\right)M_{\langle
i}M^{(1)}_{j}M^{(1)}_{k\rangle}$
$\displaystyle+\frac{52}{7}\frac{G^{2}}{c^{6}}\left(\Pi_{\varepsilon}+\frac{14\,387}{4\,095}\right)MM_{\langle
ij}^{(1)}M^{(1)}_{k\rangle}-\frac{12}{5}\frac{G^{2}}{c^{6}}\left(\Pi_{\varepsilon}+\frac{63\,421}{7560}\right)MM_{\langle
ij}^{(2)}M_{k\rangle}\,,$ (57a)
$\displaystyle\mathcal{D}S_{i|jk}^{\text{3PN}}$
$\displaystyle=-\frac{214}{105}\frac{G^{2}}{c^{6}}\left(\Pi_{\varepsilon}+\frac{54\,989}{44\,940}\right)M^{2}S_{i|jk}^{(2)}$
$\displaystyle-\frac{856}{105}\frac{G^{2}}{c^{6}}\left(\Pi_{\varepsilon}+\frac{229\,289}{44\,940}\right)MM_{\\{i|j}^{(2)}M^{(1)}_{k\\}}+\frac{106}{7}\frac{G^{2}}{c^{6}}\,MM_{\\{i|j}^{(3)}M_{k\\}}\,.$
(57b)

Similarly to the case of the 3PN mass quadrupole, we find that the proper
regularization of the non-linear interactions exactly compensates the
correction due to the dimensional regularization of the source mass octupole
and current quadrupole. As clear from (56) and (57), this is proven not only
for the pole parts but also for the finite terms beyond the poles. In other
words, for all these 3PN moments, the corrections due to dimensional
regularization exactly vanish in the radiative moments, which are the
physically relevant quantities. In conclusion, the present investigations have
confirmed the previous calculations of multipole moments at the 3PN order
presented in Refs. [10, 12, 29, 28].

##  VI The renormalized mass quadrupole moment at 4PN order

###  VI.1 Regularization of non-linear interactions at 4PN order

Having proved that Hadamard and dimensional (IR) regularizations are
equivalent at the 3PN order, we now turn to the main goal of this work, namely
the IR regularization of the radiative mass quadrupole up to 4PN.

In addition to the 3PN interaction (55) which has now to be considered up to
4PN, five additional interactions can show up: the “tail-of-memory” $M\times
M_{ij}\times M_{ij}$, the interaction $M\times M_{ij}\times S_{i|j}$ involving
the constant angular momentum $S_{i|j}$, and three interactions $M\times
M_{i}\times M_{ijk}$, $M_{i}\times M_{i}\times M_{ij}$ and $M\times
M_{i}\times S_{i|jk}$ involving the mass dipole $M_{i}$. As discussed earlier,
the difference between the regularizations of the non-linear interactions are
local and can be expressed as corrections to the canonical moments $M_{L}$ and
$S_{i|L}$.

For the sake of simplicity, we will work in the CoM frame, where $M_{i}$
vanishes. Indeed, the main application of the present work for PN templates
can be performed in such frame. Only three interactions survive, yielding the
following corrections in the canonical moment:

$\displaystyle\mathcal{D}M_{ij}^{\text{CoM, 4PN}}=\,$
$\displaystyle-\frac{214}{105}\,\frac{G^{2}M^{2}}{c^{6}}\left(\Pi_{\varepsilon}+\frac{246\,299}{44\,940}\right)\,M_{ij}^{(2)}$
(58)
$\displaystyle+\frac{G^{2}M}{c^{8}}\Biggl{[}\frac{12}{7}\left(\Pi_{\varepsilon}+\frac{3581}{7560}\right)\,M_{a\langle
i}^{(2)}\,M_{j\rangle
a}^{(2)}-\frac{24}{7}\left(\Pi_{\varepsilon}-\frac{338}{2835}\right)\,M_{a\langle
i}^{(1)}\,M_{j\rangle a}^{(3)}$
$\displaystyle\qquad\qquad-\frac{4}{7}\left(\Pi_{\varepsilon}-\frac{1447}{216}\right)\,M_{a\langle
i}\,M_{j\rangle
a}^{(4)}+\frac{4}{3}\left(\Pi_{\varepsilon}+\frac{11\,243}{7560}\right)M_{a\langle
i}^{(3)}S_{j\rangle|a}\Biggl{]}\,,$

where we recall that the pole has been “dressed-up” according to (53).

As we have discussed in Sec. V the tail-of-tail contribution given by the
first line of (58) exactly cancels the correction due to dimensional
regularization in the source moment at the 3PN order. Note that this keeps
being true up to the 4PN order because the canonical and source moments are
equivalent at the 1PN order. Now, we discover that the coefficients of the
remaining poles entering the non-linear interactions in the purely 4PN terms
of (58) exactly cancel those coming from the IR regularization of the source
mass quadrupole, given by Eq. (4.4) of the companion paper [14].

Therefore, we arrive at the main result of this paper, namely that the
radiative mass quadrupole moment is free of poles up to the 4PN level. This
important feature (but expected for a physical quantity) naturally leads us to
define a “renormalized mass quadrupole” as the sum of the properly regularized
source mass quadrupole, which was computed in [13, 14], and augmented by the
contributions from the dimensional regularization of the cubic non-linear
interactions, as displayed in Eqs. (58), including the finite corrections
beyond the poles explicitly found in (58). The interest of such quantity is
that it is finite when $\varepsilon\to 0$, so that it can be considered and
manipulated in ordinary three dimensions. This notably allows us to express it
in the CoM frame and, then, for quasi-circular orbits using the usual three-
dimensional reduction procedure.

Recalling the work conducted in [14], the renormalized quadrupole is thus
defined as

$I^{\text{renorm}}_{ij}\equiv
I_{ij}+\mathcal{D}M_{ij}=I_{ij}^{\text{Had}}+I_{ij}^{\text{non-
loc}}+\mathcal{D}I_{ij}+\delta_{\chi}I_{ij}+\mathcal{D}M_{ij}\,,$ (59)

where $I_{ij}^{\text{Had}}$ denotes the result of [13], i.e. the local source
quadrupole computed with a UV dimensional regularization and an IR Hadamard
one; $I_{ij}^{\text{non-loc}}$ is the non-local effect of the source
quadrupole, coming from the 4PN tail effect in the conservative sector
(equations of motion and Fokker Lagrangian), given in Eq. (2.14) of [14];
$\mathcal{D}I_{ij}$ is the correction due to dimensional regularization in the
source quadrupole ($\mathcal{D}I_{ij}$ is itself the sum of four contributions
in Eq. (3.1) of [14]); $\delta_{\chi}I_{ij}$ is the effect of the local part
of the IR shift applied in the conservative sector [14, 19]; and, finally,
$\mathcal{D}M_{ij}$ is the correction due to dimensional regularization of the
non-linear interactions, computed in the present work.

The non-local tail term in the source quadrupole found in [14] can be restated
as

$\displaystyle I_{ij}^{\text{non-
loc}}=I_{ij}^{\text{inst}}+\hat{I}_{ij}^{\text{non-loc}}\,.$ (60)

The first term corresponds to instantaneous (non-tail) contributions given by

$\displaystyle\hat{I}_{ij}^{\text{inst}}=\frac{24}{7}\frac{G^{2}M}{c^{8}}\left(\Pi_{\varepsilon}+\frac{74}{105}\right)\,I_{a\langle
i}\,I_{j\rangle a}^{(4)}\,,$ (61)

and which combine with (58) and the other terms in (59) to cancel the pole.
The second term is, properly speaking, the purely non-local tail contribution
which reads as

$\displaystyle\hat{I}_{ij}^{\text{non-
loc}}=\frac{24}{7}\frac{G^{2}M}{c^{8}}\,I_{k\langle
i}(t)\int_{0}^{+\infty}\mathrm{d}\tau\,\ln\left(\frac{c\tau}{2r_{0}}\right)\,I^{(5)}_{j\rangle
k}(t-\tau)\,.$ (62)

Such contribution will naturally combine with the non-localities induced by
the non-linear interactions in the radiative quadrupole. The non-local tail
term contains both conservative and dissipative effects,

$\hat{I}_{ij}^{\text{non-loc}}=\hat{I}_{ij}^{\text{non-
loc}}\bigg{|}_{\text{cons}}+\hat{I}_{ij}^{\text{non-
loc}}\bigg{|}_{\text{diss}}\,,$ (63)

where the conservative sector is given by the time-symmetrized integral

$\displaystyle\hat{I}_{ij}^{\text{non-
loc}}\bigg{|}_{\text{cons}}=\frac{12}{7}\frac{G^{2}M}{c^{8}}\,I_{k\langle
i}\int_{0}^{+\infty}\mathrm{d}\tau\,\ln\left(\frac{c\tau}{2r_{0}}\right)\left[I^{(5)}_{j\rangle
k}(t-\tau)-I^{(5)}_{j\rangle k}(t+\tau)\right]\,,$ (64)

or can also be expressed in the form of the Hadamard partie finie (Pf)
integral

$\displaystyle\hat{I}_{ij}^{\text{non-
loc}}\bigg{|}_{\text{cons}}=\frac{12}{7}\frac{G^{2}M}{c^{8}}\,I_{k\langle
i}\,\text{Pf}_{2r_{0}/c}\int_{-\infty}^{+\infty}\frac{\mathrm{d}t^{\prime}}{|t-t^{\prime}|}\,I^{(4)}_{j\rangle
k}(t^{\prime})\,,$ (65)

with associated time scale $2r_{0}/c$. The dissipative effect is naturally
encoded in the complementary time-antisymmetrized integral

$\displaystyle\hat{I}_{ij}^{\text{non-
loc}}\bigg{|}_{\text{diss}}=\frac{12}{7}\frac{G^{2}M}{c^{8}}\,I_{k\langle
i}\int_{0}^{+\infty}\mathrm{d}\tau\,\ln\left(\frac{c\tau}{2r_{0}}\right)\left[I^{(5)}_{j\rangle
k}(t-\tau)+I^{(5)}_{j\rangle k}(t+\tau)\right]\,.$ (66)

Note that the notion of renormalized mass quadrupole applies to the source
type moment, hence we appropriately called it $I^{\text{renorm}}_{ij}$ and not
$M^{\text{renorm}}_{ij}$. In order to obtain the physical radiative moment, we
still have to add to the renormalized source moment the corrections coming
from (i) the relation between the canonical quadrupole on the one hand, and
the source and gauge moments on the other hand (see [30] for discussion), and
(ii) the non-linear tails, tails-of-tails, tails-of-memory, etc. multipole
interactions. However, the interest of having defined such renormalized source
quadrupole moment $I^{\text{renorm}}_{ij}$ is that the latter corrections (i)
and (ii) can be computed in ordinary three dimensions using the standard MPM
algorithm, as all dimensional regularization contributions have already been
included. Many of the corrections (i)-(ii) are already known (see notably [37,
38]), although not all of them at the 4PN order. Notably, the involved cubic
tail-of-memory terms $M\times M_{ij}\times M_{ij}$ will still have to be
computed in 3 dimensions, which is left for future work.

###  VI.2 Center-of-mass frame and quasi-circular orbits

As, by construction, the renormalized quadrupole is a three-dimensional
quantity, one can safely reduce it to the CoM frame and, then, specialize it
to the case of quasi-circular orbits, using usual three-dimensional reduction
formulas. Unfortunately, the complete quadrupole is too long to be presented,
even when expressed in the CoM frame. Therefore, we will only present the
difference between the purely local renormalized quadrupole and the Hadamard
one $\Delta I_{ij}\equiv I^{\text{renorm}}_{ij}-\hat{I}_{ij}^{\text{non-
loc}}-I^{\text{Had}}_{ij}$ in the CoM frame.

We denote by $r=|\bm{y}_{1}-\bm{y}_{2}|$ the radial separation in harmonic
coordinates, $\bm{x}=\bm{y}_{1}-\bm{y}_{2}$ the relative distance, and
$\bm{v}=\bm{v}_{1}-\bm{v}_{2}$ the relative velocity. We recall also the
dimensionless 1PN parameter $\gamma$ defined as $\gamma=\frac{Gm}{rc^{2}}$
with $m=m_{1}+m_{2}$ the total mass and $\nu=m_{1}m_{2}/m^{2}$ the symmetric
mass ratio. The difference $\Delta I_{ij}$ is rather compact, when written in
the CoM frame. It reads

$\displaystyle\Delta I_{ij}^{\text{CoM}}=$
$\displaystyle\frac{G^{2}m^{3}\nu}{c^{8}}\Biggl{\\{}\frac{G^{2}m^{2}}{r^{2}}\left[\frac{1072}{105}-\frac{571826}{14175}\nu\right]n^{\langle
i}n^{j\rangle}+\left[\frac{7516}{735}-\frac{123944}{6615}\nu\right]v^{2}\,v^{\langle
i}v^{j\rangle}$ (67)
$\displaystyle\quad+\frac{Gm}{r}\biggl{(}\left[\frac{3004}{147}+\frac{1800322}{11025}\nu\right](n.v)^{2}\,n^{\langle
i}n^{j\rangle}+\left[-\frac{30004}{2205}-\frac{9985783}{99225}\nu\right]v^{2}\,n^{\langle
i}n^{j\rangle}$
$\displaystyle\qquad\qquad+\left[-\frac{7516}{245}-\frac{2437018}{19845}\nu\right](n.v)\,n^{\langle
i}v^{j\rangle}+\left[-\frac{15056}{2205}+\frac{161096}{6615}\nu\right]v^{\langle
i}v^{j\rangle}\biggr{)}\Biggr{\\}}\,.$

It is fully apparent that the dimensional regularization brings in a finite
contribution to the quadrupole moment at the 4PN order. In previous work on
4PN equations of motion [19, 20], adding the analogue of (67) permitted
resolving the long-standing problem of regularization ambiguities.

Finally, we give the complete end result for the renormalized quadrupole on
quasi-circular orbits. We write it in the form999Note that, after publishing
the Hadamard source quadrupole in [13], we have spotted an error in the
$d$-dimensional computation of the value of the potential $\hat{R}_{i}$ at 1PN
order, when evaluated in $\bm{y}_{1,2}$. This error induces a small change in
the value of the coefficients, namely $\delta A=-4\gamma^{4}\nu^{2}/63$ and
$\delta B=4\gamma^{3}\nu^{2}/63$, which has been taken into account here.

$I_{ij}^{\text{renorm}}=m\nu\,\left(A\,x_{\langle
i}x_{j\rangle}+B\,\frac{r^{2}}{c^{2}}v_{\langle
i}v_{j\rangle}+\frac{G^{2}m^{2}\nu}{c^{5}r}\,C\,x_{\langle
i}v_{j\rangle}\right)+\mathcal{O}\left(\frac{1}{c^{9}}\right)\,,$ (68)

where the coefficients read

$\displaystyle A$
$\displaystyle=1+\gamma\biggl{(}-\frac{1}{42}-\frac{13}{14}\nu\biggr{)}+\gamma^{2}\biggl{(}-\frac{461}{1512}-\frac{18395}{1512}\nu-\frac{241}{1512}\nu^{2}\biggr{)}$
$\displaystyle\quad+\gamma^{3}\biggl{(}\frac{395899}{13200}-\frac{428}{105}\ln\biggl{(}\frac{r}{r_{0}{}}\biggr{)}+\biggl{[}\frac{3304319}{166320}-\frac{44}{3}\ln\biggl{(}\frac{r}{r^{\prime}_{0}}\biggr{)}\biggr{]}\nu+\frac{162539}{16632}\nu^{2}+\frac{2351}{33264}\nu^{3}\biggr{)}$
$\displaystyle\quad+\gamma^{4}\biggl{(}-\frac{1067041075909}{12713500800}+\frac{31886}{2205}\ln\biggl{(}\frac{r}{r_{0}{}}\biggr{)}+\biggl{[}-\frac{85244498897}{470870400}-\frac{2783}{1792}\pi^{2}-\frac{64}{7}\ln\left(16\gamma\mathrm{e}^{2\gamma_{\text{E}}}\right)$
$\displaystyle\qquad\quad-\frac{10886}{735}\ln\biggl{(}\frac{r}{r_{0}{}}\biggr{)}+\frac{8495}{63}\ln\biggl{(}\frac{r}{r^{\prime}_{0}}\biggr{)}\biggr{]}\nu+\biggl{[}\frac{171906563}{4484480}+\frac{44909}{2688}\pi^{2}-\frac{4897}{21}\ln\biggl{(}\frac{r}{r^{\prime}_{0}}\biggr{)}\biggr{]}\nu^{2}$
$\displaystyle\qquad\quad-\frac{22063949}{5189184}\nu^{3}+\frac{71131}{314496}\nu^{4}\biggl{)}\,,$
(69a) $\displaystyle B$
$\displaystyle=\frac{11}{21}-\frac{11}{7}\nu+\gamma\biggl{(}\frac{1607}{378}-\frac{1681}{378}\nu+\frac{229}{378}\nu^{2}\biggr{)}$
$\displaystyle\quad+\gamma^{2}\biggl{(}-\frac{357761}{19800}+\frac{428}{105}\ln\biggl{(}\frac{r}{r_{0}{}}\biggr{)}-\frac{92339}{5544}\nu+\frac{35759}{924}\nu^{2}+\frac{457}{5544}\nu^{3}\biggr{)}$
$\displaystyle\quad+\gamma^{3}\biggl{(}\frac{23006898527}{1589187600}-\frac{4922}{2205}\ln\biggl{(}\frac{r}{r_{0}{}}\biggr{)}+\biggl{[}\frac{8431514969}{529729200}+\frac{143}{192}\pi^{2}-\frac{32}{7}\ln\left(16\gamma\mathrm{e}^{2\gamma_{\text{E}}}\right)$
$\displaystyle\qquad\quad-\frac{1266}{49}\ln\biggl{(}\frac{r}{r_{0}{}}\biggr{)}-\frac{968}{63}\ln\biggl{(}\frac{r}{r^{\prime}_{0}}\biggr{)}\biggr{]}\nu$
$\displaystyle\qquad\quad+\biggl{[}\frac{351838141}{5045040}-\frac{41}{24}\pi^{2}+\frac{968}{21}\ln\biggl{(}\frac{r}{r^{\prime}_{0}}\biggr{)}\biggr{]}\nu^{2}-\frac{1774615}{81081}\nu^{3}-\frac{3053}{432432}\nu^{4}\biggl{)}\,,$
(69b) $\displaystyle C$
$\displaystyle=\frac{48}{7}+\gamma\left(-\frac{4096}{315}-\frac{24512}{945}\nu\right)-\frac{32}{7}\pi\,\gamma^{3/2}\,.$
(69c)

The term $C$ corresponds to time-odd contributions. The 2.5PN and 3.5PN terms
were computed in [38], whereas the new 4PN term
$-\frac{32}{7}\pi\,\gamma^{3/2}$ comes from the dissipative contribution of
the non-local term (66). As for the conservative contribution of the tail term
(65), it yields the logarithmic dependences
$\ln(16\gamma\mathrm{e}^{2\gamma_{\text{E}}})$ in the time-even coefficients
$A$ and $B$.

Finally, the two constants $r_{0}$ and $r_{0}^{\prime}$ entering the
coefficients (69) are associated with the Partie Finie and the dimensional UV
regularization schemes, respectively (see [13]). As unphysical scales, they
are expected to be exactly compensated in the radiative moment, by virtue of
non-linear effects and application of the time derivatives. This will be
checked in future work.

###### Acknowledgements.

We thank Gabriel Luz Almeida, Laura Bernard, Stefano Foffa, Sylvain Marsat,
Rafael Porto and Riccardo Sturani for interesting discussions. 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 Check with the $d=3$ limit of the retarded integral

In three dimensions, the function $\gamma_{\frac{1-d}{2}}(z)$ defined by (3)
takes the distributional form $\gamma_{-1}(z)=\delta(z-1)$, where
$\delta(z-1)$ means the Dirac function at point $z=1$. More generally, we have
for any $\ell\in\mathbb{N}$ (see App. C of [19])

$\displaystyle\gamma_{-1-\ell}(z)=\frac{1}{(2\ell-1)!!}\sum_{i=0}^{\ell}\beta_{i}^{\ell}\,\delta^{(i)}(z-1)\,,$
(70)

where $\delta^{(i)}$ denotes the $i$-th derivative of the Dirac function. The
numerical coefficients $\beta_{i}^{\ell}\equiv
2^{i-\ell}\frac{(2\ell-i)!}{i!(\ell-i)!}$ are such that the multipolar
retarded or advanced homogeneous wave in 3 dimensions reads

$\displaystyle\hat{\partial}_{L}\\!\left(\frac{F(t\mp r)}{r}\right)$
$\displaystyle=(-)^{\ell}\frac{\hat{n}_{L}}{r^{\ell+1}}\sum_{i=0}^{\ell}\beta_{i}^{\ell}\,(\pm
r)^{i}\,F^{(i)}(t\mp r)\,,$ (71)

We compute the 3-dimensional limit of our final $d$-dimensional result (18).
In this limit the two integrations over variables $y$ and $z$ can be effected
explicitly thanks to Eq. (70), while the integration over $\tau$ follows from
the limit (14) which just amounts to replacing the source term
$\widehat{N}_{\varepsilon}$ by $N$. Thus, we obtain in 3 dimensions

$\displaystyle
h_{L}=\frac{(-)^{\ell}}{2}\,\hat{n}_{L}\sum_{i=0}^{\ell}\sum_{j=0}^{\ell}\beta_{i}^{\ell}\,\beta_{j}^{\ell}\,r^{-1-\ell+j}$
$\displaystyle\qquad\times\Biggl{\\{}\int_{0}^{+\infty}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1+i}N^{(-2\ell-1+i+j)}({r^{\prime}},t-r-{r^{\prime}})$
$\displaystyle\qquad\qquad-(-)^{j}\int_{r}^{+\infty}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1+i}N^{(-2\ell-1+i+j)}({r^{\prime}},t+r-{r^{\prime}})$
$\displaystyle\qquad\qquad-(-)^{i}\int_{0}^{r}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1+i}N^{(-2\ell-1+i+j)}({r^{\prime}},t-r+{r^{\prime}})\Biggr{\\}}\,.$
(72)

After inserting the well-known formula for the multiple time anti-derivative
(valid for the source term becoming zero when ${t^{\prime}}\to-\infty$)

$\displaystyle
N^{(-2\ell-1+i+j)}({r^{\prime}},t^{\prime})=\int_{-\infty}^{t^{\prime}-r^{\prime}}\mathrm{d}s\,\frac{(t^{\prime}-r^{\prime}-s)^{2\ell-
i-j}}{(2\ell-i-j)!}N({r^{\prime}},s+{r^{\prime}})\,,$ (73)

and inverting the summations over ${r^{\prime}}$ and $s$, we can re-express
the result (A) into the form

$\displaystyle
h_{L}=\int_{-\infty}^{t-r}\mathrm{d}s\,\Biggl{\\{}\int_{0}^{\frac{t-r-s}{2}}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1}N({r^{\prime}},s+{r^{\prime}})\,\mathcal{A}_{L}$
(74)
$\displaystyle\qquad\qquad\qquad-\int_{r}^{\frac{t+r-s}{2}}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1}N({r^{\prime}},s+{r^{\prime}})\,\mathcal{B}_{L}-\int_{0}^{r}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1}N({r^{\prime}},s+{r^{\prime}})\,\mathcal{C}_{L}\Biggr{\\}}\,.$

Here, we have posed as intermediate notation

$\displaystyle\mathcal{A}_{L}$
$\displaystyle=\frac{(-)^{\ell}}{2}\frac{\hat{n}_{L}}{r^{1+\ell}}\sum_{i=0}^{\ell}\sum_{j=0}^{\ell}\beta_{i}^{\ell}\,\beta_{j}^{\ell}\,r^{j}{r^{\prime}}^{i}\frac{(t-r-s-2{r^{\prime}})^{2\ell-
i-j}}{(2\ell-i-j)!}\,,$ (75a) $\displaystyle\mathcal{B}_{L}$
$\displaystyle=\frac{(-)^{\ell}}{2}\frac{\hat{n}_{L}}{r^{1+\ell}}\sum_{i=0}^{\ell}\sum_{j=0}^{\ell}\beta_{i}^{\ell}\,\beta_{j}^{\ell}\,(-r)^{j}{r^{\prime}}^{i}\frac{(t+r-s-2{r^{\prime}})^{2\ell-
i-j}}{(2\ell-i-j)!}\,,$ (75b) $\displaystyle\mathcal{C}_{L}$
$\displaystyle=\frac{(-)^{\ell}}{2}\frac{\hat{n}_{L}}{r^{1+\ell}}\sum_{i=0}^{\ell}\sum_{j=0}^{\ell}\beta_{i}^{\ell}\,\beta_{j}^{\ell}\,r^{j}(-{r^{\prime}})^{i}\frac{(t-r-s)^{2\ell-
i-j}}{(2\ell-i-j)!}\,.$ (75c)

All the point is now to rewrite these quantities into a more convenient form.
For this we dispose of the following formula, valid for any $\rho$ and
${r^{\prime}}$,

$\displaystyle\sum_{i=0}^{\ell}\beta_{i}^{\ell}\,\frac{{r^{\prime}}^{i}(2\rho-2{r^{\prime}})^{2\ell-
i-j}}{(2\ell-
i-j)!}=\frac{2^{\ell-j}}{\ell!}\left(\frac{\mathrm{d}}{\mathrm{d}\rho}\right)^{j}\Bigl{[}\rho^{\ell}(\rho-{r^{\prime}})^{\ell}\Bigr{]}\,,$
(76)

together with the fact that the coefficients $\beta_{i}^{\ell}$ also enter the
formula (71), which can alternatively be further written as

$\displaystyle\left(\frac{1}{r}\frac{\partial}{\partial
r}\right)^{\ell}\left[\frac{F(t\mp
r)}{r}\right]=\frac{(-)^{\ell}}{r^{2\ell+1}}\sum_{i=0}^{\ell}\beta_{i}^{\ell}\,(\pm
r)^{i}\,F^{(i)}(t\mp r)\,.$ (77)

These two facts permit proving that the three terms defined in (75) are
actually all identical:

$\mathcal{A}_{L}=\mathcal{B}_{L}=\mathcal{C}_{L}\quad\text{(for any
{\hbox{t}}, {\hbox{s}}, {\hbox{r}}, {\hbox{{r'}}})}\,,$ (78)

so that we may glue together the three integrals in (74) and arrive at

$\displaystyle
h_{L}=\int_{-\infty}^{t-r}\mathrm{d}s\,\Biggl{\\{}\int_{0}^{\frac{t-r-s}{2}}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1}N({r^{\prime}},s+{r^{\prime}})\,\mathcal{A}_{L}-\int_{0}^{\frac{t+r-s}{2}}\\!\mathrm{d}{r^{\prime}}\,{r^{\prime}}^{B-\ell+1}N({r^{\prime}},s+{r^{\prime}})\,\mathcal{B}_{L}\Biggr{\\}}\,,$
(79)

together with the following summed-up expressions

$\displaystyle\mathcal{A}_{L}$
$\displaystyle=\frac{2^{\ell-1}}{\ell!}\hat{\partial}_{L}\\!\left[\frac{1}{r}\left(\frac{t-r-s}{2}\right)^{\ell}\left(\frac{t-r-s}{2}-{r^{\prime}}\right)^{\ell}\right]\,,$
(80a) $\displaystyle\mathcal{B}_{L}$
$\displaystyle=\frac{2^{\ell-1}}{\ell!}\hat{\partial}_{L}\\!\left[\frac{1}{r}\left(\frac{t+r-s}{2}\right)^{\ell}\left(\frac{t+r-s}{2}-{r^{\prime}}\right)^{\ell}\right]\,.$
(80b)

The fact that $\mathcal{A}_{L}$ and $\mathcal{B}_{L}$ are actually identical
follows immediately from Eq. (A36) in [32]. At last, we can commute the multi-
spatial derivative operator $\hat{\partial}_{L}$ in $\mathcal{A}_{L}$ and
$\mathcal{B}_{L}$ with the integration over ${r^{\prime}}$ in (79). Indeed,
the contributions coming from the differentiation of the bound $(t\mp r-s)/2$
are seen to be zero from the structure of $\mathcal{A}_{L}$ and
$\mathcal{B}_{L}$ respectively. So, finally, posing

$\displaystyle
R(\rho,s)=\rho^{\ell}\int_{0}^{\rho}\mathrm{d}{r^{\prime}}\frac{(\rho-{r^{\prime}})^{\ell}}{\ell!}\left(\frac{2}{{r^{\prime}}}\right)^{\ell-1}\\!{r^{\prime}}^{B}N({r^{\prime}},s+{r^{\prime}})\,,$
(81)

we get our final result in 3 dimensions,

$\displaystyle
h_{L}=\int_{-\infty}^{t-r}\mathrm{d}s\,\hat{\partial}_{L}\left[\frac{R\bigl{(}\frac{t-r-s}{2},s\bigr{)}-R\bigl{(}\frac{t+r-s}{2},s\bigr{)}}{r}\right]\,.$
(82)

This is exactly the result given by Eq. (6.4) in [32].

## 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] .
  * 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._ [2021] 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 (2021), companion paper.
  * 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).
  * 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] .
  * Foffa and Sturani [2012] S. Foffa and R. Sturani, Tail terms in gravitational radiation reaction via effective field theory, Phys. Rev. D 87, 044056 (2012), arXiv:1111.5488 [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] .
  * 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._ [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._ [2021] D. Bini, T. Damour, and A. Geralico, Radiative contributions to gravitational scattering, arXiv preprint arXiv:2107.08896 (2021).
  * Blanchet [1998a] L. Blanchet, Gravitational-wave tails of tails, Class. Quant. Grav. 15, 113 (1998a), gr-qc/9710038 .
  * 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] .
  * Blanchet [1998b] L. Blanchet, On the multipole expansion of the gravitational field, Class. Quant. Grav. 15, 1971 (1998b), gr-qc/9801101 .
  * 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 [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).
  * Thorne [1980] K. Thorne, Multipole expansions of gravitational radiation, Rev. Mod. Phys. 52, 299 (1980).
  * Blanchet [1998c] L. Blanchet, Quadrupole-quadrupole gravitational waves, Class. Quant. Grav. 15, 89 (1998c), gr-qc/9710037 .
  * 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] .
  * Almeida  _et al._ [2021] G. L. Almeida, S. Foffa, and R. Sturani, Gravitational Multipole Renormalization (2021), arXiv:2107.02634 [gr-qc] .
  * 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] .