# Gravitational-Wave Tails of Memory

David Trestini<EMAIL_ADDRESS>Laboratoire Univers et Théories,
Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, F-92190
Meudon, France  $\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  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

###### Abstract

Gravitational-wave tails are linear waves that backscatter on the curvature of
space-time generated by the total mass-energy of the source. The non-linear
memory effect arises from gravitational waves sourced by the stress-energy
distribution of linear waves themselves. These two effects are due to
quadratic multipolar interactions (mass-quadrupole and quadrupole-quadrupole)
and are well known. Also known are the tails generated by tails themselves
(cubic “tails-of-tails”) and the tails generated by tails-of-tails or _vice-
versa_ (quartic “tails-of-tails-of-tails”). In this work, we focus on the
cubic “tails-of-memory” corresponding to the mass-quadrupole-quadrupole
interaction, as well as the “spin-quadrupole tails”, which are due to the
cubic interaction between the mass, the total angular momentum and the
quadrupole. The tails-of-memory and the spin-quadrupole tails contribute to
the asymptotic waveform at the fourth-post-Newtonian (4PN) order beyond
quadrupolar radiation.

###### pacs:

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

##  I Introduction

Due to the non-linear nature of general relativity, non-linear effects in the
propagation of gravitational waves from their source to a distant detector
play an important role in predictions for gravitational waves generated by
compact binary systems [1, 2, 3, 4]. The prototype of such effects is the
gravitational-wave _tail_ , namely the backscattering of linear waves on the
curvature of space-time generated by the total mass-energy of the source [5,
6, 7, 8]. The tail effect has actually already been detected by the LIGO-Virgo
observations [9]. This constitutes an interesting test of the non-linear
structure of general relativity [10, 11, 12, 13].

At leading order the tail effect is a quadratic coupling between the Arnowitt-
Deser-Misner (ADM) total mass of space-time and the mass quadrupole moment of
the source. This effect appears at one-and-a-half post-Newtonian (1.5PN) order
in the waveform beyond the Einstein quadrupole formula [14, 15, 16, 17, 18,
19]. At 3PN order there is the _tail-of-tail_ effect which is a cubic coupling
between two ADM masses and the quadrupole [20, 21]. The previous effects have
been included into the PN templates for compact binary inspiral that are
routinely used by the LIGO-Virgo detectors. Furthermore, at quartic order
there is the _tail-of-tail-of-tail_ which arises at 4.5PN order and involves
three masses and a quadrupole [22, 23].

The non-linear memory effect, _i.e._ the permanent change in the wave
amplitude after the passage of a gravitational wave burst, is due to the re-
radiation of quadrupole gravitational waves by a linear quadrupole wave and
corresponds to the coupling between two quadrupole moments [24, 25, 26, 27,
14, 28, 29, 30, 31]. It enters at 2.5PN order in the waveform. The oscillatory
(“AC”) piece of this interaction is routinely used by the LIGO-Virgo detector,
but the genuine, secular memory effect (which is almost zero frequency or
“DC”) has not yet been detected experimentally, but could be observed in the
coming years with ground-based gravitational wave detectors [32, 33, 34].

In the present paper we shall proceed one step further and compute the cubic
non-linear interaction between the constant ADM mass $\mathrm{M}$ and two
quadrupole moments $\mathrm{M}_{ij}$, which we denote as
$\mathrm{M}\times\mathrm{M}_{ij}\times\mathrm{M}_{ij}$. Physically this can be
viewed as a combination between the tails produced by the memory and the
memory associated with the tail; we coin this effect “ _tails-of-memory_ ”.
From an effective field theory (EFT) perspective, the tails-of-memory
correspond to the three Feynman diagrams shown in Fig. 1, which we consider
here for illustrative purposes; see [35, 36, 37, 38] for their precise
computational meaning within the EFT framework.

(a)

(b)

(c)

Figure 1: Feynman diagrams associated to the tails-of-memory.

Interestingly, the tails-of-memory arise at 4PN order in the waveform, _i.e._
, at the level which is our overall goal currently [39, 40, 41, 42, 43, 44],
motivated by third generation detectors onground and by LISA in space. For
this reason, we also compute the only other cubic contribution that arises at
4PN order, namely the “ _spin-quadrupole tail_ ” interaction
$\mathrm{M}\times\mathrm{S}_{i}\times\mathrm{M}_{ij}$ between the constant ADM
mass, the constant ADM angular momentum and a quadrupole moment.

Note that the contributions due to the dimensional regularization of all non-
linear interactions relevant at 4PN order (in particular the tails-of-memory)
have already been included into the definition of a “renormalized” source
quadrupole moment, see [42] for details. Therefore we shall compute the tails-
of-memory and spin-quadrupole tails purely in three dimensions; adding them to
the renormalized source quadrupole moment [42] will naturally yield the
physical radiative moment which is observed at infinity.

Since the tails-of-memory arise at 4PN order in the waveform, and that
radiation reaction adds another 2.5PN order, we expect the tails-of-memory to
contribute in the form of radiation modes to the conservative equations of
motion and Lagrangian/Hamiltonian at 6.5PN order. This is far beyond the
state-of-the-art on equations of motion which is 4PN order [45, 46, 47, 48,
49, 50] (see also the results [51, 52, 53, 54] at 5PN order and [55] at 6PN
order).

We shall compute the tails-of-memory using the multipolar-post-Minkowskian
(MPM) formalism [56, 57] but with an important variant with respect to
previous works such as [28, 20]: we shall construct the metric directly in
_radiative_ coordinates (see [58, 59, 60] for a general definition), following
the method proposed in [57], instead of the usual construction in harmonic
coordinates. The great advantage is that we shall avoid the appearance of the
far zone logarithms that plague harmonic coordinates. A minor disadvantage is
that we shall have to apply a correction to account for the different
expressions of the multipole moments in radiative and harmonic coordinates,
but this correction is already known [44].

Finally, despite the fact that we are mainly motivated by the case of quasi-
circular compact binaries, the results presented in this paper are valid for a
generic isolated source described in terms of its multipole moments; the
source could be a general hydrodynamical fluid with compact support, and it
needs not even be a post-Newtonian source.

The plan of this paper is as follows. In Sec. II, we review the radiative
version of the MPM algorithm. In Sec. III, we present our results for the
radiative metric at quadratic order, and discuss the structure of the cubic
source. In Sec IV, we present a novel integration technique for the cubic
source, and apply it as a check to the known cases of tails-of-tails and
tails. Then, in Sec. V, we present our “raw” result in the form of complicated
functionals of the multipole moments, along with a set of coefficients given
by Table 1(d). We then introduce a method of simplification, which drastically
reduces the “raw” result, and allows for a test of our integration method.
Finally, we present in Sec. VI the final cubic radiative quadrupole moment (in
its simplified form) in terms of the harmonic canonical moments. Appendix A
presents a practical method for extracting multipole moments, and Appendix B
presents a proof for the convergence of our final integrals.

##  II The radiative MPM algorithm

The vacuum field equations of general relativity are recast in the form of a
wave like equation by defining the gothic metric deviation
$h^{\mu\nu}\equiv\sqrt{-g}g^{\mu\nu}-\eta^{\mu\nu}$, where $\eta^{\mu\nu}$ is
the flat Minkowski metric (in the mostly plus signature), $g^{\mu\nu}$ is the
inverse of the usual covariant metric and $g\equiv\det(g_{\mu\nu})$ is the
determinant. In a general coordinate system, the vacuum field equations become

$\Box h^{\mu\nu}-\partial H^{\mu\nu}=\Lambda^{\mu\nu}\bigl{[}h,\partial
h,\partial^{2}h\bigr{]}\,,$ (1)

where $\Box$ is the flat d’Alembertian operator, the divergence of the metric
is $H^{\mu}\equiv\partial_{\nu}h^{\mu\nu}$ and we introduce the useful
shorthand notation $\partial\theta^{\mu\nu}\equiv
2\partial^{(\mu}\theta^{\nu)}-\eta^{\mu\nu}\partial_{\rho}\theta^{\rho}$,
valid for any vector $\theta^{\mu}$. The non-linear source term in vacuum
reads explicitly

$\Lambda^{\mu\nu}=\Lambda_{\text{harm}}^{\mu\nu}+2h^{\rho(\mu}\partial_{\rho}H^{\nu)}-\partial_{\rho}\left(h^{\mu\nu}H^{\rho}\right)\,,$
(2)

where $\Lambda_{\text{harm}}^{\mu\nu}$ is the source term in a harmonic gauge
satisfying $H^{\mu}=0$ and is given explicitly by Eq. (24) of [3]. Taking the
divergence of (1), and using the fact that $\partial_{\nu}\partial
H^{\mu\nu}=\Box H^{\mu}$, we find that the source term is identically
divergenceless in any coordinate system, _i.e._ ,
$\partial_{\nu}\Lambda^{\mu\nu}=0$, which is equivalent to the contracted
Bianchi identity.

From now on, $h^{\mu\nu}$ will denote the metric in a radiative type
coordinate system [58, 59, 60], characterized by the absence of logarithms of
the radial distance $r$ in the far retarded zone, _i.e._ , when $r\to+\infty$
with $t-r=$ const. With the MPM construction the radiative metric is formally
decomposed to any post-Minkowskian order as

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

where each of the coefficients will be a functional of two sets of STF
multipole moments $\overline{\mathrm{M}}_{L}$ and $\overline{\mathrm{S}}_{L}$
which by definition parametrize the linearized approximation $h_{1}^{\mu\nu}$
defined as

$h^{\mu\nu}_{1}\equiv
h_{\text{can}\,1}^{\mu\nu}[\overline{\mathrm{M}}_{L},\overline{\mathrm{S}}_{L}]+\partial\xi_{1}^{\mu\nu}\,,$
(4)

where the “canonical” linear metric reads explicitly [61, 62, 7, 56]111We use
the symmetric-trace-free (STF) formalism, with $L=i_{1}\cdots i_{\ell}$ a
multi-index with $\ell$ spatial indices, and
$\hat{n}_{L}=\text{STF}[n_{i_{1}}\cdots n_{i_{\ell}}]$ the usual STF harmonics
of order $\ell$. Similarly we shall denote
$\hat{\partial}_{L}=\text{STF}[\partial_{i_{1}}\cdots\partial_{i_{\ell}}]$. We
generally pose $c=1$.

$\displaystyle
h_{\text{can}\,1}^{00}[\overline{\mathrm{M}}_{L},\overline{\mathrm{S}}_{L}]$
$\displaystyle=-4\sum_{\ell\geqslant
0}\frac{(-)^{\ell}}{\ell!}\partial_{L}\left[\frac{1}{r}\,\overline{\mathrm{M}}_{L}\left(t-r\right)\right]\,,$
(5a) $\displaystyle
h_{\text{can}\,1}^{0i}[\overline{\mathrm{M}}_{L},\overline{\mathrm{S}}_{L}]$
$\displaystyle=4\sum_{\ell\geqslant
1}\frac{(-)^{\ell}}{\ell!}\left\\{\partial_{L-1}\left[\frac{1}{r}\,\overline{\mathrm{M}}_{iL-1}^{(1)}\left(t-r\right)\right]+\frac{\ell}{\ell+1}\epsilon_{iab}\partial_{aL-1}\left[\frac{1}{r}\,\overline{\mathrm{S}}_{bL-1}\left(t-r\right)\right]\right\\}\,,$
(5b) $\displaystyle
h_{\text{can}\,1}^{ij}[\overline{\mathrm{M}}_{L},\overline{\mathrm{S}}_{L}]$
$\displaystyle=-4\sum_{\ell\geqslant
2}\frac{(-)^{\ell}}{\ell!}\left\\{\partial_{L-2}\left[\frac{1}{r}\,\overline{\mathrm{M}}_{ijL-2}^{(2)}\left(t-r\right)\right]+\frac{2\ell}{\ell+1}\partial_{aL-2}\left[\frac{1}{r}\,\epsilon_{ab(i}\overline{\mathrm{S}}_{j)bL-2}^{(1)}\left(t-r\right)\right]\right\\}\,.$
(5c)

It satisfies the harmonic gauge condition,
$\partial_{\nu}h_{\text{can}\,1}^{\mu\nu}=0$. The linear gauge transformation
vector in (4) is

$\xi_{1}^{\mu}=2M\eta^{0\mu}\ln\left(\frac{r}{b_{0}}\right)\,,$ (6)

with $b_{0}$ an arbitrary constant length scale. The role of this gauge
transformation is to ensure that the retarded time $u=t-r$ is asymptotically a
null coordinate when $r\to+\infty$.

The multipole moments $\overline{\mathrm{M}}_{L}$ and
$\overline{\mathrm{S}}_{L}$ will depend specifically of the radiative MPM
construction that we shall perform. In particular, as we have shown in [44],
the moments $\overline{\mathrm{M}}_{L}$ and $\overline{\mathrm{S}}_{L}$ will
differ from the set of multipole moments $\mathrm{M}_{L}$ and $\mathrm{S}_{L}$
defined in the MPM construction in harmonic coordinates, so we shall have to
correct for that difference using the results of [44]. The moments in harmonic
coordinates $\mathrm{M}_{L}$ and $\mathrm{S}_{L}$ can be directly related the
to source _via_ asymptotic matching, and are known to 4PN order [39, 41, 42,
43]. Consistently the relation between the moments
$\\{\overline{\mathrm{M}}_{L},\overline{\mathrm{S}}_{L}\\}$ and the harmonic
coordinate moments $\\{\mathrm{M}_{L},\mathrm{S}_{L}\\}$ has been
investigated, and in the case of the quadrupole moment, determined to 4PN
order in [44]; we recall the result in Eq. (92) below. Note also that the
conserved quantities are equal in both coordinate systems, _i.e._ ,
$\mathrm{M}=\overline{\mathrm{M}}$ and
$\mathrm{S}_{i}=\overline{\mathrm{S}}_{i}$, so in these cases we will use both
notations interchangeably.

Once the linearized solution is defined, the non-linear corrections are
readily obtained by injecting the PM expansion into the field equations and
solving these equations iteratively over the PM order. At any order $n$, the
general equation to solve is

$\Box h_{n}^{\mu\nu}-\partial H^{\mu\nu}_{n}=\Lambda_{n}^{\mu\nu}\,,$ (7)

where $\Lambda_{n}^{\mu\nu}\equiv\Lambda^{\mu\nu}[h_{1},\cdots,h_{n-1}]$ is
built out from previous iterations. We first construct a particular retarded
solution to the wave equation, satisfying $\Box
u_{n}^{\mu\nu}=\Lambda_{n}^{\mu\nu}$, as

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

where $\Box^{-1}_{\text{ret}}$ is the standard retarded integral operator, and
a particular regularization with regulator $(r/r_{0})^{B}$, where $B$ is
complex and $r_{0}$ a constant scale, is introduced to cope with the
singularity of the multipole expansion when $r\to 0$. The finite part (FP)
when $B\to 0$, _i.e._ , the coefficient of the zeroth power of $B$ in the
Laurent expansion when $B\to 0$, provides the requested solution in this
context.

Next we compute its divergence
$w_{n}^{\mu}\equiv\partial_{\nu}u_{n}^{\mu\nu}$, and from that divergence we
are able to construct an homogeneous solution $v_{n}^{\mu\nu}$ to the wave
solution, which satisfies at once $\Box v_{n}^{\mu\nu}=0$ and
$\partial_{\nu}v_{n}^{\mu\nu}=-w_{n}^{\mu}$. We describe in Appendix A the
practical method which permits to obtain $v_{n}^{\mu\nu}$ starting from
$w_{n}^{\mu}$. At this point it is clear that the sum
$s_{n}^{\mu\nu}=u_{n}^{\mu\nu}+v_{n}^{\mu\nu}$ satisfies the field equation
(7), since we have $\Box s_{n}^{\mu\nu}=\Lambda_{n}^{\mu\nu}$ and
$\partial_{\nu}s_{n}^{\mu\nu}=0$. However in order to obtain the sought-for
radiative metric, we still have to correct for the light cone deviation to
order $n$, and for that we apply a linear gauge transformation parametrized by
a vector $\xi_{n}^{\mu}$, hence defining

$h_{n}^{\mu\nu}\equiv
u_{n}^{\mu\nu}+v_{n}^{\mu\nu}+\partial\xi_{n}^{\mu\nu}\,.$ (9)

The construction of this gauge vector crucially relies on the particular
structure of the source term, proved by induction over the order $n$, _i.e._ ,
an asymptotic expansion when $r\to+\infty$ with $t-r=$ const in simple powers
of $1/r$ (without logarithms) whose leading piece $1/r^{2}$ takes the
particular form

$\Lambda_{n}^{\mu\nu}=\frac{k^{\mu}k^{\nu}}{r^{2}}\sigma_{n}(u,\mathbf{n})+{\mathcal{O}}\left(\frac{1}{r^{3}}\right)\,,$
(10)

where $k^{\mu}=(1,\mathbf{n})$ is the outgoing Minkowskian vector. An
appropriate choice (but by no means unique) for the gauge vector is

$\xi_{n}^{\mu\nu}=\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}\left(\frac{r}{r_{0}}\right)^{B}\frac{k^{\mu}}{2r^{2}}\int_{0}^{+\infty}\mathrm{d}\tau\,\sigma_{n}(u-\tau,\mathbf{n})\biggr{]}\,.$
(11)

With this choice of gauge vector, one can prove by induction that the metric
is free of far-zone logarithms at all post-Minkowskian orders [57, 44], and
thus falls into the class of radiative coordinate systems for which the
extraction of physical observables at infinity is easy. The implementation of
the above construction for the case of the tails-of-memory will indeed prove
very important.

##  III Results at quadratic order

In the present paper we shall compute from first principles the following
three _cubic_ non-linear multipole interactions (i.e., $n=3$ in the PM
expansion):

$\mathrm{M}^{2}\times\overline{\mathrm{M}}_{ij}\,,\qquad\d{M}\times\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}\quad\text{and}\quad\d{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}\,.$

The first one is known as the quadrupole “tail-of-tail” [20], the second one
can be called the “spin-quadrupole tail” while the third one is the quadrupole
“tail-of-memory” interaction which is our main goal here. To reach these goals
we evidently need to control the following quadratic interactions:

$\mathrm{M}\times\mathrm{M}\,,\qquad\mathrm{M}\times\mathrm{S}_{i}\,,\qquad\mathrm{M}\times\overline{\mathrm{M}}_{ij}\,,\qquad\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}\quad\text{and}\quad\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}\,,$

which we shall compute at any distance $r$ from the source, while the cubic
interactions will just be computed asymptotically when $r\to+\infty$.

###  III.1 Linear and quadratic results

The linear interactions $\mathrm{M}$, $\mathrm{S}_{i}$ and
$\overline{\mathrm{M}}_{ij}$ are read off directly from the linear metric
(4)–(5), including the gauge transformation (6):

$\displaystyle h_{\mathrm{M}}^{00}$ $\displaystyle=-\frac{4\mathrm{M}}{r}\,,$
$\displaystyle h_{\mathrm{M}}^{0i}$
$\displaystyle=-\frac{2\mathrm{M}n^{i}}{r}\,,$ $\displaystyle
h_{\mathrm{M}}^{ij}$ $\displaystyle=0\,,$ (12a) $\displaystyle
h_{\mathrm{S}_{i}}^{00}$ $\displaystyle=0\,,$ $\displaystyle
h_{\mathrm{S}_{i}}^{0i}$
$\displaystyle=\frac{2n^{a}\mathrm{S}_{i|a}}{r^{2}}\,,$ $\displaystyle
h_{\mathrm{S}_{i}}^{ij}$ $\displaystyle=0\,,$ (12b) $\displaystyle
h_{\overline{\mathrm{M}}_{ij}}^{00}$
$\displaystyle=-2\partial_{ab}\biggl{[}\frac{\overline{\mathrm{M}}_{ij}}{r}\biggr{]}\,,$
$\displaystyle h_{\overline{\mathrm{M}}_{ij}}^{0i}$
$\displaystyle=2\partial_{a}\biggl{[}\frac{\overline{\mathrm{M}}_{ia}^{(1)}}{r}\biggr{]}\,,$
$\displaystyle h_{\overline{\mathrm{M}}_{ij}}^{ij}$
$\displaystyle=-\frac{2\overline{\mathrm{M}}_{ij}^{(2)}}{r}\,,$ (12c)

where we have henceforth posed
$\mathrm{S}_{i|j}\equiv\epsilon_{ija}\mathrm{S}_{a}$.

In the harmonic construction, the quadratic metrics that are known are
$\mathrm{M}\times\mathrm{M}$ [20], $\mathrm{M}\times\mathrm{M}_{ij}$ [14] and
$\mathrm{M}_{ij}\times\mathrm{M}_{ij}$ [28]. We also know the
$\mathrm{S}_{i}\times\mathrm{M}_{ij}$ metric asymptotically, to leading order
in $1/r$ [28].

The radiative metrics for the static $\mathrm{M}\times\mathrm{M}$ and
stationary $\mathrm{M}\times\mathrm{S}_{i}$ interactions, that are
respectively part of the Schwarzschild and Kerr metrics (in “radiative”
coordinates), read

$\displaystyle h_{\mathrm{M}\times\mathrm{M}}^{00}$
$\displaystyle=-\frac{3\mathrm{M}^{2}}{r^{2}}\,,$ $\displaystyle
h_{\mathrm{M}\times\mathrm{M}}^{0i}$ $\displaystyle=0\,,$ $\displaystyle
h_{\mathrm{M}\times\mathrm{M}}^{ij}$
$\displaystyle=-\frac{\mathrm{M}^{2}n^{ij}}{r^{2}}\,,$ (13a) $\displaystyle
h_{\mathrm{M}\times\mathrm{S}_{i}}^{00}$ $\displaystyle=0\,,$ $\displaystyle
h_{\mathrm{M}\times\mathrm{S}_{i}}^{0i}$
$\displaystyle=\frac{2\mathrm{M}}{r^{3}}n_{a}\mathrm{S}_{i|a}\,,$
$\displaystyle h_{\mathrm{M}\times\mathrm{S}_{i}}^{ij}$ $\displaystyle=0\,.$
(13b)

