# Gravitational Wave Flux and Quadrupole Modes from Quasi-Circular Non-Spinning Compact Binaries to the Fourth Post-Newtonian Order

Luc Blanchet, Guillaume Faye, Quentin Henry, François Larrouturou, David Trestini

► To cite this version:

Luc Blanchet, Guillaume Faye, Quentin Henry, François Larrouturou, David Trestini. Gravitational Wave Flux and Quadrupole Modes from Quasi-Circular Non-Spinning Compact Binaries to the Fourth Post-Newtonian Order. Physical Review D, 2023, 108 (6), pp.064041. 10.1103/Phys-RevD.108.064041. hal-04093506

HAL Id: hal-04093506

https://hal.science/hal-04093506v1

Submitted on 28 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Gravitational-wave flux and quadrupole modes from quasicircular nonspinning compact binaries to the fourth post-Newtonian order

Luc Blanchet<sup>1,∗</sup> Guillaume Faye<sup>1,2,†</sup> Quentin Henry<sup>1,3,‡</sup> François Larrouturou<sup>1,4,§</sup> and David Trestini<sup>1,5,||</sup>

<sup>1</sup>*ARCO*, Institut d'Astrophysique de Paris, UMR 7095, CNRS, Sorbonne Université,
98<sup>bis</sup> boulevard Arago, 75014 Paris, France

<sup>2</sup>Centre for Strings, Gravitation and Cosmology, Department of Physics,
Indian Institute of Technology Madras, Chennai 600036, India

<sup>3</sup>Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-14476 Potsdam, Germany

<sup>4</sup>Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany

<sup>5</sup>Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, Université Paris Cité,
CNRS, F-92190 Meudon, France

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

(Received 16 May 2023; accepted 29 June 2023; published 21 September 2023)

This article provides the details on the technical derivation of the gravitational waveform and total gravitational-wave energy flux of nonspinning compact binary systems to the fourth post-Newtonian (4PN) order beyond the Einstein quadrupole formula. In particular: (i) we overview the link between the radiative multipole moments measured at infinity and the source moments in the framework of dimensional regularization; (ii) we compute special corrections to the source moments due to “infrared” commutators arising at the 4PN order; (iii) we derive a “postadiabatic” correction needed to evaluate the tail integral with 2.5PN relative precision; (iv) we discuss the relation between the binary’s orbital frequency in quasicircular orbit and the gravitational-wave frequency measured at infinity; (v) we compute the hereditary effects at the 4PN order, including those coming from the recently derived tails-of-memory; and (vi) we describe the various tests we have performed to ensure the correctness of the results. Those results are collected in an ancillary file.

DOI: 10.1103/PhysRevD.108.064041

I. INTRODUCTION

The post-Newtonian (PN) approximation is a paramount
tool for computing the generation of gravitational waves
(GWs) by isolated sources. This perturbation technique
relies on two joint expansions, valid for systems that have
slow velocity and are self-gravitating, and which are thus
linked by weak gravitational fields. In consequence, it is
ideally suited to study the inspiral phase of compact binary
systems, made of black holes (BHs) or neutron stars (NSs).
While the derivation of lowest-order PN results, starting
from the famous Einstein quadrupole formula [1–3], is
fairly straightforward, the systematic expansion to high

PN orders entails a number of subtle issues. One is the
control of physical nonlinear effects in the propagation of
GWs from source to detector. Other more technical
difficulties are linked to the appearance of divergent
integrals in the calculations, and the subsequent need to
implement a proper regularization scheme.

There are two modern ways of carrying out the PN approximation. The first one is to embed it into the more general PN-MPM framework [4–8], consisting of a multipolar and post-Minkowskian (MPM) expansion for the field exterior to the isolated source, subsequently matched to the inner PN field of the source. This method achieved the completion of the gravitational waveform to 3.5PN order for nonspinning objects [9–18], i.e., up to the (v/c)⁷ correction beyond the Einstein quadrupole formula. A variant of the PN-MPM framework is the DIRE formalism, developed to 2PN order [19]. The second way relies on an effective field theory (EFT) [20,21], extracting the radiation sector from a derivative expansion at the level of the action [22]. This more recent framework has reached the 2PN precision [23]. We refer the interested reader to the reviews [24–27] for more details on the PN scheme.

The main output of the PN calculation is the observable
waveform (phase and amplitude modes) as an expansion

Published by the American Physical Society under the terms of
the Creative Commons Attribution 4.0 International license.
Further distribution of this work must maintain attribution to
the author(s) and the published article's title, journal citation,
and DOI. Open access publication funded by the Max Planck
Society.

<footnote>*luc.blanchet@iap.fr</footnote>

<footnote>†faye@iap.fr</footnote>

<footnote>‡quentin.henry@aei.mpg.de</footnote>

<footnote>§francois.larrouturou@desy.de</footnote>

<footnote>||david.trestini@obspm.fr</footnote>
series in terms of the PN parameter $x = (Gm\Omega/c^3)^{2/3}$, where $\Omega$ is the half of the GW frequency. This is obviously essential for the data analysis of GW detectors [28,29], notably for the current LIGO-Virgo-Kagra network [30–32], but also for future generations of detectors, such as the Laser Interferometer Space Antenna (LISA) [33] or the Einstein Telescope (ET) [34]. In addition, the PN results are of great interest for numerical relativity, as they are naturally used for comparisons with numerical results (see e.g., [35–38]), as well as for the calibration of initial data (see e.g., [39–42]). They also play a crucial role for comparisons with gravitational self-force (GSF) results, notably second-order ones [43–45].

On the other hand, the PN framework is essential when it comes to experimentally testing the theory of gravitation. By drawing precise predictions and confronting them against observations, it allows performing “agnostic” tests of general relativity (GR) in the Solar system, for binary pulsars and even for transiting exoplanets [46]. With GWs, the idea is to treat each coefficient in the PN expansion of the GW phase as an independent parameter, and fit it against GW data [47,48]. This method has been applied to data collected from LIGO-Virgo detectors, and has already permitted probing the nonlinear structure of GR by confirming the signature of GW tails to the 1.5PN order [49–51]. This *modus operandi* is very promising, notably regarding the possibility of multiband detections, combining LISA and detectors on ground [52]. In the same spirit, the prediction for the tidal effects affecting the binary NSs [53–60], which first contribute to the waveform at the 5PN order (in the spinless case), have been used at leading order to put a constraint on the equation of state of matter composing neutron stars [61], and could be used to distinguish usual BHs from “exotic” compact objects with ET. In addition to the agnostic tests, the PN framework is also used to perform “educated” tests, by computing observable quantities in specific alternative theories. It has been developed, for instance, in the context of massless scalar field theories [62–65], and used to make predictions concerning the gravitational radiation [66–68]. Similar studies have been conducted for various alternative theories [69–73].

This paper presents the last step of the completion of the gravitational waveform for nonspinning, compact binary systems evolving with 4PN precision, in the case of quasicircular orbits (and in GR). The first step was the computation of the equations of motion to 4PN order, from which one deduces in particular the binary’s gravitational binding energy. The equations of motion were obtained by several groups using different methods: (i) the Arnowitt-Deser-Misner (ADM) Hamiltonian formalism [74–77] led to complete results except for one “ambiguity” parameter, which was fixed by resorting to a comparison with gravitational self-force results [78]; (ii) the Fokker

Lagrangian formalism in harmonic coordinates [79–82] yielded for the first time a complete result without ambiguity parameter; and (iii) the effective field theory (EFT) approach [83–90] recovered the complete and unambiguous result. The second step was the computation and proper regularization, using the PN-MPM formalism, of the source mass quadrupole moment to 4PN [91–93], the current quadrupole to 3PN [94] and the mass octupole to the same order [95]. The next step, relying on previous works as well, was the control of the various nonlinear interactions between moments occuring through the 4PN order [16,17,95–99]. The prominent of these interactions at 4PN order is the so-called “tail-of-memory”, which is a cubically nonlinear effect [99]. Moreover, we were able to compute the flux for circular orbits up to the 4.5PN order, since only one nonlinear interaction plays a role for this computation, which is the so-called “tail-of-tail-of-tail”, a quartic effect [96] (see also [100]). In the present work, we complete the derivation of the gravitational flux to 4PN for generic orbits, and extract the $(\ell, m) = (2, 2)$ and $(\ell, m) = (2, 0)$ modes from the radiative mass quadrupole moment. With those results at hand, and the 4.5PN term of the flux for quasicircular orbits already computed in [96], the gravitational phasing at the 4.5PN precision is derived and reported in the companion Letter [101].

This work is organized as follows. In Sec. II, we explicitly present the relations between the radiative, canonical and source multipole moments, at the required 4PN accuracy (and 4.5PN for circular orbits). Notably, we derive the missing relation between the canonical and source quadrupole moment, extending [97] to arbitrary dimensions as required by dimensional regularization. Section III is devoted to the calculation of novel terms coined as “infrared (IR) commutators” appearing at the 4PN order in the metric and the source mass quadrupole moment. They are zero for circular orbits and do not contribute to the results reported in [101]. Then, Sec. IV lists the expressions of the required source moments on circular orbits, obtained in previous works, notably [91–95]. Section V exposes the technical toolbox that we used to perform the derivations and integrations needed to get the final results. Notably, an important physical contribution at 4PN order comes from the 1.5PN tail term computed to relative 2.5PN order, which necessitates a “postadiabatic” evaluation, described in Sec. V B. The results are presented in Sec. VI, and consist of: (i) the flux at 4PN for generic orbits, in the center-of-mass (c.m.) frame; (ii) the flux at 4.5PN for quasicircular orbits; and (iii) the $(\ell, m) = (2, 2)$ and $(\ell, m) = (2, 0)$ modes, the latter being purely due to the nonlinear memory interactions. The results are all collected in the ancillary file [102]. Appendix A contains the test of the boosted Schwarzschild solution (BSS) [103], successfully applied to the mass quadrupole moment, and Appendix B displays the contributions of IR commutators in the 4PN metric.
## II. NONLINEAR EFFECTS IN THE GW PROPAGATION

In this work,<fnref n="1" /> we compute the GW energy flux of compact binaries on generic orbits and the dominant GW mode (ℓ, m) = (2, 2) to the 4PN order, as well as the energy flux on quasicircular orbits to the 4.5PN order. Let us recall that the general multipole decomposition of the flux reads [104]

$$ \mathcal{F} = \sum_{\ell \ge 2} \frac{G}{c^{2\ell+1}} \left[ a_{\ell} U_{L}^{(1)} U_{L}^{(1)} + \frac{b_{\ell}}{c^2} V_{L}^{(1)} V_{L}^{(1)} \right], \quad (2.1) $$

where $U_L$ and $V_L$ (for $ℓ ≥ 2$) respectively denote the mass and current radiative multipole moments [which are symmetric trace-free (STF) in their indices], and the numerical coefficients are given by

$$ \begin{aligned} a_{\ell} &= \frac{(\ell+1)(\ell+2)}{(\ell-1)\ell\ell!(2\ell+1)!!} \quad \text{and} \\ b_{\ell} &= \frac{4\ell(\ell+2)}{(\ell-1)(\ell+1)!(2\ell+1)!!}. \end{aligned} \quad (2.2) $$

The asymptotic waveform, to order $1/R$ in the distance to the source, when expressed in a (Bondi type) radiative coordinate system $(T, X^i)$, with $X \equiv RN$, and written as a perturbation to the ordinary metric $g_{ab}$, reads [104]

$$ h_{ab}^{\text{TT}} = \frac{4G}{c^2 R} \perp_{abij}(N) \sum_{\ell \ge 2} \frac{1}{c^{\ell}\ell!} \left[ N_{L-2} U_{ijL-2}(u) - \frac{2\ell}{c(\ell+1)} N_{kL-2} \epsilon_{kl}(i V_j)_{L-2}(u) \right] + \mathcal{O}\left(\frac{1}{R^2}\right), \quad (2.3) $$

where $\perp_{abij}(N)$ is the usual transverse and traceless (TT) projector and the multipole moments are defined at the retarded time $u \equiv T - R/c$. The GW amplitude (ℓ, m) modes can be read off the radiative moments $U_L$ and $V_L$, and for instance the dominant mode (2, 2) can be computed directly from the mass quadrupole moment $U_{ij}$ (in the case of planar, nonprecessing, orbits; see [16,17,105]).

Therefore, the knowledge of $U_L$ and $V_L$ at the relevant PN order will lead to the desired results.

To describe the nonlinear effects in the GW propagation from the source to the observer, we connect the radiative moments to the so-called canonical moments $M_L$ and $S_L$. In turn, the canonical moments are expressed in terms of source moments $I_L$ and $J_L$, as well as “gauge” moments, $W_L$, $X_L$, $Y_L$, and $Z_L$. The source and gauge moments describe in our formalism the physics of the PN source. This section is thus devoted to the relations between the radiative moments and the source ones, at the required PN order. We refer to [8] for a more detailed review and physical discussions about those constructions.

### A. Radiative moments versus canonical moments

#### 1. Radiative moments entering the flux at 4PN order

In order to derive the energy flux (2.1) at 4PN order beyond the leading quadrupolar order, the obvious first input is the radiative quadrupole moment itself, $U_{ij}$, to 4PN order. Recalling that, at leading order, the radiative moments $U_L$ and $V_L$ reduce to the ℓth time derivatives (with respect to the retarded time $u$ of the radiative coordinates) of the canonical moments $M_L$ and $S_L$, we straightforwardly write

$$ U_{ij} = M_{ij}^{(2)} + U_{ij}^{1.5PN} + U_{ij}^{2.5PN} + U_{ij}^{3PN} \\ + U_{ij}^{3.5PN} + U_{ij}^{4PN} + \mathcal{O}\left(\frac{1}{c^9}\right), \quad (2.4) $$

with small PN corrections up to 4PN, as indicated. The leading correction at 1.5PN is due to the quadratic interaction between the static ADM mass M and the mass quadrupole moment $M_{ij}$, denoted “M × M_{ij}” and called the “tail”. It is well-known and reads explicitly [7]

$$ U_{ij}^{1.5PN} = \frac{2GM}{c^3} \int_{0}^{+\infty} d\tau M_{ij}^{(4)}(u-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \frac{11}{12} \right], \quad (2.5) $$

where the length scale $b_0$ is an arbitrary constant linked with the choice of the origin of time of the asymptotic radiative coordinates, with respect to the harmonic coordinates covering the source’s near zone. At the next 2.5PN order, there arises the nonlocal interaction involving two quadrupole moments $M_{ij} \times M_{kl}$, called the (displacement) “memory”, together with associated instantaneous quadrupole-quadrupole terms [106], and an instantaneous interaction $M_{ij} \times S_k$ sometimes dubbed “failed tail” [107],

$$ U_{ij}^{2.5PN} = \frac{G}{c^5} \left\{ -\frac{2}{7} \int_{0}^{+\infty} d\tau \left[ M_{a(i}^{(3)} M_{j)a}^{(3)} \right] (u-\tau) + \frac{1}{7} M_{a\langle i}^{(5)} M_{j\rangle a} \right. \\ \left. -\frac{5}{7} M_{a\langle i}^{(4)} M_{j)a}^{(1)} - \frac{2}{7} M_{a\langle i}^{(3)} M_{j)a}^{(2)} + \frac{1}{3} \epsilon_{ab\langle i} M_{j)a}^{(4)} S_b \right\}. \quad (2.6) $$

<footnote n="1">¹The conventions employed throughout are as follows: we work with a mostly plus signature; Greek letters denote spacetime indices and Latin letters denote purely spatial indices; we use the multi-index notations, with symmetric trace-free (STF) mass and current moments $M_L \equiv M_{i_1 i_2 \dots i_\ell}$ and $S_L \equiv S_{i_1 i_2 \dots i_\ell}$; a superscript (n) indicates n time derivatives; hats and angular brackets denote a symmetric trace-free projection, for instance $\hat{x}_L = x_{\langle L\rangle} = \text{STF}[x_L]$ with $x_L = x_{i_1} \cdots x_{i_\ell}$; the d'Alembertian operator is defined with respect to the flat Minkowski metric $\Box = \eta^{\mu\nu}\partial_\mu\partial_\nu = \Delta - c^{-2}\partial_t^2$; symmetrization and antisymmetrization are weighted, e.g., $A_{(ij)} = (A_{ij} + A_{ji})/2$ and $A_{[ij]} = (A_{ij} - A_{ji})/2$; the usual Levi-Civita symbol is denoted $e_{ijk}$ with $e_{123} = 1$; finally, as usual, we dubb "nPN" a quantity of order $\mathcal{O}(c^{-2n})$.</footnote>
At the 3PN order the first cubic interaction appears, consisting of the interplay between two ADM masses and the quadrupole moment, i.e., $M \times M \times M_{ij}$, rightly called “tail-of-tail” and given by [10,95]

$$U_{ij}^{3\text{PN}} = \frac{2G^2M^2}{c^6} \int_0^{+\infty} d\tau M_{ij}^{(5)}(u-\tau) \left[ \ln^2\left(\frac{c\tau}{2b_0}\right) + \frac{11}{6}\ln\left(\frac{c\tau}{2b_0}\right) - \frac{107}{105}\ln\left(\frac{c\tau}{2r_0}\right) + \frac{124627}{44100} \right]. \quad (2.7)$$

Note the appearance here of a new constant length scale $r_0$, which should be distinguished from the

previously introduced scale $b_0$. The scale $r_0$ is a running scale in the sense of renormalization group theory [21,108,109]. The coefficient in front of the $\ln r_0$ term in (2.7) is exactly the $\beta$-function coefficient associated to the renormalization of the mass quadrupole moment, say $\beta_2 = -\frac{214}{105}$. In the present formalism, this scale originates from the Hadamard regularization scheme. The next-order 3.5PN term has a structure similar to the 2.5PN one, i.e., with some memory type integrals and instantaneous terms. The interactions between moments, still of quadratic nature, are however more complicated,

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

At the 4PN order appears the “tail-of-memory”, which is a cubic interaction involving two time-varying quadrupole moments and the mass, i.e., $M \times M_{ij} \times M_{kl}$. It also contains a double integration over the two quadrupole moments; see the first line of (2.9). The tail-of-memory comes with many tail-looking interactions with only one integration over the quadrupole moments. Similarly to the fact that the memory came with the failed tail [see (2.6)], the tail-of-memory comes along with a cubic interaction involving the constant angular momentum, i.e., $M \times S_i \times M_{jk}$, see the last line of (2.9). We have [99]

$$U_{ij}^{4\text{PN}} = \frac{8G^2M}{7c^8} \left\{ \begin{aligned} & \int_0^{+\infty} d\rho M_{a\langle i}^{(4)}(u-\rho) \int_0^{+\infty} d\tau M_{j)a}^{(4)}(u-\rho-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) - \frac{1613}{270} \right] \\ & - \frac{5}{2} \int_0^{+\infty} d\tau \left[ M_{a\langle i}^{(3)} M_{j)a}^{(4)} \right] (u-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{3}{2} \ln\left(\frac{c\tau}{2b_0}\right) \right] - 3 \int_0^{+\infty} d\tau \left[ M_{a\langle i}^{(2)} M_{j)a}^{(5)} \right] (u-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{11}{12} \ln\left(\frac{c\tau}{2b_0}\right) \right] \\ & - \frac{5}{2} \int_0^{+\infty} d\tau \left[ M_{a\langle i}^{(1)} M_{j)a}^{(6)} \right] (u-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{3}{10} \ln\left(\frac{c\tau}{2b_0}\right) \right] - \int_0^{+\infty} d\tau \left[ M_{a\langle i} M_{j)a}^{(7)} \right] (u-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) - \frac{1}{4} \ln\left(\frac{c\tau}{2b_0}\right) \right] \\ & - 2 M_{a\langle i}^{(2)} \int_0^{+\infty} d\tau M_{j)a}^{(5)}(u-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{27521}{5040} \right] - \frac{5}{2} M_{a\langle i}^{(1)} \int_0^{+\infty} d\tau M_{j)a}^{(6)}(u-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{15511}{3150} \right] \\ & + \frac{1}{2} M_{a\langle i} \int_0^{+\infty} d\tau M_{j)a}^{(7)}(u-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) - \frac{6113}{756} \right] \\ & - \frac{2G^2M}{3c^8} S_a \epsilon_{ab\langle i} \int_0^{+\infty} d\tau M_{j\rangle b}^{(6)}(u-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + 2\ln\left(\frac{c\tau}{2r_0}\right) + \frac{1223}{1890} \right] \end{aligned} \right\}. \quad (2.9)$$
Note that the “genuine” tail-of-memory given by the first term of (2.9) can be retrieved by replacing the canonical moment $M_{ij}$ in the expression of the memory term of (2.6) by the radiative moment itself, taking into account the dominant tail effect, i.e.,

$$M_{ij}^{(2)} \rightarrow U_{ij} = M_{ij}^{(2)} + U_{ij}^{1.5PN}, \quad (2.10)$$

along with a reexpansion at cubic order and an integration by parts. This agrees with the general prediction for memory effects computed using the radiative multipole moments defined at future null infinity, see [110,111].

Besides the mass quadrupole moment, the expression of the flux (2.1) at the 4PN order requires the knowledge of moments of higher multipolarity, but evaluated at lower PN orders. At the next multipolar order, we thus need the mass

octupole $U_{ijk}$ and current quadrupole $V_{ij}$ moments with a 3PN precision,

$$U_{ijk} = M_{ijk}^{(3)} + U_{ijk}^{1.5PN} + U_{ijk}^{2.5PN} + U_{ijk}^{3PN} + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (2.11a)$$

$$V_{ij} = S_{ij}^{(2)} + V_{ij}^{1.5PN} + V_{ij}^{2.5PN} + V_{ij}^{3PN} + \mathcal{O}\left(\frac{1}{c^7}\right). \quad (2.11b)$$

