The Quadrupole Moment of Compact Binaries to the Fourth post-Newtonian Order: I. Non-Locality in Time and Infra-Red Divergencies

With the aim of providing high accuracy post-Newtonian (PN) templates for the analysis of gravitational waves generated by compact binary systems, we complete the analytical derivation of the source type mass quadrupole moment of compact binaries (without spins) at the fourth PN order of general relativity. Similarly to the case of the conservative 4PN equations of motion, we show that the quadrupole moment at that order contains a non-local (in time) contribution, arising from the tail-transported interaction entering the conservative part of the dynamics. Furthermore, we investigate the infra-red (IR) divergences of the quadrupole moment. In a previous work, this moment has been computed using a Hadamard partie finie procedure for the IR divergences, but the knowledge of the conservative equations of motion indicates that those divergences have to be dealt with by means of dimensional regularization. This work thus derives the difference between the two regularization schemes, which has to be added on top of the previous result. We show that unphysical IR poles start to appear at the 3PN order, and we determine all of these up to the 4PN order. In particular, the non-local tail term comes in along with a specific pole at the 4PN order. It will be proven in a companion paper that the poles in the source-type quadrupole are cancelled in the physical radiative type quadrupole moment measured at future null infinity.


I. INTRODUCTION
The theory of gravitational waves (GW) generated by binary systems of compact objects has been developed using perturbative methods in classical general relativity (see [1][2][3][4] for reviews). In particular the post-Newtonian (PN) approximation is a major and widely developed technique for computing analytically the dynamics and GW emission of compact binaries. The state of the art on the conservative dynamics are the equations of motion of point-particle binaries at the 4PN (fourth-post-Newtonian) order [5][6][7][8][9][10][11][12][13][14]. Recently progresses have been made so that the equations of motion are now determined up to two unknown parameters at 5PN order, and up to six unknowns at 6PN order [15,16]. Currently the field is also evolving thanks to new methods coming from effective field theory and scattering amplitudes, naturally combined with the classical post-Minkowskian approximation [17].
This paper is concerned with the GW emission aspect, which is directly related to the data analysis of GW detectors. Here, the state of the art is the 3PN approximation in the waveform, beyond the Einstein quadrupole formula [18][19][20][21][22]. Actually, the flux and orbital phase evolution due to gravitational radiation are known to the 3.5PN order [23,24], as well as the dominant mass-quadrupole mode (ℓ, m) = (2, 2) [25,26], the current-quadrupole mode (2, 1) [27] and the mass-octupole ones (3,3) and (3, 1) [28]. The gravitational flux and mass quadrupole were confirmed at the 2PN order by means of effective field theory techniques [29].
Extending the GW emission up to the 4PN (and even 4.5PN [30]) order is the target of the present program. A central part of this program is of course the control of the mass-type quadrupole moment of the system with the 4PN precision. This computation faces subtle issues regarding the choice and proper use of regularization schemes, both for the ultraviolet (UV) and infra-red (IR) divergences. Recently, a preliminary calculation of the source mass-type quadrupole moment of compact binaries (of spinless bodies) at the 4PN order has been tackled [31]. In this calculation, the UV divergences, appearing because of the pointlike structure of the source (modelling compact objects with negligible internal structure), were properly treated by means of the powerful dimensional regularization. However, the IR ones were regularized with the Hadamard partie finie (PF) regularization procedure [32]. Furthermore, in the calculation of Ref. [31], the non-local-in-time contributions, due to retarded correlations over arbitrarily large time spans in the dynamics of the source [33], were neglected.
In the present paper and the next one [34], we complete the derivation of the 4PN source mass-quadrupole moment of compact binaries. More precisely, our goals are two-fold: 1. To derive the non-local (in time) effect in the source mass quadrupole moment which is due to the radiation modes associated with propagating tails at infinity. This effect is the analogue of the one occurring in the conservative equations of motion and the Lagrangian/Hamiltonian at the 4PN order [5,7,10]; 2. To compute all the contributions to the quadrupole moment due to the IR divergences, implementing a dimensional regularization scheme rather than the Hadamard PF scheme adopted in [31]. Such procedure leads to the appearance of specific IR poles which start to arise at the 3PN order and play a crucial role at the 4PN order.
In the follow-up paper [34], we investigate the fate of these IR poles in the radiative-type quadrupole moment, which represents the actual observable moment at future null infinity, up to the 4PN order. We will find that they are exactly canceled by radiation contributions due to non-linear propagation effects (the most important ones being coined as "tails-oftails" and "tails-of-memory") and that we can therefore safely define a three-dimensional "renormalized " mass quadrupole moment at the 4PN order, which will constitute a basic ingredient in the construction of 4PN GW templates. The central object investigated in this work is thus the source-type mass-quadrupole moment, defined in d spatial dimensions by the expression [27,31,35] [3] . (1.1) We let the reader refer to [31] for a comprehensive review of the definitions and properties of the quantities entering the source quadrupole moment, as well as its computation using the IR Hadamard and UV dimensional regularizations. Let us just emphasize a few points. The main quantity over which the source quadrupole integrates is the pseudo stress-energy tensor τ µν in harmonic coordinates composed of a matter and a gravitational part, The overbar denotes the formal PN expansion; superscript parenthesis (n) denote time derivatives; the hat refers to the symmetric-trace-free (STF) product, e.g.x ija ≡ STF[x i x j x a ]; PF B=0 denotes the Hadamard partie finie with regulator (r/r 0 ) B and associated length scale r 0 (B ∈ C); the characteristic dimensional regularization length scale is ℓ 0 .
where T µν is the matter stress-energy tensor and Λ µν the non-linear gravitational source term of the Einstein field equations. It enters Eq. (1.1) through the following PN-expanded source densities: The PN expansions implicit in (1.1) read as (with Γ the Euler function) (1.4) The partie finie procedure PF B=0 comes from the matching between near zone and exterior zone. It is crucial for the proper definition of the multipole moments in 3 dimensions [36]. The derivation of the equations of motion [8,9] indicated that the PF operator should be kept even in d dimensions. It was also shown there that, in d dimensions, the limit B → 0 is finite in physical quantities (without pole 1/B); so, in the end the whole procedure is equivalent to the usual dimensional regularization. We had however to keep the PF operator as explicit in the expression of the quadrupole moment (1.1). Here, similarly to the equations of motion [8,9], we apply first the PF process when B → 0 on the d dimensional expression and, second, the usual dimensional regularization when ε ≡ d − 3 → 0. We call this mixed regularization the "Bε" regularization. Note that the last term in (1.1), proportional to both B and d − 3, will be shown to play no role with this particular regularization. This paper is organized as follows. In the next Section II, we derive the non-local in time part of the source quadrupole, relying on results derived in [8]. We then perform the proper IR dimensional regularization for all the various categories of terms composing the quadrupole in Section III. The final result is presented in Section IV in the form of a pole followed by a finite part contribution, which will be the starting point for the final renormalization of the quadrupole in the next paper [34]. Appendix A present technical formulas generalizing the Riesz formula in d dimensions. Appendix B gives the expression of the local part of the IR shift coming from the 4PN equations of motion.