Note that although $\overline{\mathrm{M}}=\mathrm{M}$ as we said, the
$\mathrm{M}\times\mathrm{M}$ radiative metric differs from the harmonic one,
see Eqs. (2.7) in [20]. In the radiative construction, the
$\mathrm{M}\times\overline{\mathrm{M}}_{ij}$ interaction is given in [44] and
we reproduce it here for completeness:

$\displaystyle h_{\mathrm{M}\times\overline{\mathrm{M}}_{ij}}^{00}$
$\displaystyle=8\mathrm{M}n_{ab}\int_{1}^{+\infty}\mathrm{d}x\,\overline{Q}_{2}(x,r)\,\overline{\mathrm{M}}^{(4)}_{ab}(t-rx)$
$\displaystyle\quad+\frac{\mathrm{M}n_{ab}}{r}\Big{(}\frac{118}{15}\overline{\mathrm{M}}^{(3)}_{ab}+\frac{23}{5}r^{-1}\overline{\mathrm{M}}^{(2)}_{ab}-\frac{117}{5}r^{-2}\overline{\mathrm{M}}^{(1)}_{ab}-21r^{-3}\overline{\mathrm{M}}_{ab}\Big{)}\,,$
(14a) $\displaystyle h_{\mathrm{M}\times\overline{\mathrm{M}}_{ij}}^{0i}$
$\displaystyle=8\mathrm{M}n_{a}\int_{1}^{+\infty}\mathrm{d}x\,\overline{Q}_{1}(x,r)\,\overline{\mathrm{M}}^{(4)}_{ai}(t-rx)$
$\displaystyle\quad+\frac{\mathrm{M}n_{iab}}{r}\bigg{[}\overline{\mathrm{M}}^{(3)}_{ab}+3r^{-1}\overline{\mathrm{M}}^{(2)}_{ab}-r^{-2}\overline{\mathrm{M}}^{(1)}_{ab}\bigg{]}$
$\displaystyle\quad+\frac{\mathrm{M}n_{a}}{r}\bigg{[}\frac{43}{15}\overline{\mathrm{M}}^{(3)}_{ai}-\frac{107}{15}r^{-1}\overline{\mathrm{M}}^{(2)}_{ai}-5r^{-2}\overline{\mathrm{M}}^{(1)}_{ai}\bigg{]}\,,$
(14b) $\displaystyle h_{\mathrm{M}\times\overline{\mathrm{M}}_{ij}}^{ij}$
$\displaystyle=8\mathrm{M}\int_{1}^{+\infty}\mathrm{d}x\,\overline{Q}_{0}(x,r)\,\overline{\mathrm{M}}^{(4)}_{ij}(t-rx)$
$\displaystyle\quad+\frac{\mathrm{M}n_{ijab}}{r}\bigg{[}-\frac{1}{2}\overline{\mathrm{M}}^{(3)}_{ab}-3r^{-1}\overline{\mathrm{M}}^{(2)}_{ab}-\frac{15}{2}r^{-2}\overline{\mathrm{M}}^{(1)}_{ab}-\frac{15}{2}r^{-3}\overline{\mathrm{M}}_{ab}\bigg{]}$
$\displaystyle\quad+\frac{\mathrm{M}\delta_{ij}n_{ab}}{r}\bigg{[}-\frac{1}{2}\overline{\mathrm{M}}^{(3)}_{ab}-2r^{-1}\overline{\mathrm{M}}^{(2)}_{ab}-\frac{1}{2}r^{-2}\overline{\mathrm{M}}^{(1)}_{ab}-\frac{1}{2}r^{-3}\overline{\mathrm{M}}_{ab}\bigg{]}$
$\displaystyle\quad+\frac{\mathrm{M}n_{a(i}}{r}\bigg{[}4\overline{\mathrm{M}}^{(3)}_{j)a}+6r^{-1}\overline{\mathrm{M}}^{(2)}_{j)a}+6r^{-2}\overline{\mathrm{M}}^{(1)}_{j)a}+6r^{-3}\overline{\mathrm{M}}_{j)a}\bigg{]}$
$\displaystyle\quad+\frac{\mathrm{M}}{r}\bigg{[}-\frac{107}{15}\overline{\mathrm{M}}^{(3)}_{ij}-4r^{-1}\overline{\mathrm{M}}^{(2)}_{ij}-r^{-2}\overline{\mathrm{M}}^{(1)}_{ij}-r^{-3}\overline{\mathrm{M}}_{ij}\bigg{]}\,.$
(14c)

The tail terms in (14) involve the modified Legendre function

$\overline{Q}_{m}(x,r)\equiv
Q_{m}(x)-\frac{1}{2}P_{m}(x)\ln\left(\frac{r}{r_{0}}\right)\,,$ (15)

which is constructed using the usual Legendre polynomial $P_{m}(x)$ and the
usual Legendre function of second kind (with branch cut on $]-\infty,1]$),
namely

$Q_{m}(x)\equiv\frac{1}{2}\int_{-1}^{1}\mathrm{d}y\,\frac{P_{m}(y)}{x-y}=\frac{1}{2}P_{m}(x)\ln\left(\frac{x+1}{x-1}\right)-\sum_{j=1}^{m}\frac{1}{m}P_{m-j}(x)P_{j-1}(x)\,.$
(16)

It was shown in [44] that the combination (15) does not produce any far-zone
logarithms in the radiative metric, although the tail term is still non-local
in time and depends on the constant scale $r_{0}$. As shown in [44] the metric
(14) differs from its counterpart in harmonic coordinates (given in Appendix B
in [14]) by a coordinate transformation and a redefinition of the quadrupole
moment [44], see (92) below.

The radiative metric for the $\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}$
interaction is new with this paper and reads

$\displaystyle h_{\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}}^{00}$
$\displaystyle=n_{ad}\,\mathrm{S}_{a|b}\left(-\frac{4}{3}r^{-1}\overline{\mathrm{M}}_{bd}^{(4)}-4r^{-2}\overline{\mathrm{M}}_{bd}^{(3)}-6r^{-3}\overline{\mathrm{M}}_{bd}^{(2)}-6r^{-4}\overline{\mathrm{M}}_{bd}^{(1)}\right)\,,$
(17a) $\displaystyle h_{\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}}^{0i}$
$\displaystyle=n_{acd}\,\mathrm{S}_{i|a}\left(-\frac{5}{6}r^{-1}\overline{\mathrm{M}}_{cd}^{(4)}-5r^{-2}\overline{\mathrm{M}}_{cd}^{(3)}-\frac{23}{2}r^{-3}\overline{\mathrm{M}}_{cd}^{(2)}-\frac{3}{2}r^{-4}\overline{\mathrm{M}}_{cd}^{(1)}-\frac{3}{2}r^{-5}\overline{\mathrm{M}}_{cd}\right)$
$\displaystyle+n_{iad}\,\mathrm{S}_{a|b}\left(-3r^{-3}\overline{\mathrm{M}}_{bd}^{(2)}-9r^{-4}\overline{\mathrm{M}}_{bd}^{(1)}-9r^{-5}\overline{\mathrm{M}}_{bd}\right)$
$\displaystyle+n_{a}\,\mathrm{S}_{a|b}\left(-\frac{4}{3}r^{-1}\overline{\mathrm{M}}_{ib}^{(4)}-\frac{4}{3}r^{-2}\overline{\mathrm{M}}_{ib}^{(3)}+r^{-3}\overline{\mathrm{M}}_{ib}^{(2)}+r^{-4}\overline{\mathrm{M}}_{ib}^{(1)}+r^{-5}\overline{\mathrm{M}}_{ib}\right)$
$\displaystyle+n_{a}\,\mathrm{S}_{i|b}\left(\frac{5}{3}r^{-2}\overline{\mathrm{M}}_{ab}^{(3)}+5r^{-3}\overline{\mathrm{M}}_{ab}^{(2)}+4r^{-4}\overline{\mathrm{M}}_{ab}^{(1)}+4r^{-5}\overline{\mathrm{M}}_{ab}\right)\,,$
(17b) $\displaystyle h_{\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}}^{ij}$
$\displaystyle=n_{ijad}\,\mathrm{S}_{a|b}\left(\frac{1}{3}r^{-1}\overline{\mathrm{M}}_{bd}^{(4)}+\frac{10}{3}r^{-2}\overline{\mathrm{M}}_{bd}^{(3)}+9r^{-3}\overline{\mathrm{M}}_{bd}^{(2)}+9r^{-4}\overline{\mathrm{M}}_{bd}^{(1)}\right)$
$\displaystyle+n_{acd(i}\,\mathrm{S}_{j)|a}\left(\frac{1}{3}r^{-1}\overline{\mathrm{M}}_{cd}^{(4)}+\frac{10}{3}r^{-2}\overline{\mathrm{M}}_{cd}^{(3)}+9r^{-3}\overline{\mathrm{M}}_{cd}^{(2)}+9r^{-4}\overline{\mathrm{M}}_{cd}^{(1)}\right)$
$\displaystyle+\delta_{ij}n_{ad}\,\mathrm{S}_{a|b}\left(r^{-1}\overline{\mathrm{M}}_{bd}^{(4)}+\frac{8}{3}r^{-2}\overline{\mathrm{M}}_{bd}^{(3)}-r^{-3}\overline{\mathrm{M}}_{bd}^{(2)}-r^{-4}\overline{\mathrm{M}}_{bd}^{(1)}\right)$
$\displaystyle+n_{c(i}\,\mathrm{S}_{j)|a}\left(-\frac{4}{3}r^{-2}\overline{\mathrm{M}}_{ac}^{(3)}-4r^{-3}\overline{\mathrm{M}}_{ac}^{(2)}-4r^{-4}\overline{\mathrm{M}}_{ac}^{(1)}\right)$
$\displaystyle+n_{a(i}\,\mathrm{S}_{a|b}\left(-\frac{8}{3}r^{-1}\overline{\mathrm{M}}_{j)b}^{(4)}-\frac{26}{3}r^{-2}\overline{\mathrm{M}}_{j)b}^{(3)}-6r^{-3}\overline{\mathrm{M}}_{j)b}^{(2)}-6r^{-4}\overline{\mathrm{M}}_{j)b}^{(1)}\right)$
$\displaystyle+n_{ac}\,\mathrm{S}_{a|(i}\left(-2r^{-1}\overline{\mathrm{M}}_{j)b}^{(4)}-\frac{20}{3}r^{-2}\overline{\mathrm{M}}_{j)b}^{(3)}-2r^{-3}\overline{\mathrm{M}}_{j)b}^{(2)}-2r^{-4}\overline{\mathrm{M}}_{j)b}^{(1)}\right)$
$\displaystyle+\,\mathrm{S}_{a|(i}\left(\frac{14}{3}r^{-2}\overline{\mathrm{M}}_{j)a}^{(3)}+2r^{-3}\overline{\mathrm{M}}_{j)a}^{(2)}+2r^{-4}\overline{\mathrm{M}}_{j)a}^{(1)}\right)\,.$
(17c)

Note that contrary to the $\mathrm{M}\times\overline{\mathrm{M}}_{ij}$
interaction in (14), the $\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}$
interaction is purely “instantaneous” or local, as it does not contain tail
integrals. Note also that in the particular case of the
$\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}$ interaction, the expressions
for the radiative and harmonic metrics as functionals of their respective
canonical moments are identical, and given by (17).

The radiative metric for the interaction
$\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}$ is too lengthy to
be presented. It is constructed following Eq. (9) as

$h_{\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}^{\mu\nu}=u_{\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}^{\mu\nu}+v_{\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}^{\mu\nu}+\partial\xi_{\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}^{\mu\nu}\,.$
(18)

The gauge vector is defined by Eqs. (10)–(11) and we have explicitly in this
case [44]

$\displaystyle\sigma_{\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}=\hat{n}_{ijab}\,\overline{\mathrm{M}}_{ij}^{(3)}\overline{\mathrm{M}}_{ab}^{(3)}-\frac{24}{7}\hat{n}_{ij}\,\overline{\mathrm{M}}_{ia}^{(3)}\overline{\mathrm{M}}_{ja}^{(3)}+\frac{4}{5}\overline{\mathrm{M}}_{ab}^{(3)}\overline{\mathrm{M}}_{ab}^{(3)}\,,$
(19a) and therefore
$\displaystyle\xi_{\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}^{0}$
$\displaystyle=\int_{-\infty}^{u}\mathrm{d}v\\!\int_{1}^{+\infty}\mathrm{d}x\biggl{\\{}-\frac{1}{2}\hat{n}_{ijab}\,Q_{4}(x)\overline{\mathrm{M}}_{ij}^{(3)}\overline{\mathrm{M}}_{ab}^{(3)}+\frac{12}{7}\hat{n}_{ij}\,Q_{2}(x)\overline{\mathrm{M}}_{ia}^{(3)}\overline{\mathrm{M}}_{ja}^{(3)}$
$\displaystyle\qquad\qquad\qquad\qquad\quad-\frac{2}{5}\,Q_{0}(x)\overline{\mathrm{M}}_{ab}^{(3)}\overline{\mathrm{M}}_{ab}^{(3)}\biggr{\\}}\,,$
(19b)
$\displaystyle\xi_{\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}^{i}$
$\displaystyle=\int_{-\infty}^{u}\mathrm{d}v\\!\int_{1}^{+\infty}\mathrm{d}x\biggl{\\{}-\frac{1}{2}\hat{n}_{iabkl}\,Q_{5}(x)\overline{\mathrm{M}}_{ab}^{(3)}\overline{\mathrm{M}}_{kl}^{(3)}+\frac{16}{9}\hat{n}_{iab}\,Q_{3}(x)\overline{\mathrm{M}}_{ak}^{(3)}\overline{\mathrm{M}}_{kb}^{(3)}$
$\displaystyle\qquad\qquad\qquad\qquad\quad-\frac{2}{9}\hat{n}_{abk}\,Q_{3}(x)\overline{\mathrm{M}}_{ai}^{(3)}\overline{\mathrm{M}}_{bk}^{(3)}-\frac{22}{35}\hat{n}_{i}\,Q_{1}(x)\overline{\mathrm{M}}_{ab}^{(3)}\overline{\mathrm{M}}_{ab}^{(3)}$
$\displaystyle\qquad\qquad\qquad\qquad\quad+\frac{24}{35}\hat{n}_{a}\,Q_{1}(x)\overline{\mathrm{M}}_{ik}^{(3)}\overline{\mathrm{M}}_{ka}^{(3)}\biggr{\\}}\,.$
(19c)

Crucially, thanks to this gauge transformation the latter quadratic metric is
free of any far-zone logarithms and so will be the cubic source built out of
it. This allows us to take the dominant asymptotic behavior of the quadratic
radiative metric as $r\to+\infty$, and extract after standard transverse-
traceless (TT) projection of the spatial metric the associated radiative
quadrupole moments (where we reintroduce $c$ and $G$):

$\displaystyle\mathcal{U}_{ij}^{\mathrm{M}\times\overline{\mathrm{M}}_{ij}}$
$\displaystyle=\frac{2G\mathrm{M}}{c^{3}}\int_{0}^{+\infty}\mathrm{d}\tau\,\bigg{[}\ln\left(\frac{\tau}{2r_{0}}\right)+\frac{107}{60}\bigg{]}\overline{\mathrm{M}}^{(4)}_{ij}(u-\tau)\,,$
(20a)
$\displaystyle\mathcal{U}_{ij}^{\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}}$
$\displaystyle=-\frac{G}{3c^{5}}\,\mathrm{S}_{a|\langle
i}\overline{\mathrm{M}}_{j\rangle a}^{(4)}\,,$ (20b)
$\displaystyle\mathcal{U}_{ij}^{\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}$
$\displaystyle=-\frac{2G}{7c^{5}}\left[\int_{0}^{+\infty}\mathrm{d}\tau\,\overline{\mathrm{M}}^{(3)}_{a\langle
i}\overline{\mathrm{M}}^{(3)}_{j\rangle
a}(u-\tau)+\overline{\mathrm{M}}^{(2)}_{a\langle
i}\overline{\mathrm{M}}^{(3)}_{j\rangle
a}+\frac{5}{2}\overline{\mathrm{M}}^{(1)}_{a\langle
i}\overline{\mathrm{M}}^{(4)}_{j\rangle
a}-\frac{1}{2}\overline{\mathrm{M}}_{a\langle
i}\overline{\mathrm{M}}^{(5)}_{j\rangle a}\right]\,.$ (20c)

The expression for the tail differs (but is physically equivalent [44]) from
its counterpart in the harmonic construction [14]; the expressions for
$\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}$ and
$\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}$ are the same in
both radiative and harmonic constructions, see the Appendix B of [28]. In
particular the first term in (20c) is the usual memory effect (which contains
both AC and DC contributions) at the dominant 2.5PN order.

###  III.2 Structure of the cubic source

The general equation we need to solve, Eq. (7), reads at cubic order

$\Box h_{3}^{\mu\nu}-\partial
H^{\mu\nu}_{3}=\Lambda^{\mu\nu}_{3}\bigl{[}h_{1},h_{2}\bigr{]}=N^{\mu\nu}\bigl{[}h_{1},h_{2}\bigr{]}+N^{\mu\nu}\bigl{[}h_{2},h_{1}\bigr{]}+M^{\mu\nu}\bigl{[}h_{1},h_{1},h_{1}\bigr{]}\,,$
(21)

where we recall the definition
$H_{3}^{\mu}\equiv\partial_{\nu}h_{3}^{\mu\nu}$. The quadratic and cubic
functionals $N^{\mu\nu}_{\text{harm}}[h,h]$ and
$M^{\mu\nu}_{\text{harm}}[h,h,h]$ in harmonic coordinates are given by
(1.3)–(1.4) of [20], and to obtain their radiative counterparts
$N^{\mu\nu}[h,h]$ and $M^{\mu\nu}[h,h,h]$ we apply (2), where the corrections
terms are purely quadratic and thus only affect $N^{\mu\nu}[h,h]$. The source
term for the cubic interaction
$\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}$
reads, with obvious notation,

$\Lambda^{\mu\nu}_{\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}=N^{\mu\nu}_{\overline{\mathrm{M}}_{ij}\times(\mathrm{M}\times\overline{\mathrm{M}}_{ij})}+N^{\mu\nu}_{\mathrm{M}\times(\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij})}+M^{\mu\nu}_{\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}\,.$
(22)

These three terms have the following structure:

  * 1.

$\overline{\mathrm{M}}_{ij}\times(\mathrm{M}\times\overline{\mathrm{M}}_{ij})$
corresponds to a linearized quadrupolar wave interacting with a quadratic
tail, as illustrated by the first Feynman diagram in Fig. 1(a). Apart from
simple instantaneous terms, this source contains terms whose structure reads

$\displaystyle\sim\frac{\hat{n}_{L}}{r^{5-p-q}}\,\mathrm{M}\,\overline{\mathrm{M}}_{ab}^{(p)}(t-r)\int_{1}^{+\infty}\mathrm{d}x\,\overline{Q}_{m}(x)\overline{\mathrm{M}}_{cd}^{(4+q)}(t-rx)\,,$
(23)

where $(p,q)\in[\\![0,2]\\!]$, and where the index structure and the possible
presence of Kronecker deltas are disregarded. These terms are novel (with
respect to ordinary tails), and their integration will be the focus of Sec.
IV.2. For integration purposes, the modified Legendre function
$\overline{Q}_{m}(x)$ will be decomposed according to (15): this will yield
similar integrals over the standard Legendre function $Q_{m}(x)$, augmented by
instantaneous terms involving $\ln\left(r/r_{0}\right)$.

  * 2.

$\mathrm{M}\times(\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij})$
corresponds to a quadrupole-quadrupole memory-like wave scattering off the
curvature generated by the total mass, and is illustrated by Fig. 1(b). This
source contains terms with the following structure:

$\displaystyle\sim\frac{\hat{n}_{L}}{r^{3-p-q}}\,\mathrm{M}\int_{1}^{+\infty}\mathrm{d}x\,\overline{Q}_{m}(x)\bigl{(}\overline{\mathrm{M}}_{ab}^{(3+p)}\,\overline{\mathrm{M}}_{cd}^{(3+q)}\bigr{)}(t-rx)\,,$
(24a) where $0\leqslant p+q\leqslant 2$, which have the same structure as
ordinary tails, as well as memory-like terms
$\displaystyle\sim\frac{\hat{n}_{L}}{r^{4}}\,\mathrm{M}\int_{0}^{+\infty}\mathrm{d}\tau\,\bigl{(}\overline{\mathrm{M}}_{ab}^{(3)}\,\overline{\mathrm{M}}_{cd}^{(3)}\bigr{)}(t-r-\tau)\,.$
(24b)

The modified Legendre function will again be decomposed according to (15), and
all these terms will be the treated in Sec. IV.3.

  * 3.

$\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}$
properly speaking, illustrated by Fig. 1(c). This piece only leads to
instantaneous terms, which are easy to integrate.

A similar analysis can be done for the other two cubic interactions
$\mathrm{M}\times\mathrm{M}\times\overline{\mathrm{M}}_{ij}$ and
$\mathrm{M}\times\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}$, which only
contain tail-like and instantaneous terms. Note that, thanks to the radiative
construction, we can see that the $r\to+\infty$ expansion of the cubic source
does not contain any logarithms of $r$. We have explicitly verified that the
novel source terms
$\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}$
and $\mathrm{M}\times\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}$ are
divergenceless, _i.e._

$\partial_{\nu}\Lambda^{\mu\nu}_{\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}=0\qquad\text{and}\qquad\partial_{\nu}\Lambda^{\mu\nu}_{\mathrm{M}\times\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}}=0\,,$
(25)

