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

This article provides the details on the technical derivation of the gravitational waveform and total gravitational-wave energy flux of non-spinning compact binary systems to the 4PN (fourth post-Newtonian) 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"post-adiabatic"correction needed to evaluate the tail integral with 2.5PN relative precision; (iv) we discuss the relation between the binary's orbital frequency in quasi-circular 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.


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][2][3], is fairly straightforward, the systematic expansion to high PN orders entails a number of subtle issues. One is the control of physical non-linear 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][5][6][7][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 non-spinning objects [9][10][11][12][13][14][15][16][17][18], i.e., up to the (v/c) 7 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][25][26][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 series in terms of the PN parameter x = ( GmΩ c 3 ) 2/3 , where Ω 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][31][32], but also for future generations of detectors, such as the Laser Interferometer Space Antenna (LISA) [33] or the

II. NON-LINEAR EFFECTS IN THE GW PROPAGATION
In this work, 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 quasi-circular orbits to the 4.5PN order. Let us recall that the general multipole decomposition of the flux reads [104] where U L and V L (for ℓ 2) respectively denote the mass and current radiative multipole moments (which are STF in their indices), and the numerical coefficients are given by a ℓ = (ℓ + 1)(ℓ + 2) (ℓ − 1)ℓ ℓ! (2ℓ + 1)!! and b ℓ = 4ℓ(ℓ + 2) (ℓ − 1) (ℓ + 1)! (2ℓ + 1)!! . (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 ≡ R N , and written as a perturbation to the ordinary metric g ab , reads [104] h TT ab = 4G where ⊥ abij (N ) is the usual transverse and traceless (TT) projector and the multipole moments are defined at the retarded time u ≡ 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, non precessing, 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 non-linear 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

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 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] ij (u − τ ) ln cτ 2b 0 + 11 12 .
(2.5) 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 and trace-free (STF) mass and current moments M L ≡ M i 1 i 2 ...i ℓ and S L ≡ S i 1 i 2 ...i ℓ ; a superscript (n) indicates n time derivatives; hats and angular brackets denote a symmetric trace-free (STF) projection, for instancex L = x L = STF[x L ] with x L = x i 1 · · · x i ℓ ; the d'Alembertian operator is defined with respect to the flat Minkowski metric ≡ η µν ∂µν = ∆ − c −2 ∂ 2 t ; 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 ǫ ijk with ǫ 123 = 1; finally, as usual, we dubb "nPN" a quantity of order O(c −2n ). 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 × M ij × 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-ofmemory comes along with a cubic interaction involving the constant angular momentum, i.e. M × S i × M jk , see the last line of (2.9). We have [99] j a (u − τ ) ln j a (u − τ ) ln j a (u − τ ) ln 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.
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: Reporting the results of [95], the above contributions read where the underlined indices within angled brackets are excluded from the STF operation, and 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 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

Radiative moments entering the quasi-circular flux at the 4.5PN order
For quasi-circular 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 quasi-circular orbits. So only non-local contributions such as tails can potentially contribute. Second, the quadratic non-local 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 quasi-circular flux. For the mass quadrupole, only the quartic M 3 × M ij , naturally dubbed "tails-of-tails-of-tails", can play a role. It has been computed in [96], and reads The full U 4.5PN ij 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 quasi-circular 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 quasi-circular flux. The only moments playing a role beyond the mass quadrupole are [95,96] 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 imposes 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/ε ≡ 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 CoM frame (see [93] for the complete expressions at 3PN in a generic frame), they are where we employ the notations of [94] for current moments, and where withq ≡ 4πe γE , γ E being the Euler constant. The constant ℓ 0 is the length-scale associated with dimensional regularization, such that the d-dimensional gravitational coupling constant G d = ℓ ε 0 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 β-coefficients of the renormalization group flows for these multipole moments [21,108,109]. 2 As explained in Sec. IV the poles present in (2.18) will actually be cancelled 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 coordinate change (i.e., fully non-linear 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, non-local terms cannot appear at that order in the relation between canonical and source/gauge moments. Indeed, in the three-dimensional case, non-local 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 quasi-circular orbits, and we only need the 4PN order in the canonical/source/gauge relations to control the 4.5PN quasi-circular flux.

Procedure in d dimensions
The solution of the linearized, d dimensional, vacuum Einstein field equations in harmonic gauge reads 20) where the canonical metric is parametrized by three types of source moments as [94] h 00 and the gauge vector ϕ µ 1 , by four types of gauge moments, The conventions, notably for the indices of J i|L and Z i|L , are still those of [94], 3 and we have shortened This shorthand satisfies lim d→3 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 App. A of [94], where the number of independent components is given by (A6b). 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 non-linear interaction terms Ω µν n and ∆ µ n as functions of h µν can n and ϕ µ n are identical to the three-dimensional case. The main, but minor, difference concerns the two key quantities X µν n and Y µν n . If their formal definition, Eqs. (3.20) of [97], does not depend of the dimension, their explicit expressions in terms of Ω µν n and ∆ µ n are slightly changed, and they read in d dimensions These quantities are endowed with a Hadamard finite part (FP) regularization when B → 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 and Y µν n will be non-zero only when the retarded integrals develop a pole ∝ 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).