II. NON-LOCALITY IN TIME OF THE SOURCE QUADRUPOLE MOMENT
Crucial to the completion of the ambiguity-free equations of motion at the 4PN order was the proper inclusion of the tail effect, i.e. the non-local in time back-scattering of emitted gravitational waves, modifying the conservative dynamics of the system at the current time [8]. This effect enters at the 4PN order in the near-zone metric and, as such, plays a key role in the computation of the 4PN source mass quadrupole by introducing a non-local term, together with a pole.
A. The tail effect in the conservative 4PN equations of motion Let us first recall how the tail effect is included in the near-zone metric and Fokker action for the conservative dynamics. In Ref. [8], it was found that the PN-expanded gravitational field in harmonic coordinates (h µν ≡ √ −gg µν − η µν , which is such that ∂ ν h µν = 0) contains the following pieces responsible for tails at the 4PN order: and introduced the short-hand notation for the logarithmic divergence with the associated pole ∝ ε −1 ≡ (d − 3) −1 : whereq ≡ 4π e γ E with γ E the Euler constant, vanishingly small terms O(ε) being neglected. In our set-up and notation, the two and only two series of multipole moments describing the radiation field generated by an isolated source are the so-called canonical moments, M L and S L . Those differ from the source moments I L and J L , e.g. (1.1) and, thus, in principle, the quadrupole moment in (2.1) should rather be viewed as the canonical moment and denoted M ij . However, we take here advantage that M ij and the source moment I ij are equivalent at Newtonian order (and even up to 2PN order).
Still in Ref. [8], we next applied a gauge transformation, at quadratic order, for the particular interaction M × I ij , so designed as to transfer all relevant tail terms in the "00ii" component of the metric. Namely, we posed h ′ µν = h µν + ∂ µ ǫ ν tail + ∂ ν ǫ µ tail − η µν ∂ ρ ǫ ρ tail , with PN-expanded gauge vector given by the following tail pieces (used in [8] but published here for the first time): In the new gauge, the 4PN tail effect is thus entirely described by the single scalar potential h ′ 00ii tail (or, equivalently, by the 00 component of the covariant metric g ′ tail 00 ), which becomes This tail piece in the metric yields the tail term in the conservative Fokker action, which is found to be manifestly symmetric under time reversal, An elegant alternative form is provided by the Hadamard partie finie (Pf) integral where only the Newtonian term in the PN series (1.4) needs to be considered, and where the ellipsis denote other terms that do not participate to the effect at 4PN order. Considering the order of appearance of the tail integrals in the metric, the only terms in the non-linear source Λ µν which may contribute to the effect at the 4PN order are Λ 00 = −h ab ∂ ab h 00 + 1 4 ∂ a h 00 ∂ a h bb and Λ ii = − 1 4 (d−2)∂ a h 00 ∂ a h bb . Now, only the components of the metric that are proportional to 1/c 8 , namely h ab tail , can subsist at this accuracy level. They are directly given by (2.1c) and are obviously trace-free. The tail sector of the effective non-linear source in (2.7) thus reduces to On the other hand, to the lowest order, we have with V being, at leading order, the Newtonian potential in d dimensions. Using also the fact that h ab tail is a function of time only, we find that the 4PN tail contribution in the quadrupole moment reads Injecting the Newtonian potential V =k G m 1 , the term (2.9) can be computed by means of the "generalized Riesz integrals" presented in Appendix A, more precisely Eq. (A2a) with y 2 = 0. Finally, applying the operator PF at B = 0 which reduces here to a simple limit when B → 0, 2 we find at 4PN order Here, the STF quadrupole moment reads I ij = m 1 y ij 1 ; we have also introduced the Newtonian moment of inertia I = m 1 y 2 C. Indirect tail contribution due to a 4PN shift and total tail effect Now, an important point for our purpose is to remark that the gauge transformation that has been applied to the 4PN equations of motion, with gauge vector (2.3), will induce a spatial shift of the particles' world-lines, which has to be taken into account when evaluating the quadrupole moment. Applying such a shift is crucial to ensure the coherence between 2 The calculation boils down to the single simple elementary integral the coordinate systems used in the derivation of the equations of motion and the multipole moments. Denoting the value on particle 1 of the gauge vector as ǫ µ tail 1 ≡ ǫ µ tail [t, y 1 (t)], the spatial shift ζ 1 (t) ≡ y 1 (t) − y ′ 1 (t) will be given in general by . It is clear from the 1/c factors of Eq. (2.3) that only the first term contributes to 4PN order; hence the shift reads 11) and the contribution brought about by this shift in the quadrupole moment at 4PN is Therefore, again denoting by I = m 1 y 2 1 + 1 ↔ 2 + O(c −2 ) the Newtonian moment of inertia, we obtain at 4PN order As the computation carried out in [31] was purely local, none of the non-local effects (neither the direct effect of the tails nor the contribution of the non-local shift) were included in this preliminary result. Combining the two new contributions of this work (2.10) and (2.13), we thus have to add a non-local in time term to the mass quadrupole moment previously computed in [31] which reads explicitly at 4PN [neglecting terms O(ε)] An interesting point is that the contribution of the moment of inertia I(t) cancels out, but there remains a pole. The next section and the following paper [34] will show how this pole finally combines with other poles coming from IR dimensional regularization to disappear from the observable radiative moment. Note finally that I non-loc ij contains both conservative and dissipative effects, together with a purely instantaneous (non-tail) piece. Those will be specified and computed on quasicircular orbits in the companion paper [34].