which is a strong test, since it would generally fail if a single coefficient
were incorrect. However, because of its length, we cannot present the cubic
source. Instead, we provide it in the form of an ancillary file in the
Supplementary Material [63].

##  IV Solution of the wave equation at cubic order

###  IV.1 General multipolar solution

To solve the problem of tails-of-memory we need to develop new techniques for
integrating the wave equation when the source term is a complicated non-local
functional of the moments such as (23). We first consider in a general way the
wave equation whose source term admits a definite multipolarity $\ell$ in STF
guise, thus takes the general form

$\Box\Psi_{L}=\hat{n}_{L}\,S(r,t-r)\,,$ (26)

where $S(r,u)$ is an arbitrary function of $r=|\mathbf{x}|$ and $u=t-r$ that
verifies straightforward smoothness properties, and tends sufficiently rapidly
to zero when $r\to 0$, namely

$S(r,u)=\mathcal{O}(r^{\ell+5})\quad\text{(when {\hbox{r\to 0}} with
{\hbox{u}} or {\hbox{t}} held fixed)}\,.$ (27)

See Theorem 6.1 of [56] for more details on the required conditions we impose.
Under these conditions, we know how to solve the wave equation (26). We first
define

$R_{\alpha}(\rho,s)\equiv\rho^{\ell}\int_{\alpha}^{\rho}\mathrm{d}\lambda\,\frac{(\rho-\lambda)^{\ell}}{\ell!}\left(\frac{2}{\lambda}\right)^{\ell-1}S(\lambda,s)\,,$
(28)

where $\alpha$ is an arbitrary constant. Then the solution of (26) can be
written as

$\Psi_{L}=\int_{-\infty}^{t-r}\mathrm{d}s\,\hat{\partial}_{L}\left[\frac{R_{\alpha}\left(\frac{t-r-s}{2},s\right)-R_{\alpha}\left(\frac{t+r-s}{2},s\right)}{r}\right]\,,$
(29)

see Eq. (6.4) in [56]. This solution is the unique retarded solution of (26),
thus satisfying the no-incoming radiation condition when $r\to+\infty$ with
$t+r=$ const. Furthermore it is independent from the constant $\alpha$. To see
this we consider separately the two terms in (29), which we decompose as
$\Psi_{L}=\Psi^{[1]}_{L}+\Psi^{[2]}_{L}$ with

$\displaystyle\Psi^{[1]}_{L}$
$\displaystyle\equiv\int_{-\infty}^{t-r}\mathrm{d}s\,\hat{\partial}_{L}\left[\frac{1}{r}R_{\alpha}\left(\frac{t-r-s}{2},s\right)\right]\,,$
(30a) $\displaystyle\Psi^{[2]}_{L}$
$\displaystyle\equiv-\int_{-\infty}^{t-r}\mathrm{d}s\,\hat{\partial}_{L}\left[\frac{1}{r}R_{\alpha}\left(\frac{t+r-s}{2},s\right)\right]\,.$
(30b)

We see that $\Psi^{[1]}_{L}$ is a retarded homogeneous solution of the wave
equation, $\Box\Psi^{[1]}_{L}=0$, but that $\Psi^{[2]}_{L}$ has a unwieldy
dependency on the advanced time $v=t+r$, and represents a particular retarded
solution of the wave equation, _i.e._ ,
$\Box\Psi^{[2]}_{L}=\hat{n}_{L}\,S(r,t-r)$. However note that the latter two
wave equations satisfied by $\Psi^{[1]}_{L}$ and $\Psi^{[2]}_{L}$ separately
are valid only in the sense of ordinary functions; only the wave equation (26)
for the full solution $\Psi_{L}$ is satisfied in the sense of distributions.

Plugging Eq. (28) into (30) we obtain

$\displaystyle\Psi^{[1]}_{L}$
$\displaystyle=\frac{1}{\ell!}\int_{-\infty}^{t-r}\mathrm{d}s\,\hat{\partial}_{L}\\!\left[\frac{1}{r}\int_{\alpha}^{\frac{t-r-s}{2}}\\!\\!\mathrm{d}\lambda\left(\frac{t-r-s}{2}\right)^{\ell}\\!\left(\frac{t-r-s}{2}-\lambda\right)^{\ell}\left(\frac{2}{\lambda}\right)^{\ell-1}\\!\\!S(\lambda,s)\right]\,,$
(31a) $\displaystyle\Psi^{[2]}_{L}$
$\displaystyle=-\frac{1}{\ell!}\int_{-\infty}^{t-r}\mathrm{d}s\,\hat{\partial}_{L}\\!\left[\frac{1}{r}\int_{\alpha}^{\frac{t+r-s}{2}}\\!\\!\mathrm{d}\lambda\left(\frac{t+r-s}{2}\right)^{\ell}\\!\left(\frac{t+r-s}{2}-\lambda\right)^{\ell}\left(\frac{2}{\lambda}\right)^{\ell-1}\\!\\!S(\lambda,s)\right]\,.$
(31b)

Next we remark that in both terms (31), one can commute the partial
differential operator $\hat{\partial}_{L}$ with the integral over $\lambda$,
since the terms coming from the differentiation of the bound of the integral,
_i.e._ $(t\pm r-s)/2$, have to be evaluated at $s=t\pm r$ and will clearly
vanish. Hence we can rewrite

$\displaystyle\Psi^{[1]}_{L}$
$\displaystyle=\frac{2^{\ell-1}}{\ell!}\int_{-\infty}^{t-r}\mathrm{d}s\\!\int_{\alpha}^{\frac{t-r-s}{2}}\\!\\!\mathrm{d}\lambda\,\lambda^{-\ell+1}S(\lambda,s)\,\hat{\partial}_{L}\\!\left[\frac{1}{r}\left(\frac{t-r-s}{2}\right)^{\ell}\,\left(\frac{t-r-s}{2}-\lambda\right)^{\ell}\right]\,,$
(32a) $\displaystyle\Psi^{[2]}_{L}$
$\displaystyle=-\frac{2^{\ell-1}}{\ell!}\int_{-\infty}^{t-r}\mathrm{d}s\\!\int_{\alpha}^{\frac{t+r-s}{2}}\\!\\!\mathrm{d}\lambda\,\lambda^{-\ell+1}S(\lambda,s)\,\hat{\partial}_{L}\\!\left[\frac{1}{r}\left(\frac{t+r-s}{2}\right)^{\ell}\\!\left(\frac{t+r-s}{2}-\lambda\right)^{\ell}\right]\,.$
(32b)

Then, using Eq. (A36) of [56] we know that

$\hat{\partial}_{L}\\!\left[\frac{1}{r}\left(\frac{t-r-s}{2}\right)^{\ell}\left(\frac{t-r-s}{2}-\lambda\right)^{\ell}\right]=\hat{\partial}_{L}\\!\left[\frac{1}{r}\left(\frac{t+r-s}{2}\right)^{\ell}\left(\frac{t+r-s}{2}-\lambda\right)^{\ell}\right]\,,$
(33)

so that actually the two contributions in (32) can be merged together in the
full solution, and we arrive at the elegant alternative form

$\Psi_{L}=-\frac{2^{\ell-1}}{\ell!}\int_{-\infty}^{t-r}\mathrm{d}s\\!\int_{\frac{t-r-s}{2}}^{\frac{t+r-s}{2}}\\!\\!\mathrm{d}\lambda\,\lambda^{-\ell+1}S(\lambda,s)\,\hat{\partial}_{L}\\!\left[\frac{1}{r}\left(\frac{t-r-s}{2}\right)^{\ell}\,\left(\frac{t-r-s}{2}-\lambda\right)^{\ell}\right]\,.$
(34)

This form is clearly independent of any choice for the constant $\alpha$. It
was also derived in a different way, using a direct multipolar expansion of
the Green function of the d’Alembertian operator, in Appendix D of [56].

By using Eq. (A35a) of [56] we can expand explicitly the last factor
containing the multipolar derivative operator. Introducing the following
coefficients, which will become ubiquitous in our practical computations,

$C_{ij}^{\ell}\equiv\frac{(\ell+i)!(\ell+j)!}{2^{i+j}i!j!(\ell-i)!(\ell-j)!(i+j)!}\\\
\qquad\text{(for {\hbox{0\leqslant i,j\leqslant\ell}})}\,,$ (35)

we conveniently write the general solution as

$\Psi_{L}=-\frac{\hat{n}_{L}}{2}\sum_{i=0}^{\ell}\sum_{j=0}^{\ell}\,\frac{(-)^{j}C_{ij}^{\ell}}{r^{i+1}}\int_{0}^{+\infty}\mathrm{d}\rho\,\rho^{i+j}\int_{\frac{\rho}{2}}^{\frac{\rho}{2}+r}\mathrm{d}\lambda\,\lambda^{-j+1}\,S(\lambda,u-\rho)\,.$
(36)

Note that although the coefficient $C_{ij}^{\ell}$ has been defined to be
symmetric in $i$ and $j$, these two indices play a different role. In
particular the index $i$ rules the behavior of the solution at infinity, when
$r\to+\infty$ with $t-r=$ const, depending on the fall-off properties at
infinity of the integrals in (36). A useful fact is that for $i=0$, the
coefficient is straightforwardly linked to the $j$-th derivative of the usual
Legendre polynomial evaluated at 1, namely

$C_{0j}^{\ell}=\frac{1}{j!}P_{\ell}^{(j)}(1)\,.$ (37)

In addition to the general formula (36), we also have an independant result
which directly provides the leading behavior of the solution at infinity, when
$r\to+\infty$ with $u=$ const, depending on the fall-off properties of the
source $S(r,u)$. Namely, suppose that $S(r,u)$ has the following asymptotic
behavior at infinity:

$S(r,u)=\mathcal{O}\left(\frac{1}{r^{3}}\right)\quad\text{(when
{\hbox{r\to+\infty}} with {\hbox{u}} fixed)}\,.$ (38)

More precisely the rigourous fall-off conditions of the source term are given
in Lemma 7.2 of [56]. Then the corresponding solution will behave dominantly
like $1/r$ at infinity, with leading term explicitly given by

$\Psi_{L}=\frac{(-)^{\ell}}{2^{\ell}}\frac{\hat{n}_{L}}{r}\int_{-\infty}^{t-r}\mathrm{d}s\,R_{\infty}^{(\ell)}\left(\frac{t-r-s}{2},s\right)+\mathcal{O}\left(\frac{1}{r^{2}}\right)\,.$
(39)

Here the function $R_{\infty}(\rho,s)$ is defined by (28) with the explicit
choice $\alpha=+\infty$, and the superscript $(\ell)$ means $\ell$ partial
derivatives with respect to $\rho$. We shall use the latter result as a
consistency check of the derivation of the leading behavior of the solution at
infinity.

###  IV.2 Application to tails-of-memory

We apply the general formalism in the previous section to the cubic iteration
and specifically to the tails-of-memory. In this case the main problem we face
is to find the solution of the wave equation (26) when the source term takes
the form (23). For ease of notation we denote the two quadrupole moments by
arbitrary time-varying functions $F(u)$ and $G(u)$, and consider the source
term (multiplying the STF harmonics $\hat{n}_{L}$):

$\underset{k,m}{\ \
S^{B}}(r,t-r)\equiv\left(\frac{r}{r_{0}}\right)^{B}r^{-k}\,G(t-r)\int_{1}^{+\infty}\mathrm{d}x\,Q_{m}(x)F(t-rx)\,,$
(40)

where $k\geqslant 1$ and $m\geqslant 0$ are integers, and the Legendre
function of the second kind $Q_{m}(x)$ is given by (16), as was obtained from
$\overline{Q}_{m}(x)$ using the relation (15). The source (40) is more complex
than the one required for the computation of tails and tails-of-tails. In the
latter case, the function $G(u)$ is actually constant (_i.e._ the ADM mass
$\mathrm{M}$), and we shall recover and extend the results found in [20, 22].
In the case where $G(u)$ and $F(u)$ are not constant the results we shall
derive are new.

In (40) we have multiplied the source by the regularization factor
$(r/r_{0})^{B}$, where $B\in\mathbb{C}$ and $r_{0}$ is an arbitrary constant,
see Eq. (8). Very importantly, this permits to ensure after analytic
continuation in $B\in\mathbb{C}$ that the fall-off condition (27) is always
satisfied. Applying the finite part, the solution we are looking for is

$\underset{k,m}{\ \Psi_{L}}\equiv\underset{B=0}{\mathrm{FP}}\underset{k,m}{\
\Psi^{B}_{L}}\,,\qquad\text{where}\qquad\underset{k,m}{\
\Psi^{B}_{L}}\equiv\Box_{\text{ret}}^{-1}\\!\underset{k,m}{\ \
S^{B}_{L}}=\Box_{\text{ret}}^{-1}\bigl{[}\hat{n}_{L}\\!\\!\underset{k,m}{\ \
S^{B}}\bigr{]}\,,$ (41)

and will be obtained by applying the formulae of the previous section. During
the computation, we shall also encounter explicit logarithms and poles
proportional to $1/B$. This will all boil down to computing

$\underset{k,m}{\
\chi_{L}}\equiv\underset{B=0}{\mathrm{FP}}\biggl{[}\frac{1}{B}\underset{k,m}{\
\Psi^{B}_{L}}\biggr{]}=\underset{B=0}{\mathrm{FP}}\left[\frac{\mathrm{d}}{\mathrm{d}B}\\!\\!\underset{k,m}{\
\
\Psi^{B}_{L}}\right]=\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\left[\ln\left(\frac{r}{r_{0}}\right)\\!\\!\underset{k,m}{\
\ S^{B}_{L}}\right]\,,$ (42)

where we have used the useful property of the regularization factor
$(r/r_{0})^{B}$, that an extra logarithm is generated by differentiating with
respect to $B$:

$\frac{\mathrm{d}}{\mathrm{d}B}\\!\underset{k,m}{\ \
S^{B}_{L}}=\ln\left(\frac{r}{r_{0}}\right)\\!\\!\underset{k,m}{\ \
S^{B}_{L}}\,.$ (43)

Applying the formulae of the previous section the solution for any
$B\in\mathbb{C}$ reads

$\displaystyle\underset{k,m}{\ \Psi^{B}_{L}}$
$\displaystyle=-\frac{\hat{n}_{L}}{2}\sum_{i=0}^{\ell}\sum_{j=0}^{\ell}\frac{(-)^{j}C_{ij}^{\ell}}{r^{i+1}}\underset{k,m}{\phi}\\!_{ij}^{B}\,,$
(44a) where $\displaystyle\underset{k,m}{\phi}\\!_{ij}^{B}$
$\displaystyle\equiv\int_{0}^{+\infty}\\!\\!\mathrm{d}\rho\,\rho^{i+j}G(u-\rho)\int_{\frac{\rho}{2}}^{\frac{\rho}{2}+r}\\!\\!\mathrm{d}\lambda\,\left(\frac{\lambda}{r_{0}}\right)^{B}\lambda^{-k-j+1}\int_{1}^{+\infty}\\!\\!\mathrm{d}x\,Q_{m}(x)F\bigl{[}u-\rho-\lambda(x-1)\bigr{]}\,.$
(44b)

Note that this solution is “exact”, valid at any radial distance $r$ except
$r=0$. But in the following, we shall mostly be interested in the asymptotic
limit when $r\to+\infty$ with $u=$ const. Permuting the integrals, introducing
the change of variable $\lambda\longrightarrow\tau=\lambda(x-1)$, and
factorizing out the expected leading behavior $1/r$ of the solution (in
anticipation of the limit $r\to+\infty$), we find that

$\displaystyle\underset{k,m}{\ \Psi^{B}_{L}}$
$\displaystyle=-\frac{\hat{n}_{L}}{2r}\int_{0}^{+\infty}\mathrm{d}\rho\,G(u-\rho)\int_{0}^{+\infty}\mathrm{d}\tau\,F(u-\rho-\tau)\underset{k,m}{\
K^{B}_{\ell}}(\rho,\tau,r)\,,$ (45a) where we have introduced for convenience
the kernel function (with $x=y+1$) $\displaystyle\underset{k,m}{\
K^{B}_{\ell}}(\rho,\tau,r)$
$\displaystyle=\left(\frac{\tau}{r_{0}}\right)^{B}\sum_{i=0}^{\ell}\frac{1}{r^{i}}\sum_{j=0}^{\ell}(-)^{j}C_{ij}^{\ell}\,\rho^{i+j}\,\tau^{-k-j+1}\int_{\frac{2\tau}{\rho+2r}}^{\frac{2\tau}{\rho}}\mathrm{d}y\,y^{k+j-2-B}Q_{m}(y+1)\,,$
$\displaystyle=\tau^{1-k}\int_{\frac{2\tau}{\rho+2r}}^{\frac{2\tau}{\rho}}\mathrm{d}y\,\left(\frac{\tau}{yr_{0}}\right)^{B}y^{k-2}Q_{m}(y+1)\,\Pi_{\ell}\left(1-\frac{\rho
y}{\tau},1+\frac{\rho}{r}\right)\,.$ (45b)

In the second line we have introduced the following symmetric bivariate
polynomial

$\Pi_{\ell}(x,y)=\sum_{i=0}^{\ell}\sum_{j=0}^{\ell}C^{\ell}_{ij}(x-1)^{i}(y-1)^{j}\,,$
(46)

which is related to the Legendre polynomial, recalling (37), by

$P_{\ell}(x)=\sum_{j=0}^{\ell}C^{\ell}_{0j}\left(x-1\right)^{j}=\Pi_{\ell}(x,1)=\Pi_{\ell}(1,x)\,.$
(47)

Similarly, the quantity ${}_{k,m}\chi_{L}^{B}$ defined in (42) reads

$\displaystyle\underset{k,m}{\ \chi^{B}_{L}}$
$\displaystyle=-\frac{\hat{n}_{L}}{2r}\int_{0}^{+\infty}\mathrm{d}\rho\,G(u-\rho)\int_{0}^{+\infty}\mathrm{d}\tau\,F(u-\rho-\tau)\underset{k,m}{\
L^{B}_{\ell}}(\rho,\tau,r)\,,$ (48a) where the kernel is defined as
$\displaystyle\underset{k,m}{\ L^{B}_{\ell}}$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}B}\biggl{[}\underset{k,m}{\
K^{B}_{\ell}}(\rho,\tau,r)\biggr{]}$
$\displaystyle=\sum_{i=0}^{\ell}\frac{1}{r^{i}}\sum_{j=0}^{\ell}(-)^{j}C_{i,j}^{\ell}\,\rho^{i+j}\,\tau^{-k-j+1}\int_{\frac{2\tau}{\rho+2r}}^{\frac{2\tau}{\rho}}\mathrm{d}y\left(\frac{\tau}{yr_{0}}\right)^{B}\,\ln\left(\frac{\tau}{r_{0}y}\right)y^{k+j-2}Q_{m}(x)\,,$
$\displaystyle=\tau^{1-k}\int_{\frac{2\tau}{\rho+2r}}^{\frac{2\tau}{\rho}}\mathrm{d}y\,\left(\frac{\tau}{yr_{0}}\right)^{B}\ln\left(\frac{\tau}{yr_{0}}\right)y^{k-2}Q_{m}(y+1)\,\Pi_{\ell}\left(1-\frac{\rho
y}{\tau},1+\frac{\rho}{r}\right)\,.$ (48b)

When $r\to+\infty$ the above kernels are dominated by the contribution $i=0$,
and the property (47) allows us to relate the asymptotic limit to the Legendre
polynomial as

$\displaystyle\underset{k,m}{\ K^{B}_{\ell}}$
$\displaystyle=\tau^{1-k}\int_{\frac{2\tau}{\rho+2r}}^{\frac{2\tau}{\rho}}\mathrm{d}y\,\left(\frac{\tau}{yr_{0}}\right)^{B}y^{k-2}Q_{m}(y+1)\,P_{\ell}\left(1-\frac{\rho
y}{\tau}\right)+o\left(1\right)\,,$ (49a) $\displaystyle\underset{k,m}{\
L^{B}_{\ell}}$
$\displaystyle=\tau^{1-k}\int_{\frac{2\tau}{\rho+2r}}^{\frac{2\tau}{\rho}}\mathrm{d}y\,\left(\frac{\tau}{yr_{0}}\right)^{B}\ln\left(\frac{\tau}{yr_{0}}\right)y^{k-2}Q_{m}(y+1)\,P_{\ell}\left(1-\frac{\rho
y}{\tau}\right)+o\left(1\right)\,.$ (49b)

We employ the Landau symbol $o$ for remainders, hence $o(1)$ means terms that
behave as $\sim\ln^{p}r/r$ with uncontrolled powers of $\ln r$ as
$r\rightarrow+\infty$.

####  IV.2.1 Cases where $k=1$ and $k=2$

In this subsection, we assume that $k\in\\{1,2\\}$. As proven in Appendix B,
${}_{k,m}\,\Psi^{B}_{L}$ and ${}_{k,m}\,\chi^{B}_{L}$ have a well-defined
limit when $B\rightarrow 0$, so we can drop the finite part prescription and
simply set $B=0$ in (45) and (48). Therefore we obtain the expression for the
kernels, only valid for $k\in\\{1,2\\}$,

$\displaystyle\underset{k,m}{\ K_{\ell}}(\rho,\tau,r)$
$\displaystyle\equiv\tau^{1-k}\int_{\frac{2\tau}{\rho+2r}}^{\frac{2\tau}{\rho}}\mathrm{d}y\,y^{k-2}Q_{m}(y+1)\,\Pi_{\ell}\left(1-\frac{\rho
y}{\tau},1+\frac{\rho}{r}\right)\,,$ (50a) $\displaystyle\underset{k,m}{\
L_{\ell}}(\rho,\tau,r)$
$\displaystyle\equiv\tau^{1-k}\int_{\frac{2\tau}{\rho+2r}}^{\frac{2\tau}{\rho}}\mathrm{d}y\,\ln\left(\frac{\tau}{yr_{0}}\right)y^{k-2}Q_{m}(y+1)\,\Pi_{\ell}\left(1-\frac{\rho
y}{\tau},1+\frac{\rho}{r}\right)\,.$ (50b)