Reporting the results of [95], the above contributions read

$$U_{ijk}^{1.5PN} = \frac{2GM}{c^3} \int_{0}^{+\infty} d\tau M_{ijk}^{(5)}(u-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \frac{97}{60} \right], \quad (2.12a)$$

$$U_{ijk}^{2.5PN} = \frac{G}{c^5} \left\{ \begin{aligned} & \int_{0}^{+\infty} d\tau \left[ -\frac{1}{3} M_{a\langle i}^{(3)} M_{jk\rangle a}^{(4)} - \frac{4}{5} \epsilon_{ab\langle i} M_{j\underline{a}}^{(3)} S_{k\rangle b}^{(3)} \right] (u-\tau) + \frac{1}{4} M_{a\langle i} M_{jk\rangle a}^{(6)} + \frac{1}{4} M_{a\langle i}^{(1)} M_{jk\rangle a}^{(5)} + \frac{1}{4} M_{a\langle i}^{(2)} M_{jk\rangle a}^{(4)} - \frac{4}{3} M_{a\langle i}^{(3)} M_{jk\rangle a}^{(3)} \\ & - \frac{9}{4} M_{a\langle i}^{(4)} M_{jk\rangle a}^{(2)} - \frac{3}{4} M_{a\langle i}^{(5)} M_{jk\rangle a}^{(1)} + \frac{1}{12} M_{a\langle i}^{(6)} M_{jk\rangle a} + \frac{12}{5} S_{i\langle} S_{jk\rangle}^{(4)} + \epsilon_{ab\langle i} \left[ \frac{9}{5} M_{j\underline{a}} S_{k\rangle b}^{(5)} + \frac{27}{5} M_{j\underline{a}}^{(1)} S_{k\rangle b}^{(4)} + \frac{8}{5} M_{j\underline{a}}^{(2)} S_{k\rangle b}^{(3)} \right. \\ & \left. + \frac{12}{5} M_{j\underline{a}}^{(3)} S_{k\rangle b}^{(2)} + \frac{3}{5} M_{j\underline{a}}^{(4)} S_{k\rangle b}^{(1)} + \frac{1}{5} M_{j\underline{a}}^{(5)} S_{k\rangle b} + \frac{9}{20} M_{jk\rangle a}^{(5)} S_b \right] \end{aligned} \right\}, \quad (2.12b)$$

$$U_{ijk}^{3PN} = \frac{2G^2 M^2}{c^6} \int_{0}^{+\infty} d\tau M_{ijk}^{(6)}(u-\tau) \left[ \ln^2\left(\frac{c\tau}{2b_0}\right) + \frac{97}{30}\ln\left(\frac{c\tau}{2b_0}\right) - \frac{13}{21}\ln\left(\frac{c\tau}{2r_0}\right) + \frac{13283}{8820} \right], \quad (2.12c)$$

where the underlined indices within angled brackets are excluded from the STF operation, and

$$V_{ij}^{1.5PN} = \frac{2GM}{c^3} \int_{0}^{+\infty} d\tau S_{ij}^{(4)}(u-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \frac{7}{6} \right], \quad (2.12d)$$

$$V_{ij}^{2.5PN} = \frac{G}{c^5} \left\{ -\frac{3}{7} M_{a\langle i} S_{j\rangle a}^{(5)} - \frac{3}{7} M_{a\langle i}^{(1)} S_{j\rangle a}^{(4)} + \frac{8}{7} M_{a\langle i}^{(2)} S_{j\rangle a}^{(3)} + \frac{4}{7} M_{a\langle i}^{(3)} S_{j\rangle a}^{(2)} + \frac{17}{7} M_{a\langle i}^{(4)} S_{j\rangle a}^{(1)} + \frac{9}{7} M_{a\langle i}^{(5)} S_{j\rangle a} - \frac{1}{28} M_{ija}^{(5)} S_a \\ + \epsilon_{ab\langle i} \left[ -\frac{15}{56} M_{j\rangle ac} M_{bc}^{(6)} - \frac{113}{112} M_{j\rangle ac}^{(1)} M_{bc}^{(5)} - \frac{353}{336} M_{j\rangle ac}^{(2)} M_{bc}^{(4)} - \frac{3}{14} M_{j\rangle ac}^{(3)} M_{bc}^{(3)} \right. \\ \left. + \frac{5}{168} M_{j\rangle ac}^{(4)} M_{bc}^{(2)} + \frac{3}{112} M_{j\rangle ac}^{(5)} M_{bc}^{(1)} - \frac{3}{112} M_{j\rangle ac}^{(6)} M_{bc} + S_{j\rangle a}^{(4)} S_b \right] \right\}, \quad (2.12e)$$

$$V_{ij}^{3PN} = \frac{2G^2 M^2}{c^6} \int_{0}^{+\infty} d\tau S_{ij}^{(5)}(u-\tau) \left[ \ln^2\left(\frac{c\tau}{2b_0}\right) + \frac{7}{3}\ln\left(\frac{c\tau}{2b_0}\right) - \frac{107}{105}\ln\left(\frac{c\tau}{2r_0}\right) - \frac{13127}{11025} \right]. \quad (2.12f)$$

Higher-order multipole moments have no cubic contributions at the required PN order. At 2PN order, we need the mass hexadecapole $U_{ijkl}$ as well as the current octupole $V_{ijk}$, which read

$$U_{ijkl} = M_{ijkl}^{(4)} + U_{ijkl}^{1.5PN} + \mathcal{O}\left(\frac{1}{c^5}\right), \quad (2.13a)$$

$$V_{ijk} = S_{ijk}^{(3)} + V_{ijk}^{1.5PN} + \mathcal{O}\left(\frac{1}{c^5}\right), \quad (2.13b)$$
where [95]

$$ U_{ijkl}^{1.5PN} = \frac{G}{c^3} \left\{ 2M \int_{0}^{+\infty} d\tau M_{ijkl}^{(6)}(u-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \frac{59}{30} \right] \right. \\ \left. + \frac{2}{5} \int_{0}^{+\infty} d\tau \left[ M_{\langle ij \rangle kl}^{(3)} M_{kl}^{(3)} \right] (u-\tau) \right. \\ \left. - \frac{21}{5} M_{\langle ij \rangle kl} M_{kl}^{(5)} - \frac{63}{5} M_{\langle ij \rangle kl}^{(1)} M_{kl}^{(4)} - \frac{102}{5} M_{\langle ij \rangle kl}^{(2)} M_{kl}^{(3)} \right\}, \quad (2.14a) $$

$$ V_{ijk}^{1.5PN} = \frac{G}{c^3} \left\{ 2M \int_{0}^{+\infty} d\tau S_{ijk}^{(5)}(u-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \frac{5}{3} \right] \right. \\ \left. - 2M_{\langle ij \rangle k}^{(4)} S_{k} \right. \\ \left. - \frac{1}{10} \epsilon_{ab\langle i} M_{ja} M_{k\rangle b}^{(5)} + \frac{1}{2} \epsilon_{ab\langle i} M_{ja}^{(1)} M_{k\rangle b}^{(4)} \right\}. \quad (2.14b) $$

Note the appearance of a memory integral at 1.5PN order in $U_{ijkl}$, in addition to the usual tail one. Finally, we need the moments $U_{ijklm}$ and $V_{ijkl}$ at 1PN, as well as $U_{ijklmn}$ and $V_{ijklm}$ at Newtonian order, which are trivially given by

$$ U_{ijklm} = M_{ijklm}^{(5)} + \mathcal{O}\left(\frac{1}{c^3}\right), \quad V_{ijkl} = S_{ijkl}^{(4)} + \mathcal{O}\left(\frac{1}{c^3}\right), \quad (2.15a) $$

$$ U_{ijklmn} = M_{ijklmn}^{(6)} + \mathcal{O}\left(\frac{1}{c^3}\right), \quad V_{ijklm} = S_{ijklm}^{(5)} + \mathcal{O}\left(\frac{1}{c^3}\right). \quad (2.15b) $$

## 2. Radiative moments entering the quasicircular flux at the 4.5PN order

For quasicircular orbits, the 4.5PN term in the flux was obtained in [96]. One could naively think that such a computation would require the complete knowledge of the relations between radiative and canonical moments, as presented above, but pushed one half PN order further. This was actually not the case, since for circular orbits, only a limited control of the relation between the radiative and canonical mass quadrupole moments was necessary. This is discussed in [96], and we only remind the key points. First, it is well-known that the contributions of instantaneous interactions entering the flux at half-integer PN order (e.g., 4.5PN order) vanish for quasicircular orbits. So only nonlocal contributions such as tails can potentially contribute. Second, the quadratic nonlocal memory interaction that enters the radiative moments [see Eqs. (2.6) and (2.8)] become instantaneous in the flux by virtue of time differentiation, so these will not contribute either. Last, the tails-of-memory and spin-quadrupole tails, which both enter at 4PN, will next contribute at 5PN but not at 4.5PN. This allows to use dimensional arguments to determine the interactions that can contribute to the 4.5PN quasicircular flux. For the mass quadrupole, only the quartic $M^3 \times M_{ij}$, naturally dubbed “tails-of-tails-of-tails”, can play a role. It has been computed in [96], and reads

$$ U_{ij}^{4.5PN}|_{ToToT} = \frac{G^3 M^3}{c^9} \int_0^{+\infty} d\tau M_{ij}^{(6)}(u-\tau) \left[ \frac{4}{3} \ln^3\left(\frac{c\tau}{2b_0}\right) + \frac{11}{3} \ln^2\left(\frac{c\tau}{2b_0}\right) - \frac{428}{105} \ln\left(\frac{c\tau}{2b_0}\right) \ln\left(\frac{c\tau}{2r_0}\right) \right. \\ \left. + \frac{124627}{11025} \ln\left(\frac{c\tau}{2b_0}\right) - \frac{1177}{315} \ln\left(\frac{c\tau}{2r_0}\right) + \frac{129268}{33075} + \frac{428}{315} \pi^2 \right]. \quad (2.16) $$

The full $U_{ij}^{4.5PN}$ also contains quadratic memory interactions, like those entering at 2.5PN and 3.5PN; see (2.6) and (2.8). If, as explained above, those are not needed to compute the flux (and phase) for quasicircular orbits, they will enter the expression of the (2, 2) mode. As they are yet undetermined, they restrict the accuracy we can reach when deriving the (2, 2) mode, which is why it is presented in Sec. VI C at 4PN and not up to 4.5PN order. In addition, radiation reaction effects at 4.5PN in the mass quadrupole should contribute and are also not under control.

Regarding other moments, the 3.5PN terms of the mass octupole $U_{ijk}$ and current quadrupole $V_{ij}$, as well as the 2.5PN terms of the mass hexadecapole $U_{ijkl}$ and current octupole $V_{ijk}$, are composed of quadratic memory integrals, which cannot contribute in the 4.5PN quasicircular flux. The only moments playing a role beyond the mass quadrupole are [95,96]

$$ U_{ijklm}^{1.5PN} = \frac{G}{c^3} \left\{ 2M \int_{0}^{+\infty} d\tau M_{ijklm}^{(7)}(u-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \frac{232}{105} \right] + \frac{20}{21} \int_{0}^{+\infty} d\tau \left[ M_{\langle ij \rangle klm}^{(3)} M_{klm}^{(4)} \right] (u-\tau) - \frac{15}{7} M_{\langle ij \rangle klm} M_{klm}^{(6)} \right. \\ \left. - \frac{41}{7} M_{\langle ij \rangle klm} M_{klm}^{(5)} - \frac{120}{7} M_{\langle ij \rangle klm} M_{klm}^{(4)} - \frac{710}{21} M_{\langle ij \rangle klm} M_{klm}^{(3)} - \frac{265}{7} M_{\langle ij \rangle klm} M_{klm}^{(4)} - \frac{155}{7} M_{\langle ij \rangle klm} M_{klm}^{(5)} - \frac{34}{7} M_{\langle ij \rangle klm} M_{klm}^{(6)} \right\}, \quad (2.17a) $$
$$V_{ijkl}^{1.5PN} = \frac{G}{c^3} \left\{ 2M \int_{0}^{+\infty} d\tau S_{ijkl}^{(6)}(u-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \frac{119}{60} \right] - \frac{11}{6} M_{\langle ij} S_{kl\rangle}^{(5)} - \frac{25}{6} M_{\langle ij}^{(1)} S_{kl\rangle}^{(4)} - \frac{25}{3} M_{\langle ij}^{(2)} S_{kl\rangle}^{(3)} - \frac{35}{3} M_{\langle ij}^{(3)} S_{kl\rangle}^{(2)} \right. \\ - \frac{65}{6} M_{\langle ij}^{(4)} S_{kl\rangle}^{(1)} - \frac{19}{6} M_{\langle ij}^{(5)} S_{kl\rangle} - \frac{11}{12} M_{\langle ijk}^{(5)} S_{l\rangle} + \epsilon_{ab\langle i} \left[ \frac{1}{12} M_{j\underline{a}} M_{kl\rangle b}^{(6)} + \frac{37}{60} M_{j\underline{a}}^{(1)} M_{kl\rangle b}^{(5)} - \frac{5}{12} M_{j\underline{a}}^{(2)} M_{kl\rangle b}^{(4)} - \frac{5}{6} M_{j\underline{a}}^{(3)} M_{kl\rangle b}^{(3)} \right. \\ \left. - \frac{11}{12} M_{j\underline{a}}^{(4)} M_{kl\rangle b}^{(2)} - \frac{1}{12} M_{j\underline{a}}^{(5)} M_{kl\rangle b}^{(1)} + \frac{3}{60} M_{j\underline{a}}^{(6)} M_{kl\rangle b} \right] \right\}. \quad (2.17b)$$

## B. Corrections due to the dimensional regularization of radiative moments

The results presented in the previous section are correct in the ordinary three-dimensional MPM algorithm. Nevertheless, the treatment of the dynamics of point masses coerces us to use a dimensional regularization scheme, starting at the 3PN order [112,113]. For consistency purposes, it is thus required to perform the MPM algorithm in d spatial dimensions, too. As proved in [21,93], this induces divergences in the expression of the radiative multipole moments in terms of the canonical ones. Those divergences

are in the form of simple poles, i.e., they scale as $1/\epsilon \equiv 1/(d-3)$. Their cancellation against counter-divergences, arising in the computation of the source multipole moments themselves, has been established in [92,93] and represents a crucial check of the method.

Let us recapitulate how the d-dimensional radiative multipole moments required for the 4PN flux differ from their three-dimensional counterparts. As the divergences hits at 3PN, only the mass quadrupole and octupole, as well as the current quadrupole, are affected. In the c.m. frame (see [93] for the complete expressions at 3PN in a generic frame), they are

$$U_{ij}^{(d)} = U_{ij} - \frac{214 G^2 M^2}{105 c^6} \left( \Pi_{\epsilon} + \frac{246299}{44940} \right) M_{ij}^{(4)} + \frac{G^2 M}{c^8} \left[ -\frac{4}{7} \left( \Pi_{\epsilon} - \frac{1447}{216} \right) M_{a\langle i} M_{j\rangle a}^{(6)} - \frac{32}{7} \left( \Pi_{\epsilon} - \frac{17783}{10080} \right) M_{a\langle i}^{(1)} M_{j\rangle a}^{(5)} \right. \\ \left. - 4 \left( \Pi_{\epsilon} - \frac{27649}{17640} \right) M_{a\langle i}^{(2)} M_{j\rangle a}^{(4)} + \frac{1921}{945} M_{a\langle i}^{(3)} M_{j\rangle a}^{(3)} + \frac{4}{3} \left( \Pi_{\epsilon} + \frac{11243}{7560} \right) M_{a\langle i}^{(5)} S_{j\rangle a} \right], \quad (2.18a)$$

$$U_{ijk}^{(d)} = U_{ijk} - \frac{21 G^2 M^2}{26 c^6} \left( \Pi_{\epsilon} + \frac{9281}{2730} \right) M_{ijk}^{(5)}, \quad (2.18b)$$

$$V_{k|ji}^{(d)} = V_{k|ji} - \frac{214 G^2 M^2}{105 c^6} \left( \Pi_{\epsilon} + \frac{4989}{44940} \right) S_{k|ji}^{(4)}, \quad (2.18c)$$

where we employ the notations of [94] for current moments, and where

$$\Pi_{\epsilon} = -\frac{1}{2\epsilon} + \ln\left(\frac{r_0\sqrt{\bar{q}}}{\ell_0}\right), \quad (2.19)$$

with $\bar{q} \equiv 4\pi e^{\gamma_E}$, $\gamma_E$ being the Euler constant. The constant $\ell_0$ is the length scale associated with dimensional regularization, such that the d-dimensional gravitational coupling constant $G_d = \ell_0^\epsilon G$ relates to the usual Newton's constant G, and $r_0$ was introduced in the previous section. As expected, the numerical constants associated with the poles at 3PN are the $\beta$-coefficients of the renormalization group flows for these multipole moments [21,108,109].<fnref n="2" />

As explained in Sec. IV the poles present in (2.18) will actually be canceled by poles coming from the d-dimensional expression of the source multipole moments. Thus, the correction terms (2.18), added to the d-dimensional source moments, permit to define a notion of finite ("renormalized") source moment, which is useful in some intermediate steps of our calculation (see more details in [93]).

## C. Relation between canonical and source moments

In the previous sections, the radiative moments have been linked to a set of canonical moments, $M_L$ and $S_L$. Nevertheless, the matching of the linearized metric, which serves as a basis for the MPM construction, to the PN source is more easily done by using a set of source moments $I_L$ and $J_L$, together with four families of gauge moments $W_L$, $X_L$, $Y_L$, and $Z_L$. We thus need to connect the canonical moments to those source and gauge moments, which is done via a coordinate transformation. This

<footnote n="2">²See [114] for the computation using the renormalization group flow of high-order logarithmic effects in the conservative energy function.</footnote>
coordinate change (i.e., fully nonlinear gauge transformation) induces a precise prescription for the relation between the canonical moments and the source moments, which has been worked out at 4PN order in three dimensions in [97]. However, since we are using full dimensional regularization in this work, we must generalize this computation to d dimensions in order to be consistent. Throughout this section, we will assume that the reader is familiar with the construction and notations of [97], and we will only highlight the main differences between the d-dimensional and the three-dimensional computations.

Regarding 4.5PN, nonlocal terms cannot appear at that order in the relation between canonical and source/gauge moments. Indeed, in the three-dimensional case, nonlocal terms do arise at cubic order, but dimensional analysis shows that there cannot be any cubic or higher-order contribution entering at exactly 4.5PN order. As for quadratic interactions entering at 4.5PN order, an analysis of the structure of the integration formulas proves that they are necessarily instantaneous. Therefore, as discussed previously, they will play no role in the flux for quasicircular orbits, and we only need the 4PN order in the canonical/source/gauge relations to control the 4.5PN quasicircular flux.

## 1. Procedure in d dimensions

The solution of the linearized, d-dimensional, vacuum Einstein field equations in harmonic gauge reads

$$h_{\text{gen}1}^{\mu\nu} = h_{\text{can}1}^{\mu\nu} + \partial^{\mu}\varphi_{1}^{\nu} + \partial^{\nu}\varphi_{1}^{\mu} - \partial_{\lambda}\varphi_{1}^{\lambda}\eta^{\mu\nu}, \quad (2.20)$$

where the canonical metric is parametrized by three types of source moments as [94]

$$h_{\text{can}1}^{00} = -\frac{4}{c^2} \sum_{\ell \ge 0} \frac{(-)^{\ell}}{\ell!} \partial_{\ell} \tilde{I}_{L}, \quad (2.21a)$$

$$h_{\text{can}1}^{0i} = \frac{4}{c^3} \sum_{\ell \ge 1} \frac{(-)^{\ell}}{\ell!} \left[ \partial_{L-1} \tilde{I}_{iL-1}^{(1)} + \frac{\ell}{\ell+1} \partial_L \tilde{J}_{i|L} \right], \quad (2.21b)$$

$$h_{\text{can}1}^{ij} = -\frac{4}{c^4} \sum_{\ell \ge 2} \frac{(-)^{\ell}}{\ell!} \left[ \partial_{L-2} \tilde{I}_{ijL-2}^{(2)} + \frac{2\ell}{\ell+1} \partial_{L-1} \tilde{J}_{(i|L-1)j}^{(1)} \right. \\ \left. + \frac{\ell-1}{\ell+1} \partial_L \tilde{K}_{ij|L} \right], \quad (2.21c)$$

and the gauge vector $\varphi_1^\mu$, by four types of gauge moments,

$$\varphi_1^0 = \frac{4}{c^3} \sum_{\ell \ge 0} \frac{(-)^{\ell}}{\ell!} \partial_L \tilde{W}_L, \quad (2.22a)$$

$$\varphi_1^i = -\frac{4}{c^4} \sum_{\ell \ge 0} \frac{(-)^{\ell}}{\ell!} \partial_{iL} \tilde{X}_L \\ -\frac{4}{c^4} \sum_{\ell \ge 1} \frac{(-)^{\ell}}{\ell!} \left[ \partial_{L-1} \tilde{Y}_{iL-1} + \frac{\ell}{\ell+1} \partial_L \tilde{Z}_{i|L} \right]. \quad (2.22b)$$

The conventions, notably for the indices of $\tilde{J}_{i|L}$ and $\tilde{Z}_{i|L}$, are still those of [94],<fnref n="3" /> and we have shortened

$$\tilde{F}(t, r) = \frac{\tilde{k}}{r^{d-2}} \int_{1}^{+\infty} dz \gamma_{\frac{1-d}{2}}(z) F(t - zr/c), \quad (2.23a)$$

with $\gamma_s(z) \equiv \frac{2\sqrt{\pi}}{\Gamma(1+s)\Gamma(-\frac{1}{2}-s)}(z^2-1)^s$