III. IR REGULARIZATION OF THE SOURCE QUADRUPOLE MOMENT
The quadrupole moment in d dimensions given by (1.1), when PN-expanded using Eqs. (1.4), and after some suitable integrations by part, can be decomposed into four different types of contributions (see App. C of [31] for the exhaustive list of all those terms): 1. Volume terms, where the integrands are made of products of derivatives of elementary potentials. The complete list of potentials is provided in Appendix A of [31]; 3 2. Compact-support terms, where the integrands are proportional to Dirac distributions modeling the compact objects, multiplying products of derivatives of elementary potentials; 3. Surface terms, where the integrands are total spatial derivatives and can be replaced by their expansions at spatial infinity; 4. The "extra" term, which is the last term of (1.1), proportional to B(d−3), and formally zero in Hadamard's sense, in 3 dimensions within a UV dimensional regularization scheme.
Since the calculations concern only the IR bound, we only need the appropriate accurate values for the potentials when expanded at spatial infinity. Interestingly, the integral giving the value of a potential at the location of particles extends up to infinity, thus this applies also to the case of compact terms, which may be affected by the change of regularization [see Eq. (3.22) below]. In Ref. [31], the applied IR regularization scheme was a pure Partie Finie (PF) one as B → 0 (with poles 1/B discarded) on the three-dimensional expression of (1.1). Instead, we resort here to the mixed "Bε" regularization which consists of computing first the limit B → 0 on the d-dimensional expression of the quadrupole (1.1), and only then apply the dimensional regularization when ε → 0, of course keeping track of all the poles 1/ε. By contrast with the pure PF regularization, we expect from the work [8,9] on equations of motion that the first limit B → 0, i.e., performed on the top of dimensional regularization, will be finite in that case (no pole 1/B), which is confirmed by our calculations below.
Each different type of terms is to be dimensionally regularized (following the Bε scheme) with specific techniques, which are exposed in the rest of the section. In a nutshell, we aim at computing where I Bε ij represents the source mass quadrupole properly regularized following the Bε prescription, I Had ij is an abuse of notation for the result of [31], i.e. the source quadrupole computed with dimensional regularization for the UV divergences and the Hadamard regularization for the IR ones. The terms in the right-hand-side (RHS) represent the differences for the four types of terms 1 to 4.