When restricting our interest to the asympotic limit $r\rightarrow+\infty$, we
must treat separately the cases $k=1$ and $k=2$.

For $k=2$, the integrands in the kernels (50) are clearly integrable at the
$y\rightarrow 0$ bound, therefore the kernels converge in the
$r\rightarrow+\infty$ limit, such that we are allowed to define the
“asymptotic” kernels:

$\displaystyle\underset{2,m}{\ \overline{K}_{\ell}}(\rho,\tau)$
$\displaystyle\equiv\lim_{r\to+\infty}\underset{2,m}{\
K_{\ell}}(\rho,\tau,r)=\tau^{-1}\int_{0}^{\frac{2\tau}{\rho}}\mathrm{d}y\,Q_{m}(y+1)P_{\ell}\Bigl{(}1-\frac{\rho
y}{\tau}\Bigr{)}\,,$ (51a) $\displaystyle\underset{2,m}{\
\overline{L}_{\ell}}(\rho,\tau)$
$\displaystyle\equiv\lim_{r\to+\infty}\underset{2,m}{\
L_{\ell}}(\rho,\tau,r)=\tau^{-1}\int_{0}^{\frac{2\tau}{\rho}}\mathrm{d}y\,\ln\left(\frac{\tau}{yr_{0}}\right)Q_{m}(y+1)P_{\ell}\Bigl{(}1-\frac{\rho
y}{\tau}\Bigr{)}\,.$ (51b)

We have verified the above asymptotic limit in the case where $k=2$ using the
general statement in Eq. (39), where the relevant function is defined by Eq.
(28) with $\alpha=+\infty$.

In the case $k=1$, the situation is more complicated since the integrands of
the kernels are no longer integrable when $y\to 0$, so it is not possible to
simply take the limit $r\to+\infty$ like in (51). Instead the kernels exhibit
a logarithmic behavior as $r\to+\infty$. Since the radiative construction will
not exhibit any logarithmic behavior, the logarithms should cancel out in the
final metric, and it is crucial to verify this by controlling the logarithmic
limit of the kernel functions. We detail the case of the “$K$”-kernel

$\underset{1,m}{\
K_{\ell}}(\rho,\tau,r)=\int_{\frac{2\tau}{\rho+2r}}^{\frac{2\tau}{\rho}}\frac{\mathrm{d}y}{y}\,Q_{m}(y+1)\,\Pi_{\ell}\left(1-\frac{\rho
y}{\tau},1+\frac{\rho}{r}\right)\,.$ (52)

To extract the logarithmic behavior, we integrate by parts. First, we
introduce the regular part of the Legendre function of second type when $x\to
1^{+}$ defined as

$R_{m}(x)\equiv Q_{m}(x)+\frac{1}{2}P_{m}(x)\ln\left(\frac{x-1}{2}\right)\,.$
(53)

A useful observation is that $R_{m}(1)=-H_{m}$, where
$H_{m}=\sum_{k=1}^{m}k^{-1}$ is the harmonic number. Substituting $Q_{m}$ by
its expression in terms of $R_{m}$ and integrating by parts, we find

$\displaystyle\underset{1,m}{\ K_{\ell}}(\rho,\tau,r)$
$\displaystyle=\Bigg{[}\ln\left(\frac{y}{2}\right)\biggl{(}R_{m}(y+1)-\frac{1}{4}\ln\left(\frac{y}{2}\right)P_{m}(y+1)\biggr{)}\Pi_{\ell}\left(1-\frac{\rho
y}{\tau},1+\frac{\rho}{r}\right)\Bigg{]}_{y=\frac{2\tau}{\rho+2r}}^{y=\frac{2\tau}{\rho}}$
$\displaystyle-\int_{\frac{2\tau}{\rho+2r}}^{\frac{2\tau}{\rho}}\mathrm{d}y\,\ln\left(\frac{y}{2}\right)\frac{\mathrm{d}}{\mathrm{d}y}\bigg{[}R_{m}(y+1)\,\Pi_{\ell}\left(1-\frac{\rho
y}{\tau},1+\frac{\rho}{r}\right)\bigg{]}$
$\displaystyle+\frac{1}{4}\int_{\frac{2\tau}{\rho+2r}}^{\frac{2\tau}{\rho}}\mathrm{d}y\,\ln^{2}\left(\frac{y}{2}\right)\frac{\mathrm{d}}{\mathrm{d}y}\bigg{[}P_{m}(y+1)\,\Pi_{\ell}\left(1-\frac{\rho
y}{\tau},1+\frac{\rho}{r}\right)\bigg{]}\,,$ (54)

where the all-integrated terms are shown in the first line. Since $P_{m}(x)$,
$R_{m}(x)$ and $\Pi_{\ell}(x,y)$ as well as all their derivatives are
perfectly integrable as $x\rightarrow 1^{+}$, we can safely take the
$r\to+\infty$ expansion, and we obtain the explicit “polylogarithmic”
structure

$\displaystyle\underset{1,m}{\ K_{\ell}}(\rho,\tau,r)$
$\displaystyle=\frac{1}{4}\ln^{2}\left(\frac{r}{r_{0}}\right)-\frac{1}{2}\ln\left(\frac{r}{r_{0}}\right)\bigg{[}\ln\left(\frac{\tau}{2r_{0}}\right)+2H_{m}\bigg{]}+\underset{1,m}{\
\overline{K}_{\ell}}(\rho,\tau)+o\left(1\right)\,,$ (55a) where
${}_{1,m}\overline{K}_{\ell}(\rho,\tau)$ does not exhibit any $r$-dependence.
We obtain its explicit expression as $\displaystyle\underset{1,m}{\
\overline{K}_{\ell}}(\rho,\tau)$
$\displaystyle=\frac{1}{4}\ln^{2}\left(\frac{\tau}{2r_{0}}\right)+H_{m}\ln\left(\frac{\tau}{2r_{0}}\right)$
$\displaystyle-\frac{(-)^{\ell}}{4}\left[\ln^{2}\left(\frac{\tau}{2r_{0}}\right)-2\ln\left(\frac{\tau}{2r_{0}}\right)\ln\left(\frac{\rho}{2r_{0}}\right)+\ln^{2}\left(\frac{\rho}{2r_{0}}\right)\right]P_{m}\left(1+\frac{2\tau}{\rho}\right)$
$\displaystyle-(-)^{\ell}\left[\ln\left(\frac{\rho}{2r_{0}}\right)-\ln\left(\frac{\tau}{2r_{0}}\right)\right]R_{m}\left(1+\frac{2\tau}{\rho}\right)$
$\displaystyle+\frac{1}{4}\int_{0}^{\frac{2\tau}{\rho}}\mathrm{d}y\,\ln^{2}\left(\frac{y}{2}\right)\frac{\mathrm{d}}{\mathrm{d}y}\bigg{[}P_{m}(y+1)\,P_{\ell}\left(1-\frac{\rho
y}{\tau}\right)\bigg{]}$
$\displaystyle-\int_{0}^{\frac{2\tau}{\rho}}\mathrm{d}y\,\ln\left(\frac{y}{2}\right)\frac{\mathrm{d}}{\mathrm{d}y}\bigg{[}R_{m}(y+1)\,P_{\ell}\left(1-\frac{\rho
y}{\tau}\right)\bigg{]}\,.$ (55b)

The same reasonning applied to the “$L$”-kernel also gives a polylogarithmic
structure but in this case with powers of the logarithm up to three:

$\displaystyle\underset{1,m}{\ L_{\ell}}(\rho,\tau,r)$
$\displaystyle=\frac{1}{6}\ln^{3}\left(\frac{r}{r_{0}}\right)-\frac{1}{4}\ln^{2}\left(\frac{r}{r_{0}}\right)\left[\ln\left(\frac{\tau}{2r_{0}}\right)+2H_{m}\right]+\underset{1,m}{\
\overline{L}_{\ell}}(\rho,\tau)+o\left(1\right)\,,$ (56a) and where
${}_{1,m}\overline{L}_{\ell}(\rho,\tau)$ is explicitly given by
$\displaystyle\underset{1,m}{\ \overline{L}_{\ell}}(\rho,\tau)$
$\displaystyle=\frac{1}{12}\ln^{3}\left(\frac{\tau}{2r_{0}}\right)+\frac{1}{2}H_{m}\ln^{2}\left(\frac{\tau}{2r_{0}}\right)$
(56b)
$\displaystyle-\frac{(-)^{\ell}}{12}\left[\ln^{3}\left(\frac{\tau}{2r_{0}}\right)-3\ln\left(\frac{\tau}{2r_{0}}\right)\ln^{2}\left(\frac{\rho}{2r_{0}}\right)+2\ln^{3}\left(\frac{\rho}{2r_{0}}\right)\right]P_{m}\left(1+\frac{2\tau}{\rho}\right)$
$\displaystyle-\frac{(-)^{\ell}}{2}\left[\ln^{2}\left(\frac{\rho}{2r_{0}}\right)-\ln^{2}\left(\frac{\tau}{2r_{0}}\right)\right]R_{m}\left(1+\frac{2\tau}{\rho}\right)$
$\displaystyle+\int_{0}^{\frac{2\tau}{\rho}}\mathrm{d}y\,\ln^{2}\left(\frac{y}{2}\right)\left[\frac{1}{4}\ln\left(\frac{\tau}{2r_{0}}\right)-\frac{1}{6}\ln\left(\frac{y}{2}\right)\right]\frac{\mathrm{d}}{\mathrm{d}y}\bigg{[}P_{m}(y+1)\,P_{\ell}\left(1-\frac{\rho
y}{\tau}\right)\bigg{]}$
$\displaystyle+\int_{0}^{\frac{2\tau}{\rho}}\mathrm{d}y\,\ln\left(\frac{y}{2}\right)\left[-\ln\left(\frac{\tau}{2r_{0}}\right)+\frac{1}{2}\ln\left(\frac{y}{2}\right)\right]\frac{\mathrm{d}}{\mathrm{d}y}\bigg{[}R_{m}(y+1)\,P_{\ell}\left(1-\frac{\rho
y}{\tau}\right)\bigg{]}\,.$

Again, a beautiful check of the radiative construction of the metric in Sec.
II, is that all these explicitely determined far zone logarithms will be
compensated by those induced by the applied gauge transformations, notably the
one described in Eqs. (19).

In practical computations, it is easier to compute the full kernels
${}_{1,m}K_{\ell}(\rho,\tau,r)$ and ${}_{1,m}L_{\ell}(\rho,\tau,r)$, and then
to remove the logarithmic dependencies to obtain the reduced asymptotic
kernels ${}_{1,m}\overline{K}_{\ell}(\rho,\tau)$ and
${}_{1,m}\overline{L}_{\ell}(\rho,\tau)$. However, computing the two later
quantities directly from their explicit expressions (55) and (56b) yields, of
course, the same result.

####  IV.2.2 Cases where $k\geqslant 3$

In the previous subsection we could integrate the source term (40) when
$k\in\\{1,2\\}$. _A priori_ , the cases $k\geqslant 3$ are more difficult and
(when $G$ is not constant) analytic closed-form expressions for the retarded
integral seem non trivial. Hence we proceed differently and prove that we can
retrieve the cases $k\geqslant 3$ from the known cases $k\in\\{1,2\\}$. This
will show that analytic closed-form expressions also exist in the cases
$k\geqslant 3$.

Given a source term with $k\geqslant 3$, endowed with the associated
regularization factor $r^{B}$ (posing here $r_{0}=1$), we reduce it to a fully
integrated part and new source terms with decreased values of $k$ by two steps
at most: $k-1$ and $k-2$, by means of the identity:

$\displaystyle\hat{n}_{L}r^{B-k}G(t-r)\int_{1}^{+\infty}\mathrm{d}xQ_{m}(x)F(t-rx)$
$\displaystyle\qquad=\Box\left[\frac{\hat{n}_{L}r^{B-k+2}}{(k+\ell-2-B)(k-\ell-3-B)}G(t-r)\int_{1}^{+\infty}\mathrm{d}x\,Q_{m}(x)F(t-rx)\right]$
$\displaystyle\qquad\quad-\frac{2(k-3-B)\hat{n}_{L}r^{B-k+1}}{(k+\ell-2-B)(k-\ell-3-B)}\bigg{(}\overset{(1)}{G}(t-r)\int_{1}^{+\infty}\mathrm{d}x\,Q_{m}(x)F(t-rx)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad+\
G(t-r)\int_{1}^{+\infty}\mathrm{d}x\,x\,Q_{m}(x)\overset{(1)}{F}(t-rx)\bigg{)}$
$\displaystyle\qquad\quad-\frac{\hat{n}_{L}r^{B-k+2}}{(k+\ell-2-B)(k-\ell-3-B)}\bigg{(}2\,\overset{(1)}{G}(t-r)\int_{1}^{+\infty}\mathrm{d}x\,(x-1)Q_{m}(x)\overset{(1)}{F}(t-rx)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad+G(t-r)\int_{1}^{+\infty}\mathrm{d}x\,(x^{2}-1)\,Q_{m}(x)\overset{(2)}{F}(t-rx)\bigg{)}\,.$
(57)

This formula, when applied iteratively, allows us to reduce any case
$k\geqslant 3$ in terms of the cases when $k=1$ and $2$, modulo a series of
all-integrated terms (_i.e._ , inside the d’Alembertian operator). Note that
we introduced a new type of integral,
$\int_{1}^{+\infty}\mathrm{d}x\,x^{n}\,Q_{m}(x)H(t-rx)$ where
$n\in\mathbb{N}$, but which can easily be recast in the previous form (_i.e._
$n=0$) by recursively applying Bonnet’s recursion formula for the Legendre
function, namely

$\displaystyle x\,Q_{m}(x)$
$\displaystyle=\frac{m+1}{2m+1}Q_{m+1}(x)+\frac{m}{2m+1}Q_{m-1}(x)$
$\displaystyle\text{if $m\geqslant 1$}\,,$ (58a) $\displaystyle x\,Q_{0}(x)$
$\displaystyle=Q_{1}(x)+1$ $\displaystyle\text{if $m=0$}\,.$ (58b)

An important point is that we keep the Hadamard regulator $r^{B}$ “alive” in
(IV.2.2). Indeed, by applying the inverse d’Alembertian operator on both sides
of (IV.2.2), we see that the first term in the right side being the
d’Alembertian of a source term containing the regulator $r^{B}$, will directly
yield that source term proportional to $r^{B}$ without any additional
homogeneous solution (since no homogeneous solution can be proportional to
$r^{B}$). Hence the iteration with $r^{B}$ can be done blindly, ignoring
homogeneous solutions, and only at the end do we apply the finite part when
$B\to 0$.

However a price we have to pay is that the $B$-dependent coefficients in
(IV.2.2) can generate a simple pole when $B\to 0$ which will compete with the
higher contribution $\propto B$ in the retarded integral, and _vice-versa_ , a
term $\propto B$ in the coefficient will be compensated by a pole coming from
the retarded integral. A useful fact to remember in this respect is that the
integrals are convergent for $k\in\\{1,2\\}$ and hence no poles $\propto 1/B$
can be generated in these cases. Furthermore, the structure of the identity
(IV.2.2) can only generate simple poles $\propto 1/B$, and no double poles
$\propto 1/B^{2}$ or any poles of higher order. The presence of these simple
poles is the main reason why we introduced the ${}_{k,m}\chi_{L}$ integral, as
is clear by its first definition in (42).

Finally, although we now have all the needed formulae to integrate the tails-
of-memory, our recursion formula (IV.2.2) coupled to the $m=0$ case of the
Bonnet formula (58b) can generate instantaneous terms in the source, which
must also be integrated. Although the integration of such terms is well known
using standard integration techniques [14, 28], its generalization to the case
where simple poles $\propto 1/B$ can appear in the source was unknown.

###  IV.3 Application to tails-of-tails

The previous observation motivates us to extend our formalism developed for
complicated tails-of-memory to the easier cases of tails-of-tails, tails and
even instantaneous terms, such that we can treat the entire problem
consistently, using one single formalism. When comparing with previous works
[20, 22] this will provide important tests of the results of Section IV.2.
Thus, we specialize the formulae of the previous section to the case where
$G\equiv 1$, since in the case of the tails-of-tails two of the moments are
just the mass $\mathrm{M}$.

When $k\geqslant 3$, we can again reduce the source to the cases $k=1$ and
$k=2$ using (IV.2.2), with now $G^{(1)}(u)=0$. Then, after lengthy
computations where we perform two of the three integrations of the tails-of-
memory (one on the $y$ variable in the kernel and one on the time variable
$\tau$), we find drastically simpler expression when $k\in\\{1,2\\}$, which
read

$\displaystyle\underset{2,m}{\ \overline{\Psi}_{\ell}}{\bigg{|}}_{G=1}$
$\displaystyle=2\,\underset{{2,m}}{\alpha_{\ell}}\overset{(-1)}{F}(u)\,,$
(59a) $\displaystyle\underset{2,m}{\ \overline{\chi}_{\ell}}{\bigg{|}}_{G=1}$
$\displaystyle=2\int_{0}^{+\infty}\mathrm{d}\tau\,\left[\underset{{2,m}}{\beta_{\ell}}\ln\left(\frac{\tau}{2r_{0}}\right)+\underset{{2,m}}{\gamma_{\ell}}\right]F(u-\tau)\,,$
(59b) $\displaystyle\underset{1,m}{\ \overline{\Psi}_{\ell}}{\bigg{|}}_{G=1}$
$\displaystyle=\frac{1}{4}\int_{0}^{+\infty}\mathrm{d}\tau\,\overset{(-1)}{F}(u-\tau)\left[\ln^{2}\left(\frac{\tau}{2r_{0}}\right)+4H_{m}\ln\left(\frac{\tau}{2r_{0}}\right)-8\underset{{1,m}}{\delta_{\ell}}\right]\,,$
(59c) $\displaystyle\underset{1,m}{\ \overline{\chi}_{\ell}}{\bigg{|}}_{G=1}$
$\displaystyle=\frac{1}{12}\int_{0}^{+\infty}\mathrm{d}\tau\,\overset{(-1)}{F}(u-\tau)\bigg{[}\ln^{3}\left(\frac{\tau}{2r_{0}}\right)+6H_{m}\ln^{2}\left(\frac{\tau}{2r_{0}}\right)-24\underset{{1,m}}{\delta_{\ell}}\ln\left(\frac{\tau}{2r_{0}}\right)+12\underset{{1,m}}{\varepsilon_{\ell}}\bigg{]}\,.$
(59d)

where ${}_{k,m}{\overline{\Psi}_{L}}\equiv-\frac{\hat{n}_{L}}{2r}\
{}_{k,m}{\overline{\Psi}_{\ell}}$ and
${}_{k,m}{\overline{\chi}_{L}}\equiv-\frac{\hat{n}_{L}}{2r}\
{}_{k,m}{\overline{\chi}_{\ell}}$. The numerical constants
${}_{2,m}\alpha_{\ell}$, ${}_{2,m}\beta_{\ell}$ and ${}_{2,m}\gamma_{\ell}$
are defined by the integrals

$\displaystyle\underset{2,m}{\ \alpha_{\ell}}$
$\displaystyle\equiv\underset{2,m}{\
\beta_{\ell}}\equiv\int_{1}^{+\infty}\mathrm{d}x\,Q_{m}(x)Q_{\ell}(x)\,,$
(60a) $\displaystyle\underset{2,m}{\ \gamma_{\ell}}$
$\displaystyle\equiv\frac{1}{2^{\ell+1}}\int_{1}^{+\infty}\mathrm{d}x\,Q_{m}(x)\int_{-1}^{1}\mathrm{d}z\,\frac{(1-z^{2})^{\ell}}{(x-y)^{\ell+1}}\left(-\ln\left(\frac{x-y}{2}\right)+H_{\ell}\right)\,.$
(60b)

Note that ${}_{2,m}\alpha_{\ell}$ is a special case of
$\,{}_{k,m}\alpha_{\ell}$ defined for $2\leqslant k\leqslant\ell+2$ in (A16)
of [20], whereas $\,{}_{2,m}\beta_{\ell}$ and $\,{}_{2,m}\gamma_{\ell}$ are an
extension of (A20) of [20], in which these constants were defined only for
$k\geqslant\ell+3$. In the case $k=2$, we find it natural to define
${}_{2,m}\alpha_{\ell}$ and ${}_{2,m}\beta_{\ell}$ to be equal.

Similarly one can define ${}_{1,m}\delta_{\ell}$ and
${}_{1,m}\varepsilon_{\ell}$ using formal integrals but at the price of
introducing a Hadamard “partie finie” (pf) to cure the bound of the integrals
at $x=1$:

$\displaystyle\underset{1,m}{\ \delta_{\ell}}$
$\displaystyle\equiv\mathrm{pf}\int_{1}^{+\infty}\mathrm{d}x\,Q_{m}(x)\frac{\mathrm{d}Q_{\ell}}{\mathrm{d}x}\,,$
(61a) $\displaystyle\underset{1,m}{\ \varepsilon_{\ell}}$
$\displaystyle\equiv\mathrm{pf}\int_{1}^{+\infty}\mathrm{d}x\,Q_{m}(x)\frac{\mathrm{d}^{2}S_{\ell}}{\mathrm{d}^{2}x}\,.$
(61b)

The Hadamard partie finie is defined in the usual way by removing the
divergent part of the integral (and using for our purpose here a Hadamard
scale equal to 1):

