The Quadrupole Moment of Compact Binaries to the Fourth post-Newtonian Order: From Source to Canonical Moment

As a crucial step towards the completion of the fourth post-Newtonian (4PN) gravitational-wave generation from compact binary systems, we obtain the expressions of the so-called"canonical"multipole moments of the source in terms of the"source"and"gauge"moments. The canonical moments describe the propagation of gravitational waves outside the source's near zone, while the source and gauge moments encode explicit information about the matter source. Those two descriptions, in terms of two sets of canonical moments or in terms of six sets of source and gauge moments, are isometric. We thus construct the non-linear diffeomorphism between them up to the third post-Minkowskian order, and we exhibit the concrete expression of the canonical mass-type quadrupole moment at the 4PN order. This computation is one of the last missing pieces for the determination of the gravitational-wave phasing of compact binary systems at 4PN order.

near zone. We refer to [23][24][25] for discussions and details about the matching procedure we employ to link the near zone to the external zone.
It is important to emphasize that the two descriptions in terms of canonical moments {M L , S L } or in terms of source and gauge moments {I L , J L , W L , X L , Y L , Z L } are physically equivalent, i.e. describe the same physical matter source if and only if the canonical moments are related in a precise way to the source and gauge moments up to arbitrary high orders [23]. Hence there is a coordinate transformation linking the two descriptions of the same source, which is a non-linear deformation of the linear gauge transformation parametrized by the gauge moments {W L , X L , Y L , Z L }.
In this paper, we complete a missing step towards the knowledge of the 4PN orbital phasing for compact binary inspiral by obtaining the relation between the canonical mass quadrupole moment M ij and the corresponding source one I ij at the 4PN order. Such relation was known previously at the leading 2.5PN order [23] and at the next-to-leading 3.5PN order [25]. Up to the 3.5PN order, the correction terms in the canonical quadrupole are quadratic in the multipole moments; the cubic corrections start at 4PN, and are thus the aim of the present computation. Let us recapitulate the complete result for the canonical quadrupole moment up to the 4PN order 1 a| i I j a + Z a| i I Here M is the constant (ADM) total mass, W ij and Y ij for instance denote the quadrupoles associated with the series of gauge moments W L and Y L , W and W i are the monopole and dipole of W L , and so on. All those quantities are evaluated at time t except when 1 Hereafter, we denote with a capital letter L a multi-index with ℓ indices, e.g. I L = I i1i2...i ℓ , with angular brackets the symmetric-trace-free (STF) projection, e.g. W a i I j a ≡ STF ij [W ai I ja ]; we systematically use STF harmonics [26][27][28][29]; in the spirit of [22], we use the shorthands J i|L ≡ ε ii ℓ a J aL−1 and Z i|L ≡ ε ii ℓ a Z aL−1 for current type moments; the superscript (n) denotes n time derivatives; (−) n stands for (−1) n and c is the speed of light and G the gravitational constant.
specified otherwise. In addition, the source is considered to be stationary in the remote past so that all time derivatives of the moments vanish before some instant −T . This means in particular that the hereditary integrals are well defined. The expression (1.1) is valid in the center-of-mass frame, for which the mass dipole moment I i vanishes. The 4PN cubic terms are new with this paper, which extends the previously computed 3.5PN quadratic relation [23,25,30]. Note that the analogous relations for the mass octupole and current quadrupole moments (currently known only at the leading 2.5PN order) do not receive 3PN corrections and can be found in [25]. Interestingly, the relation between the canonical and source/gauge moments is not local, as clear from the 4PN tail integrals appearing in the last line of (1.1). The associated scale r 0 is unphysical and should naturally disappear from any physical results, such as the 4PN flux. This will be a stringent test for our computation, which will however have to wait for the complete calculation of the 4PN radiative-type moment, directly observable at future null infinity. Notice however that the tail terms in (1.1) will be zero in the case of circular orbits; see Eqs. (5.7) in [30].
Another worthy remark is that the subtleties arising from the use of an IR dimensional regularization scheme for the mass quadrupole have been completely tackled and solved in [19,20]. We can thus safely perform the computation in three dimensions, using the standard Hadamard regularization scheme.
Finally, after the result (1.1), only one last step remains before getting the 4PN mass quadrupole: the three-dimensional computation of cubic non-linear terms called "tails-ofmemory", entering the relation between the radiative quadrupole moment and the canonical one at 4PN order. This is left to future work.
The plan of this paper is as follows. After reminders about the multipolar-post-Minkowskian (MPM) formalism in Sec. II, we describe the general method for relating the canonical moments to the source and gauge moments up to any post-Minkowskian (PM) order in Sec. III (extending earlier works in [23,25]). Finally Sec. IV is devoted to the practical implementation that led to the result (1.1), together with required formulas for retarded integrals of non-linear source terms. The essential, but technical, near-zone expansion of tail integrals is presented in details in App. A. The verification of our main result (1.1) via an alternative procedure is relegated in App. B.