A. Volume terms
The most numerous terms to regularize are the volume terms. Their regularizations consist in two parts: (i) one has to compute the general expression of the difference between the two regularization schemes Bε and Had for each term, and then (ii) inject the accurate d-dimensional values for the potentials expanded at spatial infinity.

Dimensional regularization of volume terms
Let us consider a generic volume term The function F (x, t) is some product of (derivatives of) potentials and somex L ; below, we will drop the time dependence as it plays no role in the regularization procedure. As we investigate the difference between IR regularization schemes, we can restrict the integral to r > R, where R is an arbitrary constant scale, significantly larger than the distances of the particles to the origin, |y A |, so that we do not have to consider the problem coming from the point-particle approximation, already dealt with in [31]. In Ref. [31], we have computed the pure PF regularization where F (d=3) denotes naturally the three-dimensional limit of F obtained by performing the PN iteration in 3 dimensions and given by (3.9). On the other hand, in this work we consider the mixed regularization scheme Bε, Thus, for each volume term, we are to compute the difference of regularization schemes Since, in the limit ε → 0, the complementary integrals over r < R agree with each other, this difference should be independent of the cut-off scale R. We conclude that DV is the proper quantity that we have to add to the volume terms computed in [31].
As it appears at the 4PN order, the functions F we consider admit generic far-zone (or multipolar) expansions in d dimensions when r ≡ |x| → +∞, namely 4 with coefficients that can contain poles ∝ 1/ε, i.e. of the type We have verified that no double poles ∝ 1/ε 2 appear at the 4PN order. The coefficients ψ p,q of the pole are naturally defined with no dependence upon ε, while the ψ (ε) p,q are finite when ε → 0. An important point is that, despite the poles, the 3 dimensional limits of the functions F are finite, as clear from their expressions given in the App. C of [31] and explicitly verified in our computation. This means that (for all p) q ψ (−1) p,q (n) = 0 .
(3.8a) 4 By which we really mean that, for any N ∈ N, we can write where p 0 ∈ N indicates the maximal order of the IR divergence, and q 0 , q 1 ∈ Z represent a finite range of values for q depending implicitly on N .
Furthermore, by posing we get the corresponding logarithmic expansion in 3 dimensions Note that, as a confirmation of the absence of double poles in the d-dimensional expressions, no squared logarithms appeared in the Hadamard computation. Using the relations (3.8) linking the d and 3 dimensional quantities, we obtain where dΩ d−1 denotes the (d − 1)-dimensional solid angle element. As expected, we find that, modulo O(ε)-terms, this difference does not depend on the cut-off scale R. The result (3.10) generalizes that of Eq. (2.11) in [8], to include integrands containing poles. Note that the value q = 1 has to be excluded from the sum in (3.10). This is not an artificial requirement to ensure that the formula is well-defined: it comes naturally out of the derivation. A last use of (3.8) yields the more compact equivalent form (3.11)