$$\text{and} \quad \tilde{k} \equiv \frac{\Gamma\left(\frac{d-2}{2}\right)}{\pi^{\frac{d-2}{2}}}. \quad (2.23b)$$

This shorthand satisfies $\lim_{d \to \tilde{F}(t, r) = F(t - r/c)/r$. Note the appearance of a new type of moment, $K_{ij|L}$, which is a pure artifact of working in d ≠ 3 dimensions, and will not enter our results since it has zero independent components in three dimensions, see Appendix A of [94], where the number of independent components is given by (A6b).

$$I_L = \begin{bmatrix} i_\ell & \dots & i_1 \end{bmatrix},$$

$$J_{i|L} = \begin{bmatrix} i_\ell & i_{\ell-1} & \dots & i_1 \\ i \end{bmatrix},$$

$$K_{ij|L} = \begin{bmatrix} i_\ell & i_{\ell-1} & i_{\ell-2} & \dots & i_1 \\ j & i \end{bmatrix}.$$

The current moment can be expressed using a Levi-Civita symbol in the $d \to 3$ limit. This relation is in both conventions: $J_{i|L} = \epsilon_{i\ell\alpha} J_{aL-1}$, where $J_L$ is symmetric in its indices. Note also that the conventions for $Z_{i|L}$ are the same as for $J_{i|L}$. In practice, the only difference between the two conventions is that the object $J_{(i|L-1)j}$ of (2.21c), valid in the HFB convention (where the underlined indices are excluded from the symmetrization), would read instead $J_{(i|j)L-1}$ in the BFL convention; see for example (2.5c) of [97].

<footnote n="3">³There has been some confusion [94,97,98] concerning these conventions, which we clarify here. In the HFB convention [94], the ordering of the multi-index in the multipolar moments is defined in an unusual order by $J_{i|L} \equiv J_{i|_{i' \dots i_1}}$ and $K_{ij|L} \equiv K_{ij|i_1 \dots i_r}$. In the BFL convention [97] however, the multi-index is always defined in the usual order, namely $J_{i|L} \equiv J_{i|_{i_1 \dots i_r}}$ and $K_{ij|L} \equiv K_{ij|i_1 \dots i_r}$. Note that the mass-type moment is entirely symmetric, so the ordering of its indices does not matter, namely $I_L = I_{i_1 \dots i_r} = I_{i'_1 \dots i'_r}$. In both conventions, the Young tableaux of the moments read (with the notation of [97]),</footnote>
Starting from this d-dimensional metric, we can follow the lines of Sec. III. A in [97], which are independent of the dimension. Notably, at any PM order $n$, the expressions of the nonlinear interaction terms $\Omega_n^{\mu\nu}$ and $\Delta_n^\mu$ as functions of $h_{\text{can}}^{\mu\nu}$ and $q_n^\mu$ are identical to the three-dimensional case. The main, but minor, difference concerns the two key quantities $X_n^{\mu\nu}$ and $Y_n^{\mu\nu}$. If their formal definition, Eq. (3.20) of [97], does not depend of the dimension, their explicit expressions in terms of $\Omega_n^{\mu\nu}$ and $\Delta_n^\mu$ are slightly changed, and they read in d dimensions

$$X_n^{\mu\nu} = \text{FP}_{B=0} \Box_{\text{ret}}^{-1} \left[ B \left( \frac{r}{r_0} \right)^B \left( -\frac{B+d-2}{r^2} \Omega_n^{\mu\nu} - \frac{2}{r} \partial_r \Omega_n^{\mu\nu} \right) \right], \quad (2.24a)$$

$$Y_n^{\mu\nu} = \text{FP}_{B=0} \Box_{\text{ret}}^{-1} \left[ B \left( \frac{r}{r_0} \right)^B \frac{n_k}{r} \left( -2\delta^{k(\mu} \Delta_n^{\nu)} + \eta^{\mu\nu} \Delta_n^k \right) \right]. \quad (2.24b)$$

These quantities are endowed with a Hadamard finite part (FP) regularization when $B \to 0$, here defined on top of the dimensional regularization. Note the explicit prefactor $B$ in the source of the retarded d'Alembertian operators in (2.24). The quantities $X_n^{\mu\nu}$ and $Y_n^{\mu\nu}$ will be nonzero only when the retarded integrals develop a pole $\propto B^{-1}$. The remaining of the procedure described in Sec. III of [97], and notably the extraction of the correction to the canonical moments, is left unchanged. All the difference will reside in the detailed integration formulas to compute (2.24).

## 2. Integration techniques

Let us recall that the solution at an order $n \ge 2$ is composed of some corrections to the canonical moments, extracted from $X_n^{\mu\nu}$ and $Y_n^{\mu\nu}$ given by (2.24), as well as a gauge vector, denoted $q_n^\mu$. Thanks to the prefactor $B$ in (2.24), and similarly to the three-dimensional case, the computation of $X_n^{\mu\nu}$ and $Y_n^{\mu\nu}$ only requires the near-zone ($r \to 0$) behavior of the functions $\Omega_n^{\mu\nu}$ and $\Delta_n^\mu$, as the FP regularization is introduced to cope with the singular behavior of multipole expansions when $r \to 0$. Using the expansion formulas displayed in Appendix A of [81], the computation boils down to the evaluation of integrals of the form

$$\mathcal{J}_{p,q}^b \equiv \text{FP}_{B=0} \Box_{\text{ret}}^{-1} \left[ \left( \frac{r}{r_0} \right)^B B^b \frac{\hat{n}_L}{r^{p+q\epsilon}} H(t) \right], \quad (2.25)$$

where $b=1, 2$ and $H$ can bear simple poles $\propto 1/\epsilon$ and can indifferently be a local or nonlocal function of time. In this section, we set $r_0 = \ell_0 = 1$ for simplicity. The integrals (2.25) are nothing but the d-dimensional generalization of (4.13) of [97], with $a=0$, since the case which includes a $\ln r$ is irrelevant for the present generalization.

Techniques to solve wavelike equations in d dimensions, with more general radial dependence than (2.25) but same multipolarity $\ell$, have been developed in [93], see Sec. II therein. Importantly, no particular assumption about the presence of a prefactor $B^b$ has been made in this work, as clear from its Eq. (2.6). We can therefore follow the lines of [93] to compute the integral (2.25). The solution $h$ is the sum of two contributions, $h^<$ and $h^>$, involving integrals over the respective domains $\mathcal{D} = \{|\mathbf{x}'| < r\}$ and $\mathbb{R}^3 - \mathcal{D}$, as displayed explicitly in Eqs. (2.8) and (2.10) of [93]. When a prefactor $B^b$ is present, with $b=1, 2$, the finite part of $h \equiv \mathcal{J}_{p,q}^b$ vanishes unless the integration produces a pole, which may only occur when the integral diverges for $B=0$. As the system is taken to be stationary in the past in the MPM algorithm, the integral yielding $h^>$ converges for $B=0$ so that only $h^<$ can contribute to the finite part in Eq. (2.25). This finite part is evaluated explicitly and truncated at first order in $\epsilon$, in order to cope with the case where $H(t)$ contains a pole, namely $H = \frac{1}{\epsilon}H_{-1} + H_0 + \mathcal{O}(\epsilon)$. After resummation, we obtain a retarded homogeneous solution of the wave equation [93]

$$\begin{aligned} \mathcal{J}_{p,q}^b &= \hat{\partial}_L \tilde{G}_{p,q}^b(t, r) \\ &\equiv \hat{\partial}_L \left[ \frac{\tilde{k}}{r^{d-2}} \int_1^{+\infty} dz \gamma_{\frac{1}{2}d}(z) G_{p,q}^b(t - zr/c) \right]. \end{aligned} \quad (2.26)$$

Explicitly, we find that $G_{p,q}^b$ is nonvanishing only if $b=1$, $q=1$ and $p-\ell-3=2j$, where $j \in \mathbb{N}$. Under those very restrictive conditions, we have

$$\begin{aligned} G_{p,q}^b(t) = & \frac{(-)^{\ell+1}}{(2j)!!(2j+2\ell+1)!!} \\ & \times \left[ \left(1 + \epsilon \ln \sqrt{\bar{q}} - \epsilon \sum_{k=0}^{j+\ell} \frac{1}{2k+1} \right) \right. \\ & \left. \times H^{(2j)}(t) + \mathcal{O}(\epsilon) \right], \end{aligned} \quad (2.27)$$

where we recall that $H$ can bear a pole $1/\epsilon$, and our notation $\bar{q} = 4\pi\epsilon^{\gamma_E}$ (with $\gamma_E$ the Euler constant). With the previous formula at hand, we are ready to compute the quantities $X_n^{\mu\nu}$ and $Y_n^{\mu\nu}$ and, therefore, obtain the correction to the canonical moments to the nth PM order. Since the result (2.27) is zero unless $q=1$, we are able to eliminate many terms from the source for this computation solely based on the value of $q$.

As for the gauge vector $q_n^\mu$, it will be needed to compute the source term for the next PM order $n+1$, following the algorithmic procedure of [97]. It is decomposed as $q_n^\mu = \phi_n^\mu + \psi_n^\mu$, where $\psi_n^\mu$ is also extracted from $X_n^{\mu\nu}$ and $Y_n^{\mu\nu}$, and $\phi_n^\mu$ is obtained by the direct integration
$$ \phi_n^\mu \equiv \mathrm{FP}_{B=0} \Box_{\mathrm{ret}}^{-1} \left[ \left( \frac{r}{r_0} \right)^B \Delta_n^\mu \right]. \quad (2.28) $$

In d dimensions, it is in general not possible to find a convenient, explicit expression for $\phi_n^\mu$ (see, however, Sec. II in [93]). Fortunately, we do not need here the full solution valid at any field point, but only the solution in the form of a near-zone expansion when $r \to 0$. Indeed, following the procedure of [97], it is the near-zone expansion of $\phi_n^\mu$, denoted $\bar{\phi}_n^\mu$, that needs to be inserted into the expressions of $\Omega_{n+1}^{\mu\nu}$ and $\Delta_{n+1}^\mu$ sourcing the next-order quantities $X_{n+1}^{\mu\nu}$ and $Y_{n+1}^{\mu\nu}$. The quantity $\bar{\phi}_n^\mu$ can be obtained directly from the near-zone expansion of the source term, i.e., $\bar{\Delta}_n^\mu$, plus a crucial homogeneous solution, itself in the form of a near-zone expansion, namely,

$$ \bar{\phi}_n^\mu = \bar{\phi}_{\mathrm{part} n}^\mu + \bar{\phi}_{\mathrm{hom} n}^\mu. \quad (2.29) $$

The particular solution is obtained by a formal iteration of inverse Laplace operators in d dimensions, using the “Matthieu” formula, Eq. (B.26c) in [113]

$$ \bar{\phi}_{\mathrm{part} n}^{\mu} = \mathrm{FP}_{B=0} \sum_{k=0}^{+\infty} \Delta^{-1-k} \left(\frac{\partial}{c\partial t}\right)^{2k} \left[\left(\frac{r}{r_0}\right)^B \bar{\Delta}_n^{\mu}\right]. \quad (2.30) $$

The homogeneous solution $\bar{\phi}_{\mathrm{hom} n}^\mu$ is given by Eq. (3.20) of [81], which in this case yields (with $c=1$)

$$ \bar{\phi}_{\mathrm{hom} n}^{\mu} = -\sum_{\ell \ge 0} \sum_{j=0}^{+\infty} \frac{\Delta^{-j} \hat{x}_L}{d + 2\ell - 2} \mathrm{FP}_{B=0} \int_1^{+\infty} dz \gamma_{\frac{1-d}{2}-\ell}(z) \\ \times \int_0^{+\infty} dr' r'^{-\ell+1+B} \hat{n}_L' \bar{\Delta}_{n,L}^{\mu(2j)}(t-zr', r'), \quad (2.31) $$

where we have decomposed the source in its multipolar components as

$$ \bar{\Delta}_n^\mu(t, \mathbf{x}) = \sum_{\ell \ge 0} \hat{n}_L \bar{\Delta}_{n,L}^\mu(t, r), \quad (2.32) $$

and shortened the iterated inverse Laplacians as [81,115,116]

$$ \Delta^{-j} \hat{x}_L = \frac{\Gamma\left(\frac{d}{2} + \ell\right)}{\Gamma\left(\frac{d}{2} + \ell + j\right)} \frac{r^{2j} \hat{x}_L}{2^{2j} j!} \quad (2.33) $$

In the particular case where the source term has a structure given by

$$ \bar{\Delta}_n^\mu(t, r) = \sum_{\ell, k, p, q} \frac{\hat{n}_L}{r^{p+q\ell}} \int_1^{+\infty} dy \, y^k \gamma_{-1-\ell/2}(y) F_{\ell, k, p, q}(t - yr), \quad (2.34) $$

where the $F_{\ell, k, p, q}(t)$'s can be arbitrary functions of time, explicit formulas are known for performing the integrations; see (3.24) and (D1) of [81] when $q=2$, and the discussion in Sec. IV of [93] for other values of $q$. However, we will encounter source terms that exhibit more complicated structures, and for which no analogous formula is available. This is notably the case for sources involving two nonstatic multipolar moments, such as $\Delta_{W \times I_{ij}}^\mu$. Fortunately, we have shown that those more complicated homogeneous solutions do not contribute at 4PN order, see hereafter.

### 3. The 4PN relation between canonical and source moments

At the 4PN order, only quadratic and cubic interactions are allowed by dimensional analysis (and we recall that the 4.5PN sector vanishes in the quasicircular flux). Naturally, the moments allowed to enter those interactions do not depend on the space dimension. Applying the method described above, it was found that the quadratic order is identical in d and three dimensions, i.e., we recover the same “odd” corrections at 2.5PN and 3.5PN orders in the relation between canonical and source/gauge moments. More interestingly, the results for the cubic interactions arising at 4PN order greatly differ.

Let us recall that three interactions enter the 4PN mass quadrupole moment at cubic order, namely $M \times M \times W_{ij}$ and $M \times M \times Y_{ij}$ (where only one moment is nonstatic), and most importantly $M \times W \times I_{ij}$ (where the two moments W and I$_{ij}$ are dynamical). For the first two interactions, it turns out that the source term entering the integral (2.25) bears $q \ge 2$. Indeed, it is composed of two interactions: $h_{M \times M} \times \varphi_{K_{ij}}$ and $h_M \times \varphi_{M \times K_{ij}}$, where $K_{ij}$ stands for either $W_{ij}$ or $Y_{ij}$. As both $h_{M \times M}$ and $\varphi_{M \times K_{ij}}$ bear $q=2$, the cubic sources will have $q \ge 2$, and thus cannot satisfy the $q=1$ necessary condition to have a nonvanishing function $G(t)$ in Eq. (2.27). Therefore, the $M \times M \times W_{ij}$ and $M \times M \times Y_{ij}$ interactions do not contribute in d dimensions, as we have explicitly checked, using the method presented in the previous sections. This vanishing of the $M \times M \times W_{ij}$ and $M \times M \times Y_{ij}$ interactions in d dimensions contrasts with their explicit contributions in three dimensions; see Eq. (1.1) in [97].

As for the last cubic interaction, $M \times W \times I_{ij}$, it turns out that the source term does have a $q=1$ component, arising from the interaction $h_{I_{ij}} \times \varphi_{M \times W}$. More precisely, the term with $q=1$ comes exclusively from the interaction between: (i) the homogeneous solution $\bar{\phi}_{\mathrm{hom} 2}$ at quadratic order arising from the interaction $M \times W$, i.e., $\bar{\phi}_{\mathrm{hom} M \times W}$, which bears $q=0$; and (ii) the $q=1$ piece of the linear metric $h_{I_{ij}}$; see (2.21). Note that, as it involves a quadratic interaction with only one dynamical moment, $\varphi_{M \times W}$ can be calculated using (2.34). By contrast, the homogeneous solution for the interaction $W \times I_{ij}$, i.e., $\bar{\phi}_{\mathrm{hom} W \times I_{ij}}$, does not
fulfill the $p - \ell - 3 \in 2\mathbb{N}$ condition, and as such, does not contribute to the result. This is fortunate, since it does not follow the structure (2.34), and we are unable to compute this term explicitly, at least with currently-known formulas. Another interesting observation is that the explicit computation of the $q = 1$ terms (the only ones that contribute to the final result) leads to the appearance of a pole. The relation between canonical and source/gauge moments at 4PN is thus a pure $M \times W \times I_{ij}$ interaction, which reads

$$ \delta M_{ij} = \frac{8G^2M}{c^8} \int_0^{+\infty} d\tau \left[ \frac{1}{\epsilon} - 2 \ln \left( \frac{c\tau\sqrt{\bar{q}}}{2} \right) - 1 \right] \\ \times (I_{ij}^{(1)}(t)W^{(3)}(t-\tau) - I_{ij}(t)W^{(4)}(t-\tau)) \\ + \mathcal{O}\left(\frac{1}{c^{10}}\right). \qquad (2.35) $$

This result is consistent with the corresponding one in three dimensions, where the interaction $M \times W \times I_{ij}$ gives rise to an ordinary tail integral at 4PN order: the coefficient of the pole in (2.35) matches the logarithmic structure of Eq. (1.1) in [97]. Moreover, we have checked the result (2.35), as well as the vanishing of the interactions $M \times M \times W_{ij}$ and $M \times M \times Y_{ij}$, by an independent method, using the techniques exposed in [93], to compute the difference between the $d$-dimensional and Hadamard regularization schemes in the MPM algorithm, and adding it to the three-dimensional result of [97].

The new pole in (2.35) arising from dimensional regularization is actually canceled by another pole. This stems from an extra contribution to be added to the end result of [92], which we will now discuss. Indeed, when computing the equations of motion at 4PN order, the tail effect (entering as a radiation mode in the conservative dynamics described by the Fokker action) has been implemented using the canonical quadrupole moment $M_{ij}$ [81], and is therefore valid in the associated “canonical” coordinate system. By contrast, the source quadrupole moment $I_{ij}$, as given in [92], is expressed in a so-called “generic” coordinate system, following the terminology of [97]. Note that the computation of the source quadrupole moment $I_{ij}$ at 4PN does not use the 4PN equations of motion. These two results, taken independently, are therefore perfectly correct in their respective coordinate systems. However, in this work, we need to take time derivatives of the source quadrupole moment, which crucially requires the 4PN equations of motion. This operation can only be performed if the source quadrupole moment and the equations of motion are expressed in the same coordinate system. We thus compute the expression of the source quadrupole moment in the canonical coordinate system

rather than in the generic one. This is done by performing a simple coordinate change $x_{\text{gen}}^{\mu} = x_{\text{can}}^{\mu} + \xi^{\mu}$, which transforms the gravitational field according to Eqs. (3.1)–(3.3) in [97]. A simple dimensional analysis shows that only the $M \times W$ interaction can enter this coordinate change at 4PN order. Using the methods exposed previously to calculate the contribution of a given interaction to the gauge vector, we find that

$$ \begin{aligned} \zeta^0 = \varphi_{\text{hom M\times W}}^0 = & \frac{8G^2M}{c^7} \sum_{j \ge 0} \frac{\Gamma\left(\frac{d}{2}\right)}{\Gamma\left(\frac{d}{2} + j\right)} \frac{r^{2j}}{2^{2j}j!c^{2j}} \\ & \times \int_0^{+\infty} d\tau \left[ \frac{1}{\epsilon} - 2 \ln \left( \frac{c\tau\sqrt{\bar{q}}}{2} \right) - 1 \right] W^{(3+2j)}(t-\tau), \end{aligned} \qquad (2.36a) $$

$$ \zeta^i = \varphi_{\text{hom M\times W}}^i = 0, \qquad (2.36b) $$

where only $j = 0$ is relevant at 4PN. Implementing the procedure exposed in Sec. II of [92] with the “pure gauge” metric $h^{\mu\nu} = -\partial^\mu \zeta^\nu - \partial^\nu \zeta^\mu + \eta^{\mu\nu} \partial_\lambda \zeta^\lambda$, we arrive at

$$ \begin{aligned} \delta_{\zeta} I_{ij} = & -\frac{8G^2M}{c^8} \int_0^{+\infty} d\tau \left[ \frac{1}{\epsilon} - 2 \ln \left( \frac{c\tau\sqrt{\bar{q}}}{2} \right) - 1 \right] \\ & \times (I_{ij}^{(1)}(t)W^{(3)}(t-\tau) - I_{ij}(t)W^{(4)}(t-\tau)) \\ & + \mathcal{O}\left(\frac{1}{c^{10}}\right). \end{aligned} \qquad (2.37) $$

This effect is formally a correction to the source mass quadrupole moment $I_{ij}$. Nevertheless, it is more handy to identify it as a correction in the relation between the canonical and source/gauge moments. At the 4PN order, no obstacle interferes with such identification. Thus, we will consider it as an “indirect” contribution in this relation, namely $\delta_{\text{indirect}}M_{ij} \equiv \delta_{\zeta}I_{ij}$, to be added to the “direct” one given by (2.35). We then find that, at least at 4PN order, the direct and indirect contributions exactly cancel,

$$ \delta_{\text{direct}}M_{ij} + \delta_{\text{indirect}}M_{ij} = 0. \qquad (2.38) $$

Therefore, we have proved that all cubic interactions in $d$ dimensions are vanishing at 4PN order, in clear contrast with the three-dimensional result, Eq. (1.1) of [97]. Finally, as there is no pole left in the final relation between the source and canonical quadrupole, we can reduce it directly to three dimensions. This relation, using full dimensional regularization, reads
$$M_{ij} = I_{ij} + \frac{4G}{c^5} [W^{(2)}I_{ij} - W^{(1)}I_{ij}^{(1)}] + \frac{4G}{c^7} \left\{ \frac{4}{7} W_{a\langle i}^{(1)} I_{j\rangle a}^{(3)} + \frac{6}{7} W_{a\langle i}^{(1)} I_{j\rangle a}^{(4)} - \frac{1}{7} Y_{a\langle i}^{(3)} I_{j\rangle a} - Y_{a\langle i}^{(1)} I_{j\rangle a}^{(3)} - 2XI_{ij}^{(3)} \right. \\ - \frac{5}{21} W_a^{(4)} I_{ija} + \frac{1}{63} W_a^{(3)} I_{ija}^{(1)} - \frac{25}{21} Y_a^{(3)} I_{ija} - \frac{22}{63} Y_a^{(2)} I_{ija}^{(1)} + \frac{5}{63} Y_a^{(1)} I_{ija}^{(2)} \\ + 2W^{(3)}W_{ij} + 2W^{(2)}W_{ij}^{(1)} - \frac{4}{3}W_{\langle i}W_{j\rangle}^{(3)} + 2W^{(2)}Y_{ij} - 4W_{\langle i}Y_{j\rangle}^{(2)} \\ + \epsilon_{ab\langle i} \left[ \frac{1}{3}I_{j\rangle a}Z_b^{(3)} - I_{j\rangle a}^{(3)}Z_b + \frac{4}{9}J_{j\rangle a}W_b^{(3)} - \frac{4}{9}J_{j\rangle a}Y_b^{(2)} + \frac{8}{9}J_{j\rangle a}Y_b^{(1)} \right] \left. \right\} + \mathcal{O}\left(\frac{1}{c^9}\right). \quad (2.39)$$

For the other multipole moments (essentially the mass octupole and the current quadrupole), such relations are purely quadratic and consist of odd parity terms, and so are not affected by dimensional-regularization corrections. We do not report them here, as they can be found in Sec. III. B of [95].

Since we have absorbed the correction $\delta_{\text{indirect}}M_{ij}$ into the redefinition of the canonical moment, the relation (2.39), derived in the context of dimensional regularization, now holds with exactly the same definition for the source quadrupole moment as proposed in [93]. Namely, $I_{ij}$ is the “renormalized” source quadrupole moment (in the sense of [93]) which is given for quasicircular orbits in Sec. IVA, and which contains crucial contributions due to the (IR) dimensional regularization. Adding up the nonlinear contributions described in Sec. II, we obtain the radiative moment measured at infinity. As we will see the above procedure yields the correct perturbative limit for compact binaries on circular orbits, in agreement with first-order black hole perturbation theory.

Interestingly, we found that the perturbative limit turns out also to be correct when using a simpler treatment of the IR divergences of the source quadrupole moment, based on the Hadamard regularization in ordinary three dimensions, rather than the dimensional regularization. In this case, the correct relation between canonical and source moments is given by Eq. (1.1) of [97], instead of the d-dimensional result (2.39). We find that the cubic interactions $M \times M \times W_{ij}$ and $M \times M \times Y_{ij}$ present in three dimensions do contribute to the perturbative limit for circular orbits (note that W is zero in this case), but in such a way as to cancel the contribution due to the difference between the Hadamard and dimensional regularization schemes for the mass quadrupole source moment. It turns out then that we also recover the correct perturbative limit. However, in contrast to the dimensional regularization, we do not expect the Hadamard regularization in three dimensions to lead to a well-defined and fully unambiguous result beyond this limit, so the above observation should be considered only as an interesting consistency check.

### III. CORRECTIONS TO THE SOURCE MOMENTS DUE TO INFRARED COMMUTATORS

In this section, we discuss the appearance of infrared (IR) “d’Alembertian commutators” in the expression of the

PN-expanded near-zone metric, and their effects on the expression of the source moments. This feature starts at 4PN order and does not affect the 4PN equations of motion, but should have been considered in previous works [91–93]. We show that the IR commutators enter the flux and the (2, 2) mode for generic orbits, however they do not contribute in the quasicircular limit. Since they are zero for quasicircular orbits, all the results explicitly presented in the previous works [91–93] are correct. Throughout this section, we will be based on the construction and computation of the source mass quadrupole, as described in [91].

#### A. Appearance of the commutators in the PN metric

Crucial to the derivation of the source moments at a given PN order is the expression of the near-zone metric with the same precision. Hereafter, we will present the case of the source mass quadrupole at 4PN order, but the following discussion is easily generalized to any source moment. The gothic metric perturbation $h^{\mu\nu} \equiv \sqrt{-g}g^{\mu\nu} - \eta^{\mu\nu}$ obeys the gauge-fixed Einstein field equations

$$\Box h^{\mu\nu} = \frac{16\pi G}{c^4} \tau^{\mu\nu}, \quad (3.1)$$

where the pseudo stress-energy tensor $\tau^{\mu\nu}$ is defined e.g., in Eq. (2.4) of [91]. As explained in [115], the PN metric is obtained by solving iteratively the wave equation (3.1) in the near zone, in such a way as to match the metric in the far zone. The constructed expansion of the gravitational field in the near zone, say $\bar{h}^{\mu\nu}$, is then given formally, to arbitrarily high orders, by [81,117]

$$\bar{h}^{\mu\nu} = \frac{16\pi G}{c^4} [\text{FP}\bar{\Box}_{\text{ret}}^{-1}[\bar{\tau}^{\mu\nu}]] + \sum_{\ell=0}^{+\infty} \sum_{j=0}^{+\infty} \Delta^{-j} \hat{x}_L \frac{1}{c^{2j}} f_L^{(2j)\mu\nu}(t). \quad (3.2)$$

The second term on the right-hand side of (3.2) is an homogeneous solution that is regular at the coordinate origin $r = 0$. The explicit d-dimensional expressions of the time-dependent functions $f_L^{\mu\nu}(t)$, as functionals of the
MPM expansion of $\tau^{\mu\nu}$ in the exterior zone, can be found in [117].

Of interest to us here is the first term in (3.2). It is made by the PN-expanded retarded integral $FP\bar{\Box}_{\text{ret}}^{-1}$, i.e., the inverse d'Alembertian operator in which all retardations are formally expanded in a PN fashion, and regularized by means of the finite part procedure, $FP_{B=0} \int d^d\mathbf{x}' (|\mathbf{x}'|/r_0)^B[\cdots]$. Furthermore, the operator $FP\bar{\Box}_{\text{ret}}^{-1}$ acts on the near-zone expansion of the pseudo stress-energy tensor $\bar{\tau}^{\mu\nu}$, which can be derived at any PN order from the lower order expressions of $\bar{h}^{\mu\nu}$. The IR behavior of the integrals entering the propagator is naturally affected by the choice of regularization. In general, the PN expanded propagator in (3.2) has no reason to commute with the d'Alembertian operator $\Box$. This is the root of the appearance of the IR commutators in the PN metric, that are defined as

$$C\{F\} = [FP\bar{\Box}_{\text{ret}}^{-1}, \Box]F = FP\bar{\Box}_{\text{ret}}^{-1} \left[ \left( \frac{r}{r_0} \right)^B \Box F \right] - F, \quad (3.3)$$

where $F$ is a regular (smooth) function in the near zone but typically diverges at spatial infinity, when $r \to +\infty$; see (3.11) below.

In order to simplify the resolution of the Einstein field equations (3.1), we parametrize the PN metric by means of some PN “potentials”, obeying flat space-time wave equations in $d$ dimensions

$$\Box P = S, \quad (3.4)$$

where $S$ is composed of a compact-supported sector (proportional to Dirac distributions in the case of point particles), and a noncompact one (composed of products of derivatives of lower-order PN potentials). We refer to Appendix A of [91] for complete definitions of the PN potentials. Now those potentials are constructed by means of the operator $FP\bar{\Box}_{\text{ret}}^{-1}$, i.e., they are given, by definition, as

$$P = FP\bar{\Box}_{\text{ret}}^{-1} S, \quad (3.5)$$

which accounts for the first term in the right side of (3.2). The second term in (3.2) is nonlocal in time and corresponds to radiation modes that are being accounted for separately [81].

When parametrizing the PN metric in this manner, it is very handy to recognize products of lower-order potentials. For instance, if some relevant combination of the field equations (3.1) contains terms of the form $\partial_i P \partial_i P$ at some $n$PN order,<fnref n="4" />

$$\Box \bar{h}_{(n)}^{00ii} = \cdots + \partial_i P \partial_i P + \cdots, \quad (3.6)$$

then it can be solved using the PN-expanded retarded propagator $FP\bar{\Box}_{\text{ret}}^{-1}$, as

$$\begin{align*}
\bar{h}_{(n)}^{00ii} &\equiv \cdots + FP\bar{\Box}_{\text{ret}}^{-1} \{\partial_i P \partial_i P\} + \cdots \\
&= \cdots + FP\bar{\Box}_{\text{ret}}^{-1} \left\{ \frac{1}{2} \Box(P^2) - P \Box P \right\} + \cdots \\
&= \cdots + \frac{P^2}{2} - FP\bar{\Box}_{\text{ret}}^{-1} \{PS\} + \frac{1}{2} C\{P^2\} + \cdots, \quad (3.7)
\end{align*}$$

where we have used Eq. (3.4) and neglected the higher PN order term $(\partial_i P)^2/c^2$ in the second curly brackets of the first line. Because solving $FP\bar{\Box}_{\text{ret}}^{-1}\{PS\}$ is much more tractable than the original $FP\bar{\Box}_{\text{ret}}^{-1}\{\partial_i P \partial_i P\}$, the operation described by Eq. (3.7) constitutes a real improvement. As shown later on, the commutator left in Eq. (3.7) is in general not vanishing. The expression of the 4PN metric in terms of those commutators is displayed in Appendix B.

While, formally, the commutators appear at a relative 1PN order in the metric [see Eq. (B1)], they effectively only contribute starting at 4PN order, as shown in Sec. III C. Therefore, the only source multipole moment to be affected is the 4PN source mass quadrupole. Moreover, the 4PN contribution in the metric arises as a mere function of time in $\bar{g}_{00}$ only, through the combination $\bar{h}^{00ii}$. Since it is the spatial gradient $\partial_i \bar{g}_{00}$ that contributes to the equations of motion, the corresponding contribution vanishes. More generally, following the “n + 2 method” [79] to determine the PN order of the metric needed for insertion in the Fokker Lagrangian, we see that the commutators play no role in the derivation of the conservative equations of motion at 4PN order.

## B. Expression of the commutators

We immediately see from (3.3) that the commutator must be a homogeneous solution of the wave equation. Furthermore, for the IR commutator, which depends on the far zone behavior $r \to +\infty$, that homogeneous solution must be regular in the near zone, when $r \to 0$. These facts constrain its form: given the function $F$, there exist a set $\{\hat{g}_L\}$ of (yet unspecified) functions of time only, such that

$$C\{F\} = \sum_{\ell=0}^{+\infty} \frac{(-)^{\ell}}{\ell!} \sum_{j=0}^{+\infty} \frac{1}{c^{2j}} \Delta^{-j} \hat{x}_L \hat{g}_L^{(2j)}(t), \quad (3.8)$$

where we recall that the inverse Laplacians are defined in Eq. (2.33). In order to explicitly determine the functions $\hat{g}_L$, we work out the propagator applied to a generic source $S(\mathbf{x}, t)$,

<footnote n="4">⁴For convenience we pose $h^{00ii} \equiv \frac{2}{(d-1)}[(d-2)h^{00} + h^{ii}]$ in $d$ dimensions.</footnote>
$$ \mathrm{FP}\Box_{\mathrm{ret}}^{-1}S = -\frac{\tilde{k}}{4\pi} \mathrm{FP}_{B=0} \sum_{k \ge 0} \int d^d \mathbf{x}' \left(\frac{r'}{r_0}\right)^B \frac{\Gamma\left(2-\frac{d}{2}\right) S^{(2k)}(\mathbf{x}', t)}{2^{2k} k! c^{2k} \Gamma\left(k+2-\frac{d}{2}\right) |\mathbf{x}-\mathbf{x}'|^{d-2}} \left[ \frac{\sqrt{\pi}}{2\Gamma\left(\frac{5-d}{2}\right) \Gamma\left(k+\frac{d}{2}\right)} \int_0^{+\infty} d\tau \frac{\tau^{3-d}}{c^{d-2}} S^{(2k+2)}(\mathbf{x}', t-\tau) \right]. \quad (3.9) $$

This expression comes from the near-zone expansion of the homogeneous solution of the wave equation in d dimensions. The two terms in (3.9) represent respectively the even and odd parity parts of this expansion, as given by Eqs. (A7) and (A8) in [81]. Note that the odd parity part involves a nonlocal in time integral, shown in the second term of (3.9). The PN propagator (3.9) corresponds exactly to the prescription which is required for the first term in Eq. (3.2), i.e., where all retardations of the inverse d'Alembertian operator are PN expanded. In three dimensions, this prescription recovers Eqs. (3.4) and (3.7) of [116]. In turn, the prescription ensures the correct matching of the PN metric to the MPM metric in the exterior zone.

When computing the commutator $C\{F\}$ with Eqs. (3.3) and (3.9), we may transform the d'Alembertian operator inside the source term by means of appropriate integrations by parts, which yields for the commutator a finite part integral with a global factor B appearing in the source,

$$ C\{F\} = \mathrm{FP}_{B=0} \Box_{\mathrm{ret}}^{-1} \left[ B \left( \frac{r}{r_0} \right)^B \left( -\frac{B-d}{r^2} F - \frac{2}{r} \partial_r F \right) \right]. \quad (3.10) $$

Note the similarity with the expression of the quantities $X_n^{\mu\nu}$ and $Y_n^{\mu\nu}$ in (2.24). Indeed, the latter have been formally defined as commutators [see Eq. (3.20) of [97]], but, in contrast to the IR commutators (3.3) or (3.10), $X_n^{\mu\nu}$ and $Y_n^{\mu\nu}$ are ultraviolet (UV) commutators since they depend on the near zone behavior of the source term, when $r \to 0$.

We now use the expression (3.10) of the commutator. As we said, because of the explicit factor B in the source, the commutator depends only on the IR behavior of the function

F, i.e., when $r \to +\infty$ with $t = \text{const}$ (spatial infinity). We expand the function in $d = 3 + \varepsilon$ dimensions as

$$ F(\mathbf{x}, t) = \sum_{\substack{p_0 \le p \\ q_0 \le q \le q_1}} \frac{f_{p,q}(\mathbf{n}, t)}{r^{p+q\varepsilon}}, \quad (3.11) $$

where the sums extend on all multiplicities $\ell$, and formally for all values of $p$ larger than some (generally negative) integer $p_0$, and for a finite set of values of $q$ for each given $p$. We inject this expansion into the commutator (3.10) and employ the explicit expression (3.9) for the PN expanded propagator. Consistently with the IR limit we must expand the factors $|\mathbf{x} - \mathbf{x}'|^\alpha$ in (3.9) (where $\alpha$ depends on the parity) when $|\mathbf{x}'| \to +\infty$. Actually this expansion is valid as soon as $r' = |\mathbf{x}'| > r = |\mathbf{x}|$ and reads, using the STF decomposition and the notation (2.33),<fnref n="5" />

$$ |\mathbf{x}' - \mathbf{x}|^{\alpha} = \sum_{\ell=0}^{+\infty} \frac{(-)^{\ell}}{\ell!} \sum_{j=0}^{+\infty} \Delta^{-j} \hat{x}_L \Delta^{j} \hat{\partial}_L^{\prime j} r^{\prime \alpha}, \quad (3.12) $$

the object $\Delta^{ij}\hat{\partial}_L^{\prime j}r^{\prime \alpha}$ being explicitly known as

$$ \Delta^{ij} \hat{\partial}_L^{\prime j} r^{\prime \alpha} = 2^{\ell+2j} \frac{\Gamma\left(\frac{\alpha}{2}+1\right) \Gamma\left(\frac{\alpha+d}{2}\right)}{\Gamma\left(\frac{\alpha}{2}+1-\ell-j\right) \Gamma\left(\frac{\alpha+d}{2}-j\right)} \hat{x}_L^{\prime j \alpha-2\ell-2j}. \quad (3.13) $$

With the latter formulas at hand, it is straightforward to obtain the functions $\hat{g}_L(t)$ parametrizing the commutator in (3.8) in terms of the expansion coefficients in (3.11). We find

$$ \hat{g}_L(t) = \sum_{k=0}^{+\infty} \frac{4k - 2\ell - d + 2}{2^{2k} k! c^{2k} (d-2)} \left[ \frac{2^\ell \Gamma\left(2-\frac{d}{2}\right)}{\Gamma\left(k+2-\ell-\frac{d}{2}\right)} \hat{f}_{2k-\ell,0}^{L(2k)}(t) - \frac{\sqrt{\pi}}{2^{\ell+1} \Gamma\left(\frac{5-d}{2}\right) \Gamma\left(k+\ell+\frac{d}{2}\right)} \frac{4k+2\ell+d-2}{4k-2\ell-d+2} \int_0^{+\infty} d\tau \frac{(c\tau)^{3-d}}{c^{2\ell+1}} \hat{f}_{2k+\ell+1,1}^{L(2k+2\ell+2)}(t-\tau) \right], \quad (3.14) $$

$$ |\mathbf{x}' - \mathbf{x}|^{\alpha} = r^{\prime \alpha} \sum_{\ell=0}^{+\infty} C_{\ell}^{-\alpha/2} (\mathbf{n} \cdot \mathbf{n}') \left(\frac{r}{r'}\right)^{\ell}. $$

<footnote n="5">⁵Equation (3.12) is equivalent to the more familiar (but less convenient for our purpose) expansion in terms of Gegengauer polynomials $C_\ell^\alpha(x)$,</footnote>
where we define the angular average over the $(d-1)$-dimensional sphere

$$ \hat{f}_{p,q}^{L}(t) = \int \frac{d\Omega_{d-1}}{\Omega_{d-1}} \hat{n}^{L} f_{p,q}(\mathbf{n}, t) \quad \text{with} \quad \Omega_{d-1} = \frac{2\pi^{d/2}}{\Gamma(\frac{d}{2})}. \qquad (3.15) $$

Note the nonlocal character of the relation (3.14) in $d$ dimensions, coming from the odd parity part of the PN retarded integral (3.9). Finally we see that only specific values of $q$ contribute: the $q=0$ terms induce an even sector, and the $q=1$ terms, an odd one. Finally, it is important to remark that the formula (3.14) cannot induce poles $\propto (d-3)^{-1}$ by itself. Therefore, we can safely work in the $d=3$ limit, as long as the coefficients $\hat{f}_{p,q}^L$ do not have poles.

### C. Application to the source mass quadrupole

Let us apply the formulas (3.8) and (3.14) to the case of the source mass quadrupole at 4PN order (the IR commutators do not affect the computation of the other required moments). The fact that Eq. (3.14) selects only terms with $q=0$ or $q=1$ drastically reduces the number of interactions that effectively can play a role at 4PN. To lowest orders at least, one can show that the even PN sector of the potentials bears $q \ge 1$, and the odd one, $q=0$. Therefore, the only possible interactions at lowest orders are either even-odd or odd-odd couplings, as even-even ones will have $q \ge 2$. In addition, the odd part of both the lowest-order potentials $V$ and $V_i$ (defined in the Appendix A of [91]) starts at 1.5PN (more precisely, the 0.5PN term of those two potentials is vanishing in the $d \to 3$ limit). Investigating the metric (B1), one can see that only two interactions can lead to corrections with 4PN accuracy: $V^2$ and $V\hat{W}$, where we denote $\hat{W} \equiv \hat{W}_{kk}$.<fnref n="6" /> The commutators affect the quantities $\Sigma$ and $\tilde{\mu}_1$, defined in Eqs. (2.10) and (2.17) of [91], as

$$ \Delta\Sigma = \frac{4(d-1)\sigma}{(d-2)^2 c^4} \left( C\{V^2\} + \frac{2(d-2)}{(d-1)c^2} C\{V\hat{W}\} \right) + \mathcal{O}\left(\frac{1}{c^{10}}\right), \qquad (3.16a) $$

$$ \Delta \tilde{\mu}_1 = \frac{2(d-4)}{(d-2)c^4} \left( C\{V^2\} + \frac{2(d-2)}{(d-1)c^2} C\{V\hat{W}\} \right)_1 m_1 + \mathcal{O}\left(\frac{1}{c^{10}}\right), \qquad (3.16b) $$

where $\sigma$ is the compact-supported expression defined in Eq. (2.15) of [91], $m_1$ the mass of particle 1, and where the

label 1 indicates that the quantities have to be evaluated and regularized at its location. This induces a correction to the source quadrupole

$$ \Delta I_{ij} = \frac{d(d-1)}{(d-2)^2 c^4} \int d^d \mathbf{x} \left( C\{V^2\} + \frac{2(d-2)}{(d-1)c^2} C\{V\hat{W}\} \right) \times \delta_1 m_1 \hat{y}_1^{ij} + (1 \leftrightarrow 2) + \mathcal{O}\left(\frac{1}{c^{10}}\right), \qquad (3.17) $$

where $\delta_1 = \delta^{(d)}[\mathbf{x} - \mathbf{y}_1]$ is the $d$-dimensional Dirac distribution at the position $\mathbf{y}_1$ of the particle 1.

The commutator of the interaction $V^2$ is quite easy to compute. Indeed, the potential $V$ is known in the whole space to a high PN order, well beyond the required precision. One can straightforwardly expand it as in Eq. (3.11) and read the coefficients $\hat{f}_{p,q}^L$. As discussed previously, the 0.5PN term of $V$ vanishes in the $d=3$ limit. Thus, only the (even, 0PN)–(odd, 1.5PN) coupling, having $q=1$, can enter at 4PN. Working it explicitly, we find

$$ C\{V^2\} = \frac{G^2 m}{c^4} \left[ 2\tilde{\mu}^{(2)} c^2 + \frac{1}{3} I^{(4)} \right] + \mathcal{O}(\epsilon), \qquad (3.18) $$

where $m = m_1 + m_2$ and $\tilde{\mu}(t) = \int d^3\mathbf{x}'\sigma(\mathbf{x}', t)$ reduces to a constant at leading order for arbitrary matter systems, so that $\tilde{\mu}(t) = m + \mathcal{O}(c^{-2})$ with $m$ constant, while $I = \int d^3\mathbf{x}'\sigma(\mathbf{x}', t)r'^2$ represents the effective moment of inertia. For compact binaries, those quantities read explicitly

$$ \begin{aligned} \tilde{\mu}(t) &= \tilde{\mu}_1(t) + \tilde{\mu}_2(t) \\ &= m_1 \left[ 1 + \frac{1}{c^2} \left( \frac{3}{2} v_1^2 - \frac{Gm_2}{r_{12}} \right) \right] \\ &\quad + (1 \leftrightarrow 2) + \mathcal{O}\left(\frac{1}{c^4}\right), \end{aligned} \qquad (3.19) $$

$$ I(t) = m_1 y_1^2 + m_2 y_2^2 + \mathcal{O}\left(\frac{1}{c^2}\right), \qquad (3.20) $$

with $v_1 = dy_1/dt$ and $r_{12} = |y_1 - y_2|$. On the other hand, the $d$-dimensional expression of the potential $\hat{W}$ is not known in the whole space (see Appendix C of [113] for more details on the computation). Using an integration by part, one can rewrite the wave equation it obeys, Eq. (A4d) of [91], as

$$ \begin{aligned} \Box \hat{W} &= \frac{8\pi G}{d-2} \sigma_{kk} - \frac{d-1}{2(d-2)} \partial_i V \partial_i V \\ &= \frac{8\pi G}{d-2} \sigma_{kk} - \frac{d-1}{4(d-2)} \left( \Box V^2 + \frac{2}{c^2} (\partial_t V)^2 + 8\pi G V \sigma \right), \end{aligned} \qquad (3.21) $$

where $\sigma_{kk}$ and $\sigma$ are the compact-supported expression defined in Eq. (2.15) of [91]. The compact-supported part of

<footnote n="6">⁶The near-zone potential $\hat{W}_{ij}$ and its trace $\hat{W}$ (see Appendix A of [91]) should not be confused with the $\ell$th gauge moment $W_\ell$ and W corresponding to $\ell=0$, as defined in (2.22)–(2.23).</footnote>
the source is easily integrated over the whole space, using methods described e.g., in [91]. It contains an even sector with $q = 1$ that cannot contribute (we recall that, up to 1PN, $V$ bears $q = 1$), and an odd, 0.5PN sector with $q = 0$ that will contribute, when coupled to the expansion of $V$. As for the noncompact sector, it reads

$$ \begin{aligned} \hat{W}^{\text{NC}} &= -\frac{d-1}{4(d-2)} \text{FP}\overline{\Box}_{\text{ret}}^{-1} \{\Box V^2 + \frac{2}{c^2}(\partial_t V)^2\} \\ &= -\frac{d-1}{4(d-2)} \left(V^2 + \frac{2}{c^2}(\partial_t V)^2 + C\{V^2\}\right). \end{aligned} \quad (3.22) $$

The first two terms will have $q = 2$ up to 1PN and 2PN orders, respectively, and we just proved that the third term is of order 2PN. Therefore, $\hat{W}^{\text{NC}}$ does not contribute to the IR commutator at the required order. Applying the formula (3.14), one finds the nice relation

$$ C\{V\hat{W}\} = -\frac{3c^2}{4}C\{V^2\} + \mathcal{O}(\epsilon) + \mathcal{O}\left(\frac{1}{c^4}\right). \quad (3.23) $$

With the two explicit expressions (3.18) and (3.23) at hand, one can compute the resulting correction to the source mass quadrupole $\Delta I_{ij}$ in (3.17). In the c.m. frame, it comes as

$$ \Delta I_{ij} = -\frac{4G^3 m^4 \nu^2}{r^3 c^8} \left[ \nu^2 - 3\dot{r}^2 - \frac{Gm}{r} \right] x_{\langle i}x_{j\rangle} + \mathcal{O}\left(\frac{1}{c^{10}}\right). \quad (3.24) $$

In the case of quasicircular orbits, $\nu^2 = Gm/r + \mathcal{O}(c^{-2})$ and $\dot{r} = \mathcal{O}(c^{-5})$, so this correction plays no role. Note also that the two corrections to the metric (3.18) and (3.23) are functions of the time only, and so cannot contribute to the equations of motion.

## IV. SOURCE MULTIPOLE MOMENTS FOR CIRCULAR ORBITS

Having now expressed the radiative moments in terms of the canonical moments in Sec. II A, and the canonical moments in terms of the source and gauge moments in Sec. II C, we now review the explicit expression of the source moments for compact binary systems. Note that the gauge moments [entering the relation (2.39) as well as Sec. III. B of [95]] are required only at Newtonian order, excepted for W, needed at 1PN. They are easily computed using Eq. (125) of [24], and their expressions on quasicircular orbits are displayed in Eq. (5.7) of [17]. The most crucial moment is the source quadrupole moment $I_{ij}$, needed only at 4PN order. Indeed, its 4.5PN terms, due to 2PN relative radiation reaction corrections, cannot involve nonlocal integrals, and will disappear from the flux.

### A. Source mass quadrupole at the 4PN order

A first computation, using the Hadamard regularization for the IR divergences and dimensional regularization for the UV ones, was performed in [91]. However, in order to be consistent with the equations of motion [79–82] which used dimensional regularization both in the UV and IR sectors, it should be computed with dimensional regularization also for the IR divergences. It was thus completed in [92], and the expression of the source quadrupole moment, obtained with full dimensional regularization, has the expected feature of exhibiting poles in $\epsilon = d-3$. As reviewed in Sec. II B, those poles are crucial to cancel the divergences linked with the $d$-dimensional computation of the radiative moments, obtained in [93]. Indeed, the source moments are not observables per se, but the radiative moments are. Therefore, only the latter have to be finite in the $\epsilon \to 0$ limit. However, as already mentioned, for the sake of computational simplicity, it was deemed possible to introduce the notion of a “renormalized” source quadrupole, defined by Eq. (6.2) of [93]. This renormalized quantity is nothing but the sum of the $d$-dimensional source quadrupole and the corrections (2.18a), arising from the radiative/canonical relation [we recall that we consider the corrections (2.37) as being part of the canonical/source relation]. The renormalized quadrupole, by construction, has no poles in $\epsilon$, and when injected into the three-dimensional MPM algorithm, yields the correct radiative quadrupole by definition. Last, but not least, up to 3.5PN order, the corrections due to the IR dimensional regularization exactly cancel the ones due to the renormalization, even out of the c.m. frame, as proven in [93]. Thus, up to 3.5PN order, the renormalized source quadrupole coincides with the one computed with Hadamard regularization in the IR, which is why those subtleties did not hit previous computations of the flux. However, this equivalence is no more valid at 4PN order, thus the crucial need for the implementation of dimensional regularization and the “renormalization” described above.

As the correction due to the IR commutators discussed in Sec. III vanishes on quasicircular orbits, the expression of the renormalized source quadrupole on quasicircular orbits remains identical to the one displayed in Eq. (6.11) of [93], namely

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

where, introducing the PN parameter $\gamma \equiv \frac{Gm}{rc^2}$, the coefficients are given by
$$
\begin{equation}
\begin{aligned}
A ={}& 1 + \gamma \left(-\frac{1}{42} - \frac{13}{14}\nu\right) + \gamma^2 \left(-\frac{461}{1512} - \frac{18395}{1512}\nu - \frac{241}{1512}\nu^2\right) \\
& + \gamma^3 \left(\frac{395899}{13200} - \frac{428}{105}\ln\left(\frac{r}{r_0}\right) + \left[\frac{3304319}{166320} - \frac{44}{3}\ln\left(\frac{r}{r_0}\right)\right]\nu + \frac{162539}{16632}\nu^2 + \frac{2351}{33264}\nu^3\right) \\
& + \gamma^4 \left(-\frac{1067041075909}{12713500800} + \frac{31886}{2205}\ln\left(\frac{r}{r_0}\right)\right. \\
& \qquad \left. + \left[-\frac{85244498897}{470870400} - \frac{2783}{1792}\pi^2 - \frac{64}{7}\ln(16\gamma e^{2\gamma_E}) - \frac{10886}{735}\ln\left(\frac{r}{r_0}\right) + \frac{8495}{63}\ln\left(\frac{r}{r_0}\right)\right]\nu\right. \\
& \qquad \left. + \left[\frac{171906563}{4484480} + \frac{44909}{2688}\pi^2 - \frac{4897}{21}\ln\left(\frac{r}{r_0}\right)\right]\nu^2 - \frac{22063949}{5189184}\nu^3 + \frac{71131}{314496}\nu^4\right),
\end{aligned}
\tag{4.2a}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
B ={}& \frac{11}{21} - \frac{11}{7}\nu + \gamma \left(\frac{1607}{378} - \frac{1681}{378}\nu + \frac{229}{378}\nu^2\right) + \gamma^2 \left(-\frac{357761}{19800} + \frac{428}{105}\ln\left(\frac{r}{r_0}\right) - \frac{92339}{5544}\nu + \frac{35759}{924}\nu^2 + \frac{457}{5544}\nu^3\right) \\
& + \gamma^3 \left(\frac{23006898527}{1589187600} - \frac{4922}{2205}\ln\left(\frac{r}{r_0}\right) + \left[\frac{8431514969}{529729200} + \frac{143}{192}\pi^2 - \frac{32}{7}\ln(16\gamma e^{2\gamma_E}) - \frac{1266}{49}\ln\left(\frac{r}{r_0}\right) - \frac{968}{63}\ln\left(\frac{r}{r_0}\right)\right]\nu\right. \\
& \qquad \left. + \left[\frac{351838141}{5045040} - \frac{41}{24}\pi^2 + \frac{968}{21}\ln\left(\frac{r}{r_0}\right)\right]\nu^2 - \frac{1774615}{81081}\nu^3 - \frac{3053}{432432}\nu^4\right),
\end{aligned}
\tag{4.2b}
\end{equation}
$$

$$
C = \frac{48}{7} + \gamma \left( -\frac{4096}{315} - \frac{24512}{945}\nu \right) - \frac{32}{7}\pi\gamma^{3/2}. \qquad (4.2c)
$$

The coefficients A and B represent the conservative part of the quadrupole, while C is due to the radiation reaction dissipative effects. Note that, in addition to the already-encountered scale $r_0$, the expression of the quadrupole involves a scale $r'_0$, associated with the UV regularization (see the footnote 10 of [91] for more details).

As discussed in [92], we found that the source quadrupole moment is not a local quantity at the 4PN order anymore, as it contains a nonlocal tail integral, given by Eq. (6.5) in [93]. The 4PN logarithms ln(16γe2γE) in A and B are due to the conservative part of this nonlocal tail term in the mass quadrupole, and the coefficient -32/7 πγ3/2 in C, to the corresponding dissipative part of the tail term.

Finally, as reported in Appendix A, the general expression of the source mass quadrupole moment for generic orbits has been conclusively tested in the so-called boosted Schwarzschild black hole limit.

## B. Higher-order source moments

As expected from the behavior of the associated radiative moment (2.18b), poles arise when performing the IR dimensional regularization of the source mass octupole $I_{ijk}$. Nevertheless, the renormalized source mass octupole moment exactly coincides at 3PN order with the one computed with Hadamard regularization in the IR, see [93] for discussion. It reads [95]

$$
I_{ijk} = -\nu m \Delta \left\{ D x_{\langle i} x_{j} x_{k\rangle} + E \frac{r}{c} v_{\langle i} x_{j} x_{k\rangle} + F \frac{r^2}{c^2} v_{\langle i} v_{j} x_{k\rangle} + G \frac{r^3}{c^3} v_{\langle i} v_{j} v_{k\rangle} \right\} + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (4.3)
$$

where $\Delta \equiv (m_1 - m_2)/m$, and the coefficients are

$$
D = 1 - \gamma\nu + \gamma^2 \left(-\frac{139}{330} - \frac{11923}{660}\nu - \frac{29}{110}\nu^2\right) + \gamma^3 \left(\frac{1229440}{63063} + \frac{610499}{20020}\nu + \frac{319823}{17160}\nu^2 - \frac{101}{2340}\nu^3 - \frac{26}{7}\ln\left(\frac{r}{r_0}\right) - 22\nu\ln\left(\frac{r}{r'_0}\right)\right), \quad (4.4a)
$$

$$
E = \frac{196}{15} \gamma^2 \nu, \qquad (4.4b)
$$
$$F = 1 - 2\nu + \gamma \left( \frac{1066}{165} - \frac{1433}{330}\nu + \frac{21}{55}\nu^2 \right) + \gamma^2 \left( -\frac{1130201}{48510} - \frac{989}{33}\nu + \frac{20359}{330}\nu^2 - \frac{37}{198}\nu^3 + \frac{52}{7}\ln\left(\frac{r}{r_0}\right) \right), \quad (4.4c)$$

$$G = 0. \qquad (4.4d)$$

As for the renormalized current quadrupole $J_{ij}$, it also coincides at 3PN order with the one computed with Hadamard regularization in the IR [93]. It comes [94]

$$J_{ij} = -\nu m \Delta \left[ HL^{(i}x^{j)} + K \frac{Gm}{c^3} L^{(i}v^{j)} \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \qquad (4.5)$$

where we denote $L^i \equiv \epsilon_{ijk}x^j v^k$, and where

$$H = 1 + \gamma \left( \frac{67}{28} - \frac{2}{7}\nu \right) + \gamma^2 \left( \frac{13}{9} - \frac{4651}{252}\nu - \frac{\nu^2}{168} \right) + \gamma^3 \left( \frac{2301023}{415800} - \frac{214}{105}\ln\left(\frac{r}{r_0}\right) + \left[ -\frac{243853}{9240} + \frac{123}{128}\pi^2 - 22\ln\left(\frac{r}{r_0}\right) \right]\nu \right. \\ \left. + \frac{44995}{5544}\nu^2 + \frac{599}{16632}\nu^3 \right), \qquad (4.6a)$$

$$K = \frac{188}{35}\nu\gamma. \qquad (4.6b)$$

The remaining required source moments (mass hexadecapole at 2PN, current octupole at 2PN, etc.) are easy to compute as no regularization subtleties arise. Their expressions on quasicircular orbits are displayed e.g., in Sec. 9.1 of [24]. Note that, as the first nonlocal feature cannot appear at a lower order than 4PN, all those higher-order source moments are instantaneous up to 3.5PN order. Therefore, they cannot contribute to the 4.5PN term of the quasicircular flux.

## V. TOOLBOX OF INTEGRATION FORMULAS

This technical section presents some miscellaneous material that was used to perform the explicit computations of the gravitational wave flux to 4.5PN and the (2, 2) mode to 4PN, which will be reported in Sec. VI.

### A. 4PN equations of motion for circular orbits

First, the relation between the radiative and canonical moments, as well as the expression of the flux, involve

temporal derivatives. We thus need the expression of the equations of motion for quasicircular orbits at the 4PN order, including both conservative and dissipative terms. Note that the 4.5PN piece of the equations of motion [118,119] is local and so, as already discussed, will not contribute to the 4.5PN flux. Recalling that $\gamma = \frac{Gm}{rc^2}$, the acceleration for quasicircular orbit reads [82]

$$\frac{dv^i}{dt} = -\omega^2 x^i - \frac{32 G^3 m^3 \nu}{5 r^4 c^5} \left[ 1 + \left( -\frac{743}{336} - \frac{11}{4} \nu \right) \gamma + 4\pi\gamma^{3/2} + \mathcal{O}(\gamma^2) \right] v^i. \qquad (5.1)$$

The radiation-reaction part (proportional to $v^i$) includes 1PN and 1.5PN corrections beyond the leading 2.5PN effect. In particular the 1.5PN correction is due to tails. The conservative part of the equations of motion is specified for circular orbits by the orbital frequency, which is a generalization of Kepler's law valid through 4PN,

$$\omega^2 = \frac{Gm}{r^3} \left\{ \begin{aligned} & 1 + (-3 + \nu)\gamma + \left(6 + \frac{41}{4}\nu + \nu^2\right)\gamma^2 + \left(-10 + \left[-\frac{75707}{840} + \frac{41}{64}\pi^2 + 22\ln\left(\frac{r}{r_0}\right)\right]\nu + \frac{19}{2}\nu^2 + \nu^3\right)\gamma^3 \\ & + \left(15 + \left[\frac{19644217}{33600} + \frac{163}{1024}\pi^2 + \frac{256}{5}\gamma_E + \frac{128}{5}\ln(16\gamma) - 290\ln\left(\frac{r}{r_0}\right)\right]\nu \right. \\ & \left. + \left[\frac{44329}{336} - \frac{1907}{64}\pi^2 + \frac{1168}{3}\ln\left(\frac{r}{r_0}\right)\right]\nu^2 + \frac{51}{4}\nu^3 + \nu^4\right)\gamma^4 + \mathcal{O}(\gamma^5) \end{aligned} \right\}. \qquad (5.2)$$

We recall that the scales $r_0$ and $r'_0$ are respectively associated with IR and UV regularizations. Using the flux-balance equation as well as the generalized Kepler's law, we find the secular decay of orbital quantities at 1.5PN order
$$ \dot{r} \equiv \frac{dr}{dt} = -\frac{64 G^3 m^3 \nu}{5 r^3 c^5} \left[ 1 + \left( -\frac{1751}{336} - \frac{7}{4}\nu \right) \gamma + 4\pi\gamma^{3/2} + \mathcal{O}(\gamma^2) \right], \qquad (5.3a) $$

$$ \dot{\omega} \equiv \frac{d\omega}{dt} = \frac{96 G m \nu}{5 r^3} \gamma^{5/2} \left[ 1 + \left( -\frac{2591}{336} - \frac{11}{12}\nu \right) \gamma + 4\pi\gamma^{3/2} + \mathcal{O}(\gamma^2) \right]. \qquad (5.3b) $$

The relative 2PN order is also known [118,119], but is not required for the present work. In the following, it will be important to distinguish between the orbital frequency $\omega$ defined from the equations of motion by (5.2), and the GW half-frequency $\Omega$ which is measured in the waves propagating in the far zone. We pose

$$ y \equiv \left( \frac{Gm\omega}{c^3} \right)^{2/3}, \qquad (5.4) $$

and distinguish it later from the PN parameter $x$ defined from the GW half-frequency $\Omega$, see Eq. (6.9).

## B. Postadiabatic integration of the tail effect

In order to compute the tail integrals, for instance (2.5), one needs to specify the orbit's behavior of the compact binary system in the remote past, as the effect is not localized in time, but integrates over the whole past history of the source. In the case of quasicircular orbits, and up to 3.5PN precision in the multipoles, an adiabatic approximation (considering the orbital elements $r$ and $\omega$ to be constant in time) is sufficient, and one can follow the lines of [120], together with the integrals presented in the Appendix B of [121]. However, Eq. (5.3) show that this adiabatic approximation is no longer valid at a relative 2.5PN precision. As the tail enters at 1.5PN order in the moments, the first “postadiabatic” (PA) correction will affect the moments at the 4PN order, thus we need to properly evaluate it in order to consistently derive the 4PN flux and (2, 2) mode.

As shown in [120], the tail integrals on quasicircular orbits reduce to elementary integrals of the type

$$ \mathcal{I}_{a,n}(u) = \int_{0}^{+\infty} d\tau [y(u-\tau)]^a e^{-in\phi(u-\tau)} \ln\left(\frac{\tau}{\tau_0}\right), \qquad (5.5) $$

where $n \in \mathbb{N}^*$ (the case $n = 0$ does not appear at 4PN), $a$ is usually a rational fraction, and the constant $\tau_0$ denotes either $2r_0/c$ or $2b_0/c$. The orbital phase $\phi$, the orbital frequency $\omega = \dot{\phi}$ and the parameter $y = \left(\frac{Gm\omega}{c^3}\right)^{2/3}$ are integrated over any instant $u - \tau$ in the past. Our strategy will be to notice that the integral $\mathcal{I}_{a,n}$ involves a fast-oscillating exponential, that is derivable under properly defined PA approximations. This section provides a general method valid an arbitrary PA order, which we apply to first order 1PA.

We first remark that the difference between the current phase $\phi(u)$ and the phase $\phi(u - \tau)$ in the past is of the order of the inverse of the radiation reaction scale, i.e., $\phi(u) - \phi(u - \tau)$ scales as $\mathcal{O}(c^5)$, which is also obvious from Eq. (317) of [24]. To describe the radiation-reaction scale, we introduce a dimensionless adiabatic parameter at the current time $u$, denoted by$^7$

$$ \xi(u) \equiv \frac{\dot{\omega}(u)}{\omega^2(u)} = \mathcal{O}\left(\frac{1}{c^5}\right). \qquad (5.6) $$

Using the energy balance equation we have $\xi = -\mathcal{F}(\omega)/[\omega^2 dE/d\omega]$. We want to compute the integral (5.5) in the PA limit where $\xi(u) \to 0$. To this end, let us pose

$$ \phi(u) - \phi(u - \tau) = \frac{v}{\xi(u)}, \qquad (5.7) $$

and change the integration variable from $\tau$ to $v$, so that the integral (5.5) becomes

$$ \mathcal{I}_{a,n} = \frac{Gm e^{-in\phi(u)}}{c^3 \xi(u)} \int_{0}^{+\infty} dv [y(u-\tau(v))]^{a-3/2} \ln\left(\frac{\tau(v)}{\tau_0}\right) e^{\frac{in\nu}{\xi(u)}}. \qquad (5.8) $$

Here $\tau(v)$ is obtained by inverting Eq. (5.7) for $\tau$ as a function of $v$. The main point is that, in the PA limit, $\xi(u) \to 0$, and so, as the imaginary exponential in the integrand of (5.8) oscillates very rapidly, the integral is a sum of alternatively positive and negative contributions which essentially sum up to zero. However there are no oscillations when $v = 0$, therefore the integral is essentially given by the contribution in the neighborhood of the bound at $v = 0$. In other words we are entitled to compute the integral by replacing the integrand by its formal expansion series when $v \to 0$. This will give a formal (asymptotic) series in powers of $\xi(u)$ which is the “PA” expansion of the integral.

In the PA limit we thus have $\tau \to 0$, and we can expand Eq. (5.7) to any order as

<footnote>$^7$This definition is intended to be valid at any PN order, and is different from the one adopted in previous work; see (4.8) in [120].</footnote>
$$v = \xi \sum_{n=0}^{+\infty} \frac{(-)^n}{(n+1)!} \omega^{(n)} \tau^{n+1}. \quad (5.9)$$

From now on, we pose $\xi = \xi(u)$, $\phi = \phi(u)$, $\omega = \omega(u) = \dot{\phi}(u)$, $\omega^{(n)} = (d^n\omega/du^n)(u)$ for the quantities at the current time $u$. Since $\omega^{(1)} \equiv \dot{\omega} = \xi\omega^2$ the expansion (5.9) clearly represents the PA approximation in powers of $\xi$ and arbitrary high time derivatives of $\xi$ denoted $\xi^{(n)}$. Note that each time derivative of $\xi$ adds a radiation reaction scale, hence we have $\xi^{(n)} = \mathcal{O}(\xi^{n+1}) = \mathcal{O}(c^{-5n-5})$.

To obtain the PA expansion of (5.8) we need to invert the series (5.9) and obtain $\tau$ as a power series in $v$. This is given by the Lagrange inversion theorem as

$$\tau = \sum_{p=0}^{+\infty} \frac{1}{(p+1)!} f_p \left(\frac{v}{\xi}\right)^{p+1}, \quad (5.10)$$

where the general coefficients $f_p$ read

$$\begin{aligned} f_p &= \left(\frac{d}{d\tau}\right)^p \left[ \left(\frac{\tau}{\phi(u) - \phi(u-\tau)}\right)^{p+1} \right]_{\tau=0} \\ &= \left(\frac{d}{d\tau}\right)^p \left[ \left(\sum_{n=0}^{+\infty} \frac{(-)^n}{(n+1)!} \omega^{(n)} \tau^n \right)^{-p-1} \right]_{\tau=0}, \quad (5.11) \end{aligned}$$

where $\tau = 0$ is applied after the $p$ differentiations with respect to $\tau$. For instance, up to 2PA order, we obtain

$$\tau = \frac{v}{\xi\omega} \left[ 1 + \frac{v}{2} + \left(1 - \frac{\dot{\xi}}{\xi^2\omega}\right) \frac{v^2}{6} + \mathcal{O}(v^3) \right]. \quad (5.12)$$

The previous formulas show that the method can be straightforwardly extended to any order. But as we said, to compute the tail term at 4PN order, we need only the correction of order 2.5PN, which corresponds to the 1PA approximation, hence just $\tau = \frac{v}{\xi\omega} [1 + \frac{v}{2} + \mathcal{O}(v^2)]$. To 1PA order, the integral (5.8) reads

$$\mathcal{I}_{a,n} = \frac{Gm}{c^3} \frac{y^{a-3/2} e^{-in\phi}}{\xi} \int_0^{+\infty} dv \left[ \left(1 + \left(1 - \frac{2a}{3}\right)v\right) \right. \\ \left. \times \ln\left(\frac{v}{\xi\omega\tau_0}\right) + \frac{v}{2} + \mathcal{O}(v^2) \right] e^{\frac{inv}{\xi}}. \quad (5.13)$$

The remaining integral is computed as follows. We transform the complex exponential into a real one by performing the change of variable $v = iw$ for $n > 0$, and $v = -iw$ for $n < 0$, respectively. The integration now takes place along the imaginary axis, which can be remedied by resorting to the Cauchy theorem on the closed contour made of the three following pieces, to be considered in the limit $R \to +\infty$: (i) the path from $w = 0$ to $w = R$ on the real axis; (ii) the oriented quarter of circle of radius $|w| = R$

from arg $w = 0$ to arg $w = -\pi/2$ (or arg $w = \pi/2$ if $n < 0$); and (iii) the segment of the imaginary axis going from $w = -iR$ ($w = iR$ if $n < 0$) to $w = 0$. This leads to

$$\mathcal{I}_{a,n} = \frac{Gm y^{a-3/2} e^{-in\phi}}{c^3 i\xi} \int_0^{+\infty} dw \left[ \left(-s(n) + \frac{2a-3}{3}iw\right) \right. \\ \left. \times \left( \ln\left(\frac{w}{\xi\omega\tau_0}\right) + i\frac{\pi}{2}s(n) \right) + \frac{iw}{2} + \mathcal{O}(w^2) \right] e^{-\frac{|n|w}{\xi}}, \quad (5.14)$$

where $s(n)$ is the sign function. In this form, $\mathcal{I}_{a,n}$ can be integrated, as it boils down to elementary integrals. For completeness we give the formulas needed to handle any PA approximation,

$$\begin{aligned} \int_0^{+\infty} dw \, w^j e^{-\frac{|n|w}{\xi}} &= j! \left(\frac{\xi}{|n|}\right)^{j+1}, \quad (5.15a) \\ \int_0^{+\infty} dw \, w^j \ln\left(\frac{w}{\xi\omega\tau_0}\right) e^{-\frac{|n|w}{\xi}} &= j! \left(\frac{\xi}{|n|}\right)^{j+1} [H_j - \gamma_E - \ln(|n|\omega\tau_0)], \quad (5.15b) \end{aligned}$$

where $H_j = \sum_{k=1}^j \frac{1}{k}$ denotes the harmonic number and $\gamma_E$ is the Euler constant. Finally we obtain, at desired 1PA order, the result

$$\mathcal{I}_{a,n} = \frac{Gm y^{a-3/2} e^{-in\phi}}{c^3 in} \left\{ \left(1 + \frac{2a-3}{3} \frac{\xi}{in}\right) \right. \\ \left. \times \left[ \ln(|n|\omega\tau_0) + \gamma_E - i \frac{\pi}{2} s(n) \right] \right. \\ \left. - \frac{4a-9}{6} \frac{\xi}{in} + \mathcal{O}(\xi^2) \right\}. \quad (5.16)$$

Recall that all quantities are evaluated at current time $u$ and that the adiabatic parameter $\xi = \dot{\omega}/\omega^2$ is easily computed with Eq. (5.3). With this result at hand, we are able to derive the tail integral (2.5) with 2.5PN relative precision, which impacts the computation and final results at the 4PN order.

### C. Memory effects for circular orbits up to 1.5PN relative order

Memory terms, such as the first term of Eq. (2.6), are hereditary integrals of the form $\int_0^{+\infty} d\tau F(u - \tau)G(u - \tau)$, where $F$ and $G$ represent dynamical multipole moments. Note that they enter only in mass-type moments, as clear from Sec. II A. In the case of quasicircular orbits, they reduce to a sum of terms of the form

$$\mathcal{J}_{a,n}(u) = \int_{0}^{+\infty} d\tau [y(u - \tau)]^{a} e^{-in\phi(u - \tau)}. \quad (5.17)$$
Besides the absence of logarithm, the main difference with the tail integral is the possibility of having $n = 0$, i.e., persistent or “DC” terms.

Let us first focus on oscillatory (“AC”) memory terms, having $n \neq 0$. As they only involve a simple, logarithmic-free, integration, memory terms will enter in the flux as instantaneous contributions, and can be computed as such. In particular, since they arise at the odd PN approximations 2.5PN and 3.5PN, they do not contribute to the flux for quasicircular orbits. The evaluation of the integrals of type (5.17) is thus only required for the derivation of the modes. Concerning the quadrupole moment, since the memory effect enters at 2.5PN order, as clear from (2.6), it is thus required at a relative 1.5PN precision for the mode (2, 2) and thus, one can safely compute the integral in the adiabatic approximation. As discussed in [120], for circular orbits this is equivalent to taking $y$ to be constant together with a linear phase (appropriate for the exact circular orbit), and keeping only the contribution of the integral due to the bound $\tau = 0$. One finds

$$ \mathcal{J}_{a,n} = y^a e^{-in\phi} \int_0 d\tau e^{in\omega\tau} + \mathcal{O}(\xi) = -\frac{y^a e^{-in\phi}}{in\omega} + \mathcal{O}(\xi). \quad (5.18) $$

The persistent DC terms obviously do not contribute to the flux, and they only contribute to the modes which have $m = 0$, for example, the $(\ell, m) = (2, 0)$ mode for the mass quadrupole. As is clear in the following, the absence of the fast oscillating exponential in $\mathcal{J}_{a,0}$ generates an inverse power of the adiabatic parameter $\xi$, thus degrading the precision by the radiation-reaction scale 2.5PN. This is the well-known memory effect; as it starts at 2.5PN order in the waveform, it finally enters the $(2, 0)$ mode at Newtonian order. Several methods are possible to evaluate $\mathcal{J}_{a,0}$; see e.g., [120]. In the following, we rely on a change of integration variables from $\tau$ to $y' = y(u - \tau)$:

$$ \mathcal{J}_{a,0} = \int_{0}^{+\infty} d\tau [y(u-\tau)]^a = \int_{0}^{y(u)} dy' \frac{y'^a}{\dot{y}(u-\tau)}, \quad (5.19) $$

reading $\dot{y}(u - \tau)$ as a function of $y'$ from Eq. (5.3), and supposing that $a > 4$, as is the case in practical computations. To 1.5PN relative order, we have

$$ \mathcal{J}_{a,0} = \frac{5Gm}{64c^3\nu} \frac{y^{a-4}}{a-4} \left[ 1 + \frac{a-4}{a-3} \left( \frac{743}{336} + \frac{11}{4}\nu \right)y - 8\pi \frac{a-4}{2a-5}y^{3/2} + \mathcal{O}(y^2) \right]. \quad (5.20) $$

It is interesting to observe that the tail effect in the flux directly influences the DC memory in (5.20).

Last, but not least, we need to look at the interesting case of the tails-of-memory terms (2.9). These terms can be

treated as standard tail terms (with relative Newtonian accuracy), except for the “genuine” tail-of-memory given by the first line of (2.9), namely

$$ \mathcal{K}_{ij} = \frac{8G^2 M}{7c^8} \int_0^{+\infty} d\rho M_{a\langle i \rangle}^{(4)}(u-\rho) \times \int_0^{+\infty} d\tau M_{j\rangle a}^{(4)}(u-\rho-\tau) \ln\left(\frac{\tau}{\tau_0}\right), \quad (5.21) $$

where $c\tau_0 = 2r_0e^{\frac{1613}{270}}$. We remind that this expression agrees with the tail-of-memory directly computed from the radiative quadrupole moment at infinity [110,111]. We first perform the tail-like integral over $\tau$. As we need to evaluate it at relative Newtonian order only, we can safely use the adiabatic approximation. Next we perform the integral over $\tau$, which is found to be a simple DC memory integral of the type (5.19) with $a = 13/2$. Hence we find

$$ \begin{aligned} \mathcal{K}_{ij} &= \frac{128\pi \nu^2 c^5}{7 G} \ell_{\langle i \rangle \ell_{\langle j \rangle}} \int_0^{+\infty} d\rho [y(u-\rho)]^{13/2} \\ &= \frac{4\pi}{7} m\nu c^2 y^{5/2} \ell_{\langle i \rangle \ell_{\langle j \rangle}}, \end{aligned} \quad (5.22) $$

where $\ell^i = L^i/|\mathbf{L}|$ is the constant unit vector associated with the Newtonian angular momentum and, thus, orthogonal to the orbital plane. Interestingly, this tail-of-memory result will give a contribution in the $(2, 0)$ mode that exactly cancels the one coming from the 1.5PN corrections to the ordinary memory effect, obtained in (5.20).

## VI. RESULTS

Collecting all the pieces that were discussed in previous sections, and using notably the integration techniques developed in Sec. V, we obtain our main results, namely the gravitational flux and quadrupole modes. All results displayed hereafter can be found in the ancillary file [102].

### A. Flux at 4PN order for generic orbits

For generic orbits, and in the c.m. frame, we split the gravitational energy flux (2.1) as

$$ \mathcal{F} = \mathcal{F}_{\text{can}} + \mathcal{F}_{\text{non-lin}}. \quad (6.1) $$

The first piece, $\mathcal{F}_{\text{can}}$, is the contribution due to the canonical mass and current moments, when they are fully expressed in terms of the compact binary parameters, i.e.,

$$ \mathcal{F}_{\text{can}} = \sum_{\ell \ge 2} \frac{G}{c^{2\ell+1}} [a_{\ell} M_{L}^{(\ell+1)} M_{L}^{(\ell+1)} + \frac{b_{\ell}}{c^2} S_{L}^{(\ell+1)} S_{L}^{(\ell+1)}], \quad (6.2) $$
where we recall that the numerical coefficients $a_\ell$ and $b_\ell$ are given in Eq. (2.2). We further split $\mathcal{F}_{\text{can}}$ into a local-in-time (or instantaneous) part and a nonlocal one,

$$ \mathcal{F}_{\text{can}} = \mathcal{F}_{\text{can}}^{\text{loc}} + \mathcal{F}_{\text{can}}^{\text{non-loc}}. \qquad (6.3) $$

The local piece $\mathcal{F}_{\text{can}}^{\text{loc}}$ is too long to be displayed here but can be found in the ancillary file [102]. Note the presence of the scales associated with the regularization processes, $r_0$ and $r'_0$, which is expected, as the canonical part of the flux is not a

gauge invariant quantity. As for the nonlocal piece, $\mathcal{F}_{\text{can}}^{\text{non-loc}}$, it comes from two effects: (i) the direct contribution of the nonlocal part of the mass quadrupole moment, given by Eq. (6.5) in [93]; and (ii) the tail term in the 4PN equations of motion, which contributes when time differentiating the 4PN mass quadrupole moment. [The latter tail term is given by the first line of Eq. (4.4) in [82], while the second line is included in the local part.] The nonlocal piece reads (at the required order, we are dealing with Newtonian quantities, so we can identify source and canonical moments),

$$ \mathcal{F}_{\text{can}}^{\text{non-loc}} = \frac{48 G^3 M}{5 c^{13}} M_{ij}^{(3)} \left\{ \frac{1}{7} \frac{d^3}{du^3} \left[ M_{ik}(u) \int_0^{+\infty} d\tau \ln\left(\frac{c\tau}{2r_0}\right) M_{jk}^{(5)}(u-\tau) \right] - \frac{m\nu}{15} x_{ik} \int_0^{+\infty} d\tau \ln\left(\frac{c\tau}{2r}\right) M_{jk}^{(8)}(u-\tau) \right. \\ \left. - \frac{m\nu}{15} (3x_k v_i + x_i v_k) \int_0^{+\infty} d\tau \ln\left(\frac{c\tau}{2r}\right) M_{jk}^{(7)}(u-\tau) \right\} + \mathcal{O}\left(\frac{1}{c^{14}}\right). \qquad (6.4) $$

Note the difference between the logarithmic kernel of the first integral, bearing $r_0$, and the two other ones, bearing $r$. Finally the second piece of the flux (6.1), $\mathcal{F}_{\text{non-lin}}$, is given by all the nonlinearities in the GW propagation discussed in Sec. II A. It contains all the double products between canonical moments and radiative moments, plus the square of the tail terms and the double product between the 1.5PN and 2.5PN corrections to the quadrupole. We write it as

$$ \mathcal{F}_{\text{non-lin}} = \frac{G}{5c^5} \left[ 2M_{ij}^{(3)} U_{ij}^{1.5\text{PN}(1)} + 2M_{ij}^{(3)} U_{ij}^{2.5\text{PN}(1)} + 2M_{ij}^{(3)} U_{ij}^{3\text{PN}(1)} + U_{ij}^{1.5\text{PN}(1)} U_{ij}^{1.5\text{PN}(1)} + \dots \right] \\ + \frac{G}{189c^7} \left[ 2M_{ijk}^{(4)} U_{ijk}^{1.5\text{PN}(1)} + 2M_{ijk}^{(4)} U_{ijk}^{2.5\text{PN}(1)} + \dots \right] + \dots + \mathcal{O}\left(\frac{1}{c^{14}}\right), \qquad (6.5) $$

where the ellipsis can be completed in a straightforward way from the results in Sec. II A.

## B. Flux at 4.5PN order for quasicircular orbits

Introducing the unit separation $n = x/r$ between the particles, as well as the unit vector $\lambda = \ell \times n$ such that $(n, \lambda, \ell)$ forms a direct orthonormal triad, we may write the relative velocity as $v = r\omega\lambda + \dot{r}n$, where $\dot{r}$ is given by (5.3). Then, using (5.2), we can express the 4PN flux on quasicircular orbits as a function of the orbital frequency $\omega$, or equivalently the parameter $y$ defined by (5.4).

However we find that the 4PN flux parametrized by $y$ still contains the unphysical constant $b_0$, although the other arbitrary scales $r_0$ and $r'_0$ have properly disappeared. The reason is that, starting at the 4PN order, the frequency is modified due to the propagation of tails in the wave zone. The half-phase $\psi$ of the reduced GW, i.e., half the phase of the (2, 2) mode, differs from the orbital phase $\phi$ (such that the orbital frequency is $\omega = \dot{\phi}$) by a logarithmic, tail-induced phase modulation, as

$$ \psi = \phi - \frac{2GM\omega}{c^3} \ln\left(\frac{\omega}{\omega_0}\right), \qquad (6.6) $$

where M denotes the constant ADM mass, and $\omega_0$ is related to $b_0$ by $c\omega_0^{-1} = 4b_0e^{r_E-11/12}$. We remind the reader that the constant $b_0$, which appears in the tail terms (2.5), (2.7) and (2.9), is arbitrary and parametrizes the logarithmic deviation of light cones in harmonic coordinates from the flat cones, see e.g., (2.10)–(2.11) in [99]. The constant $r_0$ is an IR scale introduced into the definition of the source multipole moments, see Eq. (2.1) in [91], while $r'_0$ is a UV scale associated with the regularization of the self-field of point particles.

This phase modulation was determined in [121,122] and it has been repeatedly argued [120,123] that it affects the waveform at the 4PN order. Now, the phase modulation (6.6) also shifts the GW half-frequency $\Omega = \dot{\psi}$ with respect to the orbital one $\omega = \dot{\phi}$. The GW half-frequency, directly measurable from the waveform at infinity, reads

$$ \Omega = \omega - \frac{2GM\dot{\omega}}{c^3} \left[ \ln\left(\frac{\omega}{\omega_0}\right) + 1 \right]. \qquad (6.7) $$

Using Eq. (5.3) for the frequency chirp at the dominant order, it explicitly comes, replacing M by $m$ at that order,
$$ \Omega = \omega \left\{ 1 - \frac{192}{5} \nu \left( \frac{Gm\omega}{c^3} \right)^{8/3} \left[ \ln\left(\frac{\omega}{\omega_0}\right) + 1 \right] + \mathcal{O}\left(\frac{1}{c^{10}}\right) \right\}, \quad (6.8) $$

where we recall that $\nu = \frac{m_1 m_2}{m^2}$ is the symmetric mass ratio. In the following, the results are expressed for quasicircular orbits in terms of the measurable GW half-frequency $\Omega$ through the PN parameter

$$ x = \left( \frac{Gm\Omega}{c^3} \right)^{2/3}. \quad (6.9) $$

Recalling also the definition (5.4) and posing $y_0 = (\frac{Gm\omega_0}{c^3})^{2/3}$ we thus have

$$ x = y \left\{ 1 - \frac{192}{5} \nu y^4 \left[ \ln\left(\frac{y}{y_0}\right) + \frac{2}{3} \right] + \mathcal{O}(y^5) \right\}. \quad (6.10) $$

showing that indeed the GW half-frequency $\Omega$ differs from the orbital one $\omega$ at the 4PN order only, hence the fact that it did not affect previous computations such as in [120,123]. Therefore, once the 4PN flux is obtained in terms of $y$, we replace $y$ in terms of $x$ using the inverse of Eq. (6.10) and reexpand to 4PN order. With this procedure, the constant $b_0$ cancels out as expected. Finally, adding also the 4.5PN piece [96], the quasicircular 4.5PN flux is

$$ \begin{aligned} \mathcal{F} = & \frac{32c^5}{5G} \nu^2 x^5 \left\{ 1 + \left( -\frac{1247}{336} - \frac{35}{12}\nu \right)x + 4\pi x^{3/2} + \left( -\frac{44711}{9072} + \frac{9271}{504}\nu + \frac{65}{18}\nu^2 \right)x^2 + \left( -\frac{8191}{672} - \frac{583}{24}\nu \right)\pi x^{5/2} \right. \\ & \left. + \left[ \frac{6643739519}{69854400} + \frac{16}{3}\pi^2 - \frac{1712}{105}\gamma_E - \frac{856}{105}\ln(16x) + \left( -\frac{134543}{7776} + \frac{41}{48}\pi^2 \right)\nu - \frac{94403}{3024}\nu^2 - \frac{775}{324}\nu^3 \right]x^3 \right. \\ & \left. + \left( -\frac{16285}{504} + \frac{214745}{1728}\nu + \frac{193385}{3024}\nu^2 \right)\pi x^{7/2} \right. \\ & \left. + \left[ -\frac{323105549467}{3178375200} + \frac{232597}{4410}\gamma_E - \frac{1369}{126}\pi^2 + \frac{39931}{294}\ln 2 - \frac{47385}{1568}\ln 3 + \frac{232597}{8820}\ln x \right. \right. \\ & \left. \left. + \left( -\frac{1452202403629}{1466942400} + \frac{41478}{245}\gamma_E - \frac{267127}{4608}\pi^2 + \frac{479062}{2205}\ln 2 + \frac{47385}{392}\ln 3 + \frac{20739}{245}\ln x \right)\nu \right. \right. \\ & \left. \left. + \left( \frac{1607125}{6804} - \frac{3157}{384}\pi^2 \right)\nu^2 + \frac{6875}{504}\nu^3 + \frac{5}{6}\nu^4 \right]x^4 \right. \\ & \left. + \left[ \frac{265978667519}{745113600} - \frac{6848}{105}\gamma_E - \frac{3424}{105}\ln(16x) + \left( \frac{2062241}{22176} + \frac{41}{12}\pi^2 \right)\nu \right. \right. \\ & \left. \left. - \frac{133112905}{290304}\nu^2 - \frac{3719141}{38016}\nu^3 \right]\pi x^{9/2} + \mathcal{O}(x^5) \right\}. \end{aligned} \quad (6.11) $$

A significant check is to observe that the leading-order terms in the test-mass limit $\nu \to 0$ perfectly agree with the results of linear black-hole perturbation theory [124–128]. Note that the flux has been confirmed by different groups up to 2PN order [19,23], and that the 4.5PN piece is in agreement with the independent work [100]. All other terms are new with the present paper. We have also explicitly verified that, at the 4PN order, Eq. (6.11) can be recovered from the gravitational modes given by Eq. (6.17) below and in [18], using

$$ \mathcal{F} = \frac{c^3}{16\pi G} \sum_{\ell \ge 0} \sum_{m=-\ell}^{\ell} |\dot{h}_{\ell m}|^2. \quad (6.12) $$

As usual, the contributions due to the absorption by the black-hole horizons should be added separately from the post-Newtonian result (6.11); see [129–134].

In the companion Letter [101] we use the quasicircular flux $\mathcal{F}(x)$ to derive the gravitational frequency chirp and phase at the 4.5PN order. Namely, we apply the energy flux-balance law

$$ \dot{E}(x) = -\mathcal{F}(x), \quad (6.13) $$

where $x$ is directly linked to the observed GW half-frequency through the definition (6.9), and is related to the orbital frequency via the relation (6.10). Here, $E(x)$ denotes the (Noetherian) binding energy of the compact binary on the quasicircular orbit, as calculated from the 4PN equations of motion, see e.g., [82]. The fact that the left-hand side of the balance equation, which concerns the motion, should be expressed in terms of the same observed GW half-frequency $x$ as the right-hand side, which concerns the radiation, and not, for instance, in
terms of the “orbital” frequency y, is worth an explanation; suppose that the compact binary system is actually a binary pulsar system. Hence, in addition to the gravitational waves generated by the orbital motion, the pulsar emits electro-magnetic (radio) waves, also received by the far away observer. Now the observer at infinity can measure the orbital frequency of the system from the instants of arrival of the radio pulses—this is the standard analysis of binary pulsars. Such frequency should be the one to be inserted in the left side of the balance equation. However, far from the system, the space-time curvature $R^{-1} \sim \sqrt{M/r^3}$ tends to zero, and therefore the geometric optics or WKB approximation applies for both the EM and gravitational waves. Thus the EM and gravitational waves follow the same geodesic, independently of their frequency, and in particular are subject to the same tail-induced phase modulation (6.6). We conclude that the distant observer measures the same frequency x from the EM radio pulses and from the gravitational wave, and this is that frequency that he inserts into both sides of the flux-balance law (6.13).<fnref n="8" />

### C. Gravitational wave modes (2, 2) and (2, 0)

Next, we can extract the $(\ell, m) = (2, 2)$ and $(2, 0)$ physical modes from the newly computed 4PN mass quadrupole radiative moment. Projecting the asymptotic metric (2.3) onto the basis of polarizations $\{+, \times\}$ and the usual basis of spin-weighted spherical harmonics, $Y_{-2}^{\ell m}$

(following the conventions of [16,17]), one can define the observable gravitational modes $h_{\ell m}$ as

$$h_{+} - ih_{\times} = \sum_{\ell=2}^{+\infty} \sum_{m=-\ell}^{\ell} h_{\ell m} Y_{-2}^{\ell m}. \quad (6.14)$$

For the sake of simplicity, we single out the (observable) GW half-phase $\psi$, and rescale the modes by the dominant contribution, defining $H_{\ell m}$ as

$$h_{\ell m} = \frac{8Gm\nu x}{c^2 R} \sqrt{\frac{\pi}{5}} H_{\ell m} e^{-im\psi}. \quad (6.15)$$

We have $H_{\ell,-m} = (-)^{\ell} \bar{H}_{\ell m}$, where the overbar denotes the complex conjugate. Using the PN-MPM framework, the current state-of-the-art is the 3.5PN accuracy for all the modes [17,18,136]. From the present completion of the 4PN mass quadrupole $U_{ij}$, one can extract the 4PN correction to the $(\ell, m) = (2, 2)$ mode as

$$H_{22} = - \frac{e^{2i\psi}}{2c^2m\nu x} \bar{m}^i \bar{m}^j U_{ij}, \quad (6.16)$$

where $m = (n + i\lambda)/\sqrt{2}$, with $\bar{m}$ its complex conjugate, and the radiative moment $U_{ij}$ is evaluated at the retarded time $u$. Explicitly, it comes

$$H_{22} = 1 + \left(-\frac{107}{42} + \frac{55}{42}\nu\right)x + 2\pi x^{3/2} + \left(-\frac{2173}{1512} - \frac{1069}{216}\nu + \frac{2047}{1512}\nu^2\right)x^2 + \left[-\frac{107\pi}{21} + \left(\frac{34\pi}{21} - 24i\right)\nu\right]x^{5/2} \\ + \left[\frac{27027409}{646800} - \frac{856}{105}\gamma_E + \frac{428i\pi}{105} + \frac{2\pi^2}{3} + \left(-\frac{278185}{33264} + \frac{41\pi^2}{96}\right)\nu - \frac{20261}{2772}\nu^2 + \frac{114635}{99792}\nu^3 - \frac{428}{105}\ln(16x)\right]x^3 \\ + \left[-\frac{2173\pi}{756} + \left(-\frac{2495\pi}{378} + \frac{14333i}{162}\right)\nu + \left(\frac{40\pi}{27} - \frac{4066i}{945}\right)\nu^2\right]x^{7/2} \\ + \left[-\frac{846557506853}{12713500800} + \frac{45796}{2205}\gamma_E - \frac{22898}{2205}i\pi - \frac{107}{63}\pi^2 + \frac{22898}{2205}\ln(16x)\right. \\ + \left. \left(-\frac{336005827477}{4237833600} + \frac{15284}{441}\gamma_E - \frac{219314}{2205}i\pi - \frac{9755}{32256}\pi^2 + \frac{7642}{441}\ln(16x)\right)\nu \right. \\ + \left. \left(\frac{256450291}{7413120} - \frac{1025}{1008}\pi^2\right)\nu^2 - \frac{81579187}{15567552}\nu^3 + \frac{26251249}{31135104}\nu^4\right]x^4 + \mathcal{O}(x^{9/2}). \quad (6.17)$$

In the test-mass limit $\nu \to 0$, this new result is also nicely in agreement with the prediction of linear black-hole perturbation theory [137].

Note that the two unphysical scales $r_0$ and $b_0$ appearing in the intermediate calculations are absent from the final

results. More precisely, these constant appear in the expressions of the source moments and in the relations between radiative and source moments, notably in (2.9). Although this is expected, the constants in both contributions cancel each other and this constitutes a highly nontrivial check of the MPM-algorithm and the robustness of the PN integration of the source moments.

Finally, we have also extracted the zero-frequency quadrupole mode $(\ell, m) = (2, 0)$ at 4PN order,

<footnote n="8">⁸See [135] for a similar argument in the context of self-forces, using an observer sitting on the particle and equipped with a flashlight.</footnote>
$$H_{20} = -\sqrt{\frac{3}{8} \frac{\ell^{\langle i} \ell^{j \rangle}}{c^2 m \nu x}} U_{ij}. \quad (6.18)$$

As discussed in Sec. V C, this mode arises from integration of the nonlinear DC memory terms over the past history of the system, assuming a model for the quasicircular evolution of the orbit in the past. Since the integration increases the effect by the inverse of the 2.5PN order, with the present 4PN formalism we are able to control the (2, 0) mode only with relative 1.5PN precision. We find

$$H_{20} = -\frac{5}{14\sqrt{6}} \left[ 1 + \left( -\frac{4075}{4032} + \frac{67}{48}\nu \right) x + \mathcal{O}(x^2) \right]. \quad (6.19)$$

Notice that the 1.5PN term of this mode vanishes. This is due to the fact that the 1.5PN correction in the model of evolution of the quasicircular orbit in the past, which results in the 1.5PN term in Eq. (5.20), exactly cancels the 1.5PN direct contribution of the “tail-of-memory” at 4PN order, and given by (5.22). That is,

$$H_{20}^{\text{ToM}} = -H_{20}^{\text{mem,1.5PN}} = -\frac{2\sqrt{2}}{7\sqrt{3}}\pi x^{3/2}. \quad (6.20)$$

The result (6.19) is in full agreement with Eq. (4.3a) of [110], obtained from the general expression of nonlinear memory terms in terms of radiative moments. Indeed, recall that the tail-of-memory integral (5.21) at 4PN order can be simply obtained from the leading memory integral (2.6) at 2.5PN order by replacing the canonical moment by the corresponding radiative moment including the tail effect (2.5) at relative 1.5PN order, see Eq. (2.10). This confirms that the tail-of-memory is adequately taken into

account in the computation of the memory using radiative moments defined at future null infinity [110,111].

## ACKNOWLEDGMENTS

It is a great pleasure to thank Laura Bernard, Alessandra Buonanno, Bala Iyer, Tanguy Marchand, Sylvain Marsat, Rafael Porto, and Adam Pound for enlightening exchanges during the project. F. L. is grateful to the Institut d'Astrophysique de Paris for its hospitality during this project. He received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 817791). G. F. thanks IIT Madras for a Visiting Faculty Fellow position under the IoE program during the completion of this work.

## APPENDIX A: TEST OF THE BOOSTED SCHWARZSCHILD SOLUTION

Among all tests that can be performed to check the expression of the mass quadrupole moment, one of the simplest is the boosted Schwarzschild limit. Despite its apparent simplicity, it is quite efficient and was crucially used to fix a remaining ambiguity constant in an early computation of the flux at 3PN order [103]. The principle is quite transparent; if we remove one of the two black holes, then our system reduces to a single Schwarzschild black hole of mass $m_1$, boosted at a (constant) speed $v_1$. The multipole moments of a boosted Schwarzschild solution (BSS) are straightforward to derive, and have been previously determined at 3PN order in [103]. Extending this work at 4PN order, we find that the quadrupole moment of a BBS of mass M, boosted with a velocity V, reads

$$I_{ij}^{\text{BSS}} = M t^2 V^{\langle i} V^{j \rangle} \left[ 1 + \frac{9}{14} \frac{V^2}{c^2} + \frac{83}{168} \frac{V^4}{c^4} + \frac{507}{1232} \frac{V^6}{c^6} + \frac{45923}{128128} \frac{V^8}{c^8} \right] \\ + \frac{4}{7} \frac{G^2 M^3}{c^6} V^{\langle i} V^{j \rangle} + \frac{10}{21} \frac{G^2 M^3}{c^8} V^2 V^{\langle i} V^{j \rangle} + \mathcal{O}\left(\frac{1}{c^{10}}\right). \quad (A1)$$

On the other hand, taking the BSS limit ($m_2, v_2$) → (0, 0) of the renormalized mass quadrupole moment defined in [93] (on generic orbits, out of the c.m. frame), it comes

$$\lim_{\text{BSS}} I_{ij} = m_1 t^2 v_1^{\langle i} v_1^{j \rangle} \left[ 1 + \frac{9}{14} \frac{v_1^2}{c^2} + \frac{83}{168} \frac{v_1^4}{c^4} + \frac{507}{1232} \frac{v_1^6}{c^6} + \frac{45923}{128128} \frac{v_1^8}{c^8} \right] \\ + \frac{4}{7} \frac{G^2 m_1^3}{c^6} v_1^{\langle i} v_1^{j \rangle} + \frac{10}{21} \frac{G^2 m_1^3}{c^8} v_1^2 v_1^{\langle i} v_1^{j \rangle} + \mathcal{O}\left(\frac{1}{c^{10}}\right). \quad (A2)$$

Both expressions coincide under the identification ($m_1, v_1$) = (M, V), therefore this test is conclusive.
APPENDIX B: CORRECTIONS TO THE METRIC DUE TO INFRARED COMMUTATORS

This appendix displays the formal contributions of the infrared commutators to the gothic metric perturbation $h^{\mu\nu} = \sqrt{-gg^{\mu\nu}} - \eta^{\mu\nu}$, up to 4PN order, as discussed in Sec. III. Those contributions are to be added to Eq. (A.2) of the work [91] in order to obtain the full metric at 4PN, and the potentials entering it are the one defined by Eq. (A.4) of [91]. It thus comes, at the required order [the commutators being defined by (3.3)]

$$ \begin{aligned} \Delta h^{00} = & -\frac{2(d-1)^2}{(d-2)^2 c^4} C\{V^2\} + \frac{1}{c^6} \left[ \frac{8(d-1)^2(d-3)}{(d-2)^3} C\{KV\} + \frac{8(d-3)}{d-2} C\{V_i V_i\} - \frac{4(d-1)}{d-2} C\{V \hat{W}\} \right. \\ & \left. - \frac{4(d-1)^3}{(d-2)^3} \left( \frac{1}{3} C\{V^3\} + V C\{V^2\} + C\{VC\{V^2\}\} \right) \right] \\ & + \frac{1}{c^8} \left[ \frac{2(d-1)^3(4-3d)}{(d-2)^4} \left( \frac{V}{3} C\{V^3\} + V C\{VC\{V^2\}\} + \frac{1}{4}(C\{V^2\})^2 - \frac{2(d-1)V^2}{4-3d} C\{V^2\} \right) \right. \\ & \left. - \frac{4(d-1)^2(d-3)(4-3d)}{(d-2)^4} \left( V C\{KV\} + \frac{K}{2} C\{V^2\} \right) - \frac{4(d-3)(4-3d)}{(d-2)^2} V C\{V_i V_i\} \right. \\ & \left. + \frac{8(d-1)(d-4)}{(d-2)^2} V_i C\{VV_i\} + \frac{2(d-1)(4-3d)}{(d-2)^2} V C\{V \hat{W}\} - \frac{4(d-1)^2}{(d-2)^2} \hat{W} C\{V^2\} \right] + \mathcal{O}\left(\frac{1}{c^{10}}\right), \end{aligned} \quad (B1a) $$

$$ \begin{aligned} \Delta h^{0i} = & \frac{4(d-1)}{(d-2)c^5} C\{VV_i\} + \frac{1}{c^7} \left[ \frac{4(d-1)^2}{(d-2)^2} \left( C\{V^2 V_i\} + V_i C\{V^2\} + V C\{VV_i\} + C\{V_i C\{V^2\}\} + C\{VC\{VV_i\}\} \right) \right. \\ & \left. - \frac{8(d-1)(d-3)}{(d-2)^2} C\{KV_i\} + 8C\{\hat{W}V_i\} - 8C\{\hat{W}_{ik}V_k\} + \frac{8(d-1)}{d-2} C\{V \hat{R}_i\} \right] + \mathcal{O}\left(\frac{1}{c^9}\right), \end{aligned} \quad (B1b) $$

$$ \begin{aligned} \Delta h^{ij} = & -\frac{8(d-1)}{(d-2)c^8} \left( V_i C\{VV_j\} + V_j C\{VV_i\} \right) \\ & + \frac{\delta_{ij}}{c^8} \left[ \frac{4(d-1)^2(d-3)}{(d-2)^3} \left( V C\{KV\} + \frac{K}{2} C\{V^2\} \right) - \frac{2(d-1)^3}{(d-2)^3} \left( \frac{V}{3} C\{V^3\} + V C\{VC\{V^2\}\} + \frac{1}{4}(C\{V^2\})^2 \right) \right. \\ & \left. + \frac{4(d-3)}{d-2} V C\{V_k V_k\} + \frac{8(d-1)}{d-2} V_k C\{VV_k\} - \frac{2(d-1)}{d-2} V C\{V \hat{W}\} \right] + \mathcal{O}\left(\frac{1}{c^{10}}\right). \end{aligned} \quad (B1c) $$

[1] A. Einstein, Über gravitationswellen, Sitzber. Preuss. Akad. Wiss. Berlin 1, 154 (1918).

[2] L. Landau and E. Lifshitz, The Classical Theory of Fields (Pergamon, Oxford, 1971).

[3] P. Peters and J. Mathews, Gravitational radiation from point masses in a Keplerian orbit, Phys. Rev. 131, 435 (1963).

[4] L. Blanchet and T. Damour, Radiative gravitational fields in general relativity. I. General structure of the field outside the source, Phil. Trans. R. Soc. A 320, 379 (1986).

[5] L. Blanchet, Radiative gravitational fields in general relativity. II. Asymptotic behaviour at future null infinity, Proc. R. Soc. A 409, 383 (1987).

[6] L. Blanchet and T. Damour, Tail-transported temporal correlations in the dynamics of a gravitating system, Phys. Rev. D 37, 1410 (1988).

[7] L. Blanchet and T. Damour, Hereditary effects in gravitational radiation, Phys. Rev. D 46, 4304 (1992).

[8] L. Blanchet, On the multipole expansion of the gravitational field, Classical Quantum Gravity 15, 1971 (1998).

[9] L. Blanchet, T. Damour, B. R. Iyer, C. M. Will, and A. G. Wiseman, Gravitational Radiation Damping of Compact Binary Systems to Second Post-Newtonian Order, Phys. Rev. Lett. 74, 3515 (1995).

[10] L. Blanchet, Gravitational-wave tails of tails, Classical Quantum Gravity 15, 113 (1998).
[11] L. Blanchet, B. R. Iyer, and B. Joguet, Gravitational waves from inspiralling compact binaries: Energy flux to third post-Newtonian order, Phys. Rev. D **65**, 064005 (2002); **71**, 129903(E) (2005).

[12] L. Blanchet, G. Faye, B. R. Iyer, and B. Joguet, Gravitational-wave inspiral of compact binary systems to 7/2 post-Newtonian order, Phys. Rev. D **65**, 061501(R) (2002); **71**, 129902(E) (2005).

[13] L. Blanchet and B. R. Iyer, Hadamard regularization of the third post-Newtonian gravitational wave generation of two point masses, Phys. Rev. D **71**, 024004 (2004).

[14] L. Blanchet, T. Damour, G. Esposito-Farèse, and B. R. Iyer, Gravitational Radiation from Inspiralling Compact Binaries Completed at the Third Post-Newtonian Order, Phys. Rev. Lett. **93**, 091101 (2004).

[15] L. Blanchet, T. Damour, G. Esposito-Farèse, and B. R. Iyer, Dimensional regularization of the third post-Newtonian gravitational wave generation of two point masses, Phys. Rev. D **71**, 124004 (2005).

[16] L. Blanchet, G. Faye, B. R. Iyer, and S. Sinha, The third post-Newtonian gravitational wave polarisations and associated spherical harmonic modes for inspiralling compact binaries in quasi-circular orbits, Classical Quantum Gravity **25**, 165003 (2008).

[17] G. Faye, S. Marsat, L. Blanchet, and B. R. Iyer, The third and a half post-Newtonian gravitational wave quadrupole mode for quasi-circular inspiralling compact binaries, Classical Quantum Gravity **29**, 175004 (2012).

[18] Q. Henry, Complete gravitational-waveform amplitude modes for quasicircular compact binaries to the 3.5PN order, Phys. Rev. D **107**, 044057 (2023).

[19] C. Will and A. Wiseman, Gravitational radiation from compact binary systems: Gravitational waveforms and energy loss to second post-Newtonian order, Phys. Rev. D **54**, 4813 (1996).

[20] W. Goldberger and I. Rothstein, Effective field theory of gravity for extended objects, Phys. Rev. D **73**, 104029 (2006).

[21] W. Goldberger and A. Ross, Gravitational radiative corrections from effective field theory, Phys. Rev. D **81**, 124015 (2010).

[22] A. Ross, Multipole expansion at the level of the action, Phys. Rev. D **85**, 125033 (2012).

[23] A. K. Leibovich, N. T. Maia, I. Z. Rothstein, and Z. Yang, Second post-Newtonian order radiative dynamics of inspiralling compact binaries in the effective field theory approach, Phys. Rev. D **101**, 084058 (2020).

[24] L. Blanchet, Gravitational radiation from post-Newtonian sources and inspiralling compact binaries, Living Rev. Relativity **17**, 2 (2014).

[25] M. Maggiore, *Gravitational Waves: Volume 1: Theory and Experiments* (Oxford University Press, New York, 2008), Vol. 1.

[26] A. Buonanno and B. S. Sathyaprakash, Sources of gravitational waves: Theory and observations, in *General Relativity and Gravitation: A Centennial Perspective*, edited by A. Ashtekar, B. K. Berger, J. Isenberg, and M. MacCallum (Cambridge University Press, Cambridge, England, 2015), pp. 287–346.

[27] R. A. Porto, The effective field theorist's approach to gravitational dynamics, Phys. Rep. **633**, 1 (2016).

[28] C. Cutler et al., The Last Three Minutes: Issues in Gravitational Wave Measurements of Coalescing Compact Binaries, Phys. Rev. Lett. **70**, 2984 (1993).

[29] C. Cutler and E. E. Flanagan, Gravitational waves from merging compact binaries: How accurately can one extract the binary's parameters from the inspiral wave form?, Phys. Rev. D **49**, 2658 (1994).

[30] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Phys. Rev. X **9**, 031040 (2019).

[31] R. Abbott et al. (LIGO Scientific and Virgo Collaborations), GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, Phys. Rev. X **11**, 021053 (2021).

[32] R. Abbott et al. (LIGO Scientific, VIRGO, and KAGRA Collaborations), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run, arXiv:2111.03606 [Phys. Rev. X (to be published)].

[33] P. Amaro-Seoane et al. (LISA Collaboration), Laser interferometer space antenna, arXiv:1702.00786.

[34] The ET Collaboration (2022), https://www.et-gw.eu/.

[35] P. Grandclement and J. Novak, Spectral methods for numerical relativity, Living Rev. Relativity **12**, 1 (2009).

[36] M. Boyle, D. A. Brown, L. E. Kidder, A. H. Mroue, H. P. Pfeiffer, M. A. Scheel, G. B. Cook, and S. A. Teukolsky, High-accuracy comparison of numerical relativity simulations with post-Newtonian expansions, Phys. Rev. D **76**, 124038 (2007).

[37] S. Bernuzzi, M. Thierfelder, and B. Bruegmann, Accuracy of numerical relativity waveforms from binary neutron star mergers and their comparison with post-Newtonian waveforms, Phys. Rev. D **85**, 104030 (2012).

[38] S. Borhanian, K. G. Arun, H. P. Pfeiffer, and B. S. Sathyaprakash, Comparison of post-Newtonian mode amplitudes with numerical relativity simulations of binary black holes, Classical Quantum Gravity **37**, 065006 (2020).

[39] W. Tichy, B. Bruegmann, M. Campanelli, and P. Diener, Binary black hole initial data for numerical general relativity based on post-Newtonian data, Phys. Rev. D **67**, 064008 (2003).

[40] G. Reifenberger and W. Tichy, Alternatives to standard puncture initial data for binary black hole evolution, Phys. Rev. D **86**, 064003 (2012).

[41] N. K. Johnson-McDaniel, N. Yunes, W. Tichy, and B. J. Owen, Conformally curved binary black hole initial data including tidal deformations and outgoing radiation, Phys. Rev. D **80**, 124039 (2009).

[42] J. Healy, C. O. Lousto, H. Nakano, and Y. Zlochower, Post-Newtonian quasicircular initial orbits for numerical relativity, Classical Quantum Gravity **34**, 145011 (2017).

[43] B. Wardell, A. Pound, N. Warburton, J. Miller, L. Durkan, and A. Le Tiec, Gravitational Waveforms for Compact Binaries from Second-Order Self-Force Theory, Phys. Rev. Lett. **130**, 241402 (2023).

[44] A. Albertini, A. Nagar, A. Pound, N. Warburton, B. Wardell, L. Durkan, and J. Miller, Comparing
second-order gravitational self-force, numerical relativity, and effective one body waveforms from inspiralling, quasicircular, and nonspinning black hole binaries, Phys. Rev. D **106**, 084061 (2022).

[45] A. Albertini, A. Nagar, A. Pound, N. Warburton, B. Wardell, L. Durkan, and J. Miller, Comparing second-order gravitational self-force and effective one body waveforms from inspiralling, quasicircular and nonspinning black hole binaries. II. The large-mass-ratio case, Phys. Rev. D **106**, 084062 (2022).

[46] L. Blanchet, G. Hébrard, and F. Larrouturou, Detecting the general relativistic orbital precession of the exoplanet HD 80606b, Astron. Astrophys. **628**, A80 (2019).

[47] L. Blanchet and B. S. Sathyaprakash, Detecting the Tail Effect in Gravitational Wave Experiments, Phys. Rev. Lett. **74**, 1067 (1995).

[48] K. Arun, B. Iyer, M. Qusailah, and B. Sathyaprakash, Testing post-Newtonian theory with gravitational wave observations, Classical Quantum Gravity **23**, L37 (2006).

[49] B. Abbott et al. (LIGO Scientific Collaboration and VIRGO Collaboration), Tests of General Relativity with GW150914, Phys. Rev. Lett. **116**, 221101 (2016).

[50] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), Tests of General Relativity with GW170817, Phys. Rev. Lett. **123**, 011102 (2019).

[51] R. Abbott et al. (LIGO Scientific, VIRGO, and KAGRA Collaborations), Tests of general relativity with GWTC-3, arXiv:2112.06861 [Phys. Rev. D (to be published)].

[52] A. Sesana, Prospects for Multiband Gravitational-Wave Astronomy after GW150914, Phys. Rev. Lett. **116**, 231102 (2016).

[53] E. E. Flanagan and T. Hinderer, Constraining neutron star tidal Love numbers with gravitational wave detectors, Phys. Rev. D **77**, 021502 (2008).

[54] T. Damour and A. Nagar, Relativistic tidal properties of neutron stars, Phys. Rev. D **80**, 084035 (2009).

[55] J. Vines, T. Hinderer, and E. Flanagan, Post-1-Newtonian tidal effects in the gravitational waveform from binary inspirals, Phys. Rev. D **83**, 084051 (2011).

[56] T. Damour, A. Nagar, and L. Villain, Measurability of the tidal polarizability of neutron stars in late-inspiral gravitational-wave signals, Phys. Rev. D **85**, 123007 (2012).

[57] D. Bini, T. Damour, and G. Faye, Effective action approach to higher-order relativistic tidal interactions in binary systems and their effective one body description, Phys. Rev. D **85**, 124034 (2012).

[58] Q. Henry, G. Faye, and L. Blanchet, Tidal effects in the equations of motion of compact binary systems to next-to-next-to-leading post-Newtonian order, Phys. Rev. D **101**, 064047 (2020).

[59] Q. Henry, G. Faye, and L. Blanchet, Tidal effects in the gravitational-wave phase evolution of compact binary systems to next-to-next-to-leading post-Newtonian order, Phys. Rev. D **102**, 044033 (2020).

[60] Q. Henry, G. Faye, and L. Blanchet, Hamiltonian for tidal interactions in compact binary systems to next-to-next-to-leading post-Newtonian order, Phys. Rev. D **102**, 124074 (2020).

[61] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), GW170817: Observation of Gravitational Waves

from a Binary Neutron Star Inspiral, Phys. Rev. Lett. **119**, 161101 (2017).

[62] S. Mirshekari and C. M. Will, Compact binary systems in scalar-tensor gravity: Equations of motion to 2.5 post-Newtonian order, Phys. Rev. D **87**, 084070 (2013).

[63] L. Bernard, Dynamics of compact binary systems in scalar-tensor theories: I. Equations of motion to the third post-Newtonian order, Phys. Rev. D **98**, 044004 (2018).

[64] L. Bernard, Dynamics of compact binary systems in scalar-tensor theories: II. Center-of-mass and conserved quantities to 3PN order, Phys. Rev. D **99**, 044047 (2019).

[65] L. Bernard, Dipolar tidal effects in scalar-tensor theories, Phys. Rev. D **101**, 021501 (2020).

[66] R. N. Lang, Compact binary systems in scalar-tensor gravity. II. Tensor gravitational waves to second post-Newtonian order, Phys. Rev. D **89**, 084014 (2014).

[67] R. N. Lang, Compact binary systems in scalar-tensor gravity. III. Scalar waves and energy flux, Phys. Rev. D **91**, 084027 (2015).

[68] L. Bernard, L. Blanchet, and D. Trestini, Gravitational waves in scalar-tensor theory to one-and-a-half post-Newtonian order, J. Cosmol. Astropart. Phys. **08** (2022) 008.

[69] T. P. Sotiriou and E. Barausse, Post-Newtonian expansion for Gauss-Bonnet gravity, Phys. Rev. D **75**, 084007 (2007).

[70] A. De Felice and S. Tsujikawa, $f(R)$ theories, Living Rev. Relativity **13**, 3 (2010).

[71] K. Yagi, L. C. Stein, N. Yunes, and T. Tanaka, Post-Newtonian, quasi-circular binary inspirals in quadratic modified gravity, Phys. Rev. D **85**, 064022 (2012); **93**, 029902(E) (2016).

[72] F.-L. Julié, Gravitational radiation from compact binary systems in Einstein-Maxwell-dilaton theories, J. Cosmol. Astropart. Phys. **10** (2018) 033.

[73] B. Shiralilou, T. Hinderer, S. M. Nissanke, N. Ortiz, and H. Witek, Post-Newtonian gravitational and scalar waves in scalar-Gauss-Bonnet gravity, Classical Quantum Gravity **39**, 035002 (2022).

[74] P. Jaranowski and G. Schäfer, Towards the fourth post-Newtonian Hamiltonian for two-point-mass systems, Phys. Rev. D **86**, 061503(R) (2012).

[75] P. Jaranowski and G. Schäfer, Dimensional regularization of local singularities in the 4th post-Newtonian two-point-mass Hamiltonian, Phys. Rev. D **87**, 081503(R) (2013).

[76] T. Damour, P. Jaranowski, and G. Schäfer, Non-local-in-time action for the fourth post-Newtonian conservative dynamics of two-body systems, Phys. Rev. D **89**, 064058 (2014).

[77] T. Damour, P. Jaranowski, and G. Schäfer, On the conservative dynamics of two-body systems at the fourth post-Newtonian approximation of general relativity, Phys. Rev. D **93**, 084014 (2016).

[78] D. Bini and T. Damour, Analytical determination of the two-body gravitational interaction potential at the fourth post-Newtonian approximation, Phys. Rev. D **87**, 121501(R) (2013).

[79] L. Bernard, L. Blanchet, A. Bohé, G. Faye, and S. Marsat, Fokker action of non-spinning compact binaries at the fourth post-Newtonian approximation, Phys. Rev. D **93**, 084037 (2016).
[80] L. Bernard, L. Blanchet, A. Bohé, G. Faye, and S. Marsat, Energy and periastron advance of compact binaries on circular orbits at the fourth post-Newtonian order, Phys. Rev. D **95**, 044026 (2017).

[81] L. Bernard, L. Blanchet, A. Bohé, G. Faye, and S. Marsat, Dimensional regularization of the IR divergences in the Fokker action of point-particle binaries at the fourth post-Newtonian order, Phys. Rev. D **96**, 104043 (2017).

[82] L. Bernard, L. Blanchet, G. Faye, and T. Marchand, Center-of-mass equations of motion and conserved integrals of compact binary systems at the fourth post-Newtonian order, Phys. Rev. D **97**, 044037 (2018).

[83] S. Foffa and R. Sturani, The dynamics of the gravitational two-body problem in the post-Newtonian approximation at quadratic order in the Newton's constant, Phys. Rev. D **87**, 064011 (2013).

[84] S. Foffa and R. Sturani, Tail terms in gravitational radiation reaction via effective field theory, Phys. Rev. D **87**, 044056 (2012).

[85] C. R. Galley, A. K. Leibovich, R. A. Porto, and A. Ross, Tail effect in gravitational radiation reaction: Time non-locality and renormalization group evolution, Phys. Rev. D **93**, 124010 (2016).

[86] S. Foffa, P. Mastrolia, R. Sturani, and C. Sturm, Effective field theory approach to the gravitational two-body dynamics at fourth post-Newtonian order and quintic in the Newton constant, Phys. Rev. D **95**, 104009 (2017).

[87] R. Porto and I. Rothstein, On the apparent ambiguities in the post-Newtonian expansion for binary systems, Phys. Rev. D **96**, 024061 (2017).

[88] S. Foffa and R. Sturani, Conservative dynamics of binary systems to fourth post-Newtonian order in the EFT approach I: Regularized Lagrangian, Phys. Rev. D **100**, 024047 (2019).

[89] S. Foffa, R. Porto, I. Rothstein, and R. Sturani, Conservative dynamics of binary systems to fourth post-Newtonian order in the EFT approach II: Renormalized Lagrangian, Phys. Rev. D **100**, 024048 (2019).

[90] J. Blümlein, A. Maier, P. Marquard, and G. Schäfer, Fourth post-Newtonian Hamiltonian dynamics of two-body systems from an effective field theory approach, Nucl. Phys. B**955**, 115041 (2020).

[91] T. Marchand, Q. Henry, F. Larrouturou, S. Marsat, G. Faye, and L. Blanchet, The mass quadrupole moment of compact binary systems at the fourth post-Newtonian order, Classical Quantum Gravity **37**, 215006 (2020).

[92] F. Larrouturou, Q. Henry, L. Blanchet, and G. Faye, The quadrupole moment of compact binaries to the fourth post-Newtonian order: I. Non-locality in time and infra-red divergencies, Classical Quantum Gravity **39**, 115007 (2022).

[93] F. Larrouturou, L. Blanchet, Q. Henry, and G. Faye, The quadrupole moment of compact binaries to the fourth post-Newtonian order: II. Dimensional regularization and renormalization, Classical Quantum Gravity **39**, 115008 (2022).

[94] Q. Henry, G. Faye, and L. Blanchet, The current-type quadrupole moment and gravitational-wave mode $(\ell, m) = (2, 1)$ of compact binary systems at the third post-Newtonian order, Classical Quantum Gravity **38**, 185004 (2021).

[95] G. Faye, L. Blanchet, and B. R. Iyer, Non-linear multipole interactions and gravitational-wave octupole modes for inspiralling compact binaries to third-and-a-half post-Newtonian order, Classical Quantum Gravity **32**, 045016 (2015).

[96] T. Marchand, L. Blanchet, and G. Faye, Gravitational-wave tail effects to quartic non-linear order, Classical Quantum Gravity **33**, 244003 (2016).

[97] L. Blanchet, G. Faye, and F. Larrouturou, The quadrupole moment of compact binaries to the fourth post-Newtonian order: From source to canonical moment, Classical Quantum Gravity **39**, 195003 (2022).

[98] D. Trestini, F. Larrouturou, and L. Blanchet, The quadrupole moment of compact binaries to the fourth post-Newtonian order: Relating the harmonic and radiative metrics, Classical Quantum Gravity **40**, 055006 (2023).

[99] D. Trestini and L. Blanchet, Gravitational-wave tails of memory, Phys. Rev. D **107**, 104048 (2023).

[100] F. Messina and A. Nagar, Parametrized-4.5PN TaylorF2 approximants and tail effects to quartic nonlinear order from the effective one body formalism, Phys. Rev. D **95**, 124001 (2017); **96**, 049907(E) (2017).

[101] L. Blanchet, G. Faye, Q. Henry, F. Larrouturou, and D. Trestini, companion Letter, Gravitational-Wave Phasing of Quasicircular Compact Binary Systems to the Fourth-and-a-Half Post-Newtonian Order, Phys. Rev. Lett. **131**, 121402 (2023).

[102] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevD.108.064041 which collects the results associated with the gravitational-wave phase at 4.5PN and the $(\ell, m) = (2, 2)$ mode at 4PN.

[103] L. Blanchet, T. Damour, and B. R. Iyer, Surface-integral expressions for the multipole moments of post-Newtonian sources and the boosted Schwarzschild solution, Classical Quantum Gravity **22**, 155 (2005).

[104] K. Thorne, Multipole expansions of gravitational radiation, Rev. Mod. Phys. **52**, 299 (1980).

[105] L. Kidder, Using full information when computing modes of post-Newtonian waveforms from inspiralling compact binaries in circular orbits, Phys. Rev. D **77**, 044016 (2008).

[106] L. Blanchet, Quadrupole-quadrupole gravitational waves, Classical Quantum Gravity **15**, 89 (1998).

[107] S. Foffa and R. Sturani, Hereditary terms at next-to-leading order in two-body gravitational dynamics, Phys. Rev. D **101**, 064033 (2020); **103**, 089901(E) (2021).

[108] W. D. Goldberger, A. Ross, and I. Z. Rothstein, Black hole mass dynamics and renormalization group evolution, Phys. Rev. D **89**, 124033 (2014).

[109] G. L. Almeida, S. Foffa, and R. Sturani, Gravitational multipole renormalization, Phys. Rev. D **104**, 084095 (2021).

[110] M. Favata, Post-Newtonian corrections to the gravitational-wave memory for quasicircular, inspiralling compact binaries, Phys. Rev. D **80**, 024002 (2009).

[111] M. Favata, The gravitational-wave memory from eccentric binaries, Phys. Rev. D **84**, 124013 (2011).

[112] T. Damour, P. Jaranowski, and G. Schäfer, Dimensional regularization of the gravitational interaction of point masses, Phys. Lett. B **513**, 147 (2001).
[113] L. Blanchet, T. Damour, and G. Esposito-Farèse, Dimensional regularization of the third post-Newtonian dynamics of point particles in harmonic coordinates, Phys. Rev. D **69**, 124007 (2004).

[114] L. Blanchet, S. Foffa, F. Larrouturou, and R. Sturani, Logarithmic tail contributions to the energy function of circular compact binaries, Phys. Rev. D **101**, 084045 (2020).

[115] O. Poujade and L. Blanchet, Post-Newtonian approximation for isolated systems calculated by matched asymptotic expansions, Phys. Rev. D **65**, 124020 (2002).

[116] L. Blanchet, G. Faye, and S. Nissanke, Structure of the post-Newtonian expansion in general relativity, Phys. Rev. D **72**, 044024 (2005).

[117] T. Marchand, L. Bernard, L. Blanchet, and G. Faye, Ambiguity-free completion of the equations of motion of compact binary systems at the fourth post-Newtonian order, Phys. Rev. D **97**, 044023 (2018).

[118] A. Gopakumar, B. R. Iyer, and S. Iyer, Second post-Newtonian gravitational radiation reaction for two-body systems: Nonspinning bodies, Phys. Rev. D **55**, 6030 (1997).

[119] A. K. Leibovich, B. A. Pardo, and Z. Yang, Radiation reaction for non-spinning bodies at 4.5PN in the effective field theory approach, arXiv:2302.11016.

[120] K. Arun, L. Blanchet, B. R. Iyer, and M. S. Qusailah, The 2.5PN gravitational wave polarisations from inspiralling compact binaries in circular orbits, Classical Quantum Gravity **21**, 3771 (2004); **22**, 3115(E) (2005).

[121] L. Blanchet and G. Schäfer, Gravitational wave tails and binary star systems, Classical Quantum Gravity **10**, 2699 (1993).

[122] A. Wiseman, Coalescing binary-systems of compact objects to (post)5/2-Newtonian order. IV. The gravitational-wave tail, Phys. Rev. D **48**, 4757 (1993).

[123] L. Blanchet, B. R. Iyer, C. M. Will, and A. G. Wiseman, Gravitational wave forms from inspiralling compact binaries to second-post-Newtonian order, Classical Quantum Gravity **13**, 575 (1996).

[124] M. Sasaki, Post-Newtonian expansion of the ingoing-wave Regge-Wheeler function, Prog. Theor. Phys. **92**, 17 (1994).

[125] H. Tagoshi and M. Sasaki, Post-Newtonian expansion of gravitational-waves from a particle in circular orbit around

a Schwarzschild black hole, Prog. Theor. Phys. **92**, 745 (1994).

[126] T. Tanaka, H. Tagoshi, and M. Sasaki, Gravitational waves by a particle in circular orbit around a Schwarzschild black hole: 5.5 Post-Newtonian formula, Prog. Theor. Phys. **96**, 1087 (1996).

[127] R. Fujita, Gravitational radiation for extreme mass ratio inspirals to the 14th post-Newtonian order, Prog. Theor. Phys. **127**, 583 (2012).

[128] R. Fujita, Gravitational waves from a particle in circular orbits around a Schwarzschild black hole to the 22nd post-Newtonian order, Prog. Theor. Phys. **128**, 971 (2012).

[129] E. Poisson and M. Sasaki, Gravitational radiation from a particle in circular orbit around a black hole. V. Black-hole absorption and tail corrections, Phys. Rev. D **51**, 5753 (1995).

[130] H. Tagoshi, S. Mano, and E. Takasugi, Post-Newtonian expansion of gravitational waves from a particle in circular orbits around a rotating black hole, Prog. Theor. Phys. **98**, 829 (1997).

[131] K. Alvi, Energy and angular momentum flow into a black hole in a Xbinary, Phys. Rev. D **64**, 104020 (2001).

[132] R. A. Porto, Absorption effects due to spin in the worldline approach to black hole dynamics, Phys. Rev. D **77**, 064026 (2008).

[133] K. Chatzioannou, E. Poisson, and N. Yunes, Tidal heating and torquing of a kerr black hole to next-to-leading order in the tidal coupling, Phys. Rev. D **87**, 044022 (2013).

[134] M. V. S. Saketh, J. Steinhoff, J. Vines, and A. Buonanno, Modeling horizon absorption in spinning binary black holes using effective worldline theory, Phys. Rev. D **107**, 084006 (2023).

[135] S. L. Detweiler, Perspective on gravitational self-force analyses, Classical Quantum Gravity **22**, S681 (2005).

[136] Q. Henry, S. Marsat, and M. Khalil, Spin contributions to the gravitational-waveform modes for spin-aligned binaries at the 3.5PN order, Phys. Rev. D **106**, 124018 (2022).

[137] R. Fujita and B. Iyer, Spherical harmonic modes of 5.5 post-Newtonian gravitational wave polarizations and associated factorized resummed waveforms for a particle in circular orbit around a Schwarzschild black hole, Phys. Rev. D **82**, 044051 (2010).