$\displaystyle\mathrm{pf}\int_{1}^{+\infty}\mathrm{d}x\,Q_{m}(x)\frac{\mathrm{d}Q_{\ell}}{\mathrm{d}x}\equiv$
$\displaystyle\lim_{\eta\to
0}\left\\{\int_{1+\eta}^{+\infty}\mathrm{d}x\,Q_{m}(x)\frac{\mathrm{d}Q_{\ell}}{\mathrm{d}x}+\frac{1}{8}\ln^{2}\left(\frac{\eta}{2}\right)+\frac{1}{2}H_{m}\ln\left(\frac{\eta}{2}\right)\right\\}\,,$
(62a)
$\displaystyle\mathrm{pf}\int_{1}^{+\infty}\mathrm{d}x\,Q_{m}(x)\frac{\mathrm{d}^{2}S_{\ell}}{\mathrm{d}^{2}x}\equiv$
$\displaystyle\lim_{\eta\to
0}\left\\{\int_{1+\eta}^{+\infty}\mathrm{d}x\,Q_{m}(x)\frac{\mathrm{d}^{2}S_{\ell}}{\mathrm{d}^{2}x}+\frac{1}{6}\ln^{3}\left(\frac{\eta}{2}\right)+\frac{1}{2}H_{m}\ln^{2}\left(\frac{\eta}{2}\right)\right\\}\,.$
(62b)

In the definition of ${}_{1,m}\varepsilon_{\ell}$ the function $S_{\ell}$ is
defined like for the Neumann formula for the Legendre function [see (16)],

$\displaystyle S_{\ell}(x)$
$\displaystyle\equiv\frac{1}{2}\int_{-1}^{1}\mathrm{d}y\,P_{\ell}(y)\ln^{2}\left(\frac{x-y}{2}\right)$
$\displaystyle=\sum_{j=0}^{\ell}\sum_{i=0}^{j}\frac{(-)^{i+1}P_{\ell}^{(j)}(1)}{2(i+1)!(j-i)!}\Bigg{\\{}(x-1)^{j+1}\left[\ln^{2}\left(\frac{x-1}{2}\right)-\frac{2}{i+1}\ln\left(\frac{x-1}{2}\right)+\frac{2}{(i+1)^{2}}\right]$
$\displaystyle\qquad+(-)^{j+\ell+1}(x+1)^{j+1}\left[\ln^{2}\left(\frac{x+1}{2}\right)-\frac{2}{i+1}\ln\left(\frac{x+1}{2}\right)+\frac{2}{(i+1)^{2}}\right]\Bigg{\\}}\,.$
(63)

The constant ${}_{1,m}\delta_{\ell}$ was used in (A7-A9) of [20] in the
special case where $m=\ell$, but was not given a name. The constants defined
in [20] were later shown in [22] to be sometimes ill-defined, _e.g._ for
$k=0$, $\ell=0$ and $m=0$, and this problem was circumvented on a case-by-case
basis. Since we restrict our attention to $k\in\\{1,2\\}$, we are assured that
our constants are always well defined.

All the coefficients we shall need can be computed analytically. The
${}_{2,m}\alpha_{\ell}$ constant (and, hence, $\,{}_{2,m}\beta_{\ell}$) has a
closed form expression given by (A18) of [20], which reads

$\underset{2,m}{\ \alpha_{\ell}}=\underset{2,m}{\
\beta_{\ell}}=\begin{cases}\displaystyle\frac{H_{m}-H_{\ell}}{(m-\ell)(m+\ell+1)}&\text{if
}m\neq\ell\,,\\\\[8.5359pt]
\displaystyle\frac{1}{2m+1}\left(\frac{\pi^{2}}{6}-\sum_{j=1}^{m}\frac{1}{j^{2}}\right)&\text{if
}m=\ell\,.\end{cases}$ (64)

The ${}_{1,m}\delta_{\ell}$ constant has a simple expression when $m=\ell$, as
pointed out in Eqs. (A7–A9) of [20], which reads
${}_{1,\ell}\delta_{\ell}=-\frac{1}{2}H_{\ell}^{2}$. When $m\neq\ell$, by
integrating by parts the integral expression of the ${}_{1,m}\delta_{\ell}$
constant, we find the simple property

$\underset{1,m}{\ \delta_{\ell}}+\underset{1,\ell}{\
\delta_{m}}=-H_{m}H_{\ell}\,.$ (65a) We can then obtain the full expression of
${}_{1,m}\delta_{\ell}$ by using recursively the differential equation
$Q^{\prime}_{n+1}-Q^{\prime}_{n-1}=(2n+1)Q_{n}$. In the case where $\ell$ and
$m$ have the same parity, _i.e._ $\ell-m\in 2\mathbb{Z}$, the differential
equation allows us to express the constant only in terms of the case
${}_{1,\ell}\delta_{\ell}$, along with many ${}_{2,m}\alpha_{\ell}$ terms,
which are easy to compute. Thanks to (65a), we can further restrict attention
to the case where $\ell\geqslant m$, and the relevant formula then reads
$\underset{1,m}{\
\delta_{\ell}}=\sum_{j=1}^{\frac{\ell-m}{2}}\bigl{(}2m+4j-1\bigr{)}\underset{2,m}{\
\alpha_{m+2j-1}}-\frac{1}{2}H_{m}^{2}\,,$ (65b) if $\ell-m\in 2\mathbb{N}$. In
the reverse case where $\ell$ and $m$ have opposite parity, _i.e._ $\ell-m\in
2\mathbb{Z}+1$, we can again use (65a) to restrict attention to the case where
$\ell$ is odd and $m$ is even. In this case, we can again use the differential
equation as well as integration by parts to express our constant solely in
terms of the explicitly known case, ${}_{1,0}\delta_{1}=\frac{\pi^{2}}{12}$,
modulo some ${}_{2,m}\alpha_{\ell}$ terms. We thus find the relevant formula,
$\displaystyle\underset{1,m}{\
\delta_{\ell}}=-H_{m}+\frac{\pi^{2}}{12}+\sum_{j=1}^{\frac{\ell-1}{2}}\bigl{(}2\ell-4j+3\bigr{)}\underset{2,m}{\
\alpha}{}_{\ell-2j+1}-\sum_{j=1}^{\frac{m}{2}}\bigl{(}2m-4j+3\bigr{)}$
$\displaystyle\underset{2,1}{\ \alpha}{}_{m-2j+1}\,,$ (65c)

if $\ell\in 2\mathbb{N}+1$ and $m\in 2\mathbb{N}$. Finally the
${}_{2,m}\gamma_{\ell}$ and ${}_{1,m}\varepsilon_{\ell}$ constants were only
needed for the following values of $m$ and $\ell$, which we calculated using
_Mathematica_ :

$\displaystyle\underset{2,0}{\ \gamma_{1}}=\frac{5}{4}\,,$
$\displaystyle\underset{2,2}{\ \gamma_{1}}=\frac{13}{32}\,,$
$\displaystyle\underset{2,4}{\ \gamma_{1}}=\frac{95}{432}\,,$
$\displaystyle\underset{1,0}{\ \varepsilon_{0}}=\frac{\zeta(3)}{2}\,,$
$\displaystyle\underset{1,2}{\
\varepsilon_{0}}=\frac{9}{2}+\frac{\zeta(3)}{2}\,,$
$\displaystyle\underset{1,4}{\
\varepsilon_{0}}=\frac{3995}{432}+\frac{\zeta(3)}{2}\,,$
$\displaystyle\underset{1,1}{\
\varepsilon_{1}}=\frac{1}{2}+\frac{\zeta(3)}{2}\,,$
$\displaystyle\underset{1,3}{\
\varepsilon_{1}}=\frac{299}{48}+\frac{\zeta(3)}{2}\,,$
$\displaystyle\underset{1,5}{\
\varepsilon_{1}}=\frac{887}{80}+\frac{\zeta(3)}{2}\,,$
$\displaystyle\underset{1,0}{\
\varepsilon_{2}}=-\frac{15}{2}+\frac{\zeta(3)}{2}\,,$
$\displaystyle\underset{1,2}{\
\varepsilon_{2}}=\frac{33}{16}+\frac{\zeta(3)}{2}\,,$
$\displaystyle\underset{1,4}{\
\varepsilon_{2}}=\frac{3425}{432}+\frac{\zeta(3)}{2}\,,$ (66)

where $\zeta$ is the Riemann function and $\zeta(3)$ is the Apéry constant. In
our computation, the contributions proportional to this constant actually
cancel out.

We end this section by a word concerning the integration of $B$-regularized
source terms that are actually instantaneous, _i.e._ just of the simple type
${}_{k}S^{B}(r,t-r)\equiv(\frac{r}{r_{0}})^{B}\,r^{-k}F(t-r)$ (multiplied by
the multipolarity factor $\hat{n}_{L}$). These source terms yield for instance
the tails at quadratic order. In line with our general formalism, if
$k\geqslant 3$, we bring ourselves to the case $k=2$ using the recursion
formula:

$\displaystyle\hat{n}_{L}r^{B-k}F(t-r)$
$\displaystyle=\Box\left[\frac{\hat{n}_{L}r^{B-k+2}}{(k+\ell-2-B)(k-\ell-3-B)}F(t-r)\right]$
$\displaystyle\quad\
-\frac{2(k-3-B)\hat{n}_{L}r^{B-k+1}}{(k+\ell-2-B)(k-\ell-3-B)}\overset{(1)}{F}(t-r)\,.$
(67)

The case $k=2$ is well known [14, 28, 20], and for the asymptotic limit reads

$\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}\left(\frac{r}{r_{0}}\right)^{B}\frac{\hat{n}_{L}}{r^{2}}F(t-r)\biggr{]}=\frac{\hat{n}_{L}}{2r}\int_{0}^{+\infty}\mathrm{d}\tau\,F(u-\tau)\Bigl{[}\ln\left(\frac{\tau}{2r}\right)+2H_{\ell}\Bigr{]}+o\left(\frac{1}{r}\right)\,,$
(68)

where we could have dropped the finite part prescription and set $B=0$,
because in this case the inverse d’Alembertian integral is convergent.
However, because of the potential appearance of single poles due to our
recursion formula (IV.3), we will also need the corresponding formula with an
extra factor $B^{-1}$ multiplying the source term. We give here the relevant
formula without proof:

$\displaystyle\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}\frac{1}{B}\left(\frac{r}{r_{0}}\right)^{B}\frac{\hat{n}_{L}}{r^{2}}F(t-r)\biggr{]}$
(69)
$\displaystyle\quad=\frac{\hat{n}_{L}}{4r}\int_{0}^{+\infty}\\!\mathrm{d}\tau\,F(u-\tau)\biggl{[}\ln^{2}\left(\frac{\tau}{2r_{0}}\right)+4H_{\ell}\ln\left(\frac{\tau}{2r_{0}}\right)+4H_{\ell}^{2}-\ln^{2}\left(\frac{r}{r_{0}}\right)\biggr{]}+o\left(\frac{1}{r}\right)\,.$

##  V Implementing the calculation of tails-of-memory

###  V.1 Explicit integration of the asymptotic kernels

Up to this point, the kernels in the asymptotic limit $r\to+\infty$ were
defined only in an integral form, given by Eqs. (51) in the relatively easy
case where $k=2$, and by the more complex forms (55) and (56b) when $k=1$.
These integrals are too complicated to seek a general explicit formula valid
for arbitrary $m$ and $\ell$. However, we can easily compute all these
integrals on a case-by-case basis. For this it suffices to insert into them
the explicit expressions of the Legendre polynomial $P_{m}(x)$ and the
Legendre function of the second kind $Q_{m}(x)$, which are simply polynomials
multiplied by some logarithms, see (16).

In this way we find that the general structure of the kernels
${}_{k,m}K_{\ell}(\rho,\tau,r)$ and ${}_{k,m}L_{\ell}(\rho,\tau,r)$ for $k=1$
and $2$, up to $o(1)$ precision, is of the type

$\left(\begin{array}[]{l}\displaystyle\underset{k,m}{\
K_{\ell}}(\rho,\tau,r)\\\\[19.91684pt] \displaystyle\underset{k,m}{\
L_{\ell}}(\rho,\tau,r)\end{array}\right)=\sum_{s,p,q}\left(\begin{array}[]{l}\displaystyle\underset{k,m}{\
\mathcal{X}}\\!\\!{}_{\,s,p,q}\\\\[19.91684pt] \displaystyle\underset{k,m}{\
\mathcal{Y}}\\!\\!{}_{\,s,p,q}\end{array}\right)\int_{\frac{2\tau}{\rho+2r}}^{\frac{2\tau}{\rho}}\mathrm{d}y\,y^{s}\ln^{p}\left(\frac{y}{2}\right)\ln^{q}\left(1+\frac{y}{2}\right)\,,$
(70)

where ${}_{k,m}\mathcal{X}_{s,p,q}$ and ${}_{k,m}\mathcal{Y}_{s,p,q}$ denote
some numerical coefficients, with $s$ an integer such that $s\geqslant k-2$,
and where we can restrict our attention only to the values $p=0,1,2,3$ and
$q=0,1$. When $s\geqslant 0$, the lower bound of the integral can be set to
$0$ in the $r\to+\infty$ limit, but not when $s=-1$ which can happen when
$k=1$. In the latter case, the integral will develop the logarithmic behavior
when $r\to+\infty$ which has already been obtained in (55a) and (56a), and we
now only consider the finite part.

When $q=0$, these integrals are of course well-known (see _e.g._ 2.721-722 in
[64]), so we now restrict our discussing to the case $q=1$. Let us first
examine the case $q=1$ and $s\geqslant 0$; we thus set the lower bound to
zero, in the limit $r\to+\infty$. We perform an integration by parts using
2.722 in [64],

$\displaystyle\int_{0}^{\frac{2\tau}{\rho}}$
$\displaystyle\mathrm{d}y\,y^{s}\ln^{p}\left(\frac{y}{2}\right)\ln\left(1+\frac{y}{2}\right)$
(71)
$\displaystyle=\sum_{i=0}^{p}\frac{(-)^{i}p!}{(p-i)!(s+1)^{i+1}}\Bigg{\\{}\left(\frac{2\tau}{\rho}\right)^{s+1}\\!\ln^{p-i}\left(\frac{\tau}{\rho}\right)\ln\left(1+\frac{\tau}{\rho}\right)-\int_{0}^{\frac{2\tau}{\rho}}\mathrm{d}y\,\frac{y^{s+1}}{y+2}\ln^{p-i}\left(\frac{y}{2}\right)\Bigg{\\}}\,.$

The rational fraction $y^{s+1}(y+2)^{-1}$ in the remaining integral can be
expanded as a polynomial, plus a function proportional to $(y+2)^{-1}$. From
this result we see that the only non trivial integrals left are of the type
$\int_{0}^{z}\mathrm{d}y\,(y+2)^{-1}\ln^{j}\left(y/2\right)$, and these are
simply related to the polylogarithm functions, for instance222The
polylogarithm, or Jonquière’s function, is defined for any $n\in\mathbb{N}$
and $z>0$ as
$\mathrm{Li}_{n}(z)\equiv\frac{(-)^{n-1}}{(n-2)!}\int_{0}^{z}\frac{\mathrm{d}s}{s}\,\ln^{n-2}\left(\frac{s}{z}\right)\ln(1-s)\,.$
In particular, the dilogarithm, or Spence’s function, is defined as
$\mathrm{Li}_{2}(z)=-\int_{0}^{z}\frac{\mathrm{d}s}{s}\,\ln(1-s)$. We also
find that $\mathrm{Li}_{1}(z)=-\ln(1-z)$.

$\displaystyle\int_{0}^{\frac{2\tau}{\rho}}\frac{\mathrm{d}y}{y+2}\ln^{2}\left(\frac{y}{2}\right)$
$\displaystyle=-2\mathrm{Li}_{3}\left(-\frac{\tau}{\rho}\right)+2\mathrm{Li}_{2}\left(-\frac{\tau}{\rho}\right)\ln\left(\frac{\tau}{\rho}\right)-\mathrm{Li}_{1}\left(-\frac{\tau}{\rho}\right)\ln^{2}\left(\frac{\tau}{\rho}\right)\,.$
(72)

Now, the only remaining cases are when $q=1$ and $s=-1$, but in these cases
too the integrals are related to polylogarithms, for instance

$\displaystyle\int_{0}^{\frac{2\tau}{\rho}}\frac{\mathrm{d}y}{y}\ln\left(\frac{y}{2}\right)\ln\left(1+\frac{y}{2}\right)=\mathrm{Li}_{3}\left(-\frac{\tau}{\rho}\right)-\,\mathrm{Li}_{2}\left(-\frac{\tau}{\rho}\right)\ln\left(\frac{\tau}{\rho}\right)\,.$
(73)

To summarize we have obtained the following general structure of the kernel
functions in the asymptotic limit $r\to+\infty$, discarding the leading
logarithmic behavior in the $k=1$ case — recall (55a) and (56a) —, as

$\left(\begin{array}[]{l}\displaystyle\underset{k,m}{\
\overline{K_{\ell}}}(\rho,\tau)\\\\[19.91684pt] \displaystyle\underset{k,m}{\
\overline{L_{\ell}}}(\rho,\tau)\end{array}\right)=\tau^{1-k}\sum_{\scriptsize{\begin{matrix}j\in\mathbb{Z}\\\
p,q,n\in\mathbb{N}\end{matrix}}}\left(\begin{array}[]{l}\displaystyle\underset{k,m}{\
\kappa}{}_{j,p,q,n}\\\\[19.91684pt] \displaystyle\underset{k,m}{\
\lambda}{}_{j,p,q,n}\end{array}\right)\left(\frac{\tau}{\rho}\right)^{j}\ln^{p}\left(\frac{\tau}{2r_{0}}\right)\ln^{q}\left(\frac{\rho}{2r_{0}}\right)\,\overline{\mathrm{Li}}_{n}\left(-\frac{\tau}{\rho}\right)\,,$
(74)

with some numerical coefficients ${}_{k,m}\kappa_{j,p,q,n}$ and
${}_{k,m}\lambda_{j,p,q,n}$ that depend on the relative integer
$j\in\mathbb{Z}$ and natural integers $p,q,n\in\mathbb{N}$, and where we have
defined for notational convenience $\overline{\mathrm{Li}}_{0}(z)\equiv 1$ and
$\overline{\mathrm{Li}}_{n}(z)\equiv\mathrm{Li}_{n}(z)$ for $n\geqslant 1$. We
also recall that $\mathrm{Li}_{1}(z)=-\ln(1-z)$.

The coefficients ${}_{k,m}\kappa_{j,p,q,n}$ and ${}_{k,m}\lambda_{j,p,q,n}$
represent six-dimensional matrices of numerical coefficients and cannot be
presented. Instead, we provide the explicit expressions of all the kernels
needed for the computation of the tails-of-memory in the Supplementary
Material [63]. For the purpose of clarity, we give two explicit examples of
such kernels, where we choose to look at $\overline{K}$ for $m=2$, $\ell=4$
and $k=1,2$:

$\displaystyle\underset{1,2}{\ \overline{K_{4}}}(\rho,\tau)=$
$\displaystyle-\frac{1}{2}\,\mathrm{Li}_{2}\left(-\frac{\tau}{\rho}\right)+\frac{7}{12}\,\mathrm{Li}_{1}\left(-\frac{\tau}{\rho}\right)-\frac{1}{4}\ln^{2}\left(\frac{\rho}{2r_{0}}\right)+\frac{1}{2}\ln\left(\frac{\rho}{2r_{0}}\right)\ln\left(\frac{\tau}{2r_{0}}\right)$
$\displaystyle-\frac{7}{12}\ln\left(\frac{\rho}{2r_{0}}\right)+\frac{25}{12}\ln\left(\frac{\tau}{2r_{0}}\right)+\frac{9}{4}-\frac{7}{12}\,\frac{\rho}{\tau}+\frac{35}{24}\left(\frac{\rho}{\tau}\right)^{2}$
$\displaystyle+\left(\frac{\rho}{\tau}\right)^{3}\left[\frac{7}{3}\,\mathrm{Li}_{1}\left(-\frac{\tau}{\rho}\right)+\frac{7}{4}\right]+\frac{7}{4}\,\left(\frac{\rho}{\tau}\right)^{4}\mathrm{Li}_{1}\left(-\frac{\tau}{\rho}\right)\,,$
(75a) $\displaystyle\underset{2,2}{\ \overline{K_{4}}}(\rho,\tau)=$
$\displaystyle\
\frac{1}{\tau}\Bigg{\\{}\frac{1}{6}-\frac{5}{6}\,\frac{\rho}{\tau}+\left(\frac{\rho}{\tau}\right)^{2}\left[-3\,\mathrm{Li}_{1}\left(-\frac{\tau}{\rho}\right)-5\right]$
$\displaystyle\qquad+\left(\frac{\rho}{\tau}\right)^{3}\left[-7\,\mathrm{Li}_{1}\left(-\frac{\tau}{\rho}\right)-4\right]-4\left(\frac{\rho}{\tau}\right)^{4}\,\mathrm{Li}_{1}\left(-\frac{\tau}{\rho}\right)\Bigg{\\}}\,.$
(75b)

Notice that despite the structure in $\rho/\tau$ of these kernels, one can
check explicitly that they are indeed integrable in a vicinity of
$(\rho,\tau)=(0,0)$.

###  V.2 Raw expression of the radiative quadrupole in the radiative
construction

We have applied the MPM construction of the metric in radiative gauge as
described in Sec. II, together with the previous integration techniques, to
the computation of the three cubic interactions
$\mathrm{M}^{2}\times\overline{\mathrm{M}}_{ij}$,
$\mathrm{M}\times\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}$ and
$\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}$.