II. THE MULTIPOLAR-POST-MINKOWSKIAN EXPANSION
The MPM formalism [29,[31][32][33] obtains the general solution of the Einstein field equations outside a matter source in the form of a post-Minkowskian expansion, with each PM coefficient expanded as a formal multipolar series. The gothic metric deviation from the Minkowski metric, h µν ≡ √ −g g µν − η µν where g µν is the inverse of the usual covariant metric, η µν , that of the Minkowski metric, and g ≡ det(g µν ), obeys the Einstein's vacuum field equations in harmonic coordinates, Here is the flat d'Alembertian operator and the gravitational source term Λ µν is at least quadratic in h and its first and second partial derivatives. In this paper, we will reason to any PM order for the general method, but in practical computations, as we are interested in the 3PM interaction, we will only use the leading quadratic ∼ h∂ 2 h + ∂h∂h and sub-leading cubic ∼ h∂h∂h pieces in the source term.

A. The generic MPM algorithm
The "generic" non-linear MPM solution of the field equations is searched in the form of a non-linear expansion in the field perturbation, labeled by the gravitational constant G. Formally, it may be represented to arbitrary high orders by the asymptotic series: The starting point is the most general solution of the linearized Einstein's vacuum equations in harmonic coordinates, h µν gen 1 = ∂ ν h µν gen 1 = 0, which can be written in terms of six sets of symmetric-trace-free (STF) multipole moments {I L , J L , W L , X L , Y L , Z L }, dispatched between a simpler linear solution called "canonical", and a linear gauge transformation, as We employ the shorthand notation ∂ϕ µν for the linear gauge transformation. Denoting functionals of multipole moments by means of capital calligraphic letters, the functional dependence of the two terms in (2.3) are For obvious reasons, the moments {I L , J L } (mass-type I L and current-type J L ) are called the source moments while the moments {W L , X L , Y L , Z L } are the gauge moments. The linear canonical solution, evaluated at field point x and at time t, reads explicitly [27][28][29]34] h 00 where r = |x| represents the radial distance to the origin located in the source. From the harmonic coordinate condition ∂ ν h µν can 1 = 0, the mass monopole I and current dipole J i must be constant, while the mass dipole I i is varying linearly with time. In applications, we choose a center-of-mass frame for which I i = 0. 2 With these exceptions, the moments are arbitrary functions of time encoding the properties of the source. The linearized gauge vector is defined by The gauge moments {W L , X L , Y L , Z L } are arbitrary functions of time without restriction. The MPM construction is defined by induction on the PM order n [29]. Suppose that for some given n 2, one has obtained the first n − 1 PM coefficients h µν gen m , ∀m n − 1. Then, the next order coefficient h µν gen n is constructed as follows. It satisfies h µν gen n = Λ µν gen n ≡ Λ µν n [h gen 1 , · · · , h gen n−1 ] , (2.7a) where the source term Λ µν gen n , being at least quadratic in h, depends only on the previous iterations as indicated. We first construct a particular retarded solution of the wave equation h µν gen n = Λ µν gen n as where −1 ret denotes the usual retarded inverse d'Alembertian operator, and the symbol FP B=0 refers to a specific operation of taking the finite part (FP) in the Laurent expansion when the complex parameter B tends to zero. This finite part involves the multiplication of the source term by the regularization factor (r/r 0 ) B , where we introduce an arbitrary constant length scale r 0 . Such FP operation is required for dealing with source terms made with multipolar expansions like in (2.5) that are singular at the origin r = 0. More generally, the regularized retarded integral operator FP −1 ret ≡ FP B=0 −1 ret (r/r 0 ) B is well defined when acting on a source term admitting an expansion when r → 0 of the form, for any N ∈ N, wheren L = STF[n i 1 n i 2 · · · n i ℓ ], with n i = x i /r, is the STF spherical harmonics of order ℓ, while the sum boundaries a min and p max are integers (depending on the PM order n). For any function in the class (2.9), we have [FP −1 ret f ] = f . In the end of the recursive process, the structure (2.9) turns out to be proved by induction.
Because of the regularization scheme, the object u µν gen n does not satisfy the harmonic gauge condition, but a simple calculation using the fact that the source term is divergenceless, ∂ ν Λ µν gen n = 0, gives w µ gen n ≡ ∂ ν u µν gen n = FP Due to the explicit factor B, this term is non zero only when the integral generates a pole ∝ 1/B in the Laurent expansion as B → 0. In turn, the pole arises only from the singular behaviour of the source term as r → 0, which is of the type (2.9). Furthermore, the coefficient of the pole is necessarily a homogeneous retarded solution of the wave equation, since the source term has no pole, hence w µ gen n = 0. At this stage, we apply the MPM "harmonicity" algorithm to construct from w µ gen n another homogeneous retarded solution, say v µν gen n ≡ V µν w gen n , (2.11) satisfying ∂ µ v µν gen n = −w µ gen n , together with v µν gen n = 0. The formulas specifying the above harmonicity algorithm w µ −→ V µν [w] are given by e.g. Eqs. (2.11)-(2.12) in [33]. Finally, the metric at the PM order n, now satisfying the full Einstein vacuum equations in harmonic coordinates at the PM order n, is naturally defined as h µν gen n = u µν gen n + v µν gen n . (2.12) This construction yields the most general solution h µν gen of the Einstein field equations in the vacuum region external to any isolated matter system [29]. It is obtained in the form of a functional of six sets of multipole moments, (2.13)

B. The canonical MPM algorithm
However, we also know that the general field of an isolated system in GR can be described by two and only two sets of STF multipole moments, which we call the "canonical" moments and denote {M L , S L }. Notably, these moments describe the GW propagation far away from the source, and they parametrize the two GW tensorial modes of GR [29]. The description of the external field of the source in terms of {M L , S L } starts at linear order just by the canonical metric (2.5) instead of the generic metric (2.3), but, in this different set up, the linear approximation is parametrized by the canonical moments M L and S L : (2.14) The canonical MPM algorithm proceeds then by induction over the PM order n 2 exactly in the same way as before, i.e. following the synthetic steps This results in a full non-linear metric, which is a functional of the canonical moments only, and represents as well the most general solution of the Einstein field equations outside the source: The next section presents the general method to obtain the relationships linking the canonical moments {M L , S L } to the set of source and gauge moments {I L , J L , · · · , Z L }, valid in principle to any PM order. We thus apply this general method to the practical case of cubic interaction to derive the 4PN relation (1.1).

A. General method
In Sec. II, we constructed two full non-linear MPM solutions, Eqs. (2.2) and (2.16), which both represent the most general solution of the Einstein field equations in the vacuum region outside the source. Requiring that these two metrics describe the exterior field of the same physical system, we now impose that they are physically equivalent, i.e., just differ by a coordinate transformation. This implies unique relations between the canonical moments {M L , S L } and the source and gauge moments {I L , J L , W L , X L , Y L , Z L }, which can be viewed as two physically equivalent sets of moments for describing the source as seen from its exterior.
We thus look for a coordinate transformation where J ≡ det(∂x ′ /∂x) is the Jacobian of the transformation. Eq. (3.1) immediately follows from the definition h µν = √ −g g µν −η µν and the law of transformation of tensors [36]. Introducing a coordinate shift ϕ µ such that x ′µ = x µ + ϕ µ (x), we can rewrite the statement (3.1), using the non-linear correction δ ϕ h µν can (x) to the metric h µν can (x) induced by the shift, as h µν It is implicit that we work perturbatively in both the vector ϕ λ and the metric h λρ can , so that we have for instance h µν gen (x ′ ) = n 0 ϕ λ 1 · · · ϕ λn ∂ λ 1 · · · ∂ λn h µν gen (x)/n!. To linear order, the correction reduces to ∂ϕ µν = ∂ µ ϕ ν + ∂ ν ϕ µ − η µν ∂ ρ ϕ ρ . Let us then pose, to any order, where Ω µν denotes a functional of ϕ λ and h λρ can , as well as their derivatives, which is at least quadratic and can be computed perturbatively up to any order using Eq. (3.1). The harmonic gauge condition satisfied by both metrics h µν gen and h µν can implies that (as a consequence of the identity ∂ ν ∂ϕ µν = ϕ µ ) Let us now look for the coordinate shift in the form of a full PM expansion series and, conjointly, for the relations between the canonical moments and the source/gauge moments in the same PM form Here M n L and S n L denote some n-th non-linear functionals of the six types of source and gauge moments I K , · · · , Z K (with, say, K = i 1 · · · i k ), which are to be determined: At linear order, as we have seen with Eq. (2.3), by definition of the two linear approximations for the two generic and canonical metrics, where the functional H µν can 1 is explicitly given by (2.5) and the gauge vector is parametrized by the gauge moments {W L , · · · , Z L } [see its expression (2.6)]. Eq. (3.8) means that the relation (3.2) is satisfied at leading order provided that (i) ϕ µ Thus, we see that the looked-for coordinate transformation will be a non-linear deformation of the linear gauge transformation associated with the shift (2.6).
Our recurrence hypothesis will be that all the ϕ µ m 's in (3.5) are known up to a given PM order n − 1, as are the functional relations (3.6) up to the corresponding order. Hence we assume that we have already determined in such a way that the equations (3.2)-(3.3) hold up to order n − 1. This means that, for all m n − 1, Recall that Ω µν m is a non-linear, at least quadratic, functional of the coordinate shift and metric, and is therefore known following our induction hypothesis. Indeed, using (3.3) and (3.8), we get Ω µν 1 ≡ 0. Crucial in our hypothesis is the assumption that the canonical metric depends on the moments so far determined to order n − 1: B. Implementation to nPM order With our recurrence hypothesis, let us see how to the next PM order ϕ µ n together with the functionals M n L and S n L are uniquely determined. We thus want to find ϕ µ n and M n L , S n L such that h µν gen n = h µν can n + ∂ϕ µν n + Ω µν n , (3.13) where h µν gen n and h µν can n are defined precisely by the two MPM constructions in Sec. II, and h µν can n is now a functional of M n L and S n L . First of all, note that Ω µν n depends on ϕ k and h can k for k n − 1 and is already known by our induction hypothesis. Second, apply the harmonic coordinate conditions on Eq. (3.13): this shows that while ϕ µ n is one of our unknowns, its d'Alembertian ∆ µ n ≡ ϕ µ n is already determined as we have [see Eq. (3.4)] 14) The explicit expressions of Ω µν n and ∆ µ n at quadratic and cubic orders are displayed in Eqs. (4.2), (4.7) and (4.8). At this stage, we must use the specific definitions of the generic and canonical metrics defined in Sec. II. Applying the d'Alembertian operator on (3.13), we find that the two source terms Λ µν gen n and Λ µν can n are related by Hence, the particular retarded solutions of the two algorithms, u µν gen n and u µν can n , defined respectively in (2.8) and (2.15), satisfy Next, we introduce a linear-looking gauge transformation with vector defined by the retarded integral of ∆ µ n , as This vector satisfies φ µ n = ∆ µ n but is not yet our looked-for vector ϕ µ n . It a priori differs from it by an homogeneous retarded solution of the wave equation. Thanks to (3.17) we can advantageously rewrite (3.16) as u µν gen n = u µν can n + ∂φ µν n + Ω µν n + X µν n + Y µν n . (3.18) The last two terms are the most interesting: they come from the non-commutation of the finite part of the retarded integral FP B=0 −1 ret with the partial derivative, due to the differentiation of the regularization factor (r/r 0 ) B therein. They are thus given by Observing that, for the class of multipole-expanded functions f we are concerned with [see Eq. (2.9)], the statement (FP −1 ret f ) = f is always correct, we see that X µν n and Y µν n represent the "commutators" of the operators that appear inside the square brackets: Since the differentiation of the regularization factor (r/r 0 ) B produces an extra factor B, the quantities X µν n and Y µν n will be non-zero only when the integral develops a pole ∼ 1/B. In that case, they are necessarily homogeneous (retarded) solutions of the wave equation: X µν n = Y µν n = 0. It is easy to figure out that if X µν n and Y µν n were actually zero, the canonical moments {M L , S L } would simply agree with their source counterparts {I L , J L }. As a result, the non trivial relations between those moments entirely follow from the evaluation of the two quantities X µν n and Y µν n . In our practical calculations, we reshuffle the commutators (3.20) and use the following expressions, exhibiting the explicit factor B in front: Next, we carry on the MPM algorithm by computing the divergence of (3.18). With Eq. (3.14), we obtain w µ gen n = w µ can n + W µ n , where we have posed U µν n ≡ X µν n + Y µν n and W µ n ≡ ∂ ν U µν n . It is in fact necessary and sufficient to apply the harmonicity algorithm V µν only to the divergence of the sum of the two commutators (3.21). We have The final step consists in remarking that H µν n ≡ U µν n + V µν n is not only divergenceless by definition of the harmonicity algorithm V µν , but that it is also a retarded homogeneous solution of the wave equation, since both U µν n and V µν n are separately such. Hence H µν n satisfies the linearized vacuum Einstein field equations, i.e. H µν n = ∂ ν H µν n = 0, to which we know the general solution. Namely, it can be decomposed in a unique way as where H µν can 1 denotes the linearized functional (2.5) of the moments, but computed with certain moments M n L and S n L , and where ∂ψ µν n is some linear-looking gauge transformation. The associated gauge "vector" ψ µ n is parametrized in a unique way by some moments {W n L , X n L , Y n L , Z n L }. The vector ψ µ n represents the homogeneous solution to be added to φ µ n as given by Eq. (3.17) in order to recover Eq. (3.13) with the shift ϕ µ n ≡ φ µ n + ψ µ n . S n−1 L } by the more accurate, newly determined set {M n L , S n L }, modulo higher-order PM terms at least of order ∝ G n+1 , which we can discard to order n. Hence we have proved that the nPM contribution to the canonical moments {M L , S L } is determined, as is the nPM piece of the coordinate shift ϕ µ , and our recurrence hypothesis is verified at the next order n. The practical implementation at quadratic and cubic orders is described in Sec. IV. As a verification, we have also followed an alternative approach, presented in App. B.