Far zone expansion of the potentials in d dimensions
From (3.11), we see that the difference in regularization schemes depends on specific orders in the far-zone expansion of the integrand of the volume terms. As those are made of products of (derivatives of) potentials (defined in App. A of [31]), we need the expansion of individual potentials in d dimensions.
The expressions of the simplest compact support potentials V , V i , K and the superpotentials (defined in Sec. III.B of [31]) in the whole d-dimensional space are already available, so that we have just to expand them when r → +∞. However, the non-compact support potentialŴ ij is also required at 1PN order, andX,R i andẐ ij at Newtonian order, none of those potentials being known in all d-dimensional space. We thus have to compute their asymptotic behaviours by iterating the propagator at infinity and, crucially, add an appropriate homogeneous solution.
The far-zone expansion (or multipolar expansion, indicated by the operator M) of a potential P with source S (assumed to be PN expanded), i.e., such that P = S = S in d dimensions, reads [8] In principle, the first term is built from the PN-expanded retarded propagator −1 R , but we trade it here for the symmetric one −1 = −1 S , since we are merely interested in even orders. The second term of the RHS is a homogeneous solution constructed out of the PN expansion of where we recall thatk ≡ π 1− d 2 Γ( d 2 − 1), and that γ 1−d 2 (z) is the kernel function entering the Green function of the d'Alembertian equation in d dimensions, defined by Eq. (3.3) in [8] and recalled in Eq. (2.3) of the follow-up paper [34]. In other words, S L ⋆ (t, r) is an elementary monopolar homogeneous solution of the d'Alembertian equation, parametrized by the moments (3.14) Note that these particular moments are chosen to be non-STF, with just x L ≡ x i 1 · · · x i ℓ . Performing explicitly the PN expansion, Eq. (3.12) becomes plus odd powers of 1/c which can be ignored for the present purpose.
Computing each of those terms requires different techniques: • Particular solutions. There are no specific issues with the multipolar expansion of the sources; so, computing the first term in (3.15) is straightforwardly done by means of iterations of the d-dimensional "Matthieu" formula [37] ∆ −1 r αn L = r α+2n L (α + d + ℓ)(α + 2 − ℓ) . (3.16) • Homogeneous solutions. The delicate part in the computation of expanded d-dimensional potentials is the control of the homogeneous solutions. Indeed, the source integral S L (u) of (3.14), of non-compact support, is composed of terms like Those can be computed either by means of the generalized Riesz formulae described in Appendix A below, or by implementing the nice method of App. A in [38], relying on the use of prolate spheroidal coordinates. In addition to the integrals (3.17), we had to deal with the cubic sector of the potentialX, which is defined by Eq. (A.4f) in [31] and whose source contains the delicate term This cubic part requires a priori the knowledge of the potentialŴ ij all over the space in d dimensions at Newtonian order, but this has not been calculated yet. We could have done it relying on the generalization of the Fock function g = ln(r 1 + r 2 + r 12 ) in d dimensions, which has been derived in [37]. However, its expression is only given in an integral form, which is not very convenient in practice. Instead, we employed the method of super-potentials [31] in order to replaceŴ ij by the expression of its source, which is now, indeed, known all over the space. Therefore, the asymptotic behaviours of required potentials have been fully determined with the appropriate accuracy.
Nevertheless, these computations are heavy 5 and some consistency checks are required. The first and most stringent one is that they were performed in a "double-blind" fashion. In addition, we investigated the three-dimensional limits of our potentials, confirming that they agree with the asymptotic expressions that were used to compute the three-dimensional surface terms in [31]. We have also verified that the harmonicity relations ∂ µ h µν = 0 hold, up to 1PN order and to the highest achievable order in 1/r. At the 1PN order, those conditions read [31] The cancellation of ∂ µ h µ0 has been checked exactly at Newtonian order, and up to O(r −6 ) at 1PN order. The fact that ∂ µ h µi vanishes has been checked up to O(r −7 ) at Newtonian order, and to O(r −5 ) at 1PN order.