The first interaction, $\mathrm{M}^{2}\times\overline{\mathrm{M}}_{ij}$, is
the tail-of-tail which enters at 3PN order and was already known in the
harmonic construction [20]. We have, on the one hand, computed the asymptotic
waveform for this interaction in the radiative construction using the
integration machinery developed in Sec. IV, and, on the other hand, used
standard integration techniques [28, 20, 21, 22] applied to the radiative
algorithm, and checked that both methods yielded identical results. We then
verified that these results, performed in the radiative construction, could be
independently retrieved from the known result in the harmonic construction
[20] solely using a moment redefinition, as explained in Ref. [44]. In the
radiative gauge, we find

$\mathcal{U}_{ij}^{\mathrm{M}^{2}\times\overline{\mathrm{M}}_{ij}}=\frac{2G^{2}\mathrm{M}^{2}}{c^{6}}\int_{0}^{+\infty}\\!\\!\mathrm{d}\tau\left[\ln^{2}\left(\frac{\tau}{2r_{0}}\right)+\frac{107}{42}\ln\left(\frac{\tau}{2r_{0}}\right)+\frac{40037}{8810}\right]\overline{\mathrm{M}}_{ij}^{(5)}(u-\tau)\,.$
(76)

For the $\mathrm{M}\times\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}$
interaction, which only enters at 4PN, the fact that the current dipole moment
$\mathrm{S}_{i}$ (or total angular momentum) is constant greatly simplifies
the computation, and _a priori_ allows for the use of standard techniques [28,
20, 22]. However, in some cases, those techniques break down, since some of
the numerical constants introduced are ill-defined, as has been mentionned in
[22]. Therefore, we have computed for the first time the full
$\mathrm{M}\times\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}$ waveform
using our new integration method, which reads in radiative gauge,

$\mathcal{U}_{ij}^{\mathrm{M}\times\mathrm{S}_{i}\times\overline{\mathrm{M}}_{ij}}=-\frac{2G^{2}\mathrm{M}}{c^{8}}\,\mathrm{S}_{a|\langle
i}\int_{0}^{+\infty}\mathrm{d}\tau\,\left[\ln\left(\frac{\tau}{2r_{0}}\right)+\frac{5381}{5670}\right]\overline{\mathrm{M}}_{j\rangle
a}^{(6)}(u-\tau)\,.$ (77)

Finally, the
$\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}$
interaction is the genuine tail-of-memory entering at 4PN; the double time-
dependence within the two associated quadrupole moments strongly complicates
the situation. We have performed this calculation only in the asymptotic limit
$r\to+\infty$. Since we are using the radiative construction [57] the waveform
is automatically free of far-zone logarithms, and can straightforwardly be
projected in a TT gauge, where the radiative quadrupole moment
$\mathcal{U}_{ij}$ is extracted in a standard way. Here we first present the
“raw” result, obtained in a direct manner from our integration formulae, in
terms of the radiative canonical moments. In Sec. V.3 we describe a very
efficient simplification method, and the resulting radiative quadrupole will
be presented in Sec. VI in terms of the harmonic canonical moments. The “raw”
result is however important and worthy to be presented, as it is the direct
result of our integration scheme.

In order to present the raw result, we define the following functionals of two
time-derivatives of quadrupole moments $F$ and $G$ for the cases $k=1,2$:

$\displaystyle\underset{k,m}{\ \overline{\Psi}_{\ell}}[F,G]$
$\displaystyle\equiv\int_{0}^{+\infty}\mathrm{d}\rho\,G(u-\rho)\int_{0}^{+\infty}\mathrm{d}\tau\,F(u-\rho-\tau)\underset{k,m}{\
\overline{K}_{\ell}}(\rho,\tau)\,,$ (78a) $\displaystyle\underset{k,m}{\
\overline{\chi}_{\ell}}[F,G]$
$\displaystyle\equiv\int_{0}^{+\infty}\mathrm{d}\rho\,G(u-\rho)\int_{0}^{+\infty}\mathrm{d}\tau\,F(u-\rho-\tau)\underset{k,m}{\
\overline{L}_{\ell}}(\rho,\tau)\,,$ (78b)

where the kernels are given by Eqs. (51) in the case $k=2$ and by Eqs. (55)
and (56b) in the case $k=1$. As explained in Sec. IV.2.2 we do not need the
cases $k\geqslant 3$.

With these definitions in hand, we decompose the cubic interaction
$\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}$
as follows: we separate terms which are purely instantaneous (_i.e._ local-in-
time), from terms whose non-locality is rather simple and looks like that for
the ordinary tails or tails-of-tail, and finally from the genuine and much
more intricated tail-of-memory (ToM) integrals. We thus write our raw result,
expressed in terms of the radiative canonical moments, as

$\mathcal{U}_{ij}^{\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}=\mathcal{U}_{ij}^{\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}{\bigg{|}}_{\text{inst}}+\mathcal{U}_{ij}^{\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}{\bigg{|}}_{\text{tail}}+\mathcal{U}_{ij}^{\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}{\bigg{|}}_{\text{ToM}}\,.$
(79)

  * 1.

The instantaneous terms (depending only on the current time $u$) are given by

$\displaystyle\mathcal{U}_{ij}^{\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}{\bigg{|}}_{\text{inst}}$
$\displaystyle=\frac{G^{2}\mathrm{M}}{c^{8}}\biggl{[}\frac{3362032}{165375}\overline{\mathrm{M}}_{a\langle
i}^{(3)}\overline{\mathrm{M}}_{j\rangle
a}^{(3)}+\frac{2014871}{165375}\overline{\mathrm{M}}_{a\langle
i}^{(2)}\overline{\mathrm{M}}_{j\rangle a}^{(4)}$
$\displaystyle\qquad\quad-\frac{5766241}{165375}\overline{\mathrm{M}}_{a\langle
i}^{(1)}\overline{\mathrm{M}}_{j\rangle
a}^{(5)}-\frac{114454}{23625}\overline{\mathrm{M}}_{a\langle
i}\overline{\mathrm{M}}_{j\rangle a}^{(6)}\biggr{]}\,.$ (80a)

  * 2.

The non-local tail like terms involve some logarithmic kernels, and naturally
the quadrupole moment of the usual tail terms is here replaced by a
combination $\propto\overline{\mathrm{M}}_{a\langle
i}^{(n)}\,\overline{\mathrm{M}}_{j\rangle a}^{(p)}$ evaluated at any time
$u-\tau$ in the past. We have

$\displaystyle\mathcal{U}_{ij}^{\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}{\bigg{|}}_{\text{tail}}=$
$\displaystyle-\frac{48}{175}\frac{G^{2}\mathrm{M}}{c^{8}}\int_{0}^{+\infty}\mathrm{d}\tau\left[\ln^{2}\left(\frac{\tau}{2r_{0}}\right)+\frac{5111}{840}\ln\left(\frac{\tau}{2r_{0}}\right)\right](\overline{\mathrm{M}}_{a\langle
i}\,\overline{\mathrm{M}}_{j\rangle a}^{(7)})(u-\tau)$
$\displaystyle-\frac{192}{175}\frac{G^{2}\mathrm{M}}{c^{8}}\int_{0}^{+\infty}\mathrm{d}\tau\left[\ln^{2}\left(\frac{\tau}{2r_{0}}\right)+\frac{475861}{30240}\ln\left(\frac{\tau}{2r_{0}}\right)\right](\overline{\mathrm{M}}_{a\langle
i}^{(1)}\,\overline{\mathrm{M}}_{j\rangle a}^{(6)})(u-\tau)$
$\displaystyle-\frac{3154}{315}\frac{G^{2}\mathrm{M}}{c^{8}}\int_{0}^{+\infty}\mathrm{d}\tau\ln\left(\frac{\tau}{2r_{0}}\right)(\overline{\mathrm{M}}_{a\langle
i}^{(2)}\,\overline{\mathrm{M}}_{j\rangle a}^{(5)})(u-\tau)$
$\displaystyle+\frac{478}{35}\frac{G^{2}\mathrm{M}}{c^{8}}\int_{0}^{+\infty}\mathrm{d}\tau\ln\left(\frac{\tau}{2r_{0}}\right)(\overline{\mathrm{M}}_{a\langle
i}^{(3)}\,\overline{\mathrm{M}}_{j\rangle a}^{(4)})(u-\tau)\,.$ (80b)

  * 3.

The genuine tail-of-memory part is fully specified by the bilinear functionals
of the two quadrupole moments defined by (78). Thus we can write

$\displaystyle\mathcal{U}_{ij}^{\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}{\bigg{|}}_{\text{ToM}}$
$\displaystyle=\frac{G^{2}\mathrm{M}}{c^{8}}\sum_{m,\ell,n}\bigg{\\{}\mathcal{A}_{m,\ell}^{n}\,\underset{1,m}{\
\overline{\Psi}_{\ell}}\bigl{[}\overline{\mathrm{M}}_{a\langle
i}^{(n)},\overline{\mathrm{M}}_{j\rangle
a}^{(8-n)}\bigr{]}+\mathcal{B}_{m,\ell}^{n}\,\underset{2,m}{\
\overline{\Psi}_{\ell}}\bigl{[}\overline{\mathrm{M}}_{a\langle
i}^{(n)},\overline{\mathrm{M}}_{j\rangle a}^{(7-n)}\bigr{]}$
$\displaystyle\qquad\quad+\mathcal{C}_{n,\ell}^{n}\,\underset{1,m}{\
\overline{\chi}_{\ell}}\bigl{[}\overline{\mathrm{M}}_{a\langle
i}^{(n)},\overline{\mathrm{M}}_{j\rangle
a}^{(8-n)}\bigr{]}+\mathcal{D}_{m,\ell}^{n}\,\underset{2,m}{\
\overline{\chi}_{\ell}}\bigl{[}\overline{\mathrm{M}}_{a\langle
i}^{(n)},\overline{\mathrm{M}}_{j\rangle a}^{(7-n)}\bigr{]}\bigg{\\}}\,,$
(80c)

where the purely numerical coefficients $\mathcal{A}_{m,\ell}^{n}$,
$\mathcal{B}_{m,\ell}^{n}$, $\mathcal{C}_{m,\ell}^{n}$ and
$\mathcal{D}_{m,\ell}^{n}$ in front of each of these integrals are provided in
the Tables 1(d), and the functionals ${}_{k,m}\overline{\Psi}_{\ell}$ and
${}_{k,m}\overline{\chi}_{\ell}$ are defined in terms of the kernels in (78).
Thus our complete results follow from these Tables together with the explicit
expressions of the kernel functions provided in the Supplementary Material
[63].

###  V.3 Simplification method

In this Section, we implement a method for simplifying the expression of the
pure _tail-of-memory_ part of the radiative quadrupole, given by (3.). The
idea is to rexpress everything as only one double integral over the two
quadrupoles and a single kernel, modulo some easy surface terms. We thus
alternatively integrate by parts the $\rho$ and $\tau$ variables of (78), so
as to transfer all the time derivatives on the quadrupole moment represented
by $F(u-\rho-\tau)$, _i.e._ the left slot in the functionals (78).

First we observe in (3.) that when $k=1$, the two quadrupoles have
respectively $n$ and $8-n$ time derivatives, while in the case $k=2$, they
instead have $n$ and $7-n$ derivatives. We first uniformize this by
transforming the case $k=2$ with the formula

$\displaystyle\int_{0}^{+\infty}\mathrm{d}\rho\,$
$\displaystyle\overline{\mathrm{M}}_{a\langle
i}^{(n)}(u-\rho)\int_{0}^{+\infty}\mathrm{d}\tau\,\overline{\mathrm{M}}_{j\rangle
a}^{(7-n)}(u-\rho-\tau)\,f(\rho,\tau)$
$\displaystyle=\int_{0}^{+\infty}\mathrm{d}\rho\,\overline{\mathrm{M}}_{a\langle
i}^{(n)}(u-\rho)\int_{0}^{+\infty}\mathrm{d}\tau\,\overline{\mathrm{M}}_{j\rangle
a}^{(8-n)}(u-\rho-\tau)\,\bigl{(}\partial_{\tau}^{-1}f\bigr{)}(\rho,\tau)\,,$
(81)

where we have introduced the $\tau$-antiderivative which vanishes at $\tau=0$,
defined for any function $f(\rho,\tau)$ with adequate regularity properties as

$\bigl{(}\partial_{\tau}^{-1}f\bigr{)}(\rho,\tau)\equiv\int_{0}^{\tau}\mathrm{d}\lambda\,f(\rho,\lambda)\,.$
(82)

After performing this operation, we are left with integrals of the moments
that only have $n$ and $8-n$ time derivatives, and the only possible cases are
$n\in[\\![0,4]\\!]$. Next we integrate by parts so as to be left with
integrals with only $0$ and $8$ time derivatives. However this operation
yields some integrals over $\rho$ that are separately divergent at the lower
bound $\rho=0$. To cure this, we introduce a regularization and replace the
$0$ at the lower bound by some small $\epsilon$, and restrict attention to the
$\epsilon\rightarrow 0$ expansion. In the end, we will verify that the final
result has a finite $\epsilon\rightarrow 0$ limit. With this caveat in mind,
we will simplify all our integrals with the formula (valid for
$n\in[\\![0,4]\\!]$)

$\displaystyle\int_{\epsilon}^{+\infty}\mathrm{d}\rho\,\overline{\mathrm{M}}_{a\langle
i}^{(n)}(u-\rho)\int_{0}^{+\infty}\mathrm{d}\tau\,\overline{\mathrm{M}}_{j\rangle
a}^{(8-n)}(u-\rho-\tau)\,f(\rho,\tau)$
$\displaystyle\qquad=\int_{\epsilon}^{+\infty}\mathrm{d}\rho\,\overline{\mathrm{M}}_{a\langle
i}(u-\rho)\int_{0}^{+\infty}\mathrm{d}\tau\,\overline{\mathrm{M}}_{j\rangle
a}^{(8)}(u-\rho-\tau)\,\mathcal{O}_{\rho,\tau}^{\,n}\bigl{[}f(\rho,\tau)\bigr{]}$
$\displaystyle\qquad\quad+\sum_{k=0}^{n-1}\overline{\mathrm{M}}_{a\langle
i}^{(k)}(u-\epsilon)\int_{0}^{+\infty}\mathrm{d}\tau\,\overline{\mathrm{M}}_{j\rangle
a}^{(7-k)}(u-\tau-\epsilon)\bigg{[}\mathcal{O}_{\rho,\tau}^{\,n-1-k}\bigl{[}f(\rho,\tau)\bigr{]}\bigg{]}_{\rho=\epsilon}\,,$
(83)

where moments involving $\epsilon$ should be Taylor-expanded when
$\epsilon\rightarrow 0$, the second term is evaluated at $\rho=\epsilon$, and
we have introduced the differential operator $\mathcal{O}_{\rho,\tau}$ defined
by

$\mathcal{O}_{\rho,\tau}\bigl{[}f(\tau,\rho)\bigr{]}\equiv\partial_{\tau}^{-1}\partial_{\rho}\bigl{[}f(\tau,\rho)\bigr{]}-f(\tau,\rho)\,,$
(84)

together with its iterations
$\mathcal{O}^{\,n}_{\rho,\tau}\equiv\mathcal{O}_{\rho,\tau}\cdots\mathcal{O}_{\rho,\tau}$.
Applying the formula (V.3) to all the terms composing the tails-of-memory
(3.), we arrive at a unique master double integral:

$\displaystyle\mathcal{U}_{ij}^{\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}{\bigg{|}}_{\text{ToM}}=\mathrm{M}\int_{\epsilon}^{+\infty}\mathrm{d}\rho\,\overline{\mathrm{M}}_{a\langle
i}(u-\rho)\int_{0}^{+\infty}\mathrm{d}\tau\,\overline{\mathrm{M}}_{j\rangle
a}^{(8)}(u-\rho-\tau)\,\Omega(\rho,\tau)+\mathcal{S}_{\epsilon}\,,$ (85)

where $\Omega(\rho,\tau)$ denotes some new kernel function and
$\mathcal{S}_{\epsilon}$ are all the surface terms coming from the second line
of (V.3), which carry at most one integral and that simplify drastically in
the $\epsilon\rightarrow 0$ expansion. Since the expression is long we do not
show here the result for $\mathcal{S}_{\epsilon}$.

At this stage, we would _a priori_ expect $\Omega(\rho,\tau)$ to have a very
complicated structure akin to (74), and in particular to involve many
polylogarithms. But instead, we find the following simple expression without
any polylogarithms:

$\displaystyle\Omega(\rho,\tau)$
$\displaystyle=\frac{7613764}{165375}-\frac{1024076}{18375}\,\frac{\tau}{\rho}-\frac{2074}{63}\left(\frac{\tau}{\rho}\right)^{2}-\frac{104}{15}\left(\frac{\tau}{\rho}\right)^{3}$
$\displaystyle+\frac{634076}{55125}\ln\left(\frac{\rho}{2r_{0}}\right)+\frac{384}{175}\,\frac{\tau}{\rho}\ln\left(\frac{\rho}{2r_{0}}\right)-\frac{144}{175}\ln\left(\frac{\rho}{2r_{0}}\right)^{2}+\frac{8}{7}\ln\left(\frac{\tau}{2r_{0}}\right)\,.$
(86)

Note that if we changed even a single coefficient in Table 1(d), the
cancellation of polylogarithms would not occur in general, and we would be
left with a much more complicated expression for $\Omega(\rho,\tau)$. When
inserting $\Omega(\rho,\tau)$ into the first term in (85), we can integrate by
parts so as to remove all the powers of $\tau/\rho$. This introduces poles in
$\epsilon$, and powers of the logarithms of $\epsilon$, but we have checked
that these poles and divergences exactly cancel when adding the surface terms
$\mathcal{S}_{\epsilon}$, namely the second term in (85). Putting all of this
together and taking the $\epsilon\rightarrow 0$ limit, we find that

$\displaystyle\mathcal{U}_{ij}^{\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}{\bigg{|}}_{\text{ToM}}$
$\displaystyle=\frac{8}{7}\,\mathrm{M}\int_{0}^{+\infty}\mathrm{d}\rho\,\overline{\mathrm{M}}_{a\langle
i}^{(4)}(u-\rho)\int_{0}^{+\infty}\mathrm{d}\tau\,\ln\left(\frac{\tau}{2r_{0}}\right)\overline{\mathrm{M}}_{j\rangle
a}^{(4)}(u-\rho-\tau)$
$\displaystyle+\frac{48}{175}\,\mathrm{M}\int_{0}^{+\infty}\mathrm{d}\tau\,\overline{\mathrm{M}}_{a\langle
i}\overline{\mathrm{M}}^{(7)}_{j\rangle
a}(u-\tau)\left[\ln^{2}\left(\frac{\tau}{2r_{0}}\right)+\frac{1243}{420}\ln\left(\frac{\tau}{2r_{0}}\right)\right]$
$\displaystyle+\frac{192}{175}\,\mathrm{M}\int_{0}^{+\infty}\mathrm{d}\tau\,\overline{\mathrm{M}}^{(1)}_{a\langle
i}\overline{\mathrm{M}}^{(6)}_{j\rangle
a}(u-\tau)\left[\ln^{2}\left(\frac{\tau}{2r_{0}}\right)+\frac{186743}{15120}\ln\left(\frac{\tau}{2r_{0}}\right)\right]$
$\displaystyle+\frac{1084}{315}\,\mathrm{M}\int_{0}^{+\infty}\mathrm{d}\tau\ln\left(\frac{\tau}{2r_{0}}\right)\overline{\mathrm{M}}^{(2)}_{a\langle
i}\overline{\mathrm{M}}^{(5)}_{j\rangle a}(u-\tau)$
$\displaystyle-\frac{104}{5}\,\mathrm{M}\int_{0}^{+\infty}\mathrm{d}\tau\ln\left(\frac{\tau}{2r_{0}}\right)\overline{\mathrm{M}}^{(3)}_{a\langle
i}\overline{\mathrm{M}}^{(4)}_{j\rangle a}(u-\tau)$
$\displaystyle+\frac{4}{7}\,\mathrm{M}\,\overline{\mathrm{M}}_{a\langle
i}\int_{0}^{+\infty}\mathrm{d}\tau\,\overline{\mathrm{M}}^{(7)}_{j\rangle
a}(u-\tau)\left[\ln\left(\frac{\tau}{2r_{0}}\right)+\frac{15667}{10500}\right]$
$\displaystyle-\frac{20}{7}\,\mathrm{M}\,\overline{\mathrm{M}}^{(1)}_{a\langle
i}\int_{0}^{+\infty}\mathrm{d}\tau\,\overline{\mathrm{M}}^{(6)}_{j\rangle
a}(u-\tau)\left[\ln\left(\frac{\tau}{2r_{0}}\right)-\frac{3590791}{472500}\right]$
$\displaystyle-\frac{16}{7}\,\mathrm{M}\,\overline{\mathrm{M}}^{(2)}_{a\langle
i}\int_{0}^{+\infty}\mathrm{d}\tau\,\overline{\mathrm{M}}^{(5)}_{j\rangle
a}(u-\tau)\left[\ln\left(\frac{\tau}{2r_{0}}\right)+\frac{981149}{94500}\right]$
$\displaystyle-\frac{3901382}{165375}\mathrm{M}\,\overline{\mathrm{M}}^{(3)}_{a\langle
i}\overline{\mathrm{M}}^{(3)}_{j\rangle a}\,.$ (87)

There is only one doubly-integrated term left, _cf._ the first line of (V.3),
which can be seen as the tails-of-memory contribution properly speaking,
whereas the structure of the other singly-integrated terms are more akin to
simpler tails-of-tails.

The terms proportional to $\ln^{2}(\frac{\tau}{2r_{0}})$ are exactly
compensated by those arising in the tail part of
$\mathcal{U}_{ij}^{\mathrm{M}\times\overline{\mathrm{M}}_{ij}\times\overline{\mathrm{M}}_{ij}}$
given by (2.). Notice that not all possible terms allowed by a dimensionality
argument are present in (V.3): for instance there is no term