C. Extraction of physical multipole moments
A straightforward procedure permits reading off the expressions of the n-th order corrections to the canonical moments {M n L , S n L } from Eq. (3.23), where the left-hand side H µν n = U µν n + V µν n follows from the algorithm (3.22). In fact, as we have seen that U µν n = X µν n + Y µν n is a retarded vacuum solution of the wave equation, we can directly, and uniquely, give the desired moments as functions of the ten sets of retarded STF moments composing U µν n . We thus pose and follow the steps (3.22), successively computing W µ n , V µν n and H µν n , which is finally put into the form (3.23) on which we read off the physical moments For completeness, we also give the corrections to the gauge moments composing the gauge vector as ψ µ n :

IV. PRACTICAL IMPLEMENTATION
Let us apply the previously described procedure to determine the quadratic and cubic corrections M 2 ij and M 3 ij to the mass-type quadrupole moment M ij at the 4PN order [see Eqs. (3.6)]. Previous investigations [23,25] focused on quadratic interactions and determined the mass quadrupole at leading order 2.5PN and sub-leading order 3.5PN; such corrections are recalled in Eq. (1.1). However, in order to obtain the 4PN correction, we need to derive the cubic interactions. A preliminary dimensional analysis shows that at 4PN order and in the center-of-mass frame (where the mass dipole I i is vanishing), the only multipole interactions between the source and gauge moments are cubic and necessarily of the three types: where M is the constant ADM mass, I ij is the source mass quadrupole moment, and where the monopole W and the two quadrupoles W ij and Y ij are gauge moments. Note that out of the center-of-mass frame, a large number of additional interactions would appear, such . , but those are not needed in concrete applications.