Integration techniques
Let us recall that the solution at an order n 2 is composed of some corrections to the canonical moments, extracted from X µν n and Y µν n given by (2.24), as well as a gauge vector, denoted ϕ µ n . Thanks to the prefactor B in (2.24), and similarly to the three-dimensional case, the computation of X µν n and Y µν n only requires the near-zone (r → 0) behavior of the functions Ω µν n and ∆ µ n , as the FP regularization is introduced to cope with the singular behavior of multipole expansions when r → 0. Using the expansion formulas displayed in App. A of [81], the computation boils down to the evaluation of integrals of the form where b = 1, 2 and H can bear simple poles ∝ 1/ε and can indifferently be a local or non-local function of time. In this section, we set r 0 = ℓ 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.
The current moment can be expressed using a Levi-Civita symbol in the d → 3 limit. This relation is in both conventions: J i|L = ǫ ii ℓ a 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−1j) 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].
Techniques to solve wave-like equations in d dimensions, with more general radial dependence than (2.25) but same multipolarity ℓ, 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 D = {|x ′ | < r} and R 3 − 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 ≡ J b p,q 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 ε, in order to cope with the case where H(t) contains a pole, namely . After resummation, we obtain a retarded homogeneous solution of the wave equation [93] Under those very restrictive conditions, we have where we recall that H can bear a pole 1/ε, and our notationq = 4πe γE (with γ E the Euler constant). With the previous formula at hand, we are ready to compute the quantities X µν n and Y µν n and, therefore, obtain the correction to the canonical moments to the n-th 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 ϕ µ n , 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 ϕ µ n = φ µ n + ψ µ n , where ψ µ n is also extracted from X µν n and Y µν n , and φ µ n is obtained by the direct integration In d dimensions, it is in general not possible to find a convenient, explicit expression for φ µ n (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 → 0. Indeed, following the procedure of [97], it is the near-zone expansion of φ µ n , denotedφ µ n , that needs to be inserted into the expressions of Ω µν n+1 and ∆ µ n+1 sourcing the next order quantities X µν n+1 and Y µν n+1 . The quantityφ µ n can be obtained directly from the near-zone expansion of the source term, i.e.∆ µ n , plus a crucial homogeneous solution, itself in the form of a near-zone expansion, namely: 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]: The homogeneous solutionφ µ hom n is given by Eq. (3.20) of [81], which in this case yields (with c = 1) where we have decomposed the source in its multipolar components as and shortened the iterated inverse Laplacians as [81,115,116] In the particular case where the source term has a structure given bȳ where the F ℓ,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 non-static multipolar moments, such as ∆ µ W×Iij . Fortunately, we have shown that those more complicated homogeneous solutions do not contribute at 4PN order, see hereafter.

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 quasi-circular 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 × M × W ij and M × M × Y ij (where only one moment is non-static), and most importantly M × W × 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 2. Indeed, it is composed of two interactions: h M×M × ϕ Kij and h M × ϕ M×Kij , where K ij stands for either W ij or Y ij . As both h M×M and ϕ M×Kij bear q = 2, the cubic sources will have q 2, and thus cannot satisfy the q = 1 necessary condition to have a non-vanishing function G(t) in Eq. (1.1) in [97].
As for the last cubic interaction, M × W × I ij , it turns out that the source term does have a q = 1 component, arising from the interaction h Iij × ϕ M×W . More precisely, the term with q = 1 comes exclusively from the interaction between: (i) the homogeneous solutionφ hom 2 at quadratic order arising from the interaction M × W, i.e.φ hom M×W , which bears q = 0; and (ii) the q = 1 piece of the linear metric h Iij ; see (2.21). Note that, as it involves a quadratic interaction with only one dynamical moment, ϕ M×W can be calculated using (2.34). By contrast, the homogeneous solution for the interaction W × I ij , i.e.φ hom W×Iij , does not fulfill the p − ℓ − 3 ∈ 2N 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 × W × I ij interaction, which reads This result is consistent with the corresponding one in three dimensions, where the interaction M×W×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 × M × W ij and M × M × 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 cancelled 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 µ gen = x µ can + ζ µ , which transforms the gravitational field according to Eqs. (3.1)-(3.3) in [97]. A simple dimensional analysis shows that only the M × 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 where only j = 0 is relevant at 4PN. Implementing the procedure exposed in Sec. II of [92] with the "pure gauge" 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 δ indirect M ij ≡ δ ζ 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: 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 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 δ 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 quasi-circular orbits in Sec. IV A, and which contains crucial contributions due to the (IR) dimensional regularization. Adding up the non-linear 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 × M × W ij and M × M × 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][92][93]. We shall show that the IR commutators enter the flux and the (2, 2) mode for generic orbits, however they do not contribute in the quasi-circular limit. Since they are zero for quasi-circular orbits, all the results explicitly presented in the previous works [91][92][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 µν ≡ √ −gg µν − η µν obeys the gauge-fixed Einstein field equations where the pseudo stress-energy tensor τ µν 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, sayh µν , is then given formally, to arbitrarily high orders, by [81,117] h µν = 16πG 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 (t), as functionals of the MPM expansion of τ µν 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 −1 ret , 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 d d x ′ (|x ′ |/r 0 ) B [· · · ]. Furthermore, the operator FP −1 ret acts on the near-zone expansion of the pseudo stress-energy tensorτ µν , which can be derived at any PN order from the lower order expressions ofh µν . 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 . This is the root of the appearance of the IR commutators in the PN metric, that are defined as where F is a regular (smooth) function in the near zone but typically diverges at spatial infinity, when r → +∞; 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 where S is composed of a compact-supported sector (proportional to Dirac distributions in the case of point particles), and a non-compact one (composed of products of derivatives of lower-order PN potentials). We refer to App. A of [91] for complete definitions of the PN potentials. Now those potentials are constructed by means of the operator FP −1 ret , i.e., they are given, by definition, as which accounts for the first term in the right side of (3.2). The second term in (3.2) is non-local 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 ∂ i P ∂ i P at some nPN order, 4 h 00ii (n) = · · · + ∂ i P ∂ i P + · · · , (3.6) then it can be solved using the PN-expanded retarded propagator FP where we have used Eq. While, formally, the commutators appear at a relative 1PN order in the metric [see Eqs. (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ḡ 00 only, through the combinationh 00ii . Since it is the spatial gradient ∂ iḡ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 → +∞, that homogeneous solution must be regular in the near zone, when r → 0. These facts constrain its form: given the function F , there exist a set {ĝ L } of (yet unspecified) functions of time only, such that where we recall that the inverse Laplacians are defined in Eq. (2.33). In order to explicitly determine the functionŝ g L , we work out the propagator applied to a generic source S(x, t): 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 non-local 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: Note the similarity with the expression of the quantities X µν n and Y µν n 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 and Y µν n are ultraviolet ("UV") commutators since they depend on the near zone behavior of the source term, when r → 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 → +∞ with t = const (spatial infinity). We expand the function in d = 3 + ε dimensions as where the sums extend on all multipolarities ℓ, 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 |x − x ′ | α in (3.9) (where α depends on the parity) when |x ′ | → +∞. Actually this expansion is valid as soon as r ′ = |x ′ | > r = |x| and reads, using the STF decomposition and the notation (2.33), 5 the object ∆ ′j∂′ L r ′α being explicitly known as With the latter formulas at hand, it is straightforward to obtain the functionsĝ L (t) parametrizing the commutator in (3.8) in terms of the expansion coefficients in (3.11). We find where we define the angular average over the (d − 1)-dimensional spherê .
(3.15) 5 Eq. (3.12) is equivalent to the more familiar (but less convenient for our purpose) expansion in terms of Gegengauer polynomials C γ ℓ (x), Note the non-local 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 ∝ (d − 3) −1 by itself. Therefore, we can safely work in the d = 3 limit, as long as the coefficientsf L p,q 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 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 2. In addition, the odd part of both the lowest order potentials V and V i (defined in the App. A of [91]) starts at 1.5PN (more precisely, the 0.5PN term of those two potentials is vanishing in the d → 3 limit). Investigating the metric (B1), one can see that only two interactions can lead to corrections with 4PN accuracy: V 2 and VŴ , where we denoteŴ ≡Ŵ kk . 6 The commutators affect the quantities Σ andμ 1 , defined in Eqs. (2.10) and (2.17) of [91], as where σ 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 Dirac distribution at the position 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 coefficientsf L p,q . 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 where m = m 1 + m 2 andμ(t) = d 3 x ′ σ(x ′ , t) reduces to a constant at leading order for arbitrary matter systems, so thatμ(t) = m + O(c −2 ) with m constant, while I = d 3 x ′ σ(x ′ , t) r ′2 represents the effective moment of inertia. For compact binaries, those quantities read explicitlỹ 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Ŵ is not known in the whole space (see App. 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

19)
where σ kk and σ are the compact-supported expression defined in Eq. (2.15) of [91]. The compact-supported part of 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 non-compact sector, it readŝ 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,Ŵ NC does not contribute to the IR commutator at the required order. Applying the formula (3.14), one finds the nice relation With the two explicit expressions (3.18) and (3.23) at hand, one can compute the resulting correction to the source mass quadrupole ∆I ij in (3.17). In the CoM frame, it comes

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 in 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 quasi-circular 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 non-local 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][80][81][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 ε = 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 ε → 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 ε, 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 CoM 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 quasi-circular orbits, the expression of the renormalized source quadrupole on quasi-circular orbits remains identical to the one displayed in Eq. (6.11) of [93], namely where, introducing the PN parameter γ ≡ Gm rc 2 , the coefficients are given by 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 non-local tail integral, given by Eq. (6.5) in [93]. The 4PN logarithms ln(16γ e 2γE ) in A and B are due to the conservative part of this non-local 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 App. 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] where ∆ ≡ (m 1 − m 2 )/m, and the coefficients are 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] where we denote L i ≡ ǫ ijk x j v k , and where 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 quasi-circular orbits are displayed e.g. in Sec. 9.1 of [24]. Note that, as the first non-local 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 quasi-circular 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 quasi-circular 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 γ = Gm rc 2 , the acceleration for quasi-circular orbit reads [82] 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: 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 ordeṙ Gmν 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 ω defined from the equations of motion by (5.2), and the GW half-frequency Ω which is measured in the waves propagating in the far zone. We pose and distinguish it later from the PN parameter x defined from the GW half-frequency Ω, see Eq. (6.9).

B. Post-adiabatic 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 quasi-circular orbits, and up to 3.5PN precision in the multipoles, an adiabatic approximation (considering the orbital elements r and ω to be constant in time) is sufficient, and one can follow the lines of [120], together with the integrals presented in the App. B of [121]. However, Eqs. (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 "post-adiabatic" (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 quasi-circular orbits reduce to elementary integrals of the type where n ∈ N ⋆ (the case n = 0 does not appear at 4PN), a is usually a rational fraction, and the constant τ 0 denotes either 2r 0 /c or 2b 0 /c. The orbital phase φ, the orbital frequency ω =φ and the parameter y = ( Gmω c 3 ) 2/3 are integrated over any instant u − τ in the past. Our strategy will be to notice that the integral 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 φ(u) and the phase φ(u − τ ) in the past is of the order of the inverse of the radiation reaction scale, i.e., φ(u) − φ(u − τ ) scales as 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 and change the integration variable from τ to v, so that the integral (5.5) becomes Here τ (v) is obtained by inverting Eq. (5.7) for τ as a function of v. The main point is that, in the PA limit, ξ(u) → 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 → 0. This will give a formal (asymptotic) series in powers of ξ(u) which is the "PA" expansion of the integral.
In the PA limit we thus have τ → 0, and we can expand Eq. (5.7) to any order as From now on, we pose ξ = ξ(u), φ = φ(u), ω = ω(u) =φ(u), ω (n) = (d n ω/du n )(u) for the quantities at the current time u. Since ω (1) ≡ω = ξ ω 2 the expansion (5.9) clearly represents the PA approximation in powers of ξ and arbitrary high time derivatives of ξ denoted ξ (n) . Note that each time derivative of ξ adds a radiation reaction scale, To obtain the PA expansion of (5.8) we need to invert the series (5.9) and obtain τ as a power series in v. This is given by the Lagrange inversion theorem as where the general coefficients f p read where τ = 0 is applied after the p differentiations with respect to τ . For instance, up to 2PA order, we obtain 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 To 1PA order, the integral (5.8) reads The remaining integral is computed as follows. We transform the complex exponential into a real one by performing the change of variable v = i w for n > 0, and v = −i w 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 → +∞: (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 = −π/2 (or arg w = π/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 where s(n) is the sign function. In this form, 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: Recall that all quantities are evaluated at current time u and that the adiabatic parameter ξ =ω/ω 2 is easily computed with Eqs. (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 , 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 quasi-circular orbits, they reduce to a sum of terms of the form 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 = 0. As they only involve a simple, logarithmicfree, 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 quasi-circular 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 τ = 0. One finds 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 (ℓ, m) = (2, 0) mode for the mass quadrupole. As is clear in the following, the absence of the fast oscillating exponential in J a,0 generates an inverse power of the adiabatic parameter ξ, 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 J a,0 ; see e.g. [120]. In the following, we rely on a change of integration variables from τ to y ′ = y(u − τ ): 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 where c τ 0 = 2r 0 e 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 τ . As we need to evaluate it at relative Newtonian order only, we can safely use the adiabatic approximation. Next we perform the integral over τ , which is found to be a simple DC memory integral of the type (5.19) with a = 13/2. Hence we find dρ y(u − ρ) 13/2 = 4π 7 mν c 2 y 5/2 ℓ i ℓ j , (5.22) where ℓ i = L i /|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 CoM frame, we split the gravitational energy flux (2.1) as (6.1) The first piece, F 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.
where we recall that the numerical coefficients a ℓ and b ℓ are given in Eq. (2.2). We further split F can into a local-in-time (or instantaneous) part and a non-local one: 3) The local piece F loc can 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 non-local piece, F non-loc can , it comes from two effects: (i) the direct contribution of the non-local part of the mass quadrupole moment, given by Eq. (6.5) in [93]; (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 non-local piece reads (at the required order, we are dealing with Newtonian quantities, so we can identify source and canonical moments): 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), F non-lin , is given by all the non-linearities 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 Introducing the unit separation n = x/r between the particles, as well as the unit vector λ = ℓ × n such that (n, λ, ℓ) forms a direct orthonormal triad, we may write the relative velocity as v = rωλ +ṙ n, whereṙ is given by (5.3). Then, using (5.2), we can express the 4PN flux on quasi-circular orbits as a function of the orbital frequency ω, 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 ψ of the reduced GW, i.e. half the phase of the (2,2) mode, differs from the orbital phase φ (such that the orbital frequency is ω =φ) by a logarithmic, tail-induced phase modulation, as where M denotes the constant ADM mass, and ω 0 is related to b 0 by c ω −1 0 = 4b 0 e γ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 Ω =ψ with respect to the orbital one ω =φ. The GW half-frequency, directly measurable from the waveform at infinity, reads Ω = ω − 2GMω c 3 ln ω ω 0 + 1 . where we recall that ν = m1m2 m 2 is the symmetric mass ratio. In the following, the results are expressed for quasi-circular orbits in terms of the measurable GW half-frequency Ω through the PN parameter showing that indeed the GW half-frequency Ω differs from the orbital one ω 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 quasi-circular 4.5PN flux is A significant check is to observe that the leading-order terms in the test-mass limit ν → 0 perfectly agree with the results of linear black-hole perturbation theory [124][125][126][127][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 F = c 3 16πG ℓ 0 ℓ m=−ℓ ḣ ℓm 2 . (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][130][131][132][133][134].
In the companion Letter [101] we use the quasi-circular flux F (x) to derive the gravitational frequency chirp and phase at the 4.5PN order. Namely, we apply the energy flux-balance laẇ 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 quasi-circular 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 electromagnetic (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 ∼ 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). 8 C. Gravitational wave modes (2, 2) and (2, 0) Next, we can extract the (ℓ, m) = (2, 2) and (2, 0) physical modes from the newly computed 4PN mass quadrupole radiative moment. Projecting the asympotic metric (2.3) onto the basis of polarizations {+, ×} and the usual basis of spin-weighted spherical harmonics, Y ℓm −2 (following the conventions of [16,17]), one can define the observable gravitational modes h ℓm as (6.14) The result (6.19) is in full agreement with Eq. (4.3a) of [110], obtained from the general expression of non-linear 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].