$\displaystyle\propto\mathrm{M}\,\overline{\mathrm{M}}_{a\langle
i}^{(4)}(u)\int_{0}^{+\infty}\mathrm{d}\tau\ln\left(\frac{\tau}{2r_{0}}\right)\overline{\mathrm{M}}_{j\rangle
a}^{(3)}(u-\tau)\,,$ (88)

which can be understood from the structure of the cubic source term (23)
corresponding to the diagram of Fig. 1(a).

###  V.4 Testing the integration method

The simplication method in Sec V.3 is also very useful to test our integration
method. Indeed, let us consider the integration of a typical tail-of-memory
source term, but which we multiply by the factor $B$. For simplicity, we
choose $k=3$, $m=0$ and $\ell=0$:

$\mathcal{I}\equiv\mathop{\mathrm{FP}}_{B=0}\Box^{-1}_{\text{ret}}\biggl{[}B\left(\frac{r}{r_{0}}\right)^{B}\frac{G(t-r)}{r^{3}}\int_{1}^{+\infty}\mathrm{d}x\,Q_{0}(x)F(t-rx)\biggr{]}\,.$
(89)

This can be computed asymptotically when $r\to+\infty$ using the arsenal of
techniques developed in Sec. IV, and we find

$\displaystyle\mathcal{I}=\frac{1}{r}\biggl{[}\Delta-\frac{1}{2}\int_{0}^{+\infty}\mathrm{d}\tau\ln\left(\frac{\tau}{2r_{0}}\right)G(u-\tau)F^{(1)}(u-\tau)\biggr{]}+o\left(\frac{1}{r}\right)\,,$
(90a) with the rather cumbersome combination of the functionals (78a):
$\displaystyle\Delta\equiv-\frac{1}{3}\,\underset{1,2}{\overline{\Psi}_{0}}[F^{(2)},G]+\frac{1}{3}\,\underset{1,0}{\overline{\Psi}_{0}}[F^{(2)},G]+\underset{1,1}{\overline{\Psi}_{0}}[F^{(1)},G^{(1)}]-\underset{1,0}{\overline{\Psi}_{0}}[F^{(1)},G^{(1)}]\,.$
(90b)

But, on the other hand, the presence of the factor $B$ has the effect of only
selecting the pole in the $B$-expansion of the integrated source. This allows
us to perform first the near-zone expansion $r\rightarrow 0$ of the source,
and then to integrate term by term this expansion using the formulae in [43,
44]. Performing the integration in this manner, we find instead

$\mathcal{I}=-\frac{G(u)}{r}\bigg{[}F(u)+\int_{0}^{+\infty}\mathrm{d}\tau\ln\left(\frac{\tau}{2r_{0}}\right)F^{(1)}(u-\tau)\biggr{]}+o\left(\frac{1}{r}\right)\,.$
(91)

The two expressions (90) and (91) must be identical, therefore we have found a
non-trivial relationship between the functionals
${}_{k,m}{\overline{\Psi}_{\ell}}[F,G]$ which must absolutely be satisfied if
our integration method is correct. Applying the simplification method
described in Sec V.3, we can prove that these two expressions are indeed
identical. This constitutes a strong test confirming simultaneously the
soundness of our integration method and of our simplification method. We have
repeated this test with many other such sources; all were successful.

##  VI Results

We present our final result, which is the radiative quadrupole moment
$\mathcal{U}_{ij}$ parametrizing the asymptotic waveform, expressed in terms
of the harmonic canonical moments up to 4.5PN order, including the tail-of-
memory contribution.

Recall that in the previous sections, we have only worked in the radiative
construction, and hence with radiative canonical moments denoted
$\overline{\mathrm{M}}_{L}$ and $\overline{\mathrm{S}}_{L}$. This is perfectly
legitimate, but until now, all the results in the literature [14, 20, 28, 65,
22] were presented in terms of the harmonic construction and hence of harmonic
canonical moments $\mathrm{M}_{L}$ and $\mathrm{S}_{L}$. Even more crucially,
it is the harmonic canonical moments that were computed in previous works in
terms of source parameters [39, 41, 42, 43]. For these reasons, we have to
express our final result in terms of the harmonic canonical moments
$\mathrm{M}_{L}$ and $\mathrm{S}_{L}$. At 4PN order, this can be done thanks
to the moment redefinition derived in [44], namely

$\displaystyle\overline{\mathrm{M}}_{ij}=\mathrm{M}_{ij}$ $\displaystyle\
+\frac{G\,\mathrm{M}}{c^{3}}\,\mathrm{M}_{ij}^{(1)}\left[-\frac{26}{15}+2\ln\left(\frac{r_{0}}{b_{0}}\right)\right]+\frac{G^{2}\mathrm{M}^{2}}{c^{6}}\,\mathrm{M}_{ij}^{(2)}\left[\frac{124}{45}-\frac{52}{15}\ln\left(\frac{r_{0}}{b_{0}}\right)+2\ln^{2}\left(\frac{r_{0}}{b_{0}}\right)\right]$
(92) $\displaystyle\
+\frac{G^{2}\mathrm{M}}{c^{8}}\bigg{[}-\frac{8}{21}\,\mathrm{M}_{a\langle
i}\mathrm{M}_{j\rangle a}^{(4)}-\frac{8}{7}\,\mathrm{M}^{(1)}_{a\langle
i}\mathrm{M}_{j\rangle a}^{(3)}-\frac{8}{9}\,\mathrm{M}^{(3)}_{a\langle
i}\mathrm{S}_{j\rangle|a}\bigg{]}+\mathcal{O}\left(\frac{1}{c^{9}}\right)\,,$

where $b_{0}$ is the arbitrary gauge constant introduced in Eq. (6), and where
$c$ is henceforth reintroduced. The moment redefinition (92) should be
inserted into the linear and quadratic contributions to the radiative
quadrupole written in terms of the radiative metric (20), so as to yield a
correction at cubic order.

This correction should be understood in the following sense: let
$\mathcal{U}_{ij}^{\text{cubic,rad}}[\overline{\mathrm{M}}_{L},\overline{\mathrm{S}}_{L}]$
be the functional expression for the cubic part of the radiative quadrupole in
terms of the radiative moments, as worked out in the present paper. Then its
counterpart in terms of the harmonic moments, simply denoted
$\mathcal{U}_{ij}^{\text{cubic}}[\mathrm{M}_{L},\mathrm{S}_{L}]$, is given by

$\mathcal{U}_{ij}^{\text{cubic}}[\mathrm{M}_{L},\mathrm{S}_{L}]=\mathcal{U}_{ij}^{\text{cubic,rad}}[\mathrm{M}_{L},\mathrm{S}_{L}]+\delta\mathcal{U}_{ij}^{\text{cubic,rad}}[\mathrm{M}_{L},\mathrm{S}_{L}]\,,$
(93)

where $\mathcal{U}_{ij}^{\text{cubic,rad}}[\mathrm{M}_{L},\mathrm{S}_{L}]$
means that we simply subsituted the radiative moments with the harmonic
moments in the radiative functional, and where
$\delta\mathcal{U}_{ij}^{\text{cubic,rad}}[\mathrm{M}_{L},\mathrm{S}_{L}]$
represents the correction to be applied due to the moment redefinition (92).
We find

$\displaystyle\delta\mathcal{U}_{ij}^{\text{cubic,rad}}$
$\displaystyle=\frac{G^{2}\mathrm{M}^{2}}{c^{6}}\Bigg{\\{}\left[-\frac{52}{15}+4\ln\left(\frac{r_{0}}{b_{0}}\right)\right]\int_{0}^{+\infty}\mathrm{d}\tau\,\ln\left(\frac{\tau}{2r_{0}}\right)\mathrm{M}_{ij}^{(5)}(u-\tau)$
$\displaystyle\qquad\qquad+\left[-\frac{257}{75}+\frac{11}{3}\ln\left(\frac{r_{0}}{b_{0}}\right)+2\ln^{2}\left(\frac{r_{0}}{b_{0}}\right)\right]\mathrm{M}_{ij}^{(4)}\Bigg{\\}}$
$\displaystyle+\frac{G^{2}\mathrm{M}}{c^{8}}\Biggl{\\{}\left[-\frac{16}{105}-\frac{8}{7}\ln\left(\frac{r_{0}}{b_{0}}\right)\right]\mathrm{M}_{a\langle
i}^{(3)}\mathrm{M}_{j\rangle
a}^{(3)}+\left[-\frac{14}{15}-2\ln\left(\frac{r_{0}}{b_{0}}\right)\right]\mathrm{M}_{a\langle
i}^{(2)}\mathrm{M}_{j\rangle a}^{(4)}$
$\displaystyle\qquad\qquad+\left[-\frac{32}{35}-\frac{8}{7}\ln\left(\frac{r_{0}}{b_{0}}\right)\right]\mathrm{M}_{a\langle
i}^{(1)}\mathrm{M}_{j\rangle
a}^{(5)}+\left[-\frac{22}{35}+\frac{2}{7}\ln\left(\frac{r_{0}}{b_{0}}\right)\right]\mathrm{M}_{a\langle
i}\mathrm{M}_{j\rangle a}^{(6)}\Biggr{\\}}$
$\displaystyle+\frac{G^{2}\mathrm{M}}{c^{8}}\mathrm{S}_{a|\langle
i}\left[\frac{22}{15}-\frac{2}{3}\ln\left(\frac{r_{0}}{b_{0}}\right)\right]\mathrm{M}_{j\rangle
a}^{(5)}\,.$ (94)

We finally present the full radiative quadrupole, parametrized by the harmonic
canonical moments $\mathrm{M}_{L}$ and $\mathrm{S}_{L}$, including all
quadratic, cubic and quartic contributions that contribute up to the 4.5PN
order. It reads

$\mathcal{U}_{ij}=\mathcal{U}_{ij}^{\text{linear}}+\,\mathcal{U}_{ij}^{\text{quadratic}}+\,\mathcal{U}_{ij}^{\text{cubic}}+\,\mathcal{U}_{ij}^{\text{quartic}}\,,$
(95)

where the linear contribution is just
$\mathcal{U}_{ij}^{\text{linear}}=\mathrm{M}_{ij}^{(2)}$, the quadratic part
$\mathcal{U}_{ij}^{\text{quadratic}}$ can be entirely found in Eqs. (4.4-5-6)
of [65], the quartic piece $\mathcal{U}_{ij}^{\text{quartic}}$ is given in Eq.
(4.8) of [22], and where the cubic contributions

$\mathcal{U}_{ij}^{\text{cubic}}=\mathcal{U}_{ij}^{\mathrm{M}^{2}\times\mathrm{M}_{ij}}+\mathcal{U}_{ij}^{\mathrm{M}\times\mathrm{S}_{i}\times\mathrm{M}_{ij}}+\mathcal{U}_{ij}^{\mathrm{M}\times\mathrm{M}_{ij}\times\mathrm{M}_{ij}}\,,$
(96)

decompose into the tail-of-tail [20, 21], and, both new with this paper, the
spin-quadrupole tail and the tail-of-memory. We have

$\displaystyle\mathcal{U}_{ij}^{\mathrm{M}^{2}\times\mathrm{M}_{ij}}$
$\displaystyle=\frac{2G^{2}\mathrm{M}^{2}}{c^{6}}\\!\\!\int_{0}^{+\infty}\\!\\!\\!\mathrm{d}\tau\,\mathrm{M}_{ij}^{(5)}(u-\tau)\\!\left[\ln^{2}\left(\frac{\tau}{2b_{0}}\right)+\frac{11}{6}\ln\left(\frac{\tau}{2b_{0}}\right)-\frac{107}{105}\ln\left(\frac{\tau}{2r_{0}}\right)+\frac{124627}{44100}\right]\\!,$
(97a)
$\displaystyle\mathcal{U}_{ij}^{\mathrm{M}\times\mathrm{S}_{i}\times\mathrm{M}_{ij}}$
$\displaystyle=-\frac{2G^{2}\mathrm{M}}{3c^{8}}\,\mathrm{S}_{a|\langle
i}\int_{0}^{+\infty}\mathrm{d}\tau\,\mathrm{M}_{j\rangle
a}^{(6)}(u-\tau)\left[\ln\left(\frac{\tau}{2b_{0}}\right)+2\ln\left(\frac{\tau}{2r_{0}}\right)+\frac{1223}{1890}\right]\,,$
(97b)
$\displaystyle\mathcal{U}_{ij}^{\mathrm{M}\times\mathrm{M}_{ij}\times\mathrm{M}_{ij}}$
$\displaystyle=\frac{8G^{2}\mathrm{M}}{7c^{8}}\Bigg{\\{}\int_{0}^{+\infty}\\!\mathrm{d}\rho\,\mathrm{M}_{a\langle
i}^{(4)}(u-\rho)\int_{0}^{+\infty}\\!\mathrm{d}\tau\,\mathrm{M}_{j\rangle
a}^{(4)}(u-\rho-\tau)\left[\ln\left(\frac{\tau}{2r_{0}}\right)-\frac{1613}{270}\right]$
$\displaystyle\quad-\frac{5}{2}\int_{0}^{+\infty}\\!\mathrm{d}\tau\,(\mathrm{M}^{(3)}_{a\langle
i}\mathrm{M}^{(4)}_{j\rangle
a})(u-\tau)\left[\ln\left(\frac{\tau}{2r_{0}}\right)+\frac{3}{2}\ln\left(\frac{\tau}{2b_{0}}\right)\right]$
$\displaystyle\quad-3\int_{0}^{+\infty}\\!\mathrm{d}\tau\,(\mathrm{M}^{(2)}_{a\langle
i}\mathrm{M}^{(5)}_{j\rangle
a})(u-\tau)\left[\ln\left(\frac{\tau}{2r_{0}}\right)+\frac{11}{12}\ln\left(\frac{\tau}{2b_{0}}\right)\right]$
$\displaystyle\quad-\frac{5}{2}\int_{0}^{+\infty}\\!\mathrm{d}\tau\,(\mathrm{M}^{(1)}_{a\langle
i}\mathrm{M}^{(6)}_{j\rangle
a})(u-\tau)\left[\ln\left(\frac{\tau}{2r_{0}}\right)+\frac{3}{10}\ln\left(\frac{\tau}{2b_{0}}\right)\right]$
$\displaystyle\quad-\int_{0}^{+\infty}\\!\mathrm{d}\tau\,(\mathrm{M}_{a\langle
i}\mathrm{M}^{(7)}_{j\rangle
a})(u-\tau)\left[\ln\left(\frac{\tau}{2r_{0}}\right)-\frac{1}{4}\ln\left(\frac{\tau}{2b_{0}}\right)\right]$
$\displaystyle\quad-2\mathrm{M}^{(2)}_{a\langle
i}\int_{0}^{+\infty}\\!\mathrm{d}\tau\,\mathrm{M}^{(5)}_{j\rangle
a}(u-\tau)\left[\ln\left(\frac{\tau}{2r_{0}}\right)+\frac{27521}{5040}\right]$
$\displaystyle\quad-\frac{5}{2}\,\mathrm{M}^{(1)}_{a\langle
i}\int_{0}^{+\infty}\\!\mathrm{d}\tau\,\mathrm{M}^{(6)}_{j\rangle
a}(u-\tau)\left[\ln\left(\frac{\tau}{2r_{0}}\right)+\frac{15511}{3150}\right]$
$\displaystyle\quad+\frac{1}{2}\,\mathrm{M}_{a\langle
i}\int_{0}^{+\infty}\\!\mathrm{d}\tau\,\mathrm{M}^{(7)}_{j\rangle
a}(u-\tau)\left[\ln\left(\frac{\tau}{2r_{0}}\right)-\frac{6113}{756}\right]\,\Bigg{\\}}\,.$
(97c)

The double integral in the first line of (97a) is the genuine tail-of-memory.
However, when computing the flux, we must actually compute the time-derivative
of the radiative quadrupole. We then find that the tail-of-memory term becomes
a simple tail term in the flux, in the same manner as the quadrupole-
quadrupole memory term becomes instantaneous in the flux [28].

We can perform two important tests of our end result by studying how the two
arbitrary constants $b_{0}$ and $r_{0}$ will vanish from physical observables.

  * •

First, the constant $b_{0}$ can be eliminated at the level of the full
radiative quadrupole at 4PN order by the shift in the time coordinate

$\displaystyle u^{\prime}=u-\frac{2G\mathrm{M}}{c^{3}}\ln b_{0}\,,$ (98)

along with a Taylor expansion of the canonical quadrupole moment. This shows
that $b_{0}$ is just associated to a choice for the origin of time in the
asymptotic radiative coordinate system and is clearly unphysical.

  * •

Second, when isolating the contribution of $\ln r_{0}$ to the radiative
quadrupole for the three cubic interactions (97), we find

$\displaystyle\mathcal{U}_{ij}^{\text{cubic}}$ $\displaystyle=\ln
r_{0}\bigg{[}\frac{214}{105}\frac{G^{2}\mathrm{M}^{2}}{c^{6}}\mathrm{M}_{ij}^{(4)}$
$\displaystyle\quad\quad+\frac{G^{2}M}{c^{8}}\Big{(}4\mathrm{M}_{a\langle
i}^{(2)}\mathrm{M}_{j\rangle a}^{(4)}+\frac{32}{7}\mathrm{M}_{a\langle
i}^{(1)}\mathrm{M}_{j\rangle a}^{(5)}+\frac{4}{7}\mathrm{M}_{a\langle
i}\mathrm{M}_{j\rangle a}^{(6)}+\frac{4}{3}\mathrm{S}_{a|\langle
i}\mathrm{M}_{j\rangle a}^{(5)}\Big{)}\bigg{]}$
$\displaystyle\quad\quad+(\text{terms independent of $r_{0}$})\,.$ (99)

This expression exactly cancels the $\ln r_{0}$ terms arising in the second
time derivative of Eq. (6.1) in [42], which accounts for the contribution of
the dimensional regularization of the cubic interactions to the renormalized
canonical quadrupole moment for compact binaries. This cancellation of the
$\ln r_{0}$ terms with those in the end result of [42] is a strong indication
of the correctness of our results (97).

###### Acknowledgements.

The authors thank François Larrouturou for useful suggestions and interesting
discussions at different stages of this project. We also thank Laura Bernard,
Guillaume Faye, Quentin Henry and Stavros Mougiakakos for interesting
discussions. L.B. acknowledges the Institut de Physique Théorique (IPhT) in
CEA/Saclay for a visiting position. D.T. thanks Laura Bernard and the Kavli
Institute for Theoretical Physics (KITP) in Santa Barbara for their invitation
to participate to the “High-Precision Gravitational Waves” program (supported
in part by the National Science Foundation under Grant No. NSF PHY-1748958).

##  Appendix A Extracting $N_{L}$, $P_{L}$, $Q_{L}$ and $R_{L}$

We describe the practical method we use to extract four sets of STF moments
$N_{L}$, $P_{L}$, $Q_{L}$ and $R_{L}$ parametrizing a vector quantity
satisfying $\Box w^{\mu}=0$. Here $w^{\mu}\equiv w_{n}^{\mu}$ represents the
divergence $w_{n}^{\mu}=\partial_{\nu}u_{n}^{\mu\nu}$ of the quantity (8)
following the MPM algorithm. Having extracted these moments we can then
construct the tensor
$v_{n}^{\mu\nu}\equiv\mathcal{V}^{\mu\nu}[N_{L},P_{L},Q_{L},R_{L}]$ using the
formulae (48) in [3], which satisfies at once $\Box v_{n}^{\mu\nu}~{}=~{}0$
and $\partial_{\nu}v_{n}^{\mu\nu}=-w_{n}^{\mu}$. By definition, see (47) in
[3],

$\displaystyle w^{0}$
$\displaystyle=\sum_{\ell=0}^{+\infty}\partial_{L}\left[r^{-1}N_{L}(t-r)\right]\,,$
(100a) $\displaystyle w^{i}$
$\displaystyle=\sum_{\ell=0}^{+\infty}\partial_{iL}\left[r^{-1}P_{L}(t-r)\right]$
$\displaystyle+\sum_{\ell=1}^{+\infty}\partial_{L-1}\left[r^{-1}Q_{iL}(t-r)\right]+\sum_{\ell=1}^{+\infty}\epsilon_{iab}\partial_{aL-1}\left[r^{-1}R_{bL-1}(t-r)\right]\,.$
(100b)

Next, we define the auxiliary quantities

$\displaystyle\widetilde{w}^{i}$ $\displaystyle\equiv
w^{i}-\sum_{\ell=0}^{+\infty}\partial_{iL}\left[r^{-1}P_{L}(t-r)\right]$
$\displaystyle=\sum_{\ell=1}^{+\infty}\partial_{L-1}\left[r^{-1}Q_{iL}(t-r)\right]+\sum_{\ell=1}^{+\infty}\epsilon_{iab}\partial_{aL-1}\left[r^{-1}R_{bL-1}(t-r)\right]\,,$
(101a) $\displaystyle\overset{\approx}{w}{}^{i}$
$\displaystyle\equiv\widetilde{w}^{i}-\sum_{\ell=1}^{+\infty}\partial_{L-1}\left[r^{-1}Q_{iL}(t-r)\right]=\sum_{\ell=1}^{+\infty}\epsilon_{iab}\partial_{aL-1}\left[r^{-1}R_{bL-1}(t-r)\right]\,.$
(101b)

Using formulae from Appendix A in [56], we express the angular integrals