A. Controlling the cubic source terms
The first step towards the practical determination of the cubic couplings is to successively construct the quadratic and cubic quantities Ω µν n and ∆ µ n that enter Eqs. (3.21) for n = 2, 3. To quadratic order we have This is valid "on-shell", as we have used the facts that h µν can 1 = ϕ µ 1 = 0, which hold at linear order. One can directly verify that [see Eq. (3.14)]  [32,37]. The quadratic couplings between source and gauge moments have been computed in [23,25]. Notably, this yields [see (3.27a)] and we have already given the 3.5PN contribution in (1.1). What remains to be computed, for insertion into the cubic quantities Ω µν 3 and ∆ µ 3 , are the quadratic couplings contributing to the quadratic coordinate shift ϕ µ 2 . As we have seen in (3.24), the shift is composed of two parts: the first one, φ µ 2 , has been defined in Eq. (3.17), and is computed by the usual techniques (see for instance the Appendix A of [33]). The second part, ψ µ 2 , is extracted from (3.23) as quadratic corrections to the four types of gauge moments; the general result has been provided in Eqs. (3.28). Working out ϕ µ 2 = φ µ 2 + ψ µ 2 for the needed quadratic interactions M × W, M × W ij , M × Y ij and I ij × W, we find Note the presence of non-local tail integrals, involving Q m (y), the Legendre functions of the second kind, here defined with a branch cut from −∞ to 1 and related to the Legendre polynomials P m (y) by We inject the expressions (4.5) together with the expressions for the canonical metric, notably h µν can 2 corresponding to the interaction M×I ij which is also non-local and given by (B3) in [32], into the cubic Ω µν 3 and ∆ µ 3 that define (3.21). For convenience, we split Ω µν 3 into quadratic-type and purely cubic interactions: Ω µν 3 ≡ Ω µν 12 + Ω µν 21 + Ω µν 111 , where Recall that our computations are done on-shell, using the wave equations satisfied at linear order, h µν can 1 = ϕ µ 1 = 0, and quadratic order, h µν can 2 = Λ µν can 2 . An important point is that Ω µν 111 only contains h can 1 × ϕ 1 × ϕ 1 and ϕ 1 × ϕ 1 × ϕ 1 terms, but there is no h can 1 × h can 1 × ϕ 1 sector. This fact allows us to discard it entirely in the practical implementation, since at the 4PN order we have only the three cubic interactions (4.1) which contain at most one gauge moment. Similarly, we have for the coordinate shift ∆ µ We have inserted ϕ µ 2 = ∆ µ 2 , which comes from (3.17), and used the fact that ψ µ 2 = 0. The latter quantities satisfy ∂ ν Ω µν 3 + ∆ µ 3 = 0, consistently with the divergencelessness of the cubic source, and in fact, even separately, At this stage, we control the cubic source terms that are required to evaluate the "commutators" X µν 3 and Y µν 3 defined by the formulas (3.21). Using canonical relations between the Legendre functions, we find that the terms to be computed fall into two and only two classes: Here, the functions F and G represent some products of (source or gauge) multipole moments, with the function G which can be constant when dealing with the ADM mass. The second type of term I tail integrates over a tail term which comes from the tail present in the quadratic metric for the interaction M×I ij , see (B3) in [32], and those which appeared in the coordinate shift [see Eqs. Thanks to the factors B or B 2 , we know that the integrals (4.10) will depend only on the behaviour of the source terms when r → 0. 3 Indeed the regularization process FP B=0 has been introduced to cope with the singular behaviour of the source when r → 0, and hence only that limit can generate poles 1/B or 1/B 2 which can compensate the factors B, B 2 and lead to a finite part when B → 0. Therefore we are entitled to restrict ourselves to a ball of finite and even infinitesimal size (say r < ǫ) and replace the source terms in (4.10) by their formal Taylor expansion when r → 0 (i.e., the near zone or PN expansion r/c → 0). Once the near zone expansion of the source is known the integration can be performed with standard techniques. The result for the simpler case I inst (with purely instantaneous source term) was already given in Eq. (A.18) of [33]. It is non-zero only when p ℓ + 3 and b = 1 (only a simple pole 1/B can appear in this case), and reads (4.11) Next, we deal with the more difficult integral I tail . By the previous argument we can replace the function G(t − r/c) by its formal Taylor expansion when r → 0. The main problem is therefore the control of the expansion series when r → 0 of the tail integral F m ≡ +∞ 1 dy Q m (y) F (t − yr/c). We devote the App. A to this not-so-easy question and state here the general result for this expansion: Here, the symbol r/c→0 = indicates the formal asymptotic expansion, the coefficients c m j and β m i are defined in (A5) and (A16), and H n denotes the usual harmonic number. Thus, we have to multiply Eq. (4.12) by the Taylor expansion of G(t − r/c) and integrate term by term. Note that all those formal expansions are convergent, due to the stationary of the source in the remote past. As we see from the structure of the expansion (4.12), the final integration will boil down to the control of just one type of term, However, the function H(t) here can be either an instantaneous function or a non-local tail integral of the type given in (4.12). Furthermore, still because of the tail term we must add the case where there is logarithm ln r in the source, see again (4.12). Thus, we consider (4.13) with a = 0 or 1.
We give a few details on the calculation of J . Writing (4.13) in ordinary three-dimensional form we have, since as we said the integration is limited to an infinitesimal ball r < ǫ, (4.14) Using the formal STF expansion of the multipolar factor when r ′ ≡ |x ′ | → 0, together with the angular integration performed using Eq. (A29a) in [29], we end up with radial integrals of the type Here the r ′ = 0 boundary vanishes by analytic continuation on B. With the factor B b in front the latter integral is zero unless there is a pole, and the pole can come only when p is of the form p = ℓ + 3 + 2j (where j ∈ N). Hence the results are immediate. When a = 0 (source without ln r) the integral is zero when b = 2 as there is only a simple pole 1/B in this case: One can check that this is perfectly consistent with the result for I inst in (4.11), in the sense that one can Taylor expand the source term of I inst when r → 0 and recover the same result by applying Eq. (4.17a) on each term of the Taylor series. By contrast, when there is a ln r the integral is also non-zero when b = 2 because of the presence of a double pole 1/B 2 : (4.17b) Finally with those results in hand, we can implement all the terms up to cubic non-linear order, using the xAct library of the Mathematica software [38]. Such computation ends up with the final result at the 4PN level already recapitulated in Eq. (1.1). In the App. B, an alternative procedure is described which permitted checking it independently.  This appendix is devoted to the determination of the near zone (or PN) expansion when r/c → 0 of the tail integral that enters the source of the integral I tail (4.10b). This will permit justifying the claim of Eq. (4.12) concerning the structure of the near zone expansion, and provide the explicit coefficients we need for the practical computation. Thus we look for the expansion of If an explicit expression of the Legendre function of the second kind Q m (y) is given in Eq. (4.6), we preferably use here the general expression of the Legendre function for generic µ ∈ C \ {−1, −2, · · · } in terms of the hypergeometric function F ≡ 2 F 1 : with |y| > 1, | arg y| < π. This representation of the Legendre function reads explicitly, for any real argument such that y > 1, It can naturally be regarded as an asymptotic expansion when y → +∞ on the real axis. However, we stress that it is valid as soon as y ∈]1, +∞[, as the series representation of the hypergeometric function is well defined in that case. Defining τ ≡ yr/c one can then rewrite (A1) as the following series where and the coefficients are just given by (A3b) in the particular case where m ∈ N, i.e., Upon the expression of F j m we perform a series of integrations by parts in order to increase the power of τ until we reach a logarithm of τ . The all-integrated terms contain the functions F (n) (t − r/c) that we replace by their Taylor expansions when r → 0. The last integral containing ln τ is split according to . In the first integral, from 0 to r/c, we are allowed to expand the integrand when τ → 0, since by definition r/c → 0 for the PN expansion. Then the second integral from 0 to +∞ just gives a tail integral in the ordinary sense. Therefore, we obtain the following asymptotic expansion when r/c → 0: Very importantly, the value i = m + 2j is to be excluded from the first summation. The last term is the tail integral, where H n denotes the usual harmonic number. Resumming on j yields We already see the type of structure claimed in Eq. (4.12). The coefficients α m i in the first sum (corresponding to instantaneous terms) are still at this stage given by an infinite series: Despite that, the result (A7) can be dealt with as it is in practical calculations. However, as it turns out the coefficient (A8) can be resummed in analytic closed form, and this will yield a more interesting and powerful expression of the near zone expansion of F m . To obtain such form of α m i we take advantage of the fact that the coefficients c m j can be generalized to any µ ∈ C \ {−1, −2, · · · } through the general definition of the Legendre function, see Eqs. (A3). Therefore, one can extend the definition of α m i to any generic value of µ by posing where the coefficients c µ j are now given by (A3b). We assume that µ is non-integral, so that the condition i = m + 2j is no longer necessary and has been dropped. Now, from the very definition of the coefficients c µ j in (A3), the fact that the hypergeometric series is absolutely convergent for |y| > 1, and using the identity +∞ 1 dy y i × y −1−µ−2j = (µ + 2j − i) −1 valid for ℜ(µ) > i, we prove that the coefficient α µ i is actually given by Furthermore, we know that [39] Hence we arrive at the closed-form expression where i k is the usual binomial coefficient. Its validity may be extended to all non-integral values of µ by analytic continuation. Still this expression is to be connected to the actual result (A8) we search for, as this result excludes the value i = m + 2j from the summation. However, posing µ = m + ε one can also substract and then re-add the contributions i = m + 2j in (A9); in this way we rewrite α m i as the following limit when ε → 0: In the case i = m + 2j an explicit pole 1/ε has to be added, and which should cancel the pole present in α m+ε m+2j so that the limit is finite. Furthermore since the coefficient c m+ε j is finite when ε → 0, with limit c m j given by (A5), we can expand it to order ε and obtain Finally, the last step is to use the closed-form expression (A12) for α m+ε i . The extra terms in (A14) are easily computed with the help of (A3b) and (A5). We verify that indeed the limit ε → 0 is finite and obtain where the new coefficient β m i reads explicitly with the convention that i k = 0 whenever i < k. Hence our final result for the near zone expansion of F m reads with the notations δU µν quad n = δX µν quad n + δY µν quad n and n−1,1 + η µν ∆ i 1,n−1 + ∆ i n−1,1 ; (B5c) the term δV µν quad n is computed from the harmonicity algorithm δV µν quad n = V µν [δW quad n ], with δW µ quad n ≡ ∂ ν δU µν quad n . The source of the coordinate shift δφ µ quad n , which involves the n − 1 order piece of the metric h µν can n−1 , is of the same undesirable type as δΛ µν quad n . However, in the current procedure, the canonical moments will be read off from the gravitational waveform, which is not sensitive to linear-looking gauge transformations. On the other hand, although the commutator contributions do contain terms proportional to h µν can n−1 or ϕ µ n−1 , their brute force integration is not required. Thus, it will be possible to determine M L n and S L n without integrating any such term. Once δh µν quad n is obtained, we compute the rest of the metric h µν rest n ≡ h µν gen n − δh µν quad n by solving the equation h µν rest n = δΛ µν rest n ≡ Λ µν can n − δΛ µν quad n , with the help of the standard MPM algorithm. By construction, the source term δΛ µν rest n does not depend on h µν can n−1 nor ϕ µ n−1 and is thus free of the difficult contributions we wanted to avoid. At this stage, the n-th order waveform may be built from h µν gen n = Gh µν gen 1 + · · · + G n h µν gen n , by taking the limit R → +∞ for constant asymptotically null time U = T − R/c, which corresponds to "radiative" coordinates (T = t − 2GM/c 3 ln(r/r 0 ), R = r) [5]. In those coordinates, the leading order contribution to the metric, H µν gen n , admits the same expression as the one in harmonic gauge, but with the logarithms ln r effectively replaced by ln r 0 , and with (t, r) replaced by their radiative counterparts (T, R); see, e.g., Ref. [21]. This yields in particular [up to terms O(R −2 )] H µν gen n = h µν rest n + ∂δφ µν quad n ln r→ln r 0 + Ω µν 1,n−1 + Ω µν n−1,1 + δU µν quad n + δV µν quad n , where ∂δφ µν quad n is the linear gauge transformation associated with δφ µ quad n . The n-th order waveform h TT n ij is then the transverse trace-free projection of the 1/R term in H µν gen n . As the TT projection of linear gauge transformations vanishes at order 1/R, the vector δφ µ quad n is actually not required. The result for h rad n ij is a certain functional of the source/gauge moments: h rad n ij = H rad n ij [I L , J L , · · · , Z L ] , and for the canonical moments, we must also have, at the same time: The expressions of M n L and S n L are finally found by guess work. We assume they are sums of terms involving the source/gauge moments, with consistent index structures and physical dimensions, but arbitrary coefficients. Those are fixed by identifying Eq. (B8a) with the outcome that ensues from inserting our ansatz into Eq. (B8b). Of course, if we wish to iterate the process to the next order n + 1, we will eventually need to tackle the difficult integrals of source terms containing h µν can n−1 , which arise in the calculation of h µν can n and δφ µ quad n . Nonetheless, we have managed to push this step to the very end. This strategy is particularly relevant to determine the canonical moments at cubic order for two reasons: (i) The latter task does not demand computing h µν can 3 nor δφ µ 3 , which means that all retarded integrals we have to consider are sourced by functions of h µν can 1 , ϕ µ 1 , or their derivatives; (ii) The corresponding integrands have the form f (t − r/c)n L /r k , with k ∈ N \ {0, 1}, whose finite part retarded integral are explicitly known [33]. In practice, we build the cubic source Λ µν 3 , subtract the "difficult" part δΛ µν 3 , and apply the MPM algorithm to the rest, which leads to h µν rest . At last, we compute the commutators δU µν quad 3 with the method developed in Sec. IV B, from which we can infer δV µν quad 3 . The cubic waveform follows from the effective metric h µν eff 3 = Gh µν gen 1 + G 2 h µν gen 2 + G 3 [h µν rest 3 + Ω µν 12 + Ω µν 21 + δU µν 3 ] .
Our final result (1.1) for the canonical quadrupole moment was obtained following the general method in Secs. III A-III B, and has then been entirely checked using this approach.