B. Compact terms
An interesting and non-trivial feature of the IR regularization scheme is that it affects the evaluation of the potentials at the location of the particles. This is due to the fact that most of the potentials have a non-compact support source: the non-linear source terms extend towards infinity, so that the value of the potential at, say, y 1 is sensitive to the IR regularization process. (Naturally, this effect does not impact the compact-support potentials, whose sources are proportional to Dirac distributions, and which are thus sensitive to the UV regularization only.) Therefore, the IR regularization scheme affects the compact terms that involve some non-compact potentials evaluated at y 1 or y 2 .
Let us consider a potential P with d-dimensional source S. As in the previous section, we are interested in the IR behaviour of S, hence we will expand it in the far zone. In practice, none of the sources we are interested in develop poles (or equivalently, logarithms in three dimensions); so, we can safely consider that p,q (n, t) , (3.20) where ϕ (ε) p,q has no pole. The source takes the three-dimensional limit Let us first deal with the Newtonian case, i.e., let us compute the difference induced by the change of IR regularization scheme in the Poisson integral evaluated in y 1 , where we have safely replaced x by y 1 in the kernel, as we integrate over the domain r ′ ≡ |x ′ | > R ≫ |y 1 |, where we are free from UV divergences. Expanding the kernel at spatial infinity according to and using the machinery developed for the treatment of volume terms, we find that the difference between the Bε and Hadamard regularization schemes turns out to be (3.24) where we recall that the volume of the (d−1)-dimensional sphere is Ω d−1 = 2π d/2 /Γ(d/2). As expected, the scale R disappears from the difference in regularizations. The q = 0 criterion is a natural consequence of the derivation. The main structural difference with the formula for the volume terms (3.11) is the sum over ℓ. Nevertheless, this sum is finite by virtue of the form of the source term (3.20), since ℓ is bounded by 2 − p 0 .
The formula (3.24), which is merely Newtonian, has to be generalized to higher PN orders. For this purpose, instead of starting from the Poisson integral (3.22), we have to use the d-dimensional propagator. Again, we can restrict ourselves to an integration region r ′ > R ≫ |y 1 | in which we are allowed to replace x by y 1 , so that where the time-dependence of the source is now crucial. By PN-expanding this propagator, it can be seen that the action of the full propagator, in a formal PN sense, up to any PN order, is equivalent to the action of the mere Poisson integral (3.22) but acting on the effective source with x ′ just playing a spectator role. The PN expansion of this effective source is given by  [8] for more details). As expected, the j = 0 term of S eff even is exactly given by S. Note that the PN-odd piece appears to be non-local in d dimensions [8]. Nevertheless, as already stated, we are merely interested in even terms 6 and, thus, will not consider it. For a potential entering the source at a given PN order, we can simply apply the Newtonian result (3.24), but using the effective source S eff even truncated at the appropriate PN order following Eqs. (3.27). The required potentials that developp a non-vanishing difference at point y 1 areŴ 2PN 1 , Z 1PN 1 ,X 1PN 1 ,T N 1 andM N 1 (with the index 1 denoting the value at particle 1 and the superscript the PN order) while the ones of the super-potentials vanish. The calculations reveal that all the non-vanishing corrections in our potentials are connected. Indeed, we found the interesting, but probably not very profound, relations None of the other potentials receives corrections at the required order. These relations are valid up to the O(ε 0 ) order only, as the O(ε) remainders do not play any role in the IR dimensional regularization of the compact-support terms. Note that only the trace of the potentialM ij at Newtonian order (i.e.M ≡M ii ) as well as the potentialX at 1PN order develop poles. Moreover, only the "scalar" sector of the potentials is affected: neitherR i norŶ i are modified at the particles' positions; regarding the tensor potentialsẐ ij andM ij , only their traces are impacted. In addition to those relations, the difference itself can be compactly written in terms of the Newtonian moment of inertia I ≡ m 1 y 2 1 + m 2 y 2 2 as where we recall thatq ≡ 4πe γ E . The regularization induced differences of those potentials yield corrections in the effective massμ 1 (see Eq. (2.17) in [31]), and in a few compactsupport terms that all enter the source mass quadrupole moment at 4PN order.