$\displaystyle\int\frac{\mathrm{d}\Omega}{4\pi}\,\hat{n}_{L}w^{0}$
$\displaystyle=\frac{\ell!}{(2\ell+1)!!}\,r^{\ell}\left(\frac{1}{r}\frac{\partial}{\partial
r}\right)^{\ell}\bigg{[}r^{-1}N_{L}(t-r)\bigg{]}\,,$ (102a)
$\displaystyle\int\frac{\mathrm{d}\Omega}{4\pi}\,\hat{n}_{iL}w^{i}$
$\displaystyle=\frac{(\ell+1)\ell!}{(2\ell+1)(2\ell+1)!!}\,r^{\ell+1}\left(\frac{1}{r}\frac{\partial}{\partial
r}\right)^{\ell+1}\bigg{[}r^{-1}P_{L}(t-r)\bigg{]}\,,$ (102b)
$\displaystyle\int\frac{\mathrm{d}\Omega}{4\pi}\,n_{i}\hat{n}_{L}\widetilde{w}^{i}$
$\displaystyle=\frac{\ell!}{(2\ell+1)!!}\,r^{\ell-1}\left(\frac{1}{r}\frac{\partial}{\partial
r}\right)^{\ell-1}\bigg{[}r^{-1}Q_{L}(t-r)\bigg{]}\,,$ (102c)
$\displaystyle\epsilon_{ab\langle
i}\\!\int\frac{\mathrm{d}\Omega}{4\pi}\,n_{L-1\rangle}n_{b}\overset{\approx}{w}{}^{a}$
$\displaystyle=\frac{(\ell+1)(\ell-1)!}{(2\ell+1)!!}\,r^{\ell}\left(\frac{1}{r}\frac{\partial}{\partial
r}\right)^{\ell}\bigg{[}r^{-1}R_{iL-1}(t-r)\bigg{]}\,.$ (102d)

If $w^{\mu}$ is known exactly to all order in $r$, we find that the multipole
moments can be computed using the near zone limit $r\to 0$ (with $t-r=$ const)
as

$\displaystyle N_{L}(t-r)$
$\displaystyle=\frac{(-)^{\ell}(2\ell+1)}{\ell!}\lim_{r\rightarrow
0}\bigg{[}r^{\ell+1}\int\frac{\mathrm{d}\Omega}{4\pi}\,\hat{n}_{L}w^{0}\bigg{]}\,,$
(103a) $\displaystyle P_{L}(t-r)$
$\displaystyle=\frac{(-)^{\ell+1}(2\ell+1)}{(\ell+1)!}\lim_{r\rightarrow
0}\bigg{[}r^{\ell+2}\int\frac{\mathrm{d}\Omega}{4\pi}\,\hat{n}_{iL}w^{i}\bigg{]}\,,$
(103b) $\displaystyle Q_{L}(t-r)$
$\displaystyle=\frac{(-)^{\ell-1}(2\ell+1)(2\ell-1)}{\ell!}\lim_{r\rightarrow
0}\bigg{[}r^{\ell}\int\frac{\mathrm{d}\Omega}{4\pi}\,n_{i}\hat{n}_{L}\widetilde{w}^{i}\bigg{]}\,,$
(103c) $\displaystyle R_{iL-1}(t-r)$
$\displaystyle=\frac{(-)^{\ell}(2\ell+1)}{(\ell+1)(\ell-1)!}\lim_{r\rightarrow
0}\bigg{[}r^{\ell+1}\int\frac{\mathrm{d}\Omega}{4\pi}\,\epsilon_{ab\langle
i}n_{L-1\rangle}n_{b}\overset{\approx}{w}{}^{a}\bigg{]}\,.$ (103d)

If we only know the leading order of the asymptotic expansion of $w^{\mu}$ as
$r\to+\infty$, the previous expressions cannot be used. In that case, we get
equivalent expressions for the time derivatives of the multipole moments:

$\displaystyle\overset{(\ell)}{N_{L}}(t-r)$
$\displaystyle=\frac{(-)^{\ell}(2\ell+1)!!}{\ell!}\lim_{r\to+\infty}\bigg{[}r\int\frac{\mathrm{d}\Omega}{4\pi}\,\hat{n}_{L}w^{0}\bigg{]}\,,$
(104a) $\displaystyle\overset{(\ell+1)}{P_{L}}(t-r)$
$\displaystyle=\frac{(-)^{\ell}(2\ell+1)(2\ell+1)!!}{(\ell+1)!}\lim_{r\to+\infty}\bigg{[}r\int\frac{\mathrm{d}\Omega}{4\pi}\,\hat{n}_{iL}w^{i}\bigg{]}\,,$
(104b) $\displaystyle\overset{(\ell-1)}{Q_{L}}(t-r)$
$\displaystyle=\frac{(-)^{\ell-1}(2\ell+1)!!}{\ell!}\lim_{r\to+\infty}\bigg{[}r\int\frac{\mathrm{d}\Omega}{4\pi}\,n_{i}\hat{n}_{L}\widetilde{w}^{i}\bigg{]}\,,$
(104c) $\displaystyle\overset{(\ell)}{R}_{iL-1}(t-r)$
$\displaystyle=\frac{(-)^{\ell}(2\ell+1)!!}{(\ell+1)(\ell-1)!}\lim_{r\to+\infty}\bigg{[}r\int\frac{\mathrm{d}\Omega}{4\pi}\,\epsilon_{ab\langle
i}n_{L-1\rangle}n_{b}\overset{\approx}{w}{}^{a}\bigg{]}\,.$ (104d)

##  Appendix B Proof of convergence when $k=1$ and $k=2$

In this Appendix, we assume that $k\in\\{1,2\\}$. We shall prove that in these
cases one can set $B=0$ in the general solution (44), for arbitrary smooth
functions $F$ and $G\in\mathcal{C}^{\infty}(\mathbb{R})$ which vanish
identically in the remote past, for $t\leqslant-\mathcal{T}$. Although it
should be possible to check this point using the expression of the three-
dimensional retarded integral of the source term (40), we present here a
detailed analysis directly based on the structure of the solution (44).

We must first do some manipulation so as to transfer the boundary conditions
of $F$ and $G$ into the bounds of the integrals. Thanks to the regularization
factor $\lambda^{B}$, we can always manipulate these integrals safely, since
they are defined by analytic continuation for any $B\in\mathbb{C}$ except at
some integer values including the value of interest $B=0$, at which we apply
the finite part in the end.

####  B.0.1 Case $k+j\leqslant 1$

In this case, we can perform the change of variable $\lambda=\mu+\rho/2$ on
Eq. (44b) and expand using the binomial formula:

$\displaystyle\phi_{ij}^{B}(u,r)$
$\displaystyle=\sum_{p=0}^{-k-j+1}\begin{pmatrix}-k-j+1\\\
p\end{pmatrix}\int_{0}^{+\infty}\mathrm{d}\rho\,\rho^{i+j}\left(\frac{\rho}{2}\right)^{-k-j+1-p}G(u-\rho)\int_{1}^{+\infty}\mathrm{d}x\,Q_{m}(x)$
$\displaystyle\qquad\qquad\qquad\times\int_{0}^{r}\mathrm{d}\mu\,\left(\frac{\mu+\rho/2}{r_{0}}\right)^{B}\mu^{p}\,F\left(u-\frac{\rho(x+1)}{2}-\mu(x-1)\right)\,.$
(105)

We then set $\tau=\mu(x-1)$ with $\mathrm{d}\tau=(x-1)\mathrm{d}\mu$, invert
the integrals, and find

$\displaystyle\phi_{ij}^{B}(u,r)$
$\displaystyle=\sum_{p=0}^{-k-j+1}2^{-k+j-1+p}\begin{pmatrix}-k-j+1\\\
p\end{pmatrix}\int_{0}^{+\infty}\mathrm{d}\tau\,\tau^{p}\int_{1+\tau/r}^{+\infty}\mathrm{d}x\,\frac{Q_{m}(x)}{(x-1)^{p+1}}\left(\frac{\frac{\tau}{x-1}+\rho/2}{r_{0}}\right)^{B}$
$\displaystyle\qquad\qquad\qquad\times\int_{0}^{+\infty}\mathrm{d}\rho\,\rho^{-k+i+1-p}G(u-\rho)F\left(u-\frac{\rho(x+1)}{2}-\tau\right)\,.$
(106)

Next we introduce

$\Upsilon_{k,i,j,p}(\tau)\equiv\tau^{p}\int_{1+\tau/r}^{+\infty}\mathrm{d}x\,\frac{Q_{m}(x)}{(x-1)^{p+1}}\int_{0}^{+\infty}\mathrm{d}\rho\,\rho^{-k+i+1-p}G(u-\rho)F\left(u-\frac{\rho(x+1)}{2}-\tau\right)\,,$
(107)

which is well defined, since $-k+i+1-p\geqslant i+j\geqslant 0$ ensures that
the integral in $\rho$ converges at $0$, and the behavior of
$Q_{m}(x)/(x-1)^{p+1}\sim x^{-m-p-2}$ when $x\to+\infty$ ensures that the
integral in $x$ converges at infinity.

Since $F$ vanishes identically in the remote past, so does
$\Upsilon_{k,i,j,p}(\tau)$, thus it is clearly integrable at
$\tau\rightarrow\infty$. We then bound this quantity with

$\left|\Upsilon_{k,i,j,p}(\tau)\right|\leqslant\tau^{p}\int_{1+\tau/r}^{+\infty}\mathrm{d}x\,\frac{Q_{m}(x)}{(x-1)^{p+1}}=\mathcal{O}_{\tau\rightarrow
0}\left(\ln^{2}\tau\right)\,,$ (108)

which proves integrability at the bound $\tau\rightarrow 0$.

####  B.0.2 Case $k+j\geqslant 2$

We first shuffle the order of the integrals in Eq. (44b) and obtain

$\displaystyle\phi_{ij}^{B}$
$\displaystyle=\int_{1}^{+\infty}\mathrm{d}x\,Q_{m}(x)\int_{0}^{+\infty}\mathrm{d}\lambda\,\left(\frac{\lambda}{r_{0}}\right)^{B}\lambda^{-k-j+1}\int_{0}^{2\lambda}\mathrm{d}\rho\,\,\rho^{i+j}G(u-\rho)F\bigl{[}u-\rho-\lambda(x-1)\bigr{]}$
$\displaystyle-\int_{1}^{+\infty}\mathrm{d}x\,Q_{m}(x)\int_{r}^{+\infty}\mathrm{d}\lambda\,\left(\frac{\lambda}{r_{0}}\right)^{B}\lambda^{-k-j+1}\int_{0}^{2(\lambda-r)}\mathrm{d}\rho\,\,\rho^{i+j}G(u-\rho)F\bigl{[}u-\rho-\lambda(x-1)\bigr{]}\,.$
(109)

Taking advantage of the conditions on $F$ and $G$, we can now restrict the
bounds of the integrals to $\lambda\leqslant\frac{u+\mathcal{T}}{x-1}$. This
yields

$\displaystyle\phi_{ij}^{B}$
$\displaystyle=\int_{1}^{+\infty}\mathrm{d}x\,Q_{m}(x)\int_{0}^{\frac{u+\mathcal{T}}{x-1}}\mathrm{d}\lambda\,\left(\frac{\lambda}{r_{0}}\right)^{B}\lambda^{-k-j+1}\int_{0}^{\min(u+\mathcal{T},\,2\lambda)}\mathrm{d}\rho\,\,\rho^{i+j}G(u-\rho)F\bigl{[}u-\rho-\lambda(x-1)\bigr{]}$
$\displaystyle-\int_{1}^{\frac{u+\mathcal{T}}{r}}\mathrm{d}x\,Q_{m}(x)\int_{r}^{\frac{u+\mathcal{T}}{x-1}}\mathrm{d}\lambda\,\left(\frac{\lambda}{r_{0}}\right)^{B}\lambda^{-k-j+1}\int_{0}^{\min(u+\mathcal{T},\,2(\lambda-r))}\mathrm{d}\rho\,\,\rho^{i+j}G(u-\rho)F\bigl{[}u-\rho-\lambda(x-1)\bigr{]}\,.$
(110)

We introduce

$\Xi_{k,i,j}(x)\equiv\int_{0}^{\frac{u+\mathcal{T}}{x-1}}\mathrm{d}\lambda\,\lambda^{-k-j+1}\int_{0}^{\min(u+\mathcal{T},\,2\lambda)}\mathrm{d}\rho\,\,\rho^{i+j}G(u-\rho)F\bigl{[}u-\rho-\lambda(x-1)\bigr{]}\,.$
(111)

Let us bound this quantity for any $x\in\ ]1,\infty]$:

$\displaystyle\big{|}\Xi_{k,i,j}(x)\big{|}$
$\displaystyle\leqslant\int_{0}^{\frac{u+\mathcal{T}}{x-1}}\mathrm{d}\lambda\,\lambda^{-k-j+1}\int_{0}^{\min(u+\mathcal{T},\,2\lambda)}\mathrm{d}\rho\,\,\rho^{i+j}\Big{|}G(u-\rho)F(u-\rho-\lambda(x-1))\Big{|}$
(112)
$\displaystyle\leqslant\mathcal{K}\int_{0}^{\frac{u+\mathcal{T}}{x-1}}\mathrm{d}\lambda\,\lambda^{-k-j+1}\int_{0}^{\min(u+\mathcal{T},\,2\lambda)}\mathrm{d}\rho\,\,\rho^{i+j}$
$\displaystyle=\frac{\mathcal{K}}{i+j+1}\int_{0}^{\frac{u+\mathcal{T}}{x-1}}\mathrm{d}\lambda\,\lambda^{-k-j+1}\Big{[}\min(u+\mathcal{T},\,2\lambda)\Big{]}^{i+j+1}$
$\displaystyle=\frac{\mathcal{K}}{i+j+1}\Bigg{\\{}\int_{0}^{\frac{u+\mathcal{T}}{2}}\mathrm{d}\lambda\,\lambda^{-k-j+1}(2\lambda)^{i+j+1}+\theta(3-x)\int_{\frac{u+\mathcal{T}}{2}}^{\frac{u+\mathcal{T}}{x-1}}\mathrm{d}\lambda\,\lambda^{-k-j+1}(u+\mathcal{T})^{i+j+1}\Bigg{\\}}$
$\displaystyle=\frac{\mathcal{K}}{i+j+1}\Bigg{\\{}\frac{2^{j+k-2}}{i+(3-k)}(u+\mathcal{T})^{i+(3-k)}+\theta(3-x)(u+\mathcal{T})^{i+j+1}\int_{\frac{u+\mathcal{T}}{2}}^{\frac{u+\mathcal{T}}{x-1}}\mathrm{d}\lambda\,\lambda^{-k-j+1}\Bigg{\\}}\,,$

where

$\mathcal{K}\equiv\sup_{(\tau,\tau^{\prime})\in[-\mathcal{T},u]^{2}}\Big{[}F(\tau)G(\tau^{\prime})\Big{]}\,.$
(113)

We now distinguish the case $j=2-k$, which yields

$\big{|}\Xi_{2,i,0}(x)\big{|}\leqslant\frac{\mathcal{K}}{i+1}\Bigg{\\{}\frac{1}{i+1}(u+\mathcal{T})^{i+1}+\theta(3-x)(u+\mathcal{T})^{i+1}\ln\left(\frac{2}{x-1}\right)\Bigg{\\}}\,,$
(114)

and the more general case $k+j\geqslant 2$, in which we find

$\big{|}\Xi_{k,i,j}(x)\big{|}\leqslant\frac{\mathcal{K}(u+\mathcal{T})^{i+(3-k)}}{i+j+1}\Bigg{\\{}\frac{2^{j+k-2}}{i+(3-k)}+\frac{\theta(3-x)}{j+(k-2)}\Big{(}2^{j+(k-2)}-(x-1)^{k+j-2}\Big{)}\Bigg{\\}}\,.$
(115)

Since $Q_{m}(x)=\mathcal{O}_{x\rightarrow
1^{+}}\left(\ln\left(x-1\right)\right)$ and
$Q_{m}(x)=\mathcal{O}_{x\to+\infty}\left(x^{-m-1}\right)$, we can now look at
the asymptotic behavior of our integrand. In the case $(k,j)=(2,0)$, we find

$Q_{m}(x)\,\Xi_{2,i,0}(x)=\begin{cases}\displaystyle{\mathcal{O}}_{x\rightarrow
1^{+}}\left(\ln(x-1)^{2}\right)\,,\\\\[8.5359pt]
\displaystyle{\mathcal{O}}_{x\to+\infty}\left(x^{-m-1}\right)\,.\end{cases}$
(116)

In the general case $k+j\geqslant 2$, we instead have

$Q_{m}(x)\,\Xi_{k,i,j}(x)=\begin{cases}\displaystyle{\mathcal{O}}_{x\rightarrow
1^{+}}\left(\ln(x-1)(x-1)^{j+(k-2)}\right)\,,\\\\[8.5359pt]
\displaystyle{\mathcal{O}}_{x\to+\infty}\left(x^{-m-1}\right)\,.\end{cases}$
(117)

It is now clear that $Q_{m}(x)\,\Xi_{k,i,j}(x)$ is integrable at both bounds
$x\rightarrow 1^{+}$ and $x\to+\infty$ as soon as $k\in\\{1,2\\}$. After doing
a similar (and much easier) analysis on the second member of Eq. (B.0.2), we
find that Eqs. (B.0.2)–(B.0.2), and consequently Eq. (44b), have a convergent
limit as $B\rightarrow 0$, _i.e._ , they do not develop poles in $1/B$.
Repeating the analysis of this section with an extra factor $\ln(r/r_{0})$ in
the integrand shows that we can also safely compute the $1/B$ coefficient in
the Laurent series.

Table 1: Coefficients $\mathcal{A}_{m,\ell}^{n}$, $\mathcal{B}_{m,\ell}^{n}$,
$\mathcal{C}_{m,\ell}^{n}$, $\mathcal{D}_{m,\ell}^{n}$ entering the radiative
quadrupole

(a) Coefficients $\mathcal{A}_{m,\ell}^{n}$ associated to ${}_{1,m}\overline{\Psi}_{\ell}[\overline{\mathrm{M}}_{a\langle i}^{(n)},\overline{\mathrm{M}}_{j\rangle a}^{(8-n)}]$ $\ell$ | 2 | 3 | 4
---|---|---|---
$m$ $n$ | 0 | 1 | 2 | 3 | 4 | 0 | 1 | 2 | 3 | 4 | 0 | 1 | 2 | 3
0 | $-\frac{136}{735}$ | $-\frac{13808}{2625}$ | $-\frac{34576}{2625}$ | $\frac{120}{7}$ | $\frac{16}{7}$ | $0$ | $\frac{832}{535}$ | $\frac{29432}{1575}$ | $\frac{88}{9}$ | -$\frac{32}{7}$ | $-\frac{928}{25725}$ | $-\frac{1168}{3675}$ | $\frac{664}{525}$ | $\frac{64}{21}$
1 | $0$ | $-\frac{3344}{6125}$ | $\frac{6184}{875}$ | $-\frac{16456}{875}$ | $-\frac{32}{7}$ | $\frac{1216}{1225}$ | $\frac{3464}{525}$ | $-\frac{11584}{525}$ | $-\frac{17032}{525}$ | $\frac{64}{7}$ | $0$ | $-\frac{9256}{8575}$ | $-\frac{6256}{1225}$ | $\frac{3776}{1225}$
2 | $\frac{8}{75}$ | $\frac{18328}{3675}$ | $\frac{22448}{3675}$ | $\frac{2528}{175}$ | $\frac{16}{7}$ | $0$ | $-\frac{3776}{735}$ | $\frac{36272}{2205}$ | $\frac{304}{45}$ | $-\frac{32}{7}$ | $-\frac{16}{343}$ | $-\frac{2896}{5145}$ | $\frac{1024}{735}$ | $\frac{4696}{735}$
3 | $0$ | $-\frac{568}{875}$ | $-\frac{1224}{875}$ | $\frac{816}{875}$ | $0$ | $-\frac{128}{175}$ | $\frac{536}{525}$ | $-\frac{5176}{525}$ | $-\frac{288}{25}$ | $0$ | $0$ | $\frac{7064}{3675}$ | $\frac{4616}{1225}$ | $\frac{1464}{1225}$
4 | $\frac{96}{275}$ | $\frac{13152}{6125}$ | $\frac{8544}{6125}$ | $0$ | $0$ | $0$ | $-\frac{4608}{1225}$ | $-\frac{3936}{1125}$ | $0$ | $0$ | $\frac{10336}{94325}$ | $\frac{4192}{8575}$ | $-\frac{1616}{1225}$ | $0$
5 | $0$ | $-\frac{832}{1225}$ | $0$ | $0$ | $0$ | $-\frac{64}{245}$ | $-\frac{32}{105}$ | $0$ | $0$ | $0$ | $0$ | $-\frac{464}{1029}$ | $0$ | $0$
6 | $-\frac{3648}{13475}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $-\frac{304}{11319}$ | $0$ | $0$ | $0$

(b) Coefficients $\mathcal{B}_{m,\ell}^{n}$ associated to $\,{}_{2,m}\overline{\Psi}_{\ell}[\overline{\mathrm{M}}_{a\langle i}^{(n)},\overline{\mathrm{M}}_{j\rangle a}^{(7-n)}]$ $\ell$ | 2 | 3 | 4
---|---|---|---
$m$ $n$ | 0 | 1 | 2 | 3 | 0 | 1 | 2 | 3 | 0 | 1 | 2 | 3
0 | $0$ | $\frac{696}{875}$ | $\frac{312}{7}$ | $\frac{48}{7}$ | $\frac{608}{175}$ | $\frac{10768}{525}$ | $-\frac{304}{21}$ | $-\frac{192}{7}$ | $0$ | $\frac{272}{147}$ | $\frac{160}{21}$ | $0$
1 | $-\frac{15336}{6125}$ | $-\frac{5632}{875}$ | $-\frac{18728}{875}$ | $\frac{64}{7}$ | $0$ | $-\frac{5952}{175}$ | $-\frac{5744}{175}$ | $\frac{176}{7}$ | $-\frac{2432}{1715}$ | $-\frac{1616}{245}$ | $\frac{1040}{49}$ | $\frac{48}{7}$