C. Surface terms
Turning now to the surface terms, we will not calculate the difference between their values obtained in the two regularization schemes, but rather directly compute them in d dimensions and show that they actually vanish. As presented in [31], surface terms appearing in the mass quadrupole are of two kinds: "Laplacian" and "divergence" terms.
• A surface term of Laplacian type reads where G is a d-dimensional function (in practice a product of potentials). As before, we have restricted ourselves to an integral over r > R. Like for the volume terms, the function G we will consider can be expanded near spatial infinity as: allowing the γ (ε) p,q to contain poles. Inserting this expansion into the integral (3.30), performing an integration by parts, using ∆(r Bx L ) = B(B + d + 2ℓ − 2)r B+ℓ−2n L , and dropping the integrated terms that are vanishing by analytic continuation in either B or ε near zero, we are led to (3.32) • A surface term of divergence type reads where H i is a product of potentials or super-potentials. We again have the far-zone expansion where the η (ε)i p,q can contain poles. A similar procedure yields The noteworthy point is that the d-dimensional values of both types of surface terms, in Eqs. (3.32) and (3.35), are non zero only in the special case q = 1. Now, as explicit in their expressions displayed in App. C of [31], the integrands of those terms are made of non-linear products of potentials. As such, the coefficients of their asymptotic expansions all bear q 2, which make them vanish.
More precisely, it is immediate to see, from the Green function ∆ −1 δ (d) (x − y 1 ) = −k r −1−ε 1 /(4π), that the expansions of compact support potentials have q = 1. Now, the sources of non-compact support potentials are made of products of compact support ones, and the iteration of the Poisson integrals cannot reduce the number of q's. The latter fact can be understood from the Matthieu formula (3.16), together with the fact that the expansion of the homogeneous solution bears q = 1, as explicitly shown by Eq. (3.15) where the second term is ∝ ∂ L r 2k+2−d ∝ r 2k−ℓ−1−ε . Thus, all individual potentials V , V i ,Ŵ ij , etc. admit asymptotic expansions with q 1. A similar argument applies to the case of super-potentials, which implies that all non-linear products of (derivatives of) potentials or super-potentials have q 2.
The conclusion is that the surface terms are vanishing in d dimensions, so that we have simply to subtract their Hadamard values from the final result: The argument developed in [31] to discard the contribution of the B(d − 3) piece of I ij [the last term in (1.1)] does not a priori hold with IR dimensional regularization; further investigation is required. We can restrict ourselves to an integration in the far zone r > R, and recast 3) is composed of two pieces: The first one, coming from the stress-energy tensor T µν , involves Dirac distributions and, thus, does not enter our computation, since we take r > R. The second piece, entailing the nonlinear (NL) gravitational interactions Λ ij , is composed, as such, of products of potentials. We expand it as Σ NL ab (x, t) = p,q ℓ qε 0 r p+qε σ (ε)ab p,q (n, t) , (3.38) where the functions σ (ε)ab p,q contain many PN orders and can develop poles. Inserting (3.38) into (3.37), we can perform the radial integration, which yields 5+2k,1 (n, t) . (3.39) Exactly as in the case of the surface terms, only the terms with q = 1 contribute. However, recalling that σ (ε)ab p,q is made of products of potentials, it cannot involve q = 1 terms. We see therefore that, even with the present IR Bε regularization scheme, the "extra" piece does not contribute:

IV. FINAL EXPRESSION OF THE SOURCE QUADRUPOLE MOMENT
Summing up the pieces computed in the previous sections, and adding the local IR shift χ i 1,2 coming from the cancellation of the remaining poles in the conservative equations of motion [8] (as described in Appendix B), we obtain the full dimensionally regularized source mass quadrupole moment at 4PN order, (4.1) Here, I Had ij stands for the end result of [31], 7 then I non-loc ij is the non-local part computed in Sec. II and given by Eq. (2.14), DI ij is the sum of the four contributions (3.1) whose calculation has been detailed above, and δ χ I ij is the contribution of the IR shift of Appendix B. Remember that the non-local sector I non-local ij is composed of both a direct tail contribution (2.10) and a shift contribution [see Eq. (2.13)].
As the complete expression for the regularized source moment is long and not very enlightening nor interesting per se, we do not display it, but rather discuss its most striking feature: the remaining poles. Indeed, the source moment I ij contains poles, which will have to be compensated by the proper treatment in dimensional regularization of the non-linear interactions (such as tails-of-tails) relating the source and radiative moments. This will be proven in the companion paper [34].
While the non-local tail I non-loc ij and shifted pieces δ χ I ij are purely 4PN contributions, the effect of the dimensional regularization DI ij starts already at the 3PN order. It can be written in a convenient form as where we have used a notation similar to (2.2) for the "dressed" pole: In Eq. (4.2), P i is the constant linear momentum and β I = − 214 105 is the coefficient associated with the renormalization of the mass quadrupole moment [39]. The appearance of this coefficient, which is known to be associated with the "tail-of-tail" interaction [23], indicates the soundness of the removal of the poles by the correct treatment in dimensional regularization of the non-linear interactions in the radiative quadrupole.
As for the poles showing up at the 4PN order, once taken the 1PN correction of Eq. (4.2) into account, they can be expressed in the center-of-mass (CoM) frame as where J i|j is the d-dimensional constant angular momentum defined in [27]. Once again, the fact that this pole can be recast as a combination of non-linear interactions (in particular the expected coupling M × I ij × I ij corresponding to "tails-of-memory") is a strong indication that it will be compensated by the proper dimensional-regularization treatment of the tail/memory effects. This will be the topic of the next paper [34].
The standard Riesz formula in arbitrary dimension d (with a, b ∈ C) reads: It is a priori valid when the integral converges, i.e., for a + d > 0, b + d > 0 and a + b + d < 0. However, the result can be extended by analytic continuation everywhere but, possibly, a countable set of parameter values. On the other hand, this formula may be generalized to similar integrals involving additional STF angular factorsx L and various r 2 . Among the most useful ones, we have with ℓ s denoting the usual binomial coefficient, and y 2 A = y 2 A , (y 1 y 2 ) = y 1 · y 2 . Let us sketch for instance the proof of the first one. The starting point consists in rewriting x L as (r 1 n 1 +y 1 ) L ≡ (r 1 n i 1 1 +y i 1 1 ) · · · (r 1 n i ℓ 1 +y i ℓ 1 ), where n i 1 is the unit vector n i 1 = (x i −y i 1 )/r 1 , and expanding the product. After permuting the integration and summation symbols, we get a sum of elementary integrals of the form for a summation index 0 s ℓ, with n L 1 ≡ n i 1 1 · · · n i ℓ 1 . For each integral, the parameters a, b and d are chosen in a domain of C 3 where the convergence is guaranteed. We can always manage to avoid the situation where one of them is an integer. We then express the factors r a+s 1n L 1 as a multi-derivative using the relation ∂ L r α 1 = (−) ℓ∂ 1L r α 1 = 2 ℓ Γ(α/2 + 1) Γ(α/2 − ℓ + 1) with ∂ 1i = ∂/∂y i 1 , valid as long as∂ 1L r α 1 does not involve any distributional contribution, which is indeed the case when α is a non-integer. Next, we commute the derivatives and the integral (taking for instance 2s < d, 2s < −a < d and a + d < −b < d to bypass any convergence issue). We find We now use the standard Riesz formula (A1), which produces a dimensional factor r a+b+d+2s
The most delicate task in the derivation is the symmetrization over the multi-indices L and K of the multi-derivative In this manner, we find that the integral on the left-hand side of Eq. (A8), after the appropriate change of indices, may be put in the form . (A12) Resorting to techniques similar to those employed to compute the previous integrals, we can perform explicitly the sums over s and r, which yields the desired result.