$$ \bar{h}^{\alpha\beta}(\mathbf{x}, t, c) = \sum_{m=-2}^{+\infty} \frac{1}{c^m} \bar{h}_m^{\alpha\beta}(\mathbf{x}, t; \ln c). \qquad (142) $$

The *m*-th post-Newtonian coefficient is naturally the factor of the *m*-th power of $1/c$. However, we know from restoring the factors *c*'s in Theorem 4 [see Eq. (87)], that the post-Newtonian expansion also involves powers of the logarithm of *c*; these are included for convenience here into the definition of the coefficients $\bar{h}_m^{\alpha\beta}$.<fnref n="34" /> For the stress-energy pseudo-tensor appearing in Eq. (118) we have the same type of expansion,

$$ \bar{\tau}^{\alpha\beta}(\mathbf{x}, t, c) = \sum_{m=-2}^{+\infty} \frac{1}{c^m} \bar{\tau}_m^{\alpha\beta}(\mathbf{x}, t; \ln c), \qquad (143) $$

with the expansion starting with a term of order $c^2$ corresponding to the rest mass-energy ($\bar{\tau}^{\alpha\beta}$ has the dimension of an energy density). As usual we shall understand the infinite sums such as (142)–(143) in the sense of formal series, i.e., merely as an ordered collection of PN coefficients. Because of our consideration of regular extended matter distributions the post-Newtonian coefficients are smooth functions of space-time, i.e., $\bar{h}_m^{\alpha\beta}(\mathbf{x}, t) \in C^{\infty}(\mathbb{R}^4)$.

Inserting the post-Newtonian ansatz (142) into the harmonic-coordinates Einstein field equation (50)–(51) and equating together the powers of $1/c$, results is an infinite set of Poisson-type equations ($\forall m \ge 2$),

$$ \Delta \bar{h}_m^{\alpha\beta} = 16\pi G \bar{\tau}_{m-4}^{\alpha\beta} + \delta_t^2 \bar{h}_m^{\alpha\beta}, \qquad (144) $$

where the second term comes from the split of the d'Alembertian operator into a Laplacian and a second time derivative: $\Box = \Delta - \frac{1}{c^2} \delta_t^2$ (this term is zero when $m=2$ and 3). We proceed by induction, i.e., we work at some given but arbitrary post-Newtonian order *m*, assume that we succeeded in constructing the sequence of previous coefficients $\bar{h}_m^{\alpha\beta}$ ($\forall p \le m-1$), and from that show how to infer the next-order coefficient $\bar{h}_m^{\alpha\beta}$.

To cure the problem of divergencies we introduce a generalized solution of the Poisson equation with non-compact support source, in the form of an appropriate *finite part* of the usual Poisson integral obtained by regularization of the bound at infinity by means of the specific FP regularization process. For any source term like $\bar{\tau}_m^{\alpha\beta}$, we multiply it by the regularization factor $\tilde{r}^B$ already extensively used in the construction of the exterior field, thus $B \in \mathbb{C}$ and $\tilde{r} = r/r_0$ as given by (72). Only then do we apply the usual Poisson integral, which therefore defines a certain function of *B*. The well-definedness of that integral heavily relies on the behaviour of the integrand at the bound at infinity. There is no problem with the vicinity of the origin inside the source because of the smoothness of the pseudo-tensor. Then the latter function of *B* generates a (unique) analytic continuation down to a neighbourhood of the value of interest $B=0$, except at $B=0$ itself, at which value

<footnote n="34">³⁴ For this argument we assume the validity of the matching equation (115) and that the post-Minkowskian series over $n=1, \dots, +\infty$ in Eq. (87) has been formally summed up.</footnote>
it admits a Laurent expansion with multiple poles up to some finite order (but growing with the post-Newtonian order $m$). Then, we consider the Laurent expansion of that function when $B \to 0$ and pick up the finite part, or coefficient of the zero-th power of $B$, of that expansion. This defines our generalized Poisson integral:

$$ \mathrm{FP}_{B=0} \Delta^{-1} [\tilde{r}_m^B \bar{\tau}_m^{\alpha\beta}] (\mathbf{x}, t) = -\frac{1}{4\pi} \mathrm{FP}_{B=0} \int \frac{d^3\mathbf{x}'}{|\mathbf{x} - \mathbf{x}'|} \tilde{r}'^B \bar{\tau}_m^{\alpha\beta}(\mathbf{x}', t). \quad (145) $$

The integral extends over all three-dimensional space but are endowed with the latter finite-part regularization. The main properties of this generalized Poisson operator is that it solves the Poisson equation,

$$ \Delta \left( \mathrm{FP}_{B=0} \Delta^{-1} [\tilde{r}_m^B \bar{\tau}_m^{\alpha\beta}] \right) = \bar{\tau}_m^{\alpha\beta}, \quad (146) $$

and that the solution $\Delta^{-1}\bar{\tau}_m^{\alpha\beta}$ owns the same properties as its source $\bar{\tau}_m^{\alpha\beta}$, i.e., the smoothness and the same type of behaviour at infinity, as given by Eq. (116). Similarly, we define the generalized iterated Poisson integral as

$$ \mathrm{FP}_{B=0} \Delta^{-k-1} [\tilde{r}_m^B \bar{\tau}_m^{\alpha\beta}] (\mathbf{x}, t) = -\frac{1}{4\pi} \mathrm{FP}_{B=0} \int d^3\mathbf{x}' \frac{|\mathbf{x} - \mathbf{x}'|^{2k-1}}{(2k)!} \tilde{r}'^B \bar{\tau}_m^{\alpha\beta}(\mathbf{x}', t). \quad (147) $$

The most general solution of the Poisson equation will be obtained by application of the previous generalized Poisson operator to the right-hand side of Eq. (144), and augmented by the most general homogeneous solution of the Poisson equation. Thus, we can write

$$ \bar{h}_m^{\alpha\beta} = 16\pi G \mathrm{FP}_{B=0} \Delta^{-1} [\tilde{r}_m^B \bar{\tau}_m^{\alpha\beta}] + \partial_t^2 \mathrm{FP}_{B=0} \Delta^{-1} [\tilde{r}_m^B \bar{h}_m^{\alpha\beta}] \\ + \sum_{\ell=0}^{+\infty} \mathcal{B}_m^{\alpha\beta}(t) \hat{x}_L. \quad (148) $$

The last term represents the most general solution of the Laplace equation that is regular at the origin $r = 0$. It can be written in STF guise as a multipolar series of terms of the type $\hat{x}_L$, and multiplied by arbitrary STF-tensorial functions of time $m$ $\mathcal{B}_m^{\alpha\beta}(t)$. These functions will be associated with the radiation reaction of the field onto the source; they will depend on which boundary conditions are to be imposed on the field at infinity from the source.

It is now trivial to iterate the process. We substitute for $\bar{h}_m^{\alpha\beta}$ in the right-hand side of Eq. (148) the same expression but with $m$ replaced by $m-2$, and similarly come down until we stop at either one of the coefficients $\bar{h}_0^{\alpha\beta} = 0$ or $\bar{h}_1^{\alpha\beta} = 0$. At this point $\bar{h}_m^{\alpha\beta}$ is expressed in terms of the previous $\bar{\tau}_m^{\alpha\beta}$'s and $\mathcal{B}_m^{\alpha\beta}$'s with $p \le m-2$. To finalize the process we introduce what we call the operator of the “instantaneous” potentials and denote it by $\Box_{\mathrm{inst}}^{-1}$. Our notation is chosen to contrast with the standard operator of the retarded potentials $\Box_{\mathrm{ret}}^{-1}$ defined by Eq. (71). However, beware of the fact that unlike $\Box_{\mathrm{ret}}^{-1}$ the operator $\Box_{\mathrm{inst}}^{-1}$ will be defined only when acting on a post-Newtonian series such as $\bar{\tau}^{\alpha\beta}$. Indeed, we pose
$$ \mathrm{FP}_{B=0} \Box_{\mathrm{inst}}^{-1} [\tilde{r}^B \bar{\tau}^{\alpha\beta}] \equiv \sum_{k=0}^{+\infty} \left(\frac{\partial}{c\hat{\partial}t}\right)^{2k} \mathrm{FP}_{B=0} \Delta^{-k-1} [\tilde{r}^B \bar{\tau}^{\alpha\beta}], \quad (149) $$

where the k-th iteration of the generalized Poisson operator is defined by Eq. (147). This operator is “instantaneous” in the sense that it does not involve any integration over time. It is readily checked that in this way we have a solution of the d’Alembertian equation,

$$ \Box \left( \mathrm{FP}_{B=0} \Box_{\mathrm{inst}}^{-1} [\tilde{r}^B \bar{\tau}^{\alpha\beta}] \right) = \bar{\tau}^{\alpha\beta}. \quad (150) $$

On the other hand, the homogeneous solution in Eq. (148) will yield by iteration an homogeneous solution of the d’Alembertian equation that is necessarily regular at the origin. Hence it should be of the *anti-symmetric* type, i.e., be made of the difference between a retarded multipolar wave and the corresponding advanced wave. We shall therefore introduce a new definition for some STF-tensorial functions $A_L^{\alpha\beta}(t)$ parametrizing those advanced-minus-retarded free waves. It is very easy to relate if necessary the post-Newtonian expansion of $A_L^{\alpha\beta}(t)$ to the functions $\hat{m}_L B_L^{\alpha\beta}(t)$ previously introduced in Eq. (148). Finally the most general post-Newtonian solution, iterated *ad infinitum* and without any divergences, is obtained into the form

$$ \bar{h}^{\alpha\beta} = \frac{16\pi G}{c^4} \mathrm{FP}_{B=0} \Box_{\mathrm{inst}}^{-1} [\tilde{r}^B \bar{\tau}^{\alpha\beta}] \\ - \frac{4G}{c^4} \sum_{\ell=0}^{+\infty} \frac{(-1)^\ell}{\ell!} \hat{\partial}_L \left\{ \frac{A_L^{\alpha\beta}(t-r/c) - A_L^{\alpha\beta}(t+r/c)}{2r} \right\}. \quad (151) $$

We shall refer to the $A_L^{\alpha\beta}(t)$'s as the *radiation-reaction* functions. If we stay at the level of the post-Newtonian iteration which is confined into the near zone we cannot do more than Eq. (151): There are no means to compute the radiation-reaction functions $A_L^{\alpha\beta}(t)$. We are here touching the second problem faced by the standard post-Newtonian approximation.

### 2.3.2 Post-Newtonian metric and radiation reaction effects

As we have understood this problem is that of the limitation to the near zone. Such limitation can be circumvented to the lowest post-Newtonian orders by considering *retarded* integrals that are formally expanded when $c \to +\infty$ as series of “instantaneous” Poisson-like integrals, see e.g., Anderson and DeCanio (1975). This procedure works well up to the 2.5PN level and has been shown to correctly fix the dominant radiation reaction term at the 2.5PN order (Ehlers 1980; Kerlick 1980a, b; Papapetrou and Linet 1981). Unfortunately the procedure assumes fundamentally that the gravitational field, after expansion of all retardations $r/c \to 0$, depends on the state of the source at a single time $t$, in keeping with the instantaneous character of the Newtonian interaction. However, we know that the post-Newtonian field (as well as the source’s dynamics) will cease at some stage to be given by a
functional of the source parameters at a single time, because of the imprint of gravitational-wave tails in the near zone field, in the form of an hereditary contribution (including a modification of the radiation reaction force) at the 4PN relative order (Blanchet and Damour 1988). Thus the formal post-Newtonian expansion of retarded Green functions is no longer valid starting at the 4PN order.

The solution of the problem resides in the matching of the near-zone field to the exterior field. We have already seen in Theorems 6 and 7 that the matching equation (115) yields the expression of the multipole expansion in the exterior domain. Now we prove that it also permits the full determination of the post-Newtonian metric in the near-zone, i.e., the radiation-reaction functions $A_L^{\alpha\beta}$ which have been left unspecified in Eq. (151). As a result we find that the functions $A_L^{\alpha\beta}$ are composed of the multipole moment functions $F_L^{\alpha\beta}$ defined by Eq. (128), which will characterize “linear-order” radiation reaction effects starting at 2.5PN order, and of an extra piece $R_L^{\alpha\beta}$, which will be due to non-linear effects in the radiation reaction (arising from radiative modes propagating at infinity) and turn out to arise at the 4PN order. Thus,

$$A_L^{\alpha\beta}(t) = F_L^{\alpha\beta}(t) + R_L^{\alpha\beta}(t). \qquad (152)$$

The piece $R_L^{\alpha\beta}$ is obtained from the multipole expansion of the pseudo-tensor $M(\tau^{\alpha\beta})$.<fnref n="35" /> Hence the radiation-reaction functions do depend on the behaviour of the field far away from the matter source (as physical intuition already told us). The explicit expression reads (Poujade and Blanchet 2002)

$$R_L^{\alpha\beta}(t) = \mathrm{FP}_{B=0} \int d^3\mathbf{x} \, \tilde{r}^B \hat{x}_L \int_1^{+\infty} dz \, \gamma_\ell(z) \mathcal{M}(\tau^{\alpha\beta}) (\mathbf{x}, t - zr/c). \qquad (153)$$

The fact that the multipolar expansion $M(\tau^{\alpha\beta})$ is the source term for the function $R_L^{\alpha\beta}$ is the consequence of the matching equation (115). The specific contributions due to $R_L^{\alpha\beta}$ in the post-Newtonian metric (151) are associated with gravitational wave tails. Notice that the FP regularization deals with the bound of the integral at $r = 0$ in Eq. (153), while in Eq. (128) it deals with the bound at $r = +\infty$. The weighting function $\gamma_\ell(z)$ therein, where $z$ extends up to infinity in contrast to the analogous function $\delta_\ell(z)$ in Eq. (128), is simply related to it by $\gamma_\ell(z) \equiv -2\delta_\ell(z)$; such definition is motivated by the fact that the integral of that function is normalized to one: $\int_1^{+\infty} dz \gamma_\ell(z) = 1$.<fnref n="36" /> The post-Newtonian metric (151) is now fully determined. However, let us now prove an interesting alternative formulation of the metric.

**Theorem 8** (Blanchet et al. 2005c) The expression of the post-Newtonian field in the near zone of a post-Newtonian source, satisfying correct boundary conditions at infinity (no incoming radiation), reads

<footnote n="35">³⁵ We mean the fully-fledge $M(\tau^{\alpha\beta})$; i.e., not the formal object $M(\bar{\tau}^{\alpha\beta})$.</footnote>

<footnote n="36">³⁶ Though the latter integral is a priori divergent, its value can be determined by invoking complex analytic continuation in $\ell \in \mathbb{C}$.</footnote>
$$ \bar{h}^{\alpha\beta} = \frac{16\pi G}{c^4} \Box_{\text{ret}}^{-1} [\bar{\tau}^{\alpha\beta}] \\ - \frac{4G}{c^4} \sum_{\ell=0}^{+\infty} \frac{(-)^\ell}{\ell!} \hat{\partial}_L \left\{ \frac{\mathcal{R}_L^{\alpha\beta}(t-r/c) - \mathcal{R}_L^{\alpha\beta}(t+r/c)}{2r} \right\}. \quad (154) $$

The first term represents a particular solution of the hierarchy of PN equations, while the second one is a homogeneous multipolar solution of the wave equation, of the “anti-symmetric” type that is regular at the origin $r = 0$ located inside the source, and parametrized by the multipole-moment functions (153).

Let us be more precise about the meaning of the first term in Eq. (154). Indeed such term is made of the formal expansion of the standard retarded integral (71) when $c \to \infty$, but acting on a post-Newtonian source term $\bar{\tau}^{\alpha\beta}$,

$$ \Box_{\text{ret}}^{-1} [\bar{\tau}^{\alpha\beta}](\mathbf{x}, t) = -\frac{1}{4\pi} \sum_{m=0}^{+\infty} \frac{(-)^m}{m!} \left(\frac{\hat{\partial}}{c\hat{\partial}t}\right)_{B=0}^m \text{FP} \quad (155) \\ \int d^3\mathbf{x}' \, \tilde{r}'^B |\mathbf{x} - \mathbf{x}'|^{m-1} \bar{\tau}^{\alpha\beta}(\mathbf{x}', t). $$

We emphasize that (155) constitutes the definition of a (formal) PN expansion, each term of which being built from the post-Newtonian expansion of the pseudo-tensor. Crucial in the present formalism, is that each of the terms is regularized by means of the FP operation in order to deal with the bound at infinity at which the post-Newtonian expansion is singular. Because of the presence of this regularization, the object (155) should carefully be distinguished from the “global” solution $\Box_{\text{ret}}^{-1}[\bar{\tau}^{\alpha\beta}]$ defined by Eq. (71), with global non-expanded pseudo-tensor $\bar{\tau}^{\alpha\beta}$.

The definition (155) is of interest because it corresponds to what one would intuitively think as the natural way of performing the PN iteration, i.e., by formally Taylor expanding the retardations in Eq. (71) as advocated by Anderson and DeCanio (1975). Moreover, each of the terms of the series (155) is mathematically well-defined thanks to the finite part operation, and can therefore be implemented in practical computations. The point is that Eq. (155) solves the wave equation in a perturbative post-Newtonian sense,

$$ \Box(\Box_{\text{ret}}^{-1}[\bar{\tau}^{\alpha\beta}]) = \bar{\tau}^{\alpha\beta}, \quad (156) $$

so constitutes a good prescription for a particular solution of the wave equation—as legitimate as the solution (149). Therefore the two solutions should differ by an homogeneous solution of the wave equation which is necessarily of the anti-symmetric type (regular inside the source). Detailed investigations yield

$$ \Box_{\text{ret}}^{-1}[\bar{\tau}^{\alpha\beta}] = \Box_{\text{inst}}^{-1}[\bar{\tau}^{\alpha\beta}] \\ - \frac{1}{4\pi} \sum_{\ell=0}^{+\infty} \frac{(-)^\ell}{\ell!} \hat{\partial}_L \left\{ \frac{\mathcal{F}_L^{\alpha\beta}(t-r/c) - \mathcal{F}_L^{\alpha\beta}(t+r/c)}{2r} \right\}, \quad (157) $$

where the homogeneous solution is parametrized by the multipole-moments $\mathcal{F}_L^{\alpha\beta}(t)$.
By combining Eqs. (154) and (157), we indeed become reconciled with the previous expression of the post-Newtonian field found in Eq. (151).

For computations limited to the 3.5PN order (level of the 1PN correction to the radiation reaction force), the first term in Eq. (154) with the “intuitive” prescription (155) is sufficient. But because of the second term in (154) there is a fundamental breakdown of this scheme at the 4PN order where it becomes necessary to take into account non-linear radiation effects associated with tails. The second term in (154) constitutes a generalization of the tail-transported radiation mode arising at 4PN order (Blanchet and Damour 1988). In particular, the tail-transported radiation reaction is required by energy conservation and the propagation of tails in the wave zone. The usual radiation reaction terms, up to 3.5PN order, are contained in the first term of Eq. (154), and are parametrized by the same multipole-moment functions $F_L^{\alpha\beta}$ as the exterior multipolar field, as Eq. (157) explicitly shows. In Sect. 2.3.4 we shall give an explicit expression of the radiation reaction force showing the usual radiation reaction terms to 3.5PN order, issued from $F_L^{\alpha\beta}$, and exhibiting the above tail-induced 4PN effect, issued from $R_L^{\alpha\beta}$.

Finally note that the post-Newtonian solution, in either (151) or (154) form, has been obtained without imposing the condition of harmonic coordinates (50) in an explicit way. We have simply matched together the PN and multipolar expansions, satisfying the “relaxed” Einstein field equations (51) in their respective domains, and found that the matching determines uniquely the solution. An important check in Poujade and Blanchet 2002; Blanchet et al. 2005c, is to verify that the harmonic coordinate condition (50) is indeed satisfied as a consequence of the conservation of the pseudo-tensor (54), so that we really grasp a solution of the full Einstein field equations.

### 2.3.3 The 3.5PN metric for general matter systems

The detailed calculations that are called for in applications necessitate having at one’s disposal some explicit expressions of the metric coefficients $g_{\alpha\beta}$, in harmonic coordinates, at the highest possible post-Newtonian order. The 3.5PN metric that we present below can be viewed as an application of the formalism of the previous section. It is expressed by means of some particular retarded-type potentials, $V, V_i, \hat{W}_{ij}, \dots$, whose main advantages are to somewhat minimize the number of terms, so that even at the 3.5PN order the metric is still tractable, and to delineate the different problems associated with the computation of different categories of terms. Of course, these potentials have no direct physical significance by themselves, but they offer a convenient parametrization of the 3.5PN metric.

The idea of this post-Newtonian parametrization comes from the following simple expression of the metric at the 1PN order.

**Theorem 9** (Blanchet and Damour 1989) The post-Newtonian metric in harmonic coordinates at the 1PN order is described by a scalar potential $V$ and a vector potential $V_i$ as
$$g_{00} = -1 + \frac{2}{c^2} V - \frac{2}{c^4} V^2 + \mathcal{O}\left(\frac{1}{c^6}\right), \qquad (158a)$$

$$g_{0i} = -\frac{4}{c^3} V_i + \mathcal{O}\left(\frac{1}{c^5}\right), \qquad (158b)$$

$$g_{ij} = \delta_{ij} \left[ 1 + \frac{2}{c^2} V \right] + \mathcal{O}\left(\frac{1}{c^4}\right). \qquad (158c)$$

The potentials represent some retarded versions of the gravito-electric (Newtonian) and gravito-magnetic potentials,

$$V = \Box_{\text{ret}}^{-1} [-4\pi G\sigma], \qquad V_i = \Box_{\text{ret}}^{-1} [-4\pi G\sigma_i], \qquad (159)$$

and the components of the matter stress-energy tensor $T^{\alpha\beta}$ are encoded into the convenient definitions [recalling Eqs. (139)],

$$\sigma \equiv \frac{T^{00} + T^{ii}}{c^2}, \qquad \sigma_i \equiv \frac{T^{0i}}{c}, \qquad \sigma_{ij} \equiv T^{ij}, \qquad (160)$$

(where $T^{ii} \equiv \delta_{ij}T^{ij}$). In particular $\sigma$ agrees up to 1PN order with the Tolman (1962) mass density valid for stationary systems.

The generalization to 3PN order—actually 3.5PN since the radiation-reaction odd parity terms are included into the definition of the potentials—reads (Blanchet and Faye 2001a)

$$g_{00} = -1 + \frac{2}{c^2} V - \frac{2}{c^4} V^2 + \frac{8}{c^6} \left( \hat{X} + V_i V_i + \frac{V^3}{6} \right) + \frac{32}{c^8} \left( \hat{T} - \frac{1}{2} V \hat{X} + \hat{R}_i V_i - \frac{1}{2} V V_i V_i - \frac{1}{48} V^4 \right) + \mathcal{O}\left(\frac{1}{c^{10}}\right), \qquad (161a)$$

$$g_{0i} = -\frac{4}{c^3} V_i - \frac{8}{c^5} \hat{R}_i - \frac{16}{c^7} \left( \hat{Y}_i + \frac{1}{2} \hat{W}_{ij} V_j + \frac{1}{2} V^2 V_i \right) + \mathcal{O}\left(\frac{1}{c^9}\right), \qquad (161b)$$

$$g_{ij} = \delta_{ij} \left[ 1 + \frac{2}{c^2} V + \frac{2}{c^4} V^2 + \frac{8}{c^6} \left( \hat{X} + V_k V_k + \frac{V^3}{6} \right) \right] + \frac{4}{c^4} \hat{W}_{ij} + \frac{16}{c^6} \left( \hat{Z}_{ij} + \frac{1}{2} V \hat{W}_{ij} - V_i V_j \right) + \mathcal{O}\left(\frac{1}{c^8}\right). \qquad (161c)$$

In addition to (159), from the 2PN order we have the potentials

$$\hat{X} = \Box_{\text{ret}}^{-1} [-4\pi G V \sigma_{ii} + \hat{W}_{ij} \partial_{ij} V + 2V_i \partial_t \partial_i V + V \partial_t^2 V + \frac{3}{2} (\partial_t V)^2 - 2 \partial_i V_j \partial_j V_i] \qquad (162a)$$
$$ \hat{R}_i = \Box_{\text{ret}}^{-1} \left[ -4\pi G(V\sigma_i - V_i\sigma) - 2\partial_k V\partial_i V_k - \frac{3}{2}\partial_t V\partial_i V \right], \quad (162\text{b}) $$

$$ \hat{W}_{ij} = \Box_{\text{ret}}^{-1} [-4\pi G(\sigma_{ij} - \delta_{ij}\sigma_{kk}) - \partial_i V \partial_j V]. \quad (162\text{c}) $$

As we see, some parts of these potentials are directly generated by compact-support matter terms, while other parts are made of non-compact-support products of V-type potentials. There exists also an important cubically non-linear term generated by the coupling between $\hat{W}_{ij}$ and V, see the second term in the $\hat{X}$-potential. Note the important point that here and below the retarded integral operator $\Box_{\text{ret}}^{-1}$ is really meant to be the one given by Eq. (155); thus it involves in principle the finite part regularization FP to deal with (IR-type) divergences occurring at high PN orders for non-compact-support integrals. For instance, such finite part regularization is important to take into account in the computation of the near zone metric at the 3PN order and for the radiation field at the 4PN order (Blanchet et al. 2010b, 2023b).

At the next level, 3PN, we have the even more complicated potentials

$$ \begin{aligned} \hat{T} = \Box_{\text{ret}}^{-1} \left[ \right. & -4\pi G \left( \frac{1}{4}\sigma_{ij}\hat{W}_{ij} + \frac{1}{2}V^2\sigma_{ii} + \sigma V_i V_i \right) + \hat{Z}_{ij}\partial_{ij}V + \hat{R}_i\partial_t V\partial_i V \\ & - 2\partial_i V_j \partial_j \hat{R}_i - \partial_i V_j \partial_t \hat{W}_{ij} + V V_i \partial_t \partial_i V + 2V_i \partial_j V_i \partial_j V \\ & \left. + \frac{3}{2}V_i \partial_t V \partial_i V + \frac{1}{2}V^2 \partial_t^2 V + \frac{3}{2}V(\partial_t V)^2 - \frac{1}{2}(\partial_t V_i)^2 \right], \end{aligned} \quad (163\text{a}) $$

$$ \begin{aligned} \hat{Y}_i = \Box_{\text{ret}}^{-1} \left[ \right. & -4\pi G \left( -\sigma \hat{R}_i - \sigma V V_i + \frac{1}{2} \sigma_k \hat{W}_{ik} + \frac{1}{2} \sigma_{ik} V_k + \frac{1}{2} \sigma_{kk} V_i \right) \\ & + \hat{W}_{kl} \partial_{kl} V_i - \partial_t \hat{W}_{ik} \partial_k V + \partial_i \hat{W}_{kl} \partial_k V_l - \partial_k \hat{W}_{il} \partial_l V_k - 2 \partial_k V \partial_i \hat{R}_k \\ & \left. - \frac{3}{2} V_k \partial_i V \partial_k V - \frac{3}{2} V \partial_t V \partial_i V - 2 V \partial_k V \partial_k V_i + V \partial_t^2 V_i + 2 V_k \partial_k \partial_t V_i \right], \end{aligned} \quad (163\text{b}) $$

$$ \begin{aligned} \hat{Z}_{ij} = \Box_{\text{ret}}^{-1} \left[ \right. & -4\pi G V (\sigma_{ij} - \delta_{ij}\sigma_{kk}) - 2\partial_{(i} V \partial_{t} V_{j)} + \partial_{i} V_{k} \partial_{j} V_{k} + \partial_{k} V_{i} \partial_{k} V_{j} \\ & - 2\partial_{(i} V_{k} \partial_{k} V_{j)} - \frac{3}{4} \delta_{ij} (\partial_{t} V)^2 - \delta_{ij} \partial_{k} V_{m} (\partial_{k} V_{m} - \partial_{m} V_{k}) \left. \right]. \end{aligned} \quad (163\text{c}) $$

These involve many types of compact-support contributions, as well as quadratic-order and cubic-order parts; but, surprisingly, there are no quartically non-linear terms. Indeed it has been possible to “integrate directly” all the quartic contributions in the 3PN metric; see the terms composed of $V^4$ and $V\hat{X}$ in the first of Eqs. (161).

Note that the 3PN metric (161) does represent the inner post-Newtonian field of an isolated system, because it contains, to this order, the correct radiation-reaction terms corresponding to outgoing radiation. These terms come from the expansions of the retardations in the retarded potentials (159)–(163); we elaborate more on radiation-reaction effects in Sect. 2.3.4.
The above potentials are not independent: They are linked together by some differential identities issued from the harmonic gauge conditions, which are equivalent, via the Bianchi identities, to the equations of motion of the matter fields; see Eq. (54). These identities read

$$ \begin{aligned} & \partial_t \left\{ V + \frac{1}{c^2} \left[ \frac{1}{2} \hat{W}_{kk} + 2V^2 \right] + \frac{4}{c^4} \left[ \hat{X} + \frac{1}{2} \hat{Z}_{kk} + \frac{1}{2} V \hat{W}_{kk} + \frac{2}{3} V^3 \right] \right\} \\ & + \partial_i \left\{ V_i + \frac{2}{c^2} [\hat{R}_i + V V_i] + \frac{4}{c^4} \left[ \hat{Y}_i - \frac{1}{2} \hat{W}_{ij} V_j + \frac{1}{2} \hat{W}_{kk} V_i + V \hat{R}_i + V^2 V_i \right] \right\} \\ & = \mathcal{O} \left( \frac{1}{c^6} \right), \end{aligned} \quad (164a) $$

$$ \begin{aligned} & \partial_t \left\{ V_i + \frac{2}{c^2} [\hat{R}_i + V V_i] \right\} + \partial_j \left\{ \hat{W}_{ij} - \frac{1}{2} \hat{W}_{kk} \delta_{ij} + \frac{4}{c^2} \left[ \hat{Z}_{ij} - \frac{1}{2} \hat{Z}_{kk} \delta_{ij} \right] \right\} \\ & = \mathcal{O} \left( \frac{1}{c^4} \right). \end{aligned} \quad (164b) $$

The above 3PN metric and differential identities have been generalized to 4PN order in the Appendix A of Marchand et al. (2020).

For later applications to systems of compact objects, let us give the geodesic equations of a particle moving in the 3.5PN metric (161). Of course the geodesic equations are appropriate for the motion of particles without spins nor internal structure; for spinning particles one has also to take into account the coupling of the spin to the space-time curvature, see the Mathisson-Papapetrou equations (574). It is convenient to write the geodesic equations in the form (where $t = x^0/c$ is the coordinate time)

$$ \frac{dP_i}{dt} = F_i, \quad (165) $$

where the “linear momentum density” $P_i$ and the “force density” $F_i$ of the particle are given by

$$ P_i = \frac{g_{i\mu}(x) v^\mu}{\sqrt{-g_{\rho\sigma}(x) \frac{v^\rho v^\sigma}{c^2}}}, \quad (166a) $$

$$ F_i = \frac{1}{2} \frac{\partial_i g_{\mu\nu}(x) v^\mu v^\nu}{\sqrt{-g_{\rho\sigma}(x) \frac{v^\rho v^\sigma}{c^2}}}, \quad (166b) $$

where $v^\mu = (c, v^j)$ with $v^j = dx^i/dt$ being the particle's ordinary coordinate velocity, and where the metric components and their gradients are taken at the location of the particle. Notice that we are viewing here the particle as moving in the fixed background metric (161). In Sect. 3, the metric will be generated by the system of particles itself, and we shall see that the same equations still hold in this case,
provided that we supplement the computation of the metric at the location of one of these particles by a suitable self-field regularization.

Further denoting $P_i = v_i + \hat{P}_i$ the equations give the ordinary acceleration $a^i = dv^i/dt$ as

$$a^i = F_i - \frac{d\hat{P}_i}{dt}, \qquad (167)$$

where the expressions of $F_i$ and $\hat{P}_i$ in terms of the non-linear potentials follow from insertion of the 3.5PN metric coefficients (161). We obtain

$$ \begin{aligned} \hat{P}_i = & \frac{1}{c^2} \left( \frac{1}{2} v^2 v^i + 3 V v^i - 4 V_i \right) \\ & + \frac{1}{c^4} \left( \frac{3}{8} v^4 v^i + \frac{7}{2} V v^2 v^i - 4 V_j v^i v^j - 2 V_i v^2 \right) \\ & + \frac{9}{2} V^2 v^i - 4 V V_i + 4 \hat{W}_{ij} v^j - 8 \hat{R}_i \right) \\ & + \frac{1}{c^6} \left( \frac{5}{16} v^6 v^i + \frac{33}{8} V v^4 v^i - \frac{3}{2} V_i v^4 - 6 V_j v^i v^j v^2 + \frac{49}{4} V^2 v^2 v^i \right) \quad (168a) \\ & + 2 \hat{W}_{ij} v^j v^2 - 4 \hat{R}_i v^2 + 2 \hat{W}_{jk} v^i v^j v^k - 10 V V_i v^2 - 20 V V_j v^i v^j \\ & - 8 \hat{R}_j v^i v^j + \frac{9}{2} V^3 v^i + 12 V_j V^i v^i + 12 \hat{W}_{ij} V v^j + 12 \hat{X} v^i + 16 \hat{Z}_{ij} v^j \\ & - 10 V^2 V_i - 8 \hat{W}_{ij} V_j - 8 V \hat{R}_i - 16 \hat{Y}_i \right) + \mathcal{O}\left(\frac{1}{c^8}\right), \end{aligned} $$

$$ \begin{aligned} F_i = & \partial_i V \\ & + \frac{1}{c^2} \left( -V \partial_i V + \frac{3}{2} \partial_i V v^2 - 4 \partial_i V_j v^j \right) \\ & + \frac{1}{c^4} \left( \frac{7}{8} \partial_i V v^4 - 2 \partial_i V_j v^j v^2 + \frac{9}{2} V \partial_i V v^2 + 2 \partial_i \hat{W}_{jk} v^j v^k \right) \\ & - 4 V_j \partial_i V v^j - 4 V \partial_i V_j v^j - 8 \partial_i \hat{R}_j v^j + \frac{1}{2} V^2 \partial_i V + 8 V_j \partial_i V_j + 4 \partial_i \hat{X} \right) \\ & + \frac{1}{c^6} \left( \frac{11}{16} v^6 \partial_i V - \frac{3}{2} \partial_i V_j v^j v^4 + \frac{49}{8} V \partial_i V v^4 + \partial_i \hat{W}_{jk} v^j v^k - 10 V_j \partial_i V v^2 v^j \right) \quad (168b) \\ & - 10 V \partial_i V_j v^2 v^j - 4 \partial_i \hat{R}_j v^j v^2 + \frac{27}{4} V^2 \partial_i V v^2 + 12 V_j \partial_i V_j v^2 \\ & + 6 \hat{W}_{jk} \partial_i V v^j v^k + 6 V \partial_i \hat{W}_{jk} v^j v^k + 6 \partial_i \hat{X} v^2 + 8 \partial_i \hat{Z}_{jk} v^j v^k \\ & - 20 V_j V \partial_i V v^j - 10 V^2 \partial_i V_j v^j - 8 V_k \partial_i \hat{W}_{jk} v^j - 8 \hat{W}_{jk} \partial_i V_k v^j \\ & - 8 \hat{R}_j \partial_i V v^j - 8 V \partial_i \hat{R}_j v^j - 16 \partial_i \hat{Y}_j v^j - \frac{1}{6} V^3 \partial_i V - 4 V_j V_j \partial_i V \\ & + 16 \hat{R}_j \partial_i V_j + 16 V_j \partial_i \hat{R}_j - 8 V V_j \partial_i V_j - 4 \hat{X} \partial_i V - 4 V \partial_i \hat{X} + 16 \partial_i \hat{T} \right) \\ & + \mathcal{O}\left(\frac{1}{c^8}\right). \end{aligned} $$
Note that all the accelerations appearing in the potentials and in the final expression of the equations of motion should be order-reduced by means of the equations of motion themselves. For instance, we see that when computing the time-derivative of $\hat{P}_i$ we shall meet some accelerations at 1PN order which are therefore to be replaced by the explicit 2.5PN equations of motion. The order-reduction is a crucial aspect of the post-Newtonian method. It is justified by the fact that the matter equations of motion, say $\nabla_\mu T^{\mu\lambda} = 0$, represent four out of the ten Einstein field equations, see Sect. 2.1.1 for discussion. In the harmonic-coordinate approach the equations of motion are equivalent to the harmonic gauge conditions $\hat{\partial}_\mu h^{\mu\lambda} = 0$. Post-Newtonian predictions based on such consistent PN order-reduction have been very successful. Note, however, that the operation of order-reduction is illicit at the level of the Lagrangian. The elimination of accelerations in a Lagrangian by substituting the equations of motion derived from that Lagrangian, results in a different Lagrangian whose equations of motion differ from those of the original Lagrangian by a gauge transformation (Schäfer 1984).

### 2.3.4 Radiation reaction potentials

We said that the metric (161) contains the correct radiation-reaction terms appropriate for an isolated system up to the 3.5PN level included. As we shall see the metric can even be generalized to include the piece of the radiation reaction which occurs at the 4PN order, and which is due to gravitational wave tails. In Sect. 1.1.2 we already discussed a particular non-harmonic coordinate system which is very convenient to describe the radiation reaction terms. In this Burke and Thorne (1970), Burke (1971) coordinate system, the radiation reaction force at the lowest 2.5PN order takes the simple form of Eq. (7), which involves only a scalar potential $V^{\text{reac}}$ so that

$$F_i^{\text{reac}} = \rho \hat{\partial}_i V^{\text{reac}} + \mathcal{O}\left(\frac{1}{c^7}\right), \qquad (169a)$$

$$V^{\text{reac}}(\mathbf{x}, t) = -\frac{G}{5c^5} x^i x^j Q_{ij}^{(5)}(t) + \mathcal{O}\left(\frac{1}{c^7}\right). \qquad (169b)$$

At such dominant 2.5PN level (“Newtonian” radiation reaction) the source quadrupole moment is given by the usual Newtonian expression (4).

Let us first emphasize that the Burke-Thorne expression (169) of the radiation reaction force is not unique as it depends on the coordinate system (Schäfer 1981, 1982). The first complete derivation of the radiation reaction force for general systems at the lowest 2.5PN order was done by Chandrasekhar and Esposito (1970)<fnref n="37" /> who obtained a different expression valid in a different coordinate system. In harmonic coordinates the radiation reaction force reads

<footnote n="37">³⁷ See also further works by Papapetrou and Linet (1981), Kerlick (1980a, b).</footnote>
$$F_i^{\text{reac}}|_{\text{harm}} = \frac{G}{c^5} \rho \left\{ \frac{3}{5} x^j Q_{ij}^{(5)} + 2 \frac{d}{dt} (\nu^j Q_{ij}^{(3)}) - Q_{jk}^{(3)} \partial_{ijk} \chi \right\} + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (170)$$

where we have introduced the “super-potential” $\chi$ of the Newtonian potential $U$, such that $\Delta\chi = 2U$ and $\Delta U = -4\pi G\rho$.<fnref n="38" /> Repeating the analysis (9)–(10) with the different expression of the radiation reaction force (170) we recover the total energy flux given by the quadrupole formula (5); the only difference is that the total time derivatives in (9) will be different, with $f$ and $g_i$ corresponding to harmonic coordinates. The different forms of the radiation reaction force were a source of confusion in the 1970s and fed the so-called radiation-reaction controversy at that time (Ehlers et al. 1976; Walker and Will 1980); see Kennefick (2007) for an historical account.

The novel feature when one extends the Newtonian radiation reaction (169) to include the 1PN corrections is that it is no longer composed of a single scalar depending on the mass-type multipole moments, but involves also a vectorial component depending in particular on the current-type quadrupole moment. This was noticed in the physically restricted case where the dominant quadrupolar radiation from the source is suppressed (Blanchet and Damour 1984). The vectorial component of the reaction force could be important in some astrophysical situations like rotating neutron stars undergoing gravitational instabilities. Here we report the results of the extension to 1PN order (and even 1.5PN) of the lowest-order Burke–Thorne (BT) scalar radiation reaction potential (169), in some appropriate coordinate system, which can be called extended BT coordinate system (Blanchet 1993, 1997b).

At that level (corresponding to 3.5PN order in the metric), and in the extended BT coordinate system, it suffices to incorporate some radiation-reaction contributions into the scalar and vectorial potentials which parametrize the 1PN metric in Eq. (9). We thus pose

$$V = V^{\text{cons}} + V^{\text{reac}}, \quad (171a)$$

$$V_i = V_i^{\text{cons}} + V_i^{\text{reac}}. \quad (171b)$$

Then the metric in extended BT coordinates, accurate to 3.5PN order regarding the radiation-reaction contributions and to 1PN order for the conservative effects—we indicate this by using the symbol $\mathcal{O}^{\text{reac}}$ for the remainders—takes the same form as the harmonic-coordinates 1PN metric (9),

$$g_{00}^{\text{BT}} = -1 + \frac{2}{c^2} V - \frac{2}{c^4} V^2 + \mathcal{O}^{\text{reac}}\left(\frac{1}{c^{11}}\right), \quad (172a)$$

$$U_{ij}(\mathbf{x}, t) = \int d^3\mathbf{x}' \rho(\mathbf{x}', t) \frac{(\mathbf{x}^i - \mathbf{x}'^i)(\mathbf{x}^j - \mathbf{x}'^j)}{|\mathbf{x} - \mathbf{x}'|^3},$$

<footnote n="38">³⁸ We have $\chi(\mathbf{x}, t) = G \int d^3\mathbf{x}' \rho(\mathbf{x}', t) |\mathbf{x} - \mathbf{x}'|$ and the super-potential is linked to the Newtonian tensor potential</footnote>

<footnote>(such that $U_{ii} = U$) by $\partial_{ij}\chi = \delta_{ij}U - U_{ij}$.</footnote>
$$g_{0i}^{\text{BT}} = -\frac{4}{c^3} V_i + \mathcal{O}^{\text{reac}}\left(\frac{1}{c^{10}}\right), \qquad (172\text{b})$$

$$g_{ij}^{\text{BT}} = \delta_{ij} \left[ 1 + \frac{2}{c^2} \mathcal{V} \right] + \mathcal{O}^{\text{reac}} \left( \frac{1}{c^9} \right). \qquad (172\text{c})$$

The other contributions, which are conservative (i.e., non radiative), are given up to 3PN order by the metric (161) in which all the potentials take the same form as in Eqs. (159)–(163), but where one neglects all the retardations, which means that the retarded integral operator $\Box_{\text{ret}}^{-1}$ should be replaced by the operator of the “instantaneous” potentials $\Box_{\text{inst}}^{-1}$ defined by Eq. (149). As such the conservative parts of the potentials (171) are defined by

$$V^{\text{cons}} = \Box_{\text{inst}}^{-1} [-4\pi G \sigma], \qquad V_i^{\text{cons}} = \Box_{\text{inst}}^{-1} [-4\pi G \sigma_i]. \qquad (173)$$

Up to 3.5PN order, the effect of all these retardations gets replaced by the effect of the radiation-reaction potentials $V^{\text{reac}}$ and $V_i^{\text{reac}}$. These are explicitly given by (Blanchet 1993, 1997b)<fnref n="39" />

$$V^{\text{reac}}(\mathbf{x}, t) = -\frac{G}{5c^5} x^{ij} I_{ij}^{(5)}(t) \\ + \frac{G}{c^7} \left[ \frac{1}{189} x^{ijk} I_{ijk}^{(7)}(t) - \frac{1}{70} r^2 x^{ij} I_{ij}^{(7)}(t) \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \qquad (174\text{a})$$

$$V_i^{\text{reac}}(\mathbf{x}, t) = \frac{G}{c^5} \left[ \frac{1}{21} \hat{x}^{ijk} I_{jk}^{(6)}(t) - \frac{4}{45} \epsilon_{ijk} x^{jl} \mathbf{J}_{kl}^{(5)}(t) \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \qquad (174\text{b})$$

where the multipole moments $I_L$ and $J_L$ denote the source multipole moments defined in a general way by Eqs. (7). The scalar potential $V^{\text{reac}}$ will obviously reproduce Eq. (169) at the dominant order. However, note that it is crucial to include in Eq. (174a) the 1PN correction in the source quadrupole moment. Following (138) the mass-type quadrupole moment $I_{ij}$ to 1PN order involves only *compact support* contributions and reads

$$I_{ij} = \int d^3\mathbf{x} \left\{ \hat{x}_{ij}\sigma + \frac{1}{14c^2} r^2 \hat{x}_{ij} \hat{\partial}_t^2 \sigma - \frac{20}{21c^2} \hat{x}_{ijk} \partial_t \sigma_k \right\} + \mathcal{O}\left(\frac{1}{c^4}\right), \qquad (175)$$

where the compact-support matter source densities $\sigma$ and $\sigma_i$ are given in Eqs. (160). Obviously to this order the other moments are Newtonian, $I_L = \int d^3\mathbf{x} \hat{x}_L \sigma + \mathcal{O}(c^{-2})$ and $J_L = \int d^3\mathbf{x} \epsilon_{ab} <i_t \hat{x}_{L-1} > a \sigma_b + \mathcal{O}(c^{-2})$.

The 3.5PN radiation reaction force in the equations of motion of compact binary systems has been derived by Iyer and Will (1993, 1995) in an arbitrary gauge, in the frame of the center of mass, assuming the validity of the flux-balance equations for energy and angular momentum at the relative 1PN order. Their result is presented in

<footnote n="39">³⁹ Recall the footnote 23 for our notation. In particular $\hat{x}^{ijk}$ in the vector potential denotes the STF combination $\hat{x}^{ijk} = x^{ijk} - \frac{c^2}{5}(x^i \delta^{jk} + x^j \delta^{ki} + x^k \delta^{ij})$ with $x^{ijk} = x^i x^j x^k$.</footnote>
terms of a set of arbitrary parameters reflecting the arbitrariness in the choice of the coordinate system. Iyer and Will (1995) have shown that the extended BT coordinate system indeed corresponds to a unique set of these parameters. Pati and Will (2002), Nissanke and Blanchet (2005), Itoh (2009) obtained the 1PN radiation in harmonic coordinates (again, corresponding to a unique set of Iyer-Will parameters) and Königsdörffer et al. (2003) obtained it in ADM coordinate (corresponding to another unique set of parameters). The method of balance equations has been extended to the 4.5PN order corresponding to 2PN relatively to lowest order by Gopakumar et al. (1997). To date the only calculation of radiation reaction from first principle to relative 2PN order was done using the EFT by Leibovich et al. (2023), however restricted to the center-of-mass frame [they match the Gopakumar et al. (1997) parameters].

The previous formalism can be extended to 4PN order (1.5PN beyond the leading radiation reaction). Indeed at this order there is a modification of the scalar radiation-reaction potential due to reaction to gravitational-wave tails propagating at infinity. The tail term is a hereditary integral characterized by a logarithmic kernel (we shall discuss thoroughly the tails in Sect. 2.4.2) and we have (Blanchet 1997b)<fnref n="40" />

$$V^{\text{react}}|_{\text{tail}} = - \frac{4G^2 M}{5c^8} x^{ij} \int_{0}^{+\infty} d\tau I_{ij}^{(7)}(t-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{11}{12} \right] + \mathcal{O}\left(\frac{1}{c^9}\right), \quad (176)$$

which has to be added to Eq. (174a) in order to complete the radiation reaction to 4PN order (in the extended BT gauge).

With the radiation-reaction potentials in hands, one can *prove* the flux-balance equation for energy up to 1.5PN order, namely

$$\begin{aligned} \frac{dE}{dt} = & -\frac{G}{5c^5} \left( I_{ij}^{(3)} + \frac{2GM}{c^3} \int_0^{+\infty} d\tau I_{ij}^{(5)}(t-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{11}{12} \right] \right)^2 \\ & -\frac{G}{c^7} \left[ \frac{1}{189} I_{ijk}^{(4)} I_{ijk}^{(4)} + \frac{16}{45} J_{ij}^{(3)} J_{ij}^{(3)} \right] + \mathcal{O}\left(\frac{1}{c^9}\right), \end{aligned} \quad (177)$$

where *E* denotes a functional of the matter in the source made of conservative terms and also, consistently with the approximation, of 2.5PN, 3.5PN and 4PN dissipative terms. The constant 11/12 depends on the gauge and corresponds to a total time-derivative in the right-hand side of (177), which can thus be transferred to the left-hand side. In the right-hand side one recognizes (minus) the positive-definite expression for the energy flux at 1.5PN order. Indeed, as we shall prove in 2.4.2 [see Eqs. (239) and (251)], the effective quadrupole moment which appears in the parenthesis of (177) agrees with the tail-modified *radiative* quadrupole moment $U_{ij}$ parametrizing the field in the far zone.

$$I_{ij}^{(7)}(t-\tau) \rightarrow \frac{1}{2} [I_{ij}^{(7)}(t-\tau) + I_{ij}^{(7)}(t+\tau)].$$

<footnote n="40">⁴⁰ Actually one must consider the time-asymmetric (dissipative) piece obtained by changing</footnote>
## 2.4 Asymptotic gravitational waveform

### 2.4.1 The radiative multipole moments

In Theorem 5 we constructed the metric in a radiative (or Bondi-type) coordinate system $(T, \mathbf{X})$, such that the retarded time $U$ becomes asymptotically a null coordinate at future null infinity (Papapetrou 1969; Madore 1970a, b). The leading-order term $1/R$ (when $R \to +\infty$) yields the operational definition of two sets of STF radiative multipole moments, mass-type $U_L(U)$ and current-type $V_L(U)$, where $U \equiv T - R/c$. The radiative moments are defined from the spatial components $ij$ of the metric in a transverse-traceless (TT) radiative coordinate system. By definition, we have (Thorne 1980)

$$h_{ij}^{\text{TT}}(U, \mathbf{X}) = -\frac{4G}{c^2 R} \perp_{ijab}(\mathbf{N}) \sum_{\ell=2}^{+\infty} \frac{1}{c^\ell \ell!} \left\{ \begin{aligned} & N_{L-2} U_{abL-2}(U) \\ & - \frac{2\ell}{c(\ell+1)} N_{cL-2} \epsilon_{cd(a} V_{b)dL-2}(U) \end{aligned} \right\} + \mathcal{O}\left(\frac{1}{R^2}\right). \quad (178)$$

The metric is written here as a perturbation of the “gothic” metric $\sqrt{-gg^{ij}}$; a change of sign is required when considering a perturbation of the ordinary metric $g_{ij}$. As before we denote for instance $N_{L-2} = N_{i_1} \cdots N_{i_{\ell-2}}$ and so on, where $N_i = (\mathbf{N})_i$ and $\mathbf{N} = \mathbf{X}/R$, and the TT algebraic projection operator $\perp_{ijab}$ is given by Eq. (3). Obviously the multipole decomposition (178) represents the generalization of the quadrupole moment formalism. Notice that the meaning of Eq. (178) is for the moment rather empty, because we do not yet know how to relate the radiative moments to the actual source parameters. Only at the Newtonian level do we know this relation, which is $U_{ij}(U) = Q_{ij}^{(2)}(U) + \mathcal{O}(c^{-2})$, where $Q_{ij}$ is the Newtonian quadrupole moment (4). We shall relate below the radiative moments to the source moments $I_L$ and $J_L$ obtained in Theorem 7.

Associated to the asymptotic waveform (178) we can compute by standard methods the total energy flux $\mathcal{F}$ and angular momentum flux $\mathcal{G}_i$ in gravitational waves (see Table 1),

$$\mathcal{F} \equiv \left(\frac{dE}{dU}\right)^{\text{GW}} = \sum_{\ell=2}^{+\infty} \frac{G}{c^{2\ell+1}} \left\{ a_{\ell} U_{L}^{(1)} U_{L}^{(1)} + \frac{b_{\ell}}{c^2} V_{L}^{(1)} V_{L}^{(1)} \right\}, \quad (179a)$$

<table><caption>**Table 1** The $\ell$-dependent coefficients in the flux formulas</caption><tr><td>Energy flux <i>F</i></td><td>a<sub>ℓ</sub> = (ℓ+1)(ℓ+2) / (ℓ-1)ℓℓ!(2ℓ+1)!!</td><td>b<sub>ℓ</sub> = 4ℓ(ℓ+2) / (ℓ-1)(ℓ+1)!(2ℓ+1)!!</td><td></td></tr><tr><td>Angular momentum flux <i>G</i><sub>i</sub></td><td>a'<sub>ℓ</sub> = (ℓ+1)(ℓ+2) / (ℓ-1)ℓ(ℓ+2)!!</td><td>b'<sub>ℓ</sub> = 4ℓ<sup>2</sup>(ℓ+2) / (ℓ-1)(ℓ+1)!(2ℓ+1)!!</td><td></td></tr><tr><td>Linear momentum flux <i>P</i><sub>i</sub></td><td>a''<sub>ℓ</sub> = 2(ℓ+2)(ℓ+3) / (ℓ(ℓ+1)!(2ℓ+3)!!)</td><td>b''<sub>ℓ</sub> = 8(ℓ+3) / (ℓ+1)!(2ℓ+3)!!</td><td>c''<sub>ℓ</sub> = 8(ℓ+2) / (ℓ-1)(ℓ+1)!(2ℓ+1)!!</td></tr><tr><td>Center of mass flux <i>C</i><sub>i</sub></td><td>a'''<sub>ℓ</sub> = 2(ℓ+2)(ℓ+3) / (ℓℓ!(2ℓ+3)!!)</td><td>b'''<sub>ℓ</sub> = 8(ℓ+3) / (ℓ(ℓ+3)!!)</td><td></td></tr></table>
$$ G_i \equiv \left( \frac{dJ_i}{dU} \right)^{\text{GW}} = \epsilon_{iab} \sum_{\ell=2}^{+\infty} \frac{G}{c^{2\ell+3}} \left\{ a'_{\ell} U_{aL-1} U_{bL-1}^{(1)} + \frac{b'_{\ell}}{c^2} V_{aL-1} V_{bL-1}^{(1)} \right\}. \quad (179b) $$

For instance the energy flux is readily obtained by plugging (178) into (109) and integrating over angles. The above fluxes satisfy the time-averaged flux-balance equations [already used in Eqs. (10)]

$$ \left\langle \frac{dE}{dt} \right\rangle = -\langle \mathcal{F} \rangle, \qquad \left\langle \frac{dJ_i}{dt} \right\rangle = -\langle G_i \rangle. \qquad (180) $$

Similarly one can compute the linear momentum flux $P_i$ which is responsible for the GW recoil of the source (see Sect. 1.1.4), as well as the flux associated to the center of mass position $G_i$. For a conservative system (for instance, stationary), there are ten Noetherian conserved invariants associated with the ten symmetries of the Poincaré group. In addition to $E$, $J_i$ and $P_i$, there is also the initial “position” of the center of mass, say $Z_i = G_i - P_i t$, where $G_i$ denotes the position of the center of mass multiplied by the total (constant) mass. The conservation of $Z_i$ is due to the invariance of the dynamics under Lorentz boosts. For self-gravitating systems, $G_i$ coincides with the mass-type dipole moment $I_i$. We have the fluxes

$$ P_i \equiv \left( \frac{dP_i}{dU} \right)^{\text{GW}} = \sum_{\ell=2}^{+\infty} \frac{G}{c^{2\ell+3}} \left\{ a''_{\ell} U_{iL}^{(1)} U_{L}^{(1)} + \frac{b''_{\ell}}{c^2} V_{iL}^{(1)} V_{L}^{(1)} + c''_{\ell} \epsilon_{iab} U_{aL-1}^{(1)} V_{bL-1}^{(1)} \right\}, \quad (181a) $$

$$ C_i \equiv \left( \frac{dG_i}{dU} \right)^{\text{GW}} = \sum_{\ell=2}^{+\infty} \frac{G}{c^{2\ell+3}} \left\{ a'''_{\ell} U_{iL} U_{L}^{(1)} + \frac{b'''_{\ell}}{c^2} V_{iL} V_{L}^{(1)} \right\}, \quad (181b) $$

which enter the following time-averaged flux-balance equations

$$ \left\langle \frac{dP_i}{dU} \right\rangle = -\langle P_i \rangle, \qquad \left\langle \frac{dG_i}{dU} \right\rangle = P_i - \langle C_i \rangle. \qquad (182) $$

While the general multipole expansions of the fluxes for energy, angular momentum and linear momentum in terms of radiative moments are well known (see notably Thorne 1980), the flux associated to the position of the center of mass has been computed more recently (Kozameh and Quirega 2016; Kozameh et al. 2018; Nichols 2018; Blanchet and Faye 2019; Compère et al. 2020).41 The coefficients in the flux formulas (179) and (181) are given in Table 1.

$$ \bar{C}_i = \frac{1}{2} \sum_{\ell=2}^{+\infty} \frac{G}{c^{2\ell+3}} \left\{ a'''_{\ell} (U_{iL} U_{L}^{(1)} - U_{iL}^{(1)} U_{L}) + \frac{b'''_{\ell}}{c^2} (V_{iL} V_{L}^{(1)} - V_{iL}^{(1)} V_{L}) \right\}, $$

which differs from $C_i$ obtained by Blanchet and Faye (2019) by a total time derivative, and is therefore physically equivalent after time averaging. An advantage of the formulation $\bar{C}_i$, as shown by Compère et al. (2020), is that it can be written by means of a covariant formula over the “celestial” sphere in Bondi gauge.

<footnote n="41">41 The result obtained by Kozameh et al. (2018), Compère et al. (2020) is</footnote>
Next we introduce two unit polarization vectors **P** and **Q**, orthogonal and transverse to the direction of propagation **N**, such that (N, **P**, **Q**) forms an oriented orthonormal triad (hence $N_iN_j + P_iP_j + Q_iQ_j = \delta_{ij}$). Our convention for the choice of **P** and **Q** is clarified in the footnote 90. Then the two “plus” and “cross” polarization states of the asymptotic waveform are defined by

$$h_+ \equiv \frac{1}{2} (P_i P_j - Q_i Q_j) h_{ij}^{TT}, \qquad (183a)$$

$$h_{\times} \equiv \frac{1}{2} (P_i Q_j + P_j Q_i) h_{ij}^{TT}. \qquad (183b)$$

Although the multipole decomposition (178) is completely general, it will also be important, having in view the comparison between the post-Newtonian and numerical results, to consider separately the various modes $(\ell, m)$ of the asymptotic waveform as defined with respect to a basis of spin-weighted spherical harmonics of weight $-2$. Those harmonics are function of the spherical angles $(\theta, \phi)$ defining the direction of propagation **N**, and given by

$$Y_{-2}^{\ell m} = \sqrt{\frac{2\ell + 1}{4\pi}} d^{\ell m}(\theta) e^{i\ell m \phi}, \qquad (184a)$$

$$d^{\ell m} = \sum_{k=k_1}^{k_2} \frac{(-)^k}{k!} e_k^{\ell m} \left(\cos \frac{\theta}{2}\right)^{2\ell+m-2k-2} \left(\sin \frac{\theta}{2}\right)^{2k-m+2}, \qquad (184b)$$

$$e_k^{\ell m} = \frac{\sqrt{(\ell + m)!(\ell - m)!(\ell + 2)!(\ell - 2)!}}{(k - m + 2)!(\ell + m - k)!(\ell - k - 2)!}, \qquad (184c)$$

where $k_1 = \max(0, m - 2)$ and $k_2 = \min(\ell + m, \ell - 2)$. We thus decompose $h_+$ and $h_\times$ onto the basis of such spin-weighted spherical harmonics, which means (Bueno et al. 2007; Kidder 2008)

$$h_+ - i h_\times = \sum_{\ell=2}^{+\infty} \sum_{m=-\ell}^{\ell} h^{\ell m} Y_{-2}^{\ell m}(\theta, \phi). \qquad (185)$$

Using the orthonormality properties of these harmonics we can invert the latter decomposition and obtain the separate modes $h^{\ell m}$ from a surface integral,

$$h^{\ell m} = \int d\Omega [h_+ - i h_\times] \bar{Y}_{-2}^{\ell m}(\theta, \phi), \qquad (186)$$

where the overline refers to the complex conjugation. On the other hand, we can also relate $h^{\ell m}$ directly to the radiative multipole moments $U_L$ and $V_L$. The result is
$$h^{\ell m} = - \frac{G}{\sqrt{2} R c^{\ell+2}} \left[ U^{\ell m} - \frac{i}{c} V^{\ell m} \right], \qquad (187)$$

where $U^{\ell m}$ and $V^{\ell m}$ denote the radiative mass and current moments in standard (non-STF) guise. These are related to the STF moments by

$$U^{\ell m} = \frac{4}{\ell!} \sqrt{\frac{(\ell+1)(\ell+2)}{2\ell(\ell-1)}} \alpha_L^{\ell m} U_L, \qquad (188a)$$

$$V^{\ell m} = -\frac{8}{\ell!} \sqrt{\frac{\ell(\ell+2)}{2(\ell+1)(\ell-1)}} \alpha_L^{\ell m} V_L. \qquad (188b)$$

Here $\alpha_L^{\ell m}$ denotes the STF tensor connecting the usual basis of spherical harmonics $Y^{\ell m}(\theta, \phi)$ to the set of STF tensors $\hat{N}_L = N_{\langle i_1 \cdots i_{\ell m} \rangle}$ with $N_i = (\sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta)$, where the brackets indicate the STF projection. Indeed both $Y^{\ell m}$ and $\hat{N}_L$ are basis of an irreducible representation of weight $\ell$ of the rotation group; the two basis are related by$^{42}$

$$\hat{N}_L(\theta, \phi) = \sum_{m=-\ell}^{\ell} \alpha_L^{\ell m} Y^{\ell m}(\theta, \phi), \qquad (189a)$$

$$Y^{\ell m}(\theta, \phi) = \frac{(2\ell + 1)!!}{4\pi\ell!} \bar{\alpha}_L^{\ell m} \hat{N}_L(\theta, \phi). \qquad (189b)$$

The STF tensorial coefficient $\alpha_L^{\ell m}$ can be computed as $\alpha_L^{\ell m} = \int d\Omega \hat{N}_L \bar{Y}^{\ell m}$. To present it in the best way we introduce the orthonormal triad $(i,j,k) = (\delta_i^1, \delta_i^2, \delta_i^3)$ such that $N = i \sin \theta \cos \phi + j \sin \theta \sin \phi + k \cos \theta$. Posing $m = (i + j)/\sqrt{2}$ we have $m \cdot m = 0$, $m \cdot \bar{m} = 1$ and $m \cdot k = 0$. The indices $L = i_1 \cdots i_\ell$ can be split into $M = i_1 \cdots i_m$ and the complementary $L - M = i_{m+1} \cdots i_\ell$, and the result reads (with explicit STF projection)$^{43}$

$$\alpha_L^{\ell m} = \ell! \sqrt{\frac{4\pi}{2\ell+1} \frac{2^m}{(\ell+m)!(\ell-m)!} \bar{m}_{\langle M} k_{L-M\rangle}}, \qquad (190)$$

for $m \ge 0$, and $\alpha_L^{\ell m} = (-)^m \bar{\alpha}_L^{\ell|m|$}$ for $m < 0$. Note that we have the normalization condition

<footnote>$^{42}$ With the usual conventions $Y^{\ell m}(\theta, \phi) = \sqrt{\frac{2\ell+1}{4\pi}} \frac{P^m}{(\ell+m)!} P^m(\cos \theta) e^{im\phi}$ where $P^m(z)$ is the associated Legendre function.</footnote>

<footnote>$^{43}$ See Eq. (4.7) of Henry et al. (2021) in which a factor $(-)^m$ should be removed. Our definition is related to that given in Eq. (2.12) of Thorne (1980) (see also Kidder 2008) by $Y_L^{\ell m} = \frac{(2\ell+1)!!}{4\pi\ell!} (-)^m \bar{\alpha}_L^{\ell m}$.</footnote>
$$ \alpha_L^{\ell m} \bar{\alpha}_L^{\ell m'} = \frac{4\pi\ell!}{(2\ell + 1)!!} \delta_{mm'} \qquad (191) $$

In Sect. 3.4.5 we shall present the dominant mode $(\ell, m) = (2, 2)$ of gravitational waves from non spinning inspiralling compact binaries up to 4PN order, and all the other sub-dominant modes to 3.5PN order.

## 2.4.2 Gravitational-wave tails and tails of tails

We learned from Theorem 5 the general method which permits the computation of the radiative multipole moments $U_L, V_L$ in terms of the source parameters. In the general case those are described by the six types of source moments $I_L, J_L$ and $W_L, X_L, Y_L, Z_L$ given by (7) and (141). However we also pointed out that an equally general description of the source is by two sets of canonical moments $M_L, S_L$. We shall now show that the relation between $U_L, V_L$ and $M_L, S_L$ includes some hereditary tail effects starting at the relative 1.5PN order. The different sets of moments $\bar{M}_L, \bar{S}_L$ used in the radiative construction of the metric are connected to the moments $M_L, S_L$ in the harmonic metric, see Eq. (5.4) of Trestini et al. (2023). All the results below are given for the moments $M_L, S_L$.

Tails are due to the back-scattering of multipolar waves off the Schwarzschild curvature generated by the total mass monopole M of the source. They correspond to the non-linear interaction between M and the multipole moments $M_L$ and $S_L$ for $\ell \ge 2$, and are given by some “hereditary” or non-local-in-time integrals, extending over the past history of the source. The effect arises at the 1.5PN order and we have (Blanchet and Damour 1992)

$$ U_L(U) = M_L^{(\ell)}(U) + \frac{2GM}{c^3} \int_0^{+\infty} d\tau M_L^{(\ell+2)}(U-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \kappa_\ell \right] + \mathcal{O}(G^2), \qquad (192a) $$

$$ V_L(U) = S_L^{(\ell)}(U) + \frac{2GM}{c^3} \int_0^{+\infty} d\tau S_L^{(\ell+2)}(U-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \pi_\ell \right] + \mathcal{O}(G^2), \qquad (192b) $$

where the constant length scale $b_0$ parametrizes the gauge transformation between radiative and harmonic coordinates, see Eq. (194), and the constants $\kappa_\ell$ and $\pi_\ell$ are given by (Blanchet 1995; Almeida et al. 2021)

$$ \kappa_{\ell} = \frac{2\ell^2 + 5\ell + 4}{\ell(\ell + 1)(\ell + 2)} + H_{\ell-2}, \qquad \pi_{\ell} = \frac{\ell - 1}{\ell(\ell + 1)} + H_{\ell-1}, \qquad (193) $$

where $H_p = \sum_{k=1}^p \frac{1}{k}$ denotes the harmonic number.

From the expression of the linearized gauge vector $\xi_{(1)}^\alpha$ in Eq. (95) the retarded time $U = T - R/c$ in radiative coordinates is related to the retarded time $u = t - r/c$ in harmonic coordinates by
$$U = u - \frac{2GM}{c^3} \ln\left(\frac{r}{b_0}\right) + \mathcal{O}(G^2). \quad (194)$$

Changing $b_0$ simply means shifting the origin of time of the radiative coordinates with respect to that of the harmonic coordinates, which has clearly no physical implication. Inserting $U$ into Eqs. (192) we obtain the radiative moments expressed in terms of “source-rooted” harmonic coordinates, i.e.

$$U_L = M_L^{(\ell)}(u) + \frac{2GM}{c^3} \int_0^{+\infty} d\tau M_L^{(\ell+2)}(u-\tau) \left[ \ln\left(\frac{c\tau}{2r}\right) + \kappa_\ell \right] + \mathcal{O}(G^2), \quad (195)$$

and similarly for $V_L$. This expression no longer depends on the arbitrary constant $b_0$, i.e., we find that $b_0$ gets replaced by $r$, so it behaves logarithmically when $r \to +\infty$. If we now replace the harmonic coordinates $(t, r)$ to some new ones, like for instance “Schwarzschild” coordinates $(t', r')$ such that $t' = t$ and $r' = r + GM/c^2$ (and $u' = u - GM/c^3$), we get

$$U_L = M_L^{(\ell)}(u') + \frac{2GM}{c^3} \int_0^{+\infty} d\tau M_L^{(\ell+2)}(u'-\tau) \left[ \ln\left(\frac{c\tau}{2r'}\right) + \kappa'_\ell \right] + \mathcal{O}(G^2), \quad (196)$$

where $\kappa'_\ell = \kappa_\ell + 1/2$. This shows that the constant $\kappa_\ell$ (and $\pi_\ell$ as well) depends on the choice of source-rooted coordinates $(t, r)$: For instance, we have $\kappa_2 = 11/12$ in harmonic coordinates from Eq. (193), but $\kappa'_2 = 17/12$ in Schwarzschild coordinates (Poisson 1993).

The tail integrals in Eqs. (192) involve all the instants from $-\infty$ in the past up to the current retarded time $U$. However, strictly speaking, they do not extend up to infinite past, since we have assumed in Eq. (60) that the metric is stationary before the date $-T$. The range of integration of the tails is therefore limited *a priori* to the time interval $[-T, U]$. But now, once we have derived the tail integrals, thanks to the latter technical assumption of stationarity in the past, we can argue that the results are in fact valid in more general situations for which the field has never been stationary. We have in mind the case of two bodies moving initially on some unbound (hyperbolic-like) orbit, and which, because of the loss of energy by gravitational radiation, will capture each other to form a gravitationally bound system around time $-T$, entering then into the inspiral regime till the final merger.

In this situation let us check, using a simple Newtonian model for the behaviour of the multipole moment $M_L(U-\tau)$ in the past, when $\tau \to +\infty$, that the tail integrals, when assumed to extend over the whole time interval $[-\infty, U]$, remain perfectly well-defined (i.e., convergent) at the integration bound $\tau = +\infty$. Indeed it was shown by Eder (1983) that the motion of initially free particles interacting gravitationally is given by $x^i(U-\tau) = V^i\tau + W^i \ln \tau + X^i + o(1)$, where $V^i, W^i$ and $X^i$ denote constant vectors, with $W^i = GmV^i/V^3$, and $o(1) \to 0$ when $\tau \to +\infty$.<fnref n="44" /> From that physical assumption we find that the multipole moments behave when $\tau \to +\infty$ like

<footnote n="44">⁴⁴ To simplify we consider the relative motion of two particles with total mass $m$.</footnote>
$$M_L(U - \tau) = A_L \tau^\ell + B_L \tau^{\ell-1} \ln \tau + C_L \tau^{\ell-1} + o(\tau^{\ell-1}), \quad (197)$$

where $A_L$, $B_L$ and $C_L$ are constant tensors. We used the fact that the moment $M_L$ will agree at the Newtonian level with the standard expression for the $\ell$-th mass multipole moment. The appropriate time derivatives of the moment appearing in Eq. (192a) are therefore dominantly like

$$M_L^{(\ell+2)}(U - \tau) = \frac{D_L}{\tau^3} + o(\tau^{-3}), \quad (198)$$

which ensures that the tail integral is convergent. This fact can be regarded as an *a posteriori* justification of our *a priori* too restrictive assumption of stationarity in the past. Thus, this assumption does not seem to yield any physical restriction on the applicability of the final formulas. However, once again, we emphasize that the past-stationarity is appropriate for real astrophysical sources of gravitational waves which have been formed and started to emit at a finite instant in the past.

The formulas (192) can be generalized to include two mass monopoles M interacting with the varying moments $M_L$ or $S_L$. We coin such cubic interactions as “tails-of-tails”. They arise at the 3PN order in the waveform. For the mass moments, including the 1.5PN tail and 3PN tail-of-tail, we have (Blanchet 1998a; Faye et al. 2015)

$$\begin{aligned} U_L(U) = M_L^{(\ell)}(U) &+ \frac{2GM}{c^3} \int_0^{+\infty} d\tau M_L^{(\ell+2)}(U-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \kappa_\ell \right] \\ &+ \frac{2G^2M^2}{c^6} \int_0^{+\infty} d\tau M_L^{(\ell+3)}(U-\tau) \left[ \ln^2\left(\frac{c\tau}{2b_0}\right) + 2\kappa_\ell \ln\left(\frac{c\tau}{2b_0}\right) \right] \\ &+ \frac{1}{2} \beta_\ell^{(m)} \ln\left(\frac{c\tau}{2r_0}\right) + \xi_\ell \right] + \mathcal{O}(G^3), \end{aligned} \quad (199)$$

The novel feature at this order is that there is a dependence not only on the pure gauge constant $b_0$, but also on the regularization scale $r_0$ which was introduced in the regulator (72) of the MPM algorithm. The constant $r_0$ cannot be removed by a coordinate transformation, but it must disappear from the radiative moments (and hence from the physical waveform) once the moments $M_L$ are explicitly related to the parameters of the matter source, e.g. the masses, trajectories and velocities of the particles in a binary system.<fnref n="45" />

Exactly as we have done for the tail term in Eq. (195), we can remove all the $b_0$'s by inserting (194) back into (199) and expanding the result keeping the necessary terms consistently. In the end we find that there remains a $r_0$-dependent term at the 3PN order, namely

$$U_L = M_L^{(\ell)}(u) + \beta_{\ell}^{(m)} \ln\left(\frac{r}{r_0}\right) \left(\frac{GM}{c^3}\right)^2 M_L^{(\ell+2)}(u), \quad (200)$$

plus terms that are independent of $r_0$. It is known from Goldberger and Ross (2010),

<footnote n="45">⁴⁵ In Eq. (199) $\kappa_\ell$ is given in (193) and the numerical constant $\xi_\ell$ is irrelevant for the discussion; but we shall obtain $\xi_2 = \frac{124627}{22050}$ in Eq. (241) below.</footnote>
Goldberger et al. (2014), Almeida et al. (2021) that the scale $r_0$ is a running scale in the sense of renormalization group theory. The coefficient in front of the logarithm $\ln(r/r_0)$ in (200) is the beta-function coefficient $\beta_{\ell}^{(m)}$, where (m) stands for “mass-type”, which is associated to the renormalization of the mass-type multipole moment. The appearance of a logarithm and its associated constant $r_0$ at the 3PN order was pointed out by Anderson et al. (1982). The beta coefficient for general $\ell$ has been obtained from a long computation of the multipole interactions $M \times M \times M_L$ by Blanchet and Damour (1988) [see Eq. (A6) there], with result

$$ \beta_{\ell}^{(m)} = -2 \frac{15\ell^4 + 30\ell^3 + 28\ell^2 + 13\ell + 24}{\ell(\ell+1)(2\ell+3)(2\ell+1)(2\ell-1)}. \quad (201) $$

In particular, $\beta_2^{(m)} = -\frac{214}{105}$ and $\beta_3^{(m)} = -\frac{26}{21}$, see Eqs. (241) and (248c). This result has also been obtained using the EFT approach by Goldberger and Ross (2010), Goldberger et al. (2014) for $\ell = 2$ and by Almeida et al. (2021) for general $\ell$. The apparent dependence of the radiative moment (200) on $r_0$ should in fine disappear. The reason is that when we explicitly compute the mass moment $M_L$ for a given matter source, we will find an extra contribution depending on $r_0$ occurring at the 3PN order which will cancel out the one in Eq. (200). This contribution is contained into the source moment $I_L$ defined in Theorem 7, where the FP regularization will produce a $\ln r_0$ at the 3PN order. This has been explicitly verified for the quadrupole moment of compact binaries, and we observe on the result (448)–(449) below the requested terms depending on $\ln r_0$. Now the canonical moments $M_L$ differ from the source moments $I_L$ by small 2.5PN and 3.5PN terms, see e.g. Equation (251), hence the statement that

$$ M_L = \mu \hat{x}_L - \beta_{\ell}^{(m)} \ln \left( \frac{r_{12}}{r_0} \right) \left( \frac{Gm}{c^3} \right)^2 (\mu \hat{x}_L)^{(2)}, \quad (202) $$

plus terms independent of $r_0$, where $\mu \hat{x}_L$ is the binary’s Newtonian multipole moment, with $r_{12}$ the orbital separation and $m$ the total mass differing from the ADM mass M by small post-Newtonian corrections; so indeed the $r_0$ cancels between (200) and (202).

Concerning current type moments $V_L$ their expression in terms of $S_L$ is the same as (199) with $\pi_{\ell}$ in place of $\kappa_{\ell}$, see Eq. (193), and with the general beta-function coefficient given by (Almeida et al. 2021)

$$ \beta_{\ell}^{(c)} = -2 \frac{60\ell^6 + 180\ell^5 + 135\ell^4 - 574\ell^3 - 368\ell^2 + 221\ell - 360}{\ell(\ell+1)(2\ell+5)(2\ell+3)(2\ell+1)(2\ell-1)(2\ell-3)}, \quad (203) $$

with (c) meaning “current-type”. We have $\beta_2^{(c)} = -\frac{214}{105}$, see Eq. (248f). The previous reasonings leading to Eqs. (200) and (202) apply as well for the current type moments. The cancellation of the renormalization scale $r_0$ has been checked in the 3PN current quadrupole moment for compact binaries (Henry et al. 2021).

The quartic tails-of-tails-of-tails terms arising at the order 4.5PN in the radiative mass quadrupole moment will be given in Eq. (246) below.
### 2.4.3 Formulas to integrate non-linearities

To obtain the explicit results like (192) and (199), we must implement the MPM algorithm presented in Sect. 2.1.3, together with some practical formulas for integrating the wave equation when the source term is a complicated (in general non-local) functional of the multipole moments. Without loss of generality we can consider the wave equation whose source term has a given multipolarity $ℓ$ in STF guise, thus takes the form

$$ \Box \Psi_L = \hat{n}_L S(r, t - r/c), \qquad (204) $$

where $S(r, u)$ is a function of $r = |\mathbf{x}|$ and $u = t - r/c$ which verifies straightforward smoothness properties, is zero when $u \le -T$ [cf. the hypothesis of stationarity in the past, Eq. (60)] and diverges along the set of scale functions $(\ln r)^q/r^p$ ($p, q \in \mathbb{N}$) when $r \to 0$. To deal with the divergence when $r \to 0$ we multiply the right side of (204) by the regulator $\tilde{r}^B \equiv (r/r_0)^B$, apply the retarded d'Alembertian operator $\Box_{\text{ret}}^{-1}$ and select the finite part (FP) in the Laurent expansion when $B \to 0$; remind Eq. (75). Hence

$$ \Psi_L = \text{FP}_{B=0} \Box_{\text{ret}}^{-1} \left[ \tilde{r}^B \hat{n}_L S(r, t - r/c) \right], \qquad (205) $$

which is the unique solution of (204), valid in the sense of distributions and stationary in the past, thus satisfying the no-incoming radiation condition when $r \to +\infty$ with $t + r/c = \text{const.}$

**Theorem 10** (Blanchet and Damour 1986) The solution (205), endowed with definite multipolarity $ℓ$, can be explicitly written as

$$ \Psi_L = \text{FP}_{B=0} \int_{-\infty}^{t-r/c} c \, ds \, \hat{\sigma}_L \left[ \frac{R_B\left(\frac{c(t-s)-r}{2}, s\right) - R_B\left(\frac{c(t-s)+r}{2}, s\right)}{r} \right], \qquad (206) $$

in which we define the auxiliary $B$-dependent function

$$ R_B(\rho, s) \equiv \rho^\ell \int_\alpha^\rho d\sigma \frac{(\rho - \sigma)^\ell}{\ell!} \left(\frac{2}{\sigma}\right)^{\ell-1} \tilde{\sigma}^B S(\sigma, s), \qquad (207) $$

where the regulator is $\tilde{\sigma}^B = (\sigma/r_0)^B$ and the analytic continuation in $B \in \mathbb{C}$ is understood throughout.

The function defined in (207) depends also on the lower bound of the integral denoted $\alpha$; however $\alpha$ can be chosen to be an arbitrary function of $t$ and $s$ because the full solution $\Psi_L$ in (206) is actually independent of $\alpha$. This can be seen from the following alternative form of the solution, which is obtained by plugging Eq. (207) into (206)$^{46}$:

$$ \hat{\sigma}_L \left\{ \frac{(\lambda - r)^p - (\lambda + r)^p}{r} \right\} = 0 \quad \text{when } p = 0, 2, \dots, 2\ell. $$

<footnote>$^{46}$ One uses the useful identity, Eq. (A36) of Blanchet and Damour (1986):</footnote>
$$ \Psi_L = -\frac{1}{\ell!} \mathrm{FP}_{B=0} \int_{-\infty}^{t-r/c} c \, ds \int_{\frac{c(t-s)-r}{2}}^{\frac{c(t-s)+r}{2}} d\sigma \left(\frac{2}{\sigma}\right)^{\ell-1} \tilde{\sigma}^B S(\sigma, s) \\ \times \hat{\mathcal{O}}_L \left[ \frac{1}{r} \left( \frac{c(t-s)-r}{2} \right)^{\ell} \left( \frac{c(t-s)-r}{2} - \sigma \right)^{\ell} \right]. \qquad (208) $$

However, the form of the solution given by (206)–(207) is generally the best in practical applications.

At the quadratic non-linear level we have to consider the simple case where the source term is of the type $S(r, u) = F(u)/r^k$. Here $F(u)$ will be a product of two time derivatives of multipole moments $M_L(u)$ and/or $S_L(u)$. The first case of interest is when the source term has $k=2$; this corresponds to the leading behaviour when $r \to \infty$, $u = \text{const}$. In this case it is easy to check that the retarded integral is convergent—see the standard form (71) of the integral—therefore the limit $B \to 0$ must be finite. Working out this limit one gets (using footnote 46)

$$ \Box_{\text{ret}}^{-1} \left[ \frac{\hat{n}_L}{r^2} F(u) \right] = \frac{(-)^\ell}{2\ell!} \int_r^{+\infty} d\lambda F(t - \lambda/c) \\ \times \hat{\mathcal{O}}_L \left\{ \frac{(\lambda - r)^\ell \ln(\lambda - r) - (\lambda + r)^\ell \ln(\lambda + r)}{r} \right\}, \qquad (209) $$

where we have removed the reference to the finite part at $B=0$, and we note that the length scale $r_0$ drops from the result. The expression (209) is clearly non-local in time, or “hereditary”, and typically describes the tails. There is an interesting connection to the Legendre function of the second kind $Q_\ell(x)$,<fnref n="47" /> since we can rewrite the previous expression as

$$ \Box_{\text{ret}}^{-1} \left[ \frac{\hat{n}_L}{r^2} F(u) \right] = - \frac{\hat{n}_L}{r} \int_r^{+\infty} d\lambda \, Q_\ell \left( \frac{\lambda}{r} \right) F(t - \lambda/c). \qquad (211) $$

We have also the following alternative expression, valid for any function $F(n, u)$ (not only with pure multipolarity $\ell$), and which can be retrieved directly from the retarded integral in its usual form (71),

$$ Q_{\ell}(x) = \frac{1}{2} \int_{-1}^{1} \frac{d\zeta P_{\ell}(\zeta)}{\zeta(x-z)} = \frac{1}{2} P_{\ell}(x) \ln\left(\frac{x+1}{x-1}\right) - \sum_{j=1}^{\ell} \frac{1}{j} P_{\ell-j}(x) P_{j-1}(x). $$

In the complex plane there is a branch cut from $-\infty$ to $1$. The first equality is the Neumann formula for the Legendre function. When $x \to 1^+$ we have $Q_\ell(x) = -\frac{1}{2} \ln\left(\frac{x-1}{2}\right) - H_\ell + o(1)$ where $H_\ell$ is the harmonic number. We also note that

$$ \hat{n}_L Q_\ell(x) = \int \frac{d\Omega'}{4\pi} \frac{\hat{n}'_L}{x - \mathbf{n} \cdot \mathbf{n}'} . \qquad (210) $$

<footnote n="47">⁴⁷ The Legendre function $Q_\ell$ is given in terms of the Legendre polynomial $P_\ell$ by</footnote>
$$ \Box_{\text{ret}}^{-1} \left[ \frac{F(n, u)}{r^2} \right] = - \int_{r}^{+\infty} d\lambda \int \frac{d\Omega'}{4\pi} \frac{F(n', t - \lambda/c)}{\lambda - r \mathbf{n} \cdot \mathbf{n}'} . \quad (212) $$

The expansion at retarded infinity $r \to +\infty$ with $u = \text{const}$ exhibits the usual logarithmic behaviour of harmonic coordinates (see Sect. 2.1.5), namely

$$ \Box_{\text{ret}}^{-1} \left[ \frac{\hat{n}_L}{r^2} F(u) \right] = \frac{\hat{n}_L}{2r} \int_0^{+\infty} d\tau F(u-\tau) \left[ \ln\left(\frac{c\tau}{2r}\right) + 2H_\ell \right] + o\left(\frac{1}{r}\right), \quad (213) $$

where the remainder is indicated as $o(1/r)$ rather than $O(1/r^2)$ because it contains powers of the logarithm of $r$; it could be more accurately written as $o(r^{\ell-2})$ for some $\ell \ll 1$. The result (213) gives directly the tail integral in the form met in Eq. (195).

Consider next the source term $F(u)/r^k$ with $k \ge 3$. Two cases must be distinguished. When $3 \le k \le \ell + 2$ (thus necessarily $\ell \ge 1$ in this case), the limit $B \to 0$ turns out to be finite and one obtains a local, “instantaneous” expression, again independent of the scale $r_0$:

$$ \Box_{\text{ret}}^{-1} \left[ \frac{\hat{n}_L}{r^k} F(u) \right] = - \frac{2^{k-3} (k-3)! (\ell+2-k)!}{(\ell+k-2)!} \hat{n}_L \sum_{j=0}^{k-3} \frac{(\ell+j)!}{2^j j!(\ell-j)!} \frac{F^{(k-3-j)}(u)}{c^{k-3-j} p^{j+1}}. \quad (214) $$

This expression can be checked by applying the d'Alembertian operator on both sides of (214). Then the case $k \ge \ell + 3$ is where the FP procedure shows its power and flexibility. In this case one obtains a non-local expression, which admits a logarithmic-free expansion at retarded infinity, but which depends on the arbitrary scale $r_0$. We have (Blanchet 1998c)

$$ \mathrm{FP}_{B=0} \Box_{\text{ret}}^{-1} \left[ \tilde{r}^B \frac{\hat{n}_L}{r^k} F(u) \right] = \sum_{i=0}^{k-3} \alpha_{k\ell}^i \hat{n}_L \frac{F^{(i)}(u)}{c^i r^{k-i-2}} + \frac{(-)^{k+\ell+1} 2^{k-2} (k-3)!}{(k+\ell-2)!(k-\ell-3)!} \frac{\hat{n}_L}{c^{k-2} r} \int_r^{+\infty} d\lambda \bar{Q}_\ell(\frac{\lambda}{r}, r) F^{(k-2)}(t-\lambda/c). \quad (215) $$

In the first term the coefficients are given by$^{48}$

$$ \alpha_{k\ell}^i = \frac{(-)^{k+\ell} 2^i (k-3)!(k+\ell-i-3)!}{(k+\ell-2)!(k-\ell-3)!(k-i-3)!} \beta_{k\ell}^i, \quad (216a) $$

$$ \beta_{k\ell}^i = \begin{cases} (-)^{k+\ell+i} (k-\ell-i-4)! & \text{when } 0 \le i \le k-\ell-4, \\ \frac{H_{k\ell}^i}{(\ell-k+i+3)!} & \text{when } k-\ell-3 \le i \le k-3, \end{cases} \quad (216b) $$

while in the second term we have introduced the modification of the Legendre

<footnote>$^{48}$ Together with the combination of harmonic numbers:</footnote>

<footnote>$H_{k\ell}^i \equiv H_{k-\ell-3} - H_{\ell+i-k+3} - H_{k-3} + H_{k-3-i} + H_{k+\ell-2} - H_{k+\ell-i-3}$,</footnote>
function defined by $\overline{Q}_m(x, r) \equiv Q_m(x) - \frac{1}{2}P_m(x) \ln(\frac{r}{r_0})$, and which explicitly depends on the scale $r_0$. The interest of this particular combination is that it does not produce any far zone logarithms when $r \to +\infty$ with $u = \text{const}$ (Trestini et al. 2023); hence the solution (215) admits a far zone expansion in simple powers of $1/r$, although it is still non-local in time and depends on the scale $r_0$. This feature of the solution is very important in the construction of the metric in radiative coordinates, see Sect. 2.1.5.

The above formulas permit the control of any quadratic non-linear interactions (e. g., tails). To cubic order (e.g., tails-of-tails) one needs some more involved formulas. In particular we have to integrate the wave equation when the source term $S(r, u)$ is itself a non-local integral of the type (211). Let us show the most difficult of the integration formulas needed for the calculation of tails-of-tails<fnref n="49" />:

$$ \mathrm{FP}_{B=0} \Box_{\mathrm{ret}}^{-1} \left[ \tilde{r}^B \frac{\hat{n}_L}{r^2} \int_r^{+\infty} d\lambda \, Q_m\left(\frac{\lambda}{r}\right) F(t - \lambda/c) \right] = \frac{\hat{n}_L}{r} \int_r^{+\infty} d\lambda \, F^{(-1)}(t - \lambda/c) \\ \times \left\{ Q_\ell\left(\frac{\lambda}{r}\right) \int_1^{\lambda/r} dx \, Q_m(x) \frac{dP_\ell}{dx}(x) + P_\ell\left(\frac{\lambda}{r}\right) \int_{\lambda/r}^{+\infty} dx \, Q_m(x) \frac{dQ_\ell}{dx}(x) \right\}, \quad (217) $$

More formulas, even more involved, are requested for the computation of the quartic tails-of-tails-of-tails and cubic tails-of-memory interactions; these can be found in Marchand et al. (2016), Trestini and Blanchet (2023).

### 2.4.4 The non-linear memory effect

The non-linear memory effect has its physical origin in the re-emission of gravitational waves by linearly generated waves (Blanchet 1990; Christodoulou 1991; Wiseman and Will 1991; Thorne 1992; Blanchet and Damour 1992).<fnref n="50" /> It implies a permanent change in the wave amplitude from before to after the passage of a gravitational wave strain. The effect can be interpreted as due to a finite accumulation of the (mass type) radiative moments during the emission of gravitational waves, producing a transition between the initial and final asymptotic Bondi et al. (1962), Sachs (1962) (BMS) frames, which turn out to be related by a supertranslation (Newman and Penrose 1966; Strominger and Zhiboedov 2016).

During the construction of the radiative metric in Sect. 2.1.5, we found that the non-linear source term for the metric takes to leading order $1/R^2$ in the distance the form of the stress-energy tensor for the massless gravitons, cf. (107). Therefore the

$$ \frac{1}{2} \int_{-1}^{1} \frac{dz P_{\ell}(z)}{\sqrt{(xy-z)^2 - (x^2-1)(y^2-1)}} = Q_{\ell}(x)P_{\ell}(y), $$

where $1 \le y < x$ is assumed; see Appendix A of Blanchet (1998a) for the proof. This relation constitutes a generalization of the Neumann formula (see the footnote 47).

<footnote n="49">⁴⁹ The formula (217) has been obtained using a not so well known mathematical relation between the Legendre functions and polynomials:</footnote>

<footnote n="50">⁵⁰ The effect was initially computed using post-Newtonian theory by Blanchet (1990), pointing out that it is dominantly of order 2.5PN, in contrast to tails at 1.5PN order. This was published later, together with physical interpretation, by Blanchet and Damour (1992).</footnote>
radiative metric will contain a piece which is the retarded integral of the stress-energy tensor of gravitational waves, see Eqs. (110)–(111). This piece is at the origin of the memory-effect, hence we pose

$$h_{\text{mem}}^{\alpha\beta}(T, \mathbf{X}) = \Box_{\text{ret}}^{-1} \left[ \frac{K^{\alpha}K^{\beta}}{R^2} \sigma(U, N) \right], \quad (218)$$

where $U = T - R/c$, $K^\alpha = (1, N)$ and we recall that the energy density in the waves is given by Eq. (108) and is related to the energy loss per steradian, see (109). Following the integration formula (212) we can write the memory term into the more explicit form

$$h_{\text{mem}}^{\alpha\beta} = - \int_{-\infty}^{U} dv \int \frac{d\Omega(\mathbf{n})}{4\pi} \frac{k^\alpha k^\beta}{T - v - RN \cdot \mathbf{n}} \sigma(\mathbf{n}, v), \quad (219)$$

which involves the spherical average with respect to the “dummy” angles $\mathbf{n}$ and we denote $k^\alpha = (1, \mathbf{n})$, while $N$ refers to the unit vector of the field point $X = RN$. As we know from Sect. 2.4.3, the asymptotic expansion when $R \to +\infty$, $U = \text{const}$ of the retarded integral of a source with radial dependence $R^{-2}$ involves some far-zone logarithms of $R$. But we proved in Sect. 2.1.5 that the logarithms can be gauged way, which permits to compute all observables in the asymptotic waveform (radiative moments). Here we shall treat the term (219) directly, following the method of Thorne (1992), Wiseman and Will (1991), which consists of applying first the TT projection operator onto (219). Because the TT projection kills any (linear) gauge term in the $1/R$ part of the metric, this method shortcuts the need of a transformation to radiative coordinates. However, one must be cautious in taking the limit $R \to +\infty$, $U = \text{const}$ using (219), as it is not allowed to work out a leading $1/R$ term from (219) because this term would involve a divergent integral, in accordance with the fact that the leading term is actually $\ln R/R$. But, the divergent parts of the integral cancel out after application of the TT projection, and in the end one recovers the correct result. Thus, we compute<fnref n="51" />

$$(h_{\text{mem}}^{ij})^{\text{TT}} = - \perp_{ijab} (N) \int_{-\infty}^{U} dv \int \frac{d\Omega}{4\pi} \frac{n^a n^b}{T - v - RN \cdot \mathbf{n}} \sigma(\mathbf{n}, v). \quad (220)$$

Next we insert the STF multipole decomposition $\sigma(\mathbf{n}, v) = \sum_{\ell \ge 0} n_\ell \sigma_L(v)$ [with STF coefficients $\sigma_L(v)$] and obtain

$$(h_{\text{mem}}^{ij})^{\text{TT}} = - \perp_{ijab} (N) \sum_{\ell=0}^{+\infty} \int_{-\infty}^{U} dv \sigma_L(v) \int \frac{d\Omega}{4\pi} \frac{n_{abL}}{T - v - RN \cdot \mathbf{n}}. \quad (221)$$

We decompose the product of unit vectors $n_{abL} = n_a n_b n_L$ on the basis of STF tensors, we drop the terms having zero TT-projection, and we express the remaining terms using the Legendre function of the second kind $Q_\ell$ using (210). Restoring the traces on the STF tensors, and dropping further terms having zero TT-projection, we obtain

<footnote n="51">⁵¹ Note that the TT projection as defined here is purely algebraic. Strictly speaking it agrees with the true TT projection only when acting on the leading $1/R$ term.</footnote>
$$ (h_{\text{mem}}^{ij})^{\text{TT}} = -\frac{1}{R} \perp_{ijab} \sum_{\ell=2}^{+\infty} N_{L-2} \int_{-\infty}^{U} \mathrm{d}v \, \sigma_{abL-2}(v) Z_{\ell} \left( \frac{T-v}{R} \right), \quad (222) $$

where $Z_\ell$ denotes the combination of Legendre functions given by

$$ \frac{Z_{\ell}}{\ell(\ell-1)} = \frac{Q_{\ell+2}}{(2\ell+3)(2\ell+1)} - \frac{2Q_{\ell}}{(2\ell+3)(2\ell-1)} + \frac{Q_{\ell-2}}{(2\ell+1)(2\ell-1)}. \quad (223) $$

The limit at future null infinity can now be applied. Indeed it suffices to insert the expansion of $Q_\ell(x)$ when $x \to 1^+$ as given in the footnote 47. As expected the limit is finite, because the terms $\ln(x-1)/2$ in the expansions of the $Q_\ell$'s cancel out. We obtain $Z_\ell = \frac{2}{(\ell+1)(\ell+2)}$ in the limit and this yields (Blanchet and Damour 1992; Blanchet 1998c)

$$ (h_{\text{mem}}^{ij})^{\text{TT}} = -\frac{2}{R} \perp_{ijab} \sum_{\ell=2}^{+\infty} \frac{N_{L-2}}{(\ell+1)(\ell+2)} \int_{-\infty}^{U} \mathrm{d}v \, \sigma_{abL-2}(v) + \mathcal{O}\left(\frac{1}{R^2}\right). \quad (224) $$

Comparing with the TT waveform (178) we see that the memory effect enters as a modification of the mass-type radiative multipole moments

$$ U_{L}^{\text{mem}}(U) = \frac{\ell!}{2(\ell+1)(\ell+2)} \frac{c^{\ell+3}}{G} \int_{-\infty}^{U} \mathrm{d}v \, \sigma_{L}(v), \quad (225) $$

but does not contribute to the current-type moments, i.e. $V_L^{\text{mem}}(U) = 0$. It is to be remarked that the result (225) is exact; it gives the full memory contribution in non-linear GR without approximation (Favata 2009, 2011c).

The convergence of the integral (225) is ensured by the fact that the time derivatives of the moments are zero in the remote past, before the instant $-\mathcal{T}$. The memory property is evident from (225). Even after the system has ceased to emit radiation, in the sense that the $(\ell+1)$-th time derivative of the source moments tend to zero at late time, the emission of gravitational waves will produce a constant (DC) contribution to the gravitational wave amplitude. More explicitly, one reads off for the total GW emission (Blanchet and Damour 1992)

$$ U_L(+\infty) - U_L(-\infty) = M_L^{(\ell)}(+\infty) - M_L^{(\ell)}(-\infty) + \frac{\ell!}{2(\ell+1)(\ell+2)} \frac{c^{\ell+3}}{G} \int_{-\infty}^{+\infty} \mathrm{d}v \, \sigma_L(v), \quad (226) $$

where the first term represents the possible contribution of the “linear” memory which has been discussed by Zel'dovich and Polnarev (1974), Braginsky and Thorne (1987b). Generically the linear memory is present if the source contains freely moving masses in its initial or final state (scattering situation).

Let us now derive an interesting alternative expression for the non-linear memory effect (224). We replace the $\ell$-th order STF coefficient $\sigma_L$ by its expression in terms of the energy distribution using Eq. (109),
$$ \sigma_L(U) = \frac{4(2\ell + 1)!!}{\ell!} \int d\Omega \, \hat{n}_L \frac{dE}{dU d\Omega}(U, \mathbf{n}). \qquad (227) $$

Hence we rewrite (224) in the following form, where the time integration has been explicitly performed,

$$ (h_{\text{mem}}^{ij})^{\text{TT}} = -\frac{8}{R} \perp_{ijab} \sum_{\ell=2}^{+\infty} \frac{(2\ell + 1)!!}{(\ell + 2)!} N_{L-2} \int d\Omega \, \hat{n}_{abL-2} \frac{dE}{d\Omega}(U, \mathbf{n}). \qquad (228) $$

We shall now show that the infinite multipole series in the formula (228) can be summed up in closed analytic form. Indeed consider the following TT projection with respect to the unit vector $N = (N_i)$:

$$ A_{ij}^{\text{TT}} \equiv \perp_{ijab}(\mathbf{N}) \frac{n_a n_b}{1 - \mathbf{N} \cdot \mathbf{n}}, \qquad (229) $$

where the TT projector is given by (3). We expand the denominator in (229) as a power series in $N \cdot n$ (which is convergent as soon as $N \cdot n < 1$):

$$ A_{ij}^{\text{TT}} = \perp_{ijab}(\mathbf{N}) \sum_{\ell=2}^{+\infty} N_{L-2} n_{abL-2}. \qquad (230) $$

The multi-index $L-2$ contains $\ell-2$ indices, namely $i_1 \cdots i_{\ell-2}$. Using Eq. (A.21a) of Blanchet and Damour (1986), we transform the ordinary product of unit vectors $n_{abL-2}$ into a sum of STF products

$$ A_{ij}^{\text{TT}} = \perp_{ijab}(\mathbf{N}) \sum_{\ell=2}^{+\infty} N_{L-2} \sum_{k=0}^{[\frac{\ell}{2}]} \alpha_k^{\ell} \delta_{\{2K} \hat{n}_{L-2K\}} , \qquad (231a) $$

$$ \text{with } \alpha_k^\ell \equiv \frac{(2\ell - 4k + 1)!!}{(2\ell - 2k + 1)!!}. \qquad (231b) $$

The operation over indices $\{\}$ is defined as the un-normalized sum over the smallest set of permutations of $i_1 \cdots i_\ell$ which makes the object symmetrical in $L = i_1 \cdots i_\ell$. The object $\delta_{\{2K} \hat{n}'_{L-2K\}}$ contains $\ell$ indices $L = abL-2$ (with, say, $a = i_\ell$ and $b = i_{\ell-1}$) of which $2k$ are displayed onto the product of $k$ Kronecker symbols denoted $\delta_{2K}$. As an example we have $\delta_{\{ij}n_k\} \equiv \delta_{ij}n_k + \delta_{jk}n_i + \delta_{ki}n_j$.

The point is that, among all the terms composing $\delta_{\{2K} \hat{n}_{L-2K\}}$, we can discard all those which contain either $\delta_{ab}$, $\delta_{ai_p}$ or $\delta_{bi_q}$, since such terms will be cancelled by the TT projection. This is obvious for $\delta_{ab}$; in the two other cases, this results from the fact that, after multiplication by $N_{L-2}$ in Eq. (231a), $\delta_{ai_p}$ or $\delta_{bi_q}$ will yield some $N_i$ or $N_j$. So, we can rewrite the expression (231a) by excluding the indices $ab$ from the operation $\{\}$, which we indicate by underlining the two indices $ab$:
$$A_{ij}^{\text{TT}} = \perp_{ijab}(\mathbf{N}) \sum_{\ell=2}^{+\infty} N_{L-2} \sum_{k=0}^{[\frac{\ell-2}{2}]} \alpha_k^{\ell} \delta_{\{2K} \hat{n}_{abL-2-2K\}}. \quad (232)$$

The next step is to notice from eq. (A.19) of Blanchet and Damour (1986) that the number of terms composing the object $\delta_{\{2K} \hat{n}_{abL-2-2K\}}$ is $\frac{(\ell-2)!}{2^k k! (\ell-2-2k)!}$. When contracted with $N_{L-2}$, all these terms will merge into a single one for each values of $\ell$ and $k$. Thus, we have

$$A_{ij}^{\text{TT}} = \perp_{ijab}(\mathbf{N}) \sum_{\ell=2}^{+\infty} \sum_{k=0}^{[\frac{\ell-2}{2}]} \alpha_k^{\ell} \frac{(\ell-2)!}{2^k k!(\ell-2-2k)!} N_{L-2-2K} \hat{n}_{abL-2-2K}. \quad (233)$$

We change $\ell$ into $\ell + 2k$ and rewrite the previous expression as

$$\begin{aligned} A_{ij}^{\text{TT}} &= \perp_{ijab}(\mathbf{N}) \sum_{\ell=2}^{+\infty} \sum_{k=0}^{+\infty} \alpha_k^{\ell+2k} \frac{(\ell+2k-2)!}{2^k k!(\ell-2)!} N_{L-2} \hat{n}_{abL-2} \\ &= \perp_{ijab}(\mathbf{N}) \sum_{\ell=2}^{+\infty} S_{\ell} N_{L-2} \hat{n}_{abL-2}, \end{aligned} \quad (234)$$

introducing the coefficient $S_\ell$ which is given as an infinite series over all integer values of $k$. However, we find, using the expression of the coefficients $\alpha_k^{\ell+2k}$ deduced from (231b), that this series can actually be re-summed in closed analytic form:

$$S_{\ell} = \frac{(2\ell + 1)!!}{(\ell - 2)!} \sum_{k=0}^{+\infty} \frac{(\ell + 2k - 2)!}{2^k k! (2\ell + 2k + 1)!!} = 2^{\ell+2} \frac{\Gamma(\ell + \frac{3}{2})}{\sqrt{\pi} \Gamma(\ell + 3)} = 2 \frac{(2\ell + 1)!!}{(\ell + 2)!}. \quad (235)$$

Hence we obtain the simple result

$$A_{ij}^{\text{TT}} = 2 \perp_{ijab}(\mathbf{N}) \sum_{\ell=2}^{+\infty} \frac{(2\ell + 1)!!}{(\ell + 2)!} N_{L-2} \hat{n}_{abL-2}, \quad (236)$$

where we recognize the infinite multipolar sum in Eq. (228). Inserting thus the resummed form (229) into (228) we obtain

$$(h_{\text{mem}}^{ij})^{\text{TT}} = -\frac{4}{R} \perp_{ijab}(\mathbf{N}) \int d\Omega \frac{n_a n_b}{1 - \mathbf{N} \cdot \mathbf{n}} \frac{dE}{d\Omega}(U, \mathbf{n}), \quad (237)$$

where the factor $1/(1 - \mathbf{N} \cdot \mathbf{n})$ is reminiscent of the Liénard-Wiechert potentials in the case of massless gravitons. The result (237), which exactly reproduces the result of Thorne (1992), constitutes the best interpretation of the non-linear memory effect as due to the re-radiation of gravitational waves by massless gravitons. Note that Eq. (237) describes the gravitational memory produced not only by a burst of gravitons, but also by a burst of any other type of massless particles, e.g. massless neutrinos (Epstein 1978).
### 2.4.5 Radiative versus canonical moments

We first give the result for the radiative quadrupole moment $U_{ij}$ expressed as a functional of the intermediate canonical moments $M_L, S_L$ up to 4PN order, plus a crucial hereditary contribution to 4.5PN order. Recalling that, at leading order, the radiative moments reduce to the $ℓ$-th time derivatives of the canonical moments, we straightforwardly write

$$U_{ij} = M_{ij}^{(2)}(U) + \sum_{n=3}^{9} U_{ij}^{\text{gPN}}(U) + \mathcal{O}\left(\frac{1}{c^{10}}\right), \qquad (238)$$

with small PN corrections up to 4.5PN, as indicated. Here $U \equiv T - R/c$ is the retarded time in harmonic coordinates. The leading correction at 1.5PN order is the quadratic quadrupolar tail $M \times M_{ij}$ (see Sect. 2.4.2) given by

$$U_{ij}^{1.5\text{PN}} = \frac{2GM}{c^3} \int_{0}^{+\infty} d\tau M_{ij}^{(4)}(U-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \frac{11}{12} \right]. \qquad (239)$$

The length scale $b_0$ is arbitrary and defined by Eq. (194). At the next 2.5PN order arises the leading non-local memory contribution (see Sect. 2.4.4) due to the interaction between two quadrupole moments $M_{ij} \times M_{kl}$, together with associated instantaneous terms (Blanchet 1998c):

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

At the 3PN order appears the first cubic correction, called the tail-of-tail and consisting of the interaction between two masses M and the quadrupole moment $M_{ij}$:

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

We have discussed in Sect. 2.4.2 the occurrence of the length scale $r_0$ (to be distinguished from the gauge constant $b_0$), which plays the role of renormalization scale in the EFT approach (Goldberger and Ross 2010) and originates from the Hadamard regularization scheme in the traditional post-Newtonian calculation (Blanchet 1998a).

The next-order 3.5PN term has a structure similar to the 2.5PN one, i.e. with some memory type integrals and instantaneous terms. The interactions between moments, still of quadratic nature, are however more complicated:
$$
\begin{equation}
\begin{aligned}
U_{ij}^{3.5PN} ={}& \frac{G}{c^7} \left\{ \int_0^{+\infty} d\tau \left[ -\frac{5}{756} M_{ab}^{(4)} M_{ijab}^{(4)} - \frac{32}{63} S_{a\langle i}^{(3)} S_{j\rangle a}^{(3)} \right] (U-\tau) \right. \\
& \left. -\frac{1}{432} M_{ab} M_{ijab}^{(7)} + \frac{1}{432} M_{ab}^{(1)} M_{ijab}^{(6)} - \frac{5}{756} M_{ab}^{(2)} M_{ijab}^{(5)} \right. \\
& \left. +\frac{19}{648} M_{ab}^{(3)} M_{ijab}^{(4)} + \frac{1957}{3024} M_{ab}^{(4)} M_{ijab}^{(3)} + \frac{1685}{1008} M_{ab}^{(5)} M_{ijab}^{(2)} \right. \\
& \left. +\frac{41}{28} M_{ab}^{(6)} M_{ijab}^{(1)} + \frac{91}{216} M_{ab}^{(7)} M_{ijab} - \frac{5}{252} M_{ab\langle i} M_{j\rangle ab}^{(7)} \right. \\
& \left. +\frac{5}{189} M_{ab\langle i}^{(1)} M_{j\rangle ab}^{(6)} + \frac{5}{126} M_{ab\langle i}^{(2)} M_{j\rangle ab}^{(5)} + \frac{5}{2268} M_{ab\langle i}^{(3)} M_{j\rangle ab}^{(4)} \right. \\
& \left. +\frac{5}{42} S_a S_{ija}^{(5)} + \frac{80}{63} S_{a\langle i} S_{j\rangle a}^{(5)} + \frac{16}{63} S_{a\langle i}^{(1)} S_{j\rangle a}^{(4)} - \frac{64}{63} S_{a\langle i}^{(2)} S_{j\rangle a}^{(3)} \right. \\
& \left. + \epsilon_{ac\langle i} \left( \int_0^{+\infty} d\tau \left[ \frac{5}{42} S_{j\rangle cb}^{(4)} M_{ab}^{(3)} - \frac{20}{189} M_{j\rangle cb}^{(4)} S_{ab}^{(3)} \right] (U-\tau) \right. \right. \\
& \left. \left. + \frac{1}{168} S_{j\rangle bc}^{(6)} M_{ab} + \frac{1}{24} S_{j\rangle bc}^{(5)} M_{ab}^{(1)} + \frac{1}{28} S_{j\rangle bc}^{(4)} M_{ab}^{(2)} - \frac{1}{6} S_{j\rangle bc}^{(3)} M_{ab}^{(3)} \right. \right. \\
& \left. + \frac{3}{56} S_{j\rangle bc}^{(2)} M_{ab}^{(4)} + \frac{187}{168} S_{j\rangle bc}^{(1)} M_{ab}^{(5)} + \frac{65}{84} S_{j\rangle bc} M_{ab}^{(6)} + \frac{1}{189} M_{j\rangle bc}^{(6)} S_{ab} \right. \\
& \left. - \frac{1}{189} M_{j\rangle bc}^{(5)} S_{ab}^{(1)} + \frac{10}{189} M_{j\rangle bc}^{(4)} S_{ab}^{(2)} + \frac{32}{189} M_{j\rangle bc}^{(3)} S_{ab}^{(3)} + \frac{65}{189} M_{j\rangle bc}^{(2)} S_{ab}^{(4)} \right. \\
& \left. - \frac{5}{189} M_{j\rangle bc}^{(1)} S_{ab}^{(5)} - \frac{10}{63} M_{j\rangle bc} S_{ab}^{(6)} \right) \right\} .
\end{aligned}
\tag{242}
\end{equation}
$$

At the 4PN order appear the tail-of-memory, made of a cubic interaction involving
two time-varying quadrupole moments and the mass: M × M<sub>ij</sub> × M<sub>kl</sub>, see Fig. 5. In
principle it comes with a double time integration over the two varying quadrupole
moments, but involves also many tail-looking interactions with only one integration
over the quadrupole moments. In addition there is also a cubic interaction involving
the constant angular momentum, i.e. M × S<sub>i</sub> × M<sub>jk</sub>, see the last line of (243). We
have (Trestini and Blanchet 2023)

<figure><img src="image_6.png" /><figcaption>**Fig. 5** From an effective field theory perspective, the tails-of-memory correspond to the three Feynman diagrams shown, which we consider here for illustrative purposes; see Foffa and Sturani (2014), Porto (2016), Levi (2020) for their computational meaning within the EFT</figcaption></figure>
$$
\begin{equation}
\begin{aligned}
U_{ij}^{4\text{PN}} = \frac{G^2 M}{c^8} & \left\{ \frac{8}{7} \int_0^{+\infty} d\rho M_{a\langle i}^{(4)}(U-\rho) \int_0^{+\infty} d\tau M_{j\rangle a}^{(4)}(U-\rho-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) - \frac{1613}{270} \right] \right. \\
& - \frac{20}{7} \int_0^{+\infty} d\tau \left[ M_{a\langle i}^{(3)} M_{j\rangle a}^{(4)} \right] (U-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{3}{2} \ln\left(\frac{c\tau}{2b_0}\right) \right] \\
& - \frac{24}{7} \int_0^{+\infty} d\tau \left[ M_{a\langle i}^{(2)} M_{j\rangle a}^{(5)} \right] (U-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{11}{12} \ln\left(\frac{c\tau}{2b_0}\right) \right] \\
& - \frac{20}{7} \int_0^{+\infty} d\tau \left[ M_{a\langle i}^{(1)} M_{j\rangle a}^{(6)} \right] (U-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{3}{10} \ln\left(\frac{c\tau}{2b_0}\right) \right] \\
& - \frac{8}{7} \int_0^{+\infty} d\tau \left[ M_{a\langle i}^{(7)} M_{j\rangle a}^{(7)} \right] (U-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) - \frac{1}{4} \ln\left(\frac{c\tau}{2b_0}\right) \right] \\
& - \frac{16}{7} M_{a\langle i}^{(2)} \int_0^{+\infty} d\tau M_{j\rangle a}^{(5)} (U-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{27521}{5040} \right] \\
& - \frac{20}{7} M_{a\langle i}^{(1)} \int_0^{+\infty} d\tau M_{j\rangle a}^{(6)} (U-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) + \frac{15511}{3150} \right] \\
& + \frac{4}{7} M_{a\langle i} \int_0^{+\infty} d\tau M_{j\rangle a}^{(7)} (U-\tau) \left[ \ln\left(\frac{c\tau}{2r_0}\right) - \frac{6113}{756} \right] \\
& - \frac{2}{3} S_a \epsilon_{ab\langle i} \int_0^{+\infty} d\tau M_{j\rangle b}^{(6)} (U-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + 2 \ln\left(\frac{c\tau}{2r_0}\right) + \frac{1223}{1890} \right] \right\}.
\end{aligned}
\tag{243}
\end{equation}
$$

It is important to verify that the first line of (243), which corresponds to the “gen-
uine” tail-of-memory, with two integrations over the quadrupole moments, can be
retrieved from the general formalism of Sect. 2.4.4. Indeed from Eq. (225) together
with the general result $\sigma = \frac{1}{2} z_{ij}^{TT} z_{ij}^{TT}$ where $z_{ij}^{TT}$ is nothing but the $1/R$ coefficient in
the asymptotic waveform, see Eqs. (108) and (109), we have

$$
U_{ij}^{\text{mem}}(U) = - \frac{2G}{7c^5} \int_{-\infty}^{U} d\nu U_{a\langle i}^{(1)}(\nu) U_{j\rangle a}^{(1)}(\nu). \quad (244)
$$

To lowest order $U_{ij}$ agrees with $M_{ij}^{(2)}$ hence (244) recovers the leading memory term displayed in (240). To the next order $U_{ij}$ will be modified by the 1.5PN tail contribution, hence we have $U_{ij} = M_{ij}^{(2)} + U_{ij}^{1.5PN}$ where $U_{ij}^{1.5PN}$ is given by (239). Inserting this into (244) and keeping only the leading 4PN correction we obtain

$$
\begin{equation}
\begin{aligned}
U_{ij}^{\text{mem}} = & -\frac{2G}{7c^5} \int_{-\infty}^{U} d\nu M_{a\langle i}^{(3)}(\nu) M_{j\rangle a}^{(3)}(\nu) \\
& -\frac{8G^2 M}{7c^8} \int_{-\infty}^{U} d\nu M_{a\langle i}^{(3)}(\nu) \int_{0}^{+\infty} d\tau M_{j\rangle a}^{(5)}(\nu-\tau) \ln\left(\frac{c\tau}{2}\right),
\end{aligned}
\tag{245}
\end{equation}
$$

where we have discarded the constant $b_0$ and the 11/12 in Eq. (239), which are
irrelevant here. It is straightforward to check that (after an integration by parts) the
4PN term in (245) agrees with the genuine tail-of-memory term given in the first line of (243).

At the next 4.5PN order only the hereditary, non-instantaneous quartic interaction $M^3 \times M_{ij}$, naturally dubbed “tail-of-tail-of-tail”, is known, but it is sufficient to control the full 4.5PN term in the energy flux for compact binaries on quasi-circular orbits. The main reason is that the contributions of instantaneous interactions entering the flux at half-integer PN order (e.g. 4.5PN order) vanish for quasi-circular orbits. The tail-of-tail-of-tail term in the radiative quadrupole moment is given by (Marchand et al. 2016)

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

We shall use this expression to extend the flux at 4.5PN order for circular orbits in Sect. 3.4.4. However note that the full $U_{ij}^{4.5PN}$ also contains quadratic memory interactions, like those entering at 2.5PN and 3.5PN, see (240) and (242). If those are not needed to compute the flux for quasi-circular orbits, they will enter the expression of the gravitational wave (2, 2) mode, therefore restricting the accuracy we can reach when deriving the (2, 2) mode; this is why it will be presented in Sect. 3.4.5 at 4PN order and not 4.5PN order.

Next we present the best known results concerning moments of higher multipolarity. Notably for the flux at 4PN order we need the mass octupole $U_{ijk}$ and current quadrupole $V_{ij}$ moments with a 3PN precision:

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

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

At the 1.5PN order there is the standard tail, at 2.5PN the memory, and at 3PN order the tail-of-tail contribution. We have (Faye et al. 2015)

$$ U_{ijk}^{1.5PN} = \frac{2GM}{c^3} \int_0^{+\infty} d\tau M_{ijk}^{(5)}(U-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \frac{97}{60} \right], \quad (248a) $$
$$
\begin{equation}
\begin{aligned}
U_{ijk}^{2.5\text{PN}} = \frac{G}{c^5} & \left\{ \int_0^{+\infty} d\tau \left[ -\frac{1}{3} M_{a\langle i}^{(3)} M_{jk\rangle a}^{(4)} - \frac{4}{5} \epsilon_{ab\langle i} M_{j\underline{a}}^{(3)} S_{k\rangle b}^{(3)} \right] (U-\tau) \right. \\
& + \frac{1}{4} M_{a\langle i} M_{jk\rangle a}^{(6)} + \frac{1}{4} M_{a\langle i}^{(1)} M_{jk\rangle a}^{(5)} + \frac{1}{4} M_{a\langle i}^{(2)} M_{jk\rangle a}^{(4)} - \frac{4}{3} M_{a\langle i}^{(3)} M_{jk\rangle a}^{(3)} \\
& - \frac{9}{4} M_{a\langle i}^{(4)} M_{jk\rangle a}^{(2)} - \frac{3}{4} M_{a\langle i}^{(5)} M_{jk\rangle a}^{(1)} + \frac{1}{12} M_{a\langle i}^{(6)} M_{jk\rangle a} + \frac{12}{5} S_{\langle i} S_{jk\rangle}^{(4)} \right. \\
& \left. + \epsilon_{ab\langle i} \left[ \frac{9}{5} M_{j\underline{a}} S_{k\rangle b}^{(5)} + \frac{27}{5} M_{j\underline{a}}^{(1)} S_{k\rangle b}^{(4)} + \frac{8}{5} M_{j\underline{a}}^{(2)} S_{k\rangle b}^{(3)} + \frac{12}{5} M_{j\underline{a}}^{(3)} S_{k\rangle b}^{(2)} \right. \right. \\
& \left. \left. + \frac{3}{5} M_{j\underline{a}}^{(4)} S_{k\rangle b}^{(1)} + \frac{1}{5} M_{j\underline{a}}^{(5)} S_{k\rangle b} + \frac{9}{20} M_{jk\rangle a}^{(5)} S_b \right] \right\},
\end{aligned}
\tag{248b}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
U_{ijk}^{3\text{PN}} = \frac{2G^2 M^2}{c^6} & \int_0^{+\infty} d\tau M_{ijk}^{(6)} (U-\tau) \left[ \ln^2\left(\frac{c\tau}{2b_0}\right) + \frac{97}{30} \ln\left(\frac{c\tau}{2b_0}\right) \right. \\
& \left. - \frac{13}{21} \ln\left(\frac{c\tau}{2r_0}\right) + \frac{13283}{8820} \right],
\end{aligned}
\tag{248c}
\end{equation}
$$

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

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

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

$$
\begin{equation}
\begin{aligned}
V_{ij}^{3\text{PN}} = \frac{2G^2 M^2}{c^6} & \int_0^{+\infty} d\tau S_{ij}^{(5)} (U-\tau) \left[ \ln^2\left(\frac{c\tau}{2b_0}\right) + \frac{7}{3} \ln\left(\frac{c\tau}{2b_0}\right) \right. \\
& \left. - \frac{107}{105} \ln\left(\frac{c\tau}{2r_0}\right) - \frac{13127}{11025} \right].
\end{aligned}
\tag{248f}
\end{equation}
$$

Finally we list all known results concerning higher order multipole moments, taken
from Faye et al. (2015). For mass type moments:
$$ U_{ijkl}^{1.5PN} = \frac{G}{c^3} \left\{ 2M \int_{0}^{+\infty} d\tau M_{ijkl}^{(6)}(U-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \frac{59}{30} \right] \right. \\ \left. + \frac{2}{5} \int_{0}^{+\infty} d\tau \left[ M_{\langle ij \rangle}^{(3)} M_{kl\rangle}^{(3)} \right] (U-\tau) \right. \\ \left. - \frac{21}{5} M_{\langle ij \rangle} M_{kl\rangle}^{(5)} - \frac{63}{5} M_{\langle ij \rangle}^{(1)} M_{kl\rangle}^{(4)} - \frac{102}{5} M_{\langle ij \rangle}^{(2)} M_{kl\rangle}^{(3)} \right\}, \quad (249a) $$

$$ U_{ijkl}^{2.5PN} = \frac{G}{c^5} \left\{ \begin{aligned} & \int_0^{+\infty} d\tau \left[ \frac{12}{55} M_{a\langle i}^{(4)} M_{jkl\rangle a}^{(4)} - \frac{14}{99} M_{a\langle ij}^{(4)} M_{kl\rangle a}^{(4)} + \frac{32}{45} S_{\langle ij}^{(3)} S_{kl\rangle}^{(3)} \right. \\ & \left. + \epsilon_{ab\langle i} \left( -\frac{4}{5} M_{j\underline{a}}^{(3)} S_{kl\rangle b}^{(4)} + \frac{32}{45} S_{j\underline{a}}^{(3)} M_{kl\rangle b}^{(4)} \right) \right] (U-\tau) \\ & + \frac{7}{55} M_{a\langle i} M_{jkl\rangle a}^{(7)} + \frac{7}{55} M_{a\langle i}^{(1)} M_{jkl\rangle a}^{(6)} + \frac{1}{25} M_{a\langle i}^{(2)} M_{jkl\rangle a}^{(5)} \\ & - \frac{28}{11} M_{a\langle i}^{(3)} M_{jkl\rangle a}^{(4)} - \frac{273}{55} M_{a\langle i}^{(4)} M_{jkl\rangle a}^{(3)} - \frac{203}{55} M_{a\langle i}^{(5)} M_{jkl\rangle a}^{(2)} \\ & - \frac{49}{55} M_{a\langle i}^{(6)} M_{jkl\rangle a}^{(1)} + \frac{14}{275} M_{a\langle i}^{(7)} M_{jkl\rangle a} + \frac{14}{33} M_{a\langle ij} M_{kl\rangle a}^{(7)} \\ & + \frac{37}{33} M_{a\langle ij}^{(1)} M_{kl\rangle a}^{(6)} + \frac{9}{11} M_{a\langle ij}^{(2)} M_{kl\rangle a}^{(5)} + \frac{8}{33} M_{a\langle ij}^{(3)} M_{kl\rangle a}^{(4)} + \frac{9}{5} S_{\langle ij} S_{jkl\rangle}^{(5)} \\ & + \frac{16}{5} S_{\langle ij} S_{kl\rangle}^{(5)} + \frac{48}{5} S_{\langle ij}^{(1)} S_{kl\rangle}^{(4)} + \frac{32}{5} S_{\langle ij}^{(2)} S_{kl\rangle}^{(3)} \\ & + \epsilon_{ab\langle i} \left( -\frac{3}{5} M_{j\underline{a}} S_{kl\rangle b}^{(6)} - \frac{63}{25} M_{j\underline{a}}^{(1)} S_{kl\rangle b}^{(5)} + \frac{3}{5} M_{j\underline{a}}^{(2)} S_{kl\rangle b}^{(4)} + \frac{18}{5} M_{j\underline{a}}^{(3)} S_{kl\rangle b}^{(3)} \right. \\ & \left. + \frac{9}{5} M_{j\underline{a}}^{(4)} S_{kl\rangle b}^{(2)} + \frac{3}{5} M_{j\underline{a}}^{(5)} S_{kl\rangle b}^{(1)} + \frac{3}{25} M_{j\underline{a}}^{(6)} S_{kl\rangle b} - \frac{8}{15} S_{j\underline{a}} M_{kl\rangle b}^{(6)} \right. \\ & \left. - \frac{24}{25} S_{j\underline{a}}^{(1)} M_{kl\rangle b}^{(5)} - \frac{8}{5} S_{j\underline{a}}^{(2)} M_{kl\rangle b}^{(4)} + \frac{16}{3} S_{j\underline{a}}^{(3)} M_{kl\rangle b}^{(3)} + \frac{72}{5} S_{j\underline{a}}^{(4)} M_{kl\rangle b}^{(2)} \right. \\ & \left. + \frac{56}{5} S_{j\underline{a}}^{(5)} M_{kl\rangle b}^{(1)} + \frac{232}{75} S_{j\underline{a}}^{(6)} M_{kl\rangle b} + \frac{29}{75} M_{jkl\rangle a}^{(6)} S_b \right\}, \quad (249b) $$
$$ U_{ijklm}^{1.5PN} = \frac{G}{c^3} \left\{ 2M \int_0^{+\infty} d\tau M_{ijklm}^{(7)}(U-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \frac{232}{105} \right] \right. \\ \left. + \frac{20}{21} \int_0^{+\infty} d\tau \left[ M_{\langle ij}^{(3)} M_{klm\rangle}^{(4)} \right] (U-\tau) \right. \\ \left. - \frac{15}{7} M_{\langle ij} M_{klm\rangle}^{(6)} - \frac{41}{7} M_{\langle ij}^{(1)} M_{klm\rangle}^{(5)} - \frac{120}{7} M_{\langle ij}^{(2)} M_{klm\rangle}^{(4)} - \frac{710}{21} M_{\langle ij}^{(3)} M_{klm\rangle}^{(3)} \right. \\ \left. - \frac{265}{7} M_{\langle ij}^{(4)} M_{klm\rangle}^{(2)} - \frac{155}{7} M_{\langle ij}^{(5)} M_{klm\rangle}^{(1)} - \frac{34}{7} M_{\langle ij}^{(6)} M_{klm\rangle}^{(4)} \right\}, \quad (249c) $$

$$ U_{ijklmn}^{1.5PN} = \frac{G}{c^3} \left\{ 2M \int_0^{+\infty} d\tau M_{ijklmn}^{(8)}(U-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \frac{403}{168} \right] \right. \\ \left. + \int_0^{+\infty} d\tau \left[ \frac{5}{7} M_{\langle ijk}^{(4)} M_{lmn\rangle}^{(4)} - \frac{15}{14} M_{\langle ij}^{(4)} M_{klmn\rangle}^{(4)} \right] (U-\tau) \right. \\ \left. - \frac{45}{28} M_{\langle ij} M_{klmn\rangle}^{(7)} - \frac{111}{28} M_{\langle ij}^{(1)} M_{klmn\rangle}^{(6)} - \frac{561}{28} M_{\langle ij}^{(2)} M_{klmn\rangle}^{(5)} \right. \\ \left. - \frac{1595}{28} M_{\langle ij}^{(3)} M_{klmn\rangle}^{(4)} - \frac{2505}{28} M_{\langle ij}^{(4)} M_{klmn\rangle}^{(3)} - \frac{2115}{28} M_{\langle ij}^{(5)} M_{klmn\rangle}^{(2)} \right. \\ \left. - \frac{909}{28} M_{\langle ij}^{(6)} M_{klmn\rangle}^{(1)} - \frac{159}{28} M_{\langle ij}^{(7)} M_{klmn\rangle}^{(4)} - \frac{15}{7} M_{\langle ijk} M_{lmn\rangle}^{(7)} \right. \\ \left. - \frac{75}{7} M_{\langle ijk}^{(1)} M_{lmn\rangle}^{(6)} - \frac{135}{7} M_{\langle ijk}^{(2)} M_{lmn\rangle}^{(5)} - \frac{505}{21} M_{\langle ijk}^{(3)} M_{lmn\rangle}^{(4)} \right\}. \quad (249d) $$

For current type moments:

$$ V_{ijk}^{1.5PN} = \frac{G}{c^3} \left\{ 2M \int_0^{+\infty} d\tau S_{ijk}^{(5)}(U-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \frac{5}{3} \right] \right. \\ \left. - 2M_{\langle ij}^{(4)} S_{k\rangle} - \frac{1}{10} \epsilon_{ab\langle i} M_{j\underline{a}} M_{k\rangle b}^{(5)} + \frac{1}{2} \epsilon_{ab\langle i} M_{j\underline{a}}^{(1)} M_{k\rangle b}^{(4)} \right\}, \quad (249e) $$
$$
V_{ijk}^{2.5PN} = \frac{G}{c^5} \left\{ \begin{aligned}
& \frac{1}{12} M_{a\langle i} S_{jk\rangle a}^{(6)} + \frac{1}{12} M_{a\langle i}^{(1)} S_{jk\rangle a}^{(5)} + \frac{5}{12} M_{a\langle i}^{(2)} S_{jk\rangle a}^{(4)} + \frac{35}{12} M_{a\langle i}^{(4)} S_{jk\rangle a}^{(2)} \\
& + \frac{49}{12} M_{a\langle i}^{(5)} S_{jk\rangle a}^{(1)} + \frac{19}{12} M_{a\langle i}^{(6)} S_{jk\rangle a} + \frac{2}{27} S_{a\langle i} M_{jk\rangle a}^{(6)} + \frac{10}{27} S_{a\langle i}^{(1)} M_{jk\rangle a}^{(5)} \\
& + \frac{2}{27} S_{a\langle i}^{(2)} M_{jk\rangle a}^{(4)} + \frac{8}{9} S_{a\langle i}^{(3)} M_{jk\rangle a}^{(3)} - \frac{10}{27} S_{a\langle i}^{(4)} M_{jk\rangle a}^{(2)} \\
& - \frac{38}{27} S_{a\langle i}^{(5)} M_{jk\rangle a}^{(1)} - \frac{2}{3} S_{a\langle i}^{(6)} M_{jk\rangle a} - \frac{1}{60} S_a M_{ijk}^{(6)} \\
& + \epsilon_{ab\langle i} \left( -\frac{1}{180} M_{jk\rangle ac}^{(7)} M_{bc} + \frac{11}{900} M_{jk\rangle ac}^{(6)} M_{bc}^{(1)} + \frac{7}{300} M_{jk\rangle ac}^{(5)} M_{bc}^{(2)} \right) \quad (249f) \\
& - \frac{37}{270} M_{jk\rangle ac}^{(4)} M_{bc}^{(3)} - \frac{191}{180} M_{jk\rangle ac}^{(3)} M_{bc}^{(4)} - \frac{65}{36} M_{jk\rangle ac}^{(2)} M_{bc}^{(5)} \\
& - \frac{367}{300} M_{jk\rangle ac}^{(1)} M_{bc}^{(6)} - \frac{133}{450} M_{jk\rangle ac} M_{bc}^{(7)} + \frac{1}{27} M_{j\underline{a}c} M_{k\rangle bc}^{(7)} \\
& + \frac{5}{162} M_{j\underline{a}c}^{(1)} M_{k\rangle bc}^{(6)} - \frac{5}{162} M_{j\underline{a}c}^{(2)} M_{k\rangle bc}^{(5)} - \frac{1}{81} M_{j\underline{a}c}^{(3)} M_{k\rangle bc}^{(4)} \\
& - \frac{11}{20} S_{jk\rangle b}^{(5)} S_a - \frac{88}{45} S_{j\underline{a}} S_{k\rangle b}^{(5)} - \frac{40}{9} S_{j\underline{a}}^{(1)} S_{k\rangle b}^{(4)} - \frac{32}{9} S_{j\underline{a}}^{(2)} S_{k\rangle b}^{(3)} \right\}
\end{aligned} \right.
$$

$$
V_{ijkl}^{1.5PN} = \frac{G}{c^3} \left\{ \begin{aligned}
& 2M \int_{0}^{+\infty} d\tau S_{ijkl}^{(6)}(U-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \frac{119}{60} \right] \\
& - \frac{11}{6} M_{\langle ij} S_{kl\rangle}^{(5)} - \frac{25}{6} M_{\langle ij}^{(1)} S_{kl\rangle}^{(4)} - \frac{25}{3} M_{\langle ij}^{(2)} S_{kl\rangle}^{(3)} - \frac{35}{3} M_{\langle ij}^{(3)} S_{kl\rangle}^{(2)} \\
& - \frac{65}{6} M_{\langle ij}^{(4)} S_{kl\rangle}^{(1)} - \frac{19}{6} M_{\langle ij}^{(5)} S_{kl\rangle}^{(2)} - \frac{11}{12} M_{\langle ij}^{(5)} S_{l\rangle}^{(3)} \\
& + \epsilon_{ab\langle i} \left[ \frac{1}{12} M_{j\underline{a}} M_{kl\rangle b}^{(6)} + \frac{37}{60} M_{j\underline{a}}^{(1)} M_{kl\rangle b}^{(5)} - \frac{5}{12} M_{j\underline{a}}^{(2)} M_{kl\rangle b}^{(4)} - \frac{5}{6} M_{j\underline{a}}^{(3)} M_{kl\rangle b}^{(3)} \right. \\
& \qquad \left. - \frac{11}{12} M_{j\underline{a}}^{(4)} M_{kl\rangle b}^{(2)} - \frac{1}{12} M_{j\underline{a}}^{(5)} M_{kl\rangle b}^{(1)} + \frac{3}{60} M_{j\underline{a}}^{(6)} M_{kl\rangle b}^{(3)} \right]
\end{aligned} \right\} \quad (249g)
$$
$$
\begin{equation}
\begin{aligned}
V_{ijklm}^{1.5PN} ={}& \frac{G}{c^3} \left\{ 2M \int_0^{+\infty} d\tau S_{ijklm}^{(7)}(U-\tau) \left[ \ln\left(\frac{c\tau}{2b_0}\right) + \frac{133}{60} \right] \right. \\
& \left. - \frac{3}{2} M_{\langle ij} S_{klm\rangle}^{(6)} - \frac{33}{10} M_{\langle ij}^{(1)} S_{klm\rangle}^{(5)} - 12 M_{\langle ij}^{(2)} S_{klm\rangle}^{(4)} - 27 M_{\langle ij}^{(3)} S_{klm\rangle}^{(3)} \right. \\
& \left. - \frac{69}{2} M_{\langle ij}^{(4)} S_{klm\rangle}^{(2)} - \frac{39}{2} M_{\langle ij}^{(5)} S_{klm\rangle}^{(1)} - \frac{21}{5} M_{\langle ij}^{(6)} S_{klm\rangle} - \frac{4}{3} S_{\langle ij} M_{klm\rangle}^{(6)} \right. \\
& \left. - \frac{76}{15} S_{\langle ij}^{(1)} M_{klm\rangle}^{(5)} - \frac{16}{3} S_{\langle ij}^{(2)} M_{klm\rangle}^{(4)} - 8 S_{\langle ij}^{(3)} M_{klm\rangle}^{(3)} - \frac{28}{3} S_{\langle ij}^{(4)} M_{klm\rangle}^{(2)} \right. \\
& \left. - \frac{20}{3} S_{\langle ij}^{(5)} M_{klm\rangle}^{(1)} - \frac{8}{5} S_{\langle ij}^{(6)} M_{klm\rangle} - \frac{3}{5} S_{\langle ij} M_{jklm}^{(6)} \right. \\
& \left. + \epsilon_{ab\langle i} \left( \frac{1}{14} M_{ja} M_{klm\rangle b}^{(7)} + \frac{1}{2} M_{ja}^{(1)} M_{klm\rangle b}^{(6)} - \frac{3}{5} M_{ja}^{(2)} M_{klm\rangle b}^{(5)} \right. \right. \\
& \left. - \frac{4}{3} M_{ja}^{(3)} M_{klm\rangle b}^{(4)} - \frac{3}{2} M_{ja}^{(4)} M_{klm\rangle b}^{(3)} - \frac{1}{2} M_{ja}^{(5)} M_{klm\rangle b}^{(2)} + \frac{1}{35} M_{ja}^{(7)} M_{klm\rangle b} \right. \\
& \left. + \frac{1}{7} M_{jka} M_{lm\rangle b}^{(7)} + \frac{2}{3} M_{jka}^{(1)} M_{lm\rangle b}^{(6)} + \frac{4}{3} M_{jka}^{(2)} M_{lm\rangle b}^{(5)} + \frac{1}{3} M_{jka}^{(3)} M_{lm\rangle b}^{(4)} \right) \right\}.
\end{aligned}
\tag{249h}
\end{equation}
$$

For all the other multipole moments, we can only state that there is agreement between the radiative and canonical moments up to 1.5PN order, namely

$$
U_L(U) = M_L^{(\ell)}(U) + \mathcal{O}\left(\frac{1}{c^3}\right), \qquad (250a)
$$

$$
V_L(U) = S_L^{(\ell)}(U) + \mathcal{O}\left(\frac{1}{c^3}\right). \qquad (250b)
$$

2.4.6 Canonical versus source moments

In a second stage of the formalism, we must express the canonical moments
${M_L, S_L}$ in terms of the six types of source moments ${I_L, J_L, W_L, X_L, Y_L, Z_L}$.
For the control of the flux and (2, 2) mode in the waveform up to 4PN order, we need to
relate the canonical quadrupole moment $M_{ij}$ to the corresponding source quadrupole
moment $I_{ij}$ up to that accuracy. The relation was known at 3.5PN order (Faye et al.
2012) and has been extended at 4PN order (Blanchet et al. 2022, 2023b) with result$^{52}$

<footnote>$^{52}$ This result consists of purely quadratic interactions. However, in a first calculation of this relation, performed in ordinary 3 dimensions using the Hadamard regularization scheme (for the IR divergences), some cubic interactions were found at the 4PN order. See Eq. (1.1) in Blanchet et al. (2022) where these interactions are of the types $M \times M \times W_{ij}$, $M \times M \times Y_{ij}$ and $M \times W \times I_{ij}$. However further work (Blanchet et al. 2023b) showed that within the dimensional regularization scheme such interactions are cancelled, so that the relation (251), derived in the context of dimensional regularization, holds with exactly the same definition for the “renormalized” source quadrupole moment $I_{ij}$ as obtained by Larrouturou et al. (2022a, b).</footnote>
$$
\begin{equation}
\begin{aligned}
M_{ij} = I_{ij} &+ \frac{4G}{c^5} [W^{(2)}I_{ij} - W^{(1)}I_{ij}^{(1)}] \\
&+ \frac{4G}{c^7} \left\{ \begin{aligned}[t]
& \frac{4}{7} W_{a\langle i \, j\rangle a}^{(1)} I_{ij}^{(3)} + \frac{6}{7} W_{a\langle i \, j\rangle a}^{(4)} - \frac{1}{7} Y_{a\langle i \, j\rangle a}^{(3)} - Y_{a\langle i \, j\rangle a} - 2X I_{ij}^{(3)} \\
& - \frac{5}{21} W_a^{(4)} I_{ija} + \frac{1}{63} W_a^{(3)} I_{ija}^{(1)} - \frac{25}{21} Y_a^{(3)} I_{ija} - \frac{22}{63} Y_a^{(2)} I_{ija}^{(1)} + \frac{5}{63} Y_a^{(1)} I_{ija}^{(2)} \\
& + 2W^{(3)} W_{ij} + 2W^{(2)} W_{ij}^{(1)} - \frac{4}{3} W_{\langle i} W_{j\rangle}^{(3)} + 2W^{(2)} Y_{ij} - 4W_{\langle i} Y_{j\rangle}^{(2)} \\
& + \epsilon_{ab\langle i} \left[ \frac{1}{3} I_{j\rangle a} Z_b^{(3)} - I_{j\rangle a}^{(3)} Z_b + \frac{4}{9} J_{j\rangle a} W_b^{(3)} - \frac{4}{9} J_{j\rangle a} Y_b^{(2)} + \frac{8}{9} J_{j\rangle a} Y_b^{(1)} \right]
\end{aligned} \right\} \\
&+ \mathcal{O}\left(\frac{1}{c^9}\right).
\end{aligned}
\tag{251}
\end{equation}
$$

Here, for instance, W denotes the monopole moment associated with the moment W_L, and Y_i is the dipole moment corresponding to Y_L. Notice that the difference between the canonical and source moments starts at the relatively high 2.5PN order.

We need also the relations for the higher-order moments, notably for the control of the full waveform, which contain similarly some correction terms starting at the 2.5PN order. Here we present the most up-to-date results Blanchet et al. (2008), Faye et al. (2012):

$$
M_{ijk} = I_{ijk} + \frac{4G}{c^5} \left[ W^{(2)} I_{ijk} - W^{(1)} I_{ijk}^{(1)} + 3 I_{\langle ij \rangle k} Y_{\langle k \rangle}^{(1)} \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (252a)
$$

$$
M_{ijkl} = I_{ijkl} + \frac{4G}{c^5} \left[ -W^{(1)} I_{ijkl} + W^{(2)} I_{ijkl} + 4Y_{\langle i}^{(1)} I_{jkl} \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (252\text{b})
$$

$$
S_{ij} = J_{ij} + \frac{2G}{c^5} \left[ \epsilon_{ab\langle i} \left( -I_{j\rangle b}^{(3)} W_a - 2I_{j\rangle b} Y_a^{(2)} + I_{j\rangle b}^{(1)} Y_a^{(1)} \right) \right. \\
\left. \qquad + 3J_{\langle i} Y_{j\rangle}^{(1)} - 2J_{ij}^{(1)} W^{(1)} \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \qquad (252c)
$$

$$
S_{ijk} = J_{ijk} + \frac{4G}{c^5} \left[ \epsilon_{ab\langle i} \left( -\frac{1}{3} I_{jk\rangle a}^{(1)} Y_b^{(1)} + I_{jk\rangle a} Y_b^{(2)} + I_{ja}^{(3)} W_{k\rangle b} \right) \right. \\
\left. - W^{(1)} J_{ijk}^{(1)} + \frac{8}{3} Y_{\langle i}^{(1)} J_{jk\rangle} \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \qquad (252\text{d})
$$

For all the other moments M_L, S_L we can just say that they agree with their source counterparts I_L, J_L up to 2.5PN order,

$$
M_L = I_L + \mathcal{O}\left(\frac{1}{c^5}\right), \quad S_L = J_L + \mathcal{O}\left(\frac{1}{c^5}\right). \qquad (253\text{a})
$$

Finally, combining the previous formulas (252) together with all the results of
Sect. 2.4.5 we control the full waveform in terms of the basic source and gauge moments up to order 3.5PN, with the dominant mode controlled to 4PN.

With all the formulas in Sects. 2.4.5 and 2.4.6 we have related the radiative moments {$U_L, V_L$} parametrizing the waveform at future null infinity (178) to the six types of source multipole moments {$I_L, J_L, W_L, X_L, Y_L, Z_L$}. The general expression of the source moments in terms of the source parameters (actually, the PN expansion of the matter + gravitation pseudo-tensor) has been obtained in Eqs. (7) and (141). Hence at this stage the problem of the wave generation by an arbitrary post-Newtonian source is solved. What is missing is the explicit implementation of the formalism for an actual model of the source. We come to grips with the important physical case of compact binary systems in the next section.

# 3 Compact binary systems

As reviewed in Sect. 1 the appropriate theoretical description of the early inspiral of compact binaries is by two structureless point-particles, characterized solely by their masses $m_1$ and $m_2$ and possibly their spins. Furthermore, when approaching the merger some extra parameters describing the internal structure of the bodies (and defined within the effective field theory) become necessary and in particular characterize the tidal interaction between the two objects (see Sect. 1.1.3). Here we compute the equations of motion and radiation of compact binary systems within this model for the compact objects.

## 3.1 Regularization of the field of point particles

One aim is to compute the metric (and its gradient needed in the equations of motion) at some high post-Newtonian order for a system of two point-like particles. A priori one is not allowed to use directly some metric expressions like in Eqs. (161), which have been derived under the assumption of a continuous (smooth) matter distribution. Applying them to a system of point particles, we find that the metric becomes divergent at the location of the particles, i.e., when $x \rightarrow y_1(t)$ or $y_2(t)$, where $y_a(t)$ denote the two trajectories. Consequently, we must supplement the calculation by a prescription for how to remove the infinite self-field of point particles. At this stage different choices for a “self-field” regularization are possible:

1. The Hadamard self-field regularization has long been used and proved to be very convenient for doing practical computations (in particular, by computer). However it suffers from the important drawback of yielding some physical ambiguity parameters, which cannot be determined within this regularization, starting essentially at the 3PN order;

2. The dimensional self-field regularization is an extremely powerful regularization which is free of ambiguities (at least up to the 4PN level). In particular it has been used to uniquely fix the values of the ambiguity parameters coming from Hadamard's regularization. However, in high post-Newtonian approximations, dimensional regularization cannot be implemented for an arbitrary space
dimension $d \in \mathbb{C}$, but has to be applied in the limit $\epsilon = d - 3 \to 0$, in which limit it is in fact connected in a precise way to the Hadamard regularization.

We expect the results to be unique and independent of the employed regularization, in agreement with the physical intuition. This can be justified from the effacing principle of general relativity (Damour 1983a)—namely that the internal structure of the compact bodies makes a contribution only at the formal 5PN approximation. However note that all PN results derived below have been obtained using dimensional regularization (or have been shown to be equivalent to those of dimensional regularization).

In principle the use of a self-field regularization, be it dimensional or based on Hadamard's partie finie, signals a somewhat unsatisfactory situation on the physical point of view, because, ideally, we would like to perform a complete calculation valid for extended bodies, taking into account the details of the internal structure of the bodies (energy density, pressure, internal velocity field, etc.). By considering the limit where the radii of the objects tend to zero, one should recover the same result as obtained by means of the point-mass regularization. This would demonstrate the suitability of the regularization. This program was undertaken at the 2PN order by Kopeikin (1985), Grishchuk and Kopeikin (1986) who derived the equations of motion of two extended fluid balls, and obtained equations of motion depending only on the two masses $m_1$ and $m_2$ of the compact bodies.<fnref n="53" /> At the 3PN or even 4PN order we expect that the extended-body program should give exactly the same result as given by dimensional regularization. Ideally, this should also be confirmed by independent and more physical methods like those of Thorne and Hartle (1985), Kopeikin (1988), Damour et al. (1991).

Another body of work, in a sense more physical than the formal use of a self-field regularization, is the one of Itoh et al. (2000), Itoh et al. (2001), Itoh and Futamase (2003), Itoh (2004, 2009), where the equations of motion are derived in harmonic coordinates by means of a particular variant of the famous “surface-integral” method (i.e., *à la* Einstein et al. 1938). The aim is to describe extended relativistic compact binary systems in the strong-field point particle limit (Futamase 1987). This approach is interesting because it is based on the physical notion of extended compact bodies in general relativity, and by-passes the problem of the self-field regularization. The end result of Itoh and Futamase (2003); Itoh (2004) is in agreement with the harmonic coordinates 3PN equations of motion (Blanchet and Faye 2000a, 2001a; Blanchet et al. 2004a).

### 3.1.1 Hadamard's regularization

In many practical computations it is useful to employ the Hadamard regularization (Hadamard 1932; Schwartz 1978; Sellier 1994), in particular it recovers in a simple way the $d \to 3$ limit of dimensional regularization in the absence of poles. See Blanchet and Faye (2000b), Schäfer and Jaranowski (2024) for reviews; Blanchet and Faye (2000b) proposed an “extended” Hadamard regularization based on a

<footnote n="53">⁵³ See some comments on the latter works by Damour (1987b), pp. 168–169.</footnote>
generalization of the theory of Schwartz (1978) distributions with generalized singular functions or pseudo-functions.

With Hadamard's regularization the PN iteration of the metric and equations of motion for a system of compact objects (point-particles) proceeds along the line of Sect. 2.3.3, i.e. computing the hierarchy of post-Newtonian potentials in three dimensions, Eqs. (159) and (163). At each step one has typically to solve a Poisson equation of the type $ΔP = F$, where the function $F(\mathbf{x})$ is smooth ($C^\infty$) on $\mathbb{R}^3$ except for the two points $\mathbf{y}_a$ at which it becomes singular (UV type divergence; for simplicity we consider only two point particles). We assume that $F$ admits a power-like singular expansion around 1 of the general structure (with $p \in \mathbb{Z}$)$^{54}$

$$ \forall \mathcal{N} \in \mathbb{N}, \quad F(\mathbf{x}) = \sum_{p_0 \le p \le \mathcal{N}} r_1^p f_p(\mathbf{n}_1) + o(r_1^\mathcal{N}), \qquad (254) $$

and similarly for the other point 2. Here $r_1 = |\mathbf{x} - \mathbf{y}_1| \to 0$, and the coefficients $f_p$ of the various powers of $r_1$ depend on the unit direction $\mathbf{n}_1 = (\mathbf{x} - \mathbf{y}_1)/r_1$ of approach to the singular point. The powers $p$ of $r_1$ are bounded from below ($p_0 \le p$ with $p_0 \in \mathbb{Z}$). The coefficients $f_p$ (and $2f_p$) for which $p \le -1$ are referred to as the singular coefficients of $F$. If both $F$ and $G$ admit expansions of the type (254) so does the ordinary product $FG$, as well as the ordinary gradient $\partial_i F$. We define the Hadamard partie finie of $F$ at the location of the point 1 where it is singular as

$$ (F)_1 = \int \frac{d\Omega_1}{4\pi} f_0(\mathbf{n}_1), \qquad (255) $$

where $d\Omega_1 = d\Omega(\mathbf{n}_1)$ denotes the solid angle element centered on $\mathbf{y}_1$ and of direction $\mathbf{n}_1$. Notice that because of the angular integration in Eq. (255), the Hadamard partie finie is “non-distributive” in the sense that

$$ (FG)_1 \neq (F)_1(G)_1 \quad \text{in general.} \qquad (256) $$

The non-distributivity of Hadamard's partie finie is the main source of the appearance of ambiguity parameters with this regularization.

The second notion of Hadamard partie finie (Pf) concerns that of the integral $\int d^3\mathbf{x} F$, which is generically divergent at the location of the two singular points $\mathbf{y}_1$ and $\mathbf{y}_2$ (we assume that the integral converges at infinity). It is defined by removing the singular parts of the expansion near the singularities:

$$ \mathrm{Pf}_{s_1, s_2} \int d^3 \mathbf{x} F(\mathbf{x}) = \lim_{s \to 0} \left\{ \int_{S(s)} d^3 \mathbf{x} F + 4\pi \sum_{p+3 \le -1} \frac{s^{p+3}}{p+3} \left(\frac{F}{r_1^p}\right)_1 + 4\pi \ln\left(\frac{s}{s_1}\right) (r_1^3 F)_1 + 1 \leftrightarrow 2 \right\}. \qquad (257) $$

The first term integrates over a domain $S(s)$ defined as $\mathbb{R}^3$ from which the two

<footnote>$^{54}$ The function $F(\mathbf{x})$ depends also on (coordinate) time $t$, through for instance its dependence on the velocities $v_1(t)$ and $v_2(t)$, but the time $t$ is purely “spectator” in the regularization process, and thus will not be indicated. See the footnote 26 for the definition of the Landau symbol $o$ for remainders.</footnote>
spherical balls $r_1 < s$ and $r_2 < s$ of radius $s$ and centered on the two singularities, i.e. $B_1(s)$ and $B_2(s)$, are excised: $S(s) = \mathbb{R}^3 \setminus B_1(s) \cup B_2(s)$. The other terms, where the value of a function at point 1 takes the meaning (255), are precisely such that they cancel out the divergent part of the first term in the limit where $s \to 0$ (the symbol $1 \leftrightarrow 2$ means the same terms but corresponding to point 2). The Hadamard partie-finie integral depends on two strictly positive constants $s_1$ and $s_2$, associated with the logarithms present in Eq. (257). The appearance of these logarithms and arbitrary constants signals in fact the ambiguousness and incompleteness of the regularization.

The Hadamard partie-finie integral in the form given by (257) is rather difficult to evaluate, because of the integration over the complicated volume $S(s)$. Fortunately, there exist several alternative expressions of the Hadamard partie finie, much better suited for practical computations. One is based on a double analytic continuation, with two complex parameters $\alpha, \beta \in \mathbb{C}$,

$$ \underset{s_1,s_2}{\mathrm{Pf}} \int d^3\mathbf{x} F = \underset{\alpha \to 0}{\mathrm{FP}} \int d^3\mathbf{x} \left(\frac{r_1}{s_1}\right)^\alpha \left(\frac{r_2}{s_2}\right)^\beta F, \qquad (258) $$
$$ \beta \to 0 $$

where the constants $s_1$ and $s_2$ are the same as in the definition (257). The operation FP means taking the finite parts in the Laurent expansions when $\alpha \to 0$ and $\beta \to 0$ successively; furthermore one can apply the two limits in whatever order. See Blanchet and Faye (2000b) for other expressions of the partie-finie integral (257).

Within the Hadamard regularization the definitions (257) or (258) are applied to Poisson integrals, solving the Poisson equation $\Delta P = F$. In this case we consider the partie-finie integral<fnref n="55" />

$$ P(\mathbf{x}') = -\frac{1}{4\pi} \underset{s_1,s_2}{\mathrm{Pf}} \int \frac{d^3\mathbf{x}}{|\mathbf{x} - \mathbf{x}'|} F(\mathbf{x}), \qquad (259) $$

where the field point $\mathbf{x}'$ is located outside the singular points $\mathbf{y}_1$ and $\mathbf{y}_2$, and is just “spectator” in the regularization process. However we want to evaluate the value of the regularized Poisson integral (259) at the singular point $\mathbf{y}_1$. Clearly the limit will not be continuous at this point, because the Poisson integral admits a singular expansion when $\mathbf{x}' \to \mathbf{y}_1$. Because of the integration in (258) this expansion will be more general than in Eq. (254), typically involving a logarithm of $r'_1 \equiv |\mathbf{x}' - \mathbf{y}_1| \to 0$. Therefore the Hadamard partie finie of $P$ at 1 in the sense of Eq. (255) is actually infinite since it a priori contains $\ln r'_1 = -\infty$. Nevertheless, in Hadamard’s regularization we must keep the “constant” $\ln r'_1$ in all calculations; such drawback will be solved when we add the corrections due to the dimensional regularization in Sect. 3.1.2. Finally, the Hadamard partie finie at the point 1 in the sense of Eq. (255), of the Poisson integral (259) is given by (Blanchet and Faye 2000b)

<footnote n="55">⁵⁵ To high post-Newtonian order, the Poisson integral will also diverge at infinity. In addition to the UV regularization, one must also apply an IR regularization, which will be given in the context of Hadamard’s regularization by Eq. (145).</footnote>
$$ (P)_1 = -\frac{1}{4\pi} \mathrm{Pf} \int \frac{d^3\mathbf{x}}{r_1} F(\mathbf{x}) + \left[ \ln\left(\frac{r'_1}{s_1}\right) - 1 \right] (r_1^2 F)_1, \quad (260) $$

with the “infinite” constant $r'_1 = |\mathbf{x}' - \mathbf{y}_1|$ appearing in the second term. Actually the constant $s_1$ cancels out from the two terms in the right side of (260) so the partie finie depends only on $r'_1$ and $s_2$.

As we said the Hadamard regularization yields some ambiguous results for the computation of certain integrals at the 3PN order. This was noticed by Jaranowski and Schäfer (1998, 1999, 2000) in their computation of the equations of motion within the ADM-Hamiltonian formulation of general relativity. They showed that there are two and only two types of ambiguous terms in the 3PN Hamiltonian, which were then parametrized by two unknown numerical coefficients. Similarly Blanchet and Faye (2000a, 2001a) obtained the 3PN equations of motion in harmonic coordinates complete except for one and only one unknown numerical constant. More precisely they showed that the infinite “constant” $r'_1$ in Eq. (260), together with $r'_2$ for the other particle, are unphysical, in the sense that they can be removed by a coordinate transformation. However the constant $s_2$ which remains in the result (260), and the $s_1$ for the other particle, are the source for the appearance of a single physical ambiguity parameter (they called $\lambda$) which was defined from the unknown ratios $r'_1/s_1$ and $r'_2/s_2$; see (7.9) of Blanchet and Faye (2001a). Concerning the binary’s radiation field, the same phenomenon occurs, with the appearance of some Hadamard regularization ambiguity parameters starting at the 3PN order in the binary’s mass quadrupole moment (Blanchet et al. 2002b), and therefore ambiguity parameters appearing in the 3PN flux and phase evolution as well (Blanchet et al. 2002a). All these problems, which concern UV type divergences associated with the model of point masses, have been resolved, and the ambiguity parameters uniquely fixed, thanks to the dimensional regularization.

### 3.1.2 Dimensional regularization of the equations of motion

Dimensional regularization was invented as a means to preserve the gauge symmetry of perturbative quantum field theories (’t Hooft and Veltman 1972; Bollini and Giambiagi 1972; Breitenlohner and Maison 1977; Collins 1984). Our basic problem here is to respect the gauge symmetry associated with the diffeomorphism invariance of the classical general relativistic description of interacting point masses. Hence, we use dimensional regularization not merely as a trick to compute some particular integrals which would otherwise be divergent, but as a powerful tool for solving in a consistent way the Einstein field equations with singular point-mass sources, while preserving its crucial symmetries. In particular, the dimensional regularization is able to correctly keep track of the global Lorentz-Poincaré invariance of the gravitational field of isolated systems. And of course, dimensional regularization is also a crucial ingredient of the EFT approach to equations of motion and gravitational radiation (Goldberger and Rothstein 2006; Galley et al. 2016; Foffa et al. 2019).

The Einstein field equations in $D = d + 1$ space-time dimensions take exactly the same form as in 4 dimensions, say, in harmonic coordinates $H^\alpha = 0$,
$$ \Box h^{\alpha\beta} = \frac{16\pi G^{(d)}}{c^4} \tau^{\alpha\beta} = \frac{16\pi G^{(d)}}{c^4} |g| T^{\alpha\beta} + A_{\text{harm}}^{\alpha\beta}. \qquad (261) $$

The differences are: (i) the box operator $\Box$ now denotes the flat space-time d'Alembertian operator in $D$ dimensions with signature $(-1, 1, 1, \dots)$; (ii) the gravitational constant $G^{(d)}$ is related to the usual three-dimensional Newton's constant $G$ by

$$ G^{(d)} = G \ell_0^{d-3}, \qquad (262) $$

where $\ell_0$ denotes an arbitrary length scale; (iii) the gravitational source term in harmonic coordinates acquires a $d$-dependent coefficient:

$$ \begin{aligned} A_{\text{harm}}^{\alpha\beta} = & -h^{\mu\nu} \partial_{\mu\nu}^2 h^{\alpha\beta} + \partial_{\mu} h^{\alpha\nu} \partial_{\nu} h^{\beta\mu} + \frac{1}{2} g^{\alpha\beta} g_{\mu\nu} \partial_{\lambda} h^{\mu\tau} \partial_{\tau} h^{\nu\lambda} \\ & - g^{\alpha\mu} g_{\nu\tau} \partial_{\lambda} h^{\beta\tau} \partial_{\mu} h^{\nu\lambda} - g^{\beta\mu} g_{\nu\tau} \partial_{\lambda} h^{\alpha\tau} \partial_{\mu} h^{\nu\lambda} + g_{\mu\nu} g^{\lambda\tau} \partial_{\lambda} h^{\alpha\mu} \partial_{\tau} h^{\beta\nu} \\ & + \frac{1}{4} (2g^{\alpha\mu} g^{\beta\nu} - g^{\alpha\beta} g^{\mu\nu}) (g_{\lambda\tau} g_{\epsilon\pi} - \frac{1}{d-1} g_{\tau\epsilon} g_{\lambda\pi}) \partial_{\mu} h^{\lambda\pi} \partial_{\nu} h^{\epsilon\nu}. \end{aligned} \qquad (263) $$

When $d=3$ we recover Eq. (56). In the following we assume, as required in dimensional regularization, that the dimension of space is a complex number, $d \in \mathbb{C}$, and prove many results by invoking complex analytic continuation in $d$. We often pose $\epsilon \equiv d - 3 \in \mathbb{C}$ and systematically consider the limit $\epsilon \to 0$.

We start by writing in a general way the retarded solution of the flat scalar wave equation in $D$ space-time dimensions (thus, with $\mathbf{x} \in \mathbb{R}^d$),

$$ \Box h(\mathbf{x}, t) = S(\mathbf{x}, t), \qquad (264) $$

where the source will represent a generic term in the matter tensor or the gravitational source term (263). The retarded Green's function of the scalar wave equation, satisfying $\Box G_{\text{ret}}(\mathbf{x}, t) = \delta(ct) \delta^{(d)}(\mathbf{x})$ where $\delta^{(d)}$ is the Dirac function in $d$ dimensions, reads in $(\mathbf{x}, t)$ space$^{56}$

$$ G_{\text{ret}}(\mathbf{x}, t) = -\frac{\tilde{k}}{4\pi} \frac{\theta(ct - r)}{r^{d-1}} \gamma_{\frac{1-d}{2}}\left(\frac{ct}{r}\right), \qquad (265) $$

$$ G_{\text{ret}}(\mathbf{x}, t) = -\frac{\theta(t)}{(2\pi)^{d/2}} \int_0^{+\infty} dk \left(\frac{k}{r}\right)^{\frac{d}{2}-1} \sin(ckt) J_{\frac{d}{2}-1}(kr), $$

where $J_{\frac{d}{2}-1}(kr)$ is the usual Bessel function.

<footnote>$^{56}$ In the Fourier domain the retarded Green's function reads (see Cardoso et al. 2003)</footnote>
where $\theta(ct - r)$ is the usual Heaviside step function, and $\tilde{k}$ is a constant related to the usual Eulerian $\Gamma$-function by$^{57}$

$$ \tilde{k} = \frac{\Gamma\left(\frac{d-2}{2}\right)}{\pi^{\frac{d-2}{2}}}. \qquad (266) $$

The corresponding advanced Green's function $G_{\text{adv}}(\mathbf{x}, t)$ is given by the same expression but with $\theta(-ct - r)$ in place of $\theta(ct - r)$. Notice that the Green's function is in fact a function of $t$ and $r = |\mathbf{x}|$ only. In Eq. (265) we have introduced the function $\gamma_s(z)$ defined for any $s \in \mathbb{C}$ and $|z| \ge 1$ by

$$ \gamma_s(z) = \frac{2\sqrt{\pi}}{\Gamma(s+1)\Gamma(-s-\frac{1}{2})} (z^2 - 1)^s, \qquad (267) $$

where the normalisation is such that $\int_1^{+\infty} dz \gamma_s(z) = 1$. In the Green's function this expression gives back the usual Minkowskian time like interval since $z^2 - 1 \propto c^2t^2 - r^2$ for $z = \frac{c}{r}$. With those notations the retarded solution of the wave equation (264) is given by

$$ h = \Box_{\text{ret}}^{-1} S = -\frac{\tilde{k}}{4\pi} \int_{1}^{+\infty} dz \gamma_{\frac{1}{2}d}(z) \int d^d \mathbf{x}' \frac{S(\mathbf{x}', t - z|\mathbf{x} - \mathbf{x}'|/c)}{|\mathbf{x} - \mathbf{x}'|^{d-2}}. \qquad (268) $$

With dimensional regularization the metric and equations of motion will be parametrized by means of some retarded potentials $V$, $V_i$, $\hat{W}_{ij}$, ..., which are straightforward $d$-dimensional generalizations of the potentials used in three dimensions and which were defined in Sect. 2.3.3. Those are obtained by post-Newtonian iteration of the $d$-dimensional field equations, starting from appropriate definitions of matter source densities generalizing Eqs. (160):

$$ \sigma = \frac{2}{d-1} \frac{(d-2)T^{00} + T^{ii}}{c^2}, \qquad \sigma_i = \frac{T^{0i}}{c}, \qquad \sigma_{ij} = T^{ij}. \qquad (269) $$

As a result all the expressions of Sect. 2.3.3 acquire some explicit $d$-dependent coefficients. For instance we find (Blanchet et al. 2004a)

$$ V = \Box_{\text{ret}}^{-1} [-4\pi G^{(d)} \sigma], \qquad (270a) $$

$$ \hat{W}_{ij} = \Box_{\text{ret}}^{-1} \left[ -4\pi G^{(d)} \left( \sigma_{ij} - \delta_{ij} \frac{\sigma_{kk}}{d-2} \right) - \frac{d-1}{2(d-2)} \partial_i V \partial_j V \right], \qquad (270b) $$

where $\Box_{\text{ret}}^{-1}$ is the retarded integral operator (268). Here again, in higher

$$ \tilde{k} \Omega_{d-1} = \frac{4\pi}{d-2}. $$

<footnote>$^{57}$ We have $\lim_{d \to 3} \tilde{k} = 1$. Notice that $\tilde{k}$ is closely linked to the volume $\Omega_{d-1}$ of the sphere with $d-1$ dimensions (i.e. embedded into Euclidean $d$-dimensional space):</footnote>
approximations there will be also IR divergences and one should really employ the d-dimensional version of Eq. (155).

As reviewed in Sect. 3.1.1, the generic functions $F(\mathbf{x})$ we have to deal with in 3 dimensions, are smooth on $\mathbb{R}^3$ except at $\mathbf{y}_1$ and $\mathbf{y}_2$, around which they admit singular Laurent-type expansions given by Eq. (254). In $d$ spatial dimensions, there is an analogue of the function $F$, which results from the same PN iteration process performed in $d$ dimensions as we just outlined. Let us call this function $F^{(d)}(\mathbf{x})$, where $\mathbf{x} \in \mathbb{R}^d$. When $r_1 \to 0$ the function $F^{(d)}$ admits a singular expansion which is more involved than in 3 dimensions:

$$F^{(d)}(\mathbf{x}) = \sum_{\substack{p_0 \le p \le \mathcal{N} \\ q_0 \le q \le q_1}} r_1^{p+q_\varepsilon} f_{p,q}^{(\varepsilon)}(\mathbf{n}_1) + o(r_1^\mathcal{N}). \quad (271)$$

The coefficients $f_{1,p,q}^{(\varepsilon)}(\mathbf{n}_1)$ depend on $\varepsilon = d-3$, and the powers of $r_1$ involve the relative integers $p$ and $q$ whose values are limited by some $p_0$, $q_0$, $q_1 \in \mathbb{Z}$ as indicated. Here we will be interested in functions $F^{(d)}(\mathbf{x})$ which have no poles as $\varepsilon \to 0$ (this will always be the case up to 4PN order). Therefore, we can deduce from the fact that $F^{(d)}(\mathbf{x})$ is continuous at $d=3$ the constraint

$$\sum_{q=q_0}^{q_1} f_{1,p,q}^{(\varepsilon=0)}(\mathbf{n}_1) = f_p(\mathbf{n}_1). \quad (272)$$

For the problem of equations of motion, we have to deal with the regularization of Poisson integrals, or iterated Poisson integrals (and their gradients), of the generic function $F^{(d)}$. The Poisson integral of $F^{(d)}$, in $d$ dimensions, is given by the Green's function for the Laplace operator,<fnref n="58" />

$$P^{(d)}(\mathbf{x}') = \Delta^{-1}[F^{(d)}(\mathbf{x})] = -\frac{\tilde{k}}{4\pi} \int \frac{d^d\mathbf{x}}{|\mathbf{x} - \mathbf{x}'|^{d-2}} F^{(d)}(\mathbf{x}). \quad (273)$$

We want to evaluate the Poisson integral at the point $\mathbf{x}' = \mathbf{y}_1$ where it is singular; this is quite easy in dimensional regularization, because the nice properties of analytic continuation in $d$ allow simply to get $[P^{(d)}(\mathbf{x}')]_{\mathbf{x}'=\mathbf{y}_1}$ by replacing $\mathbf{x}'$ by $\mathbf{y}_1$ in the integrand of (273). So we simply have

$$P^{(d)}(\mathbf{y}_1) = -\frac{\tilde{k}}{4\pi} \int \frac{d^d\mathbf{x}}{r_1^{d-2}} F^{(d)}(\mathbf{x}). \quad (274)$$

When using dimensional regularization, remind that it is not possible in general to obtain exact expressions valid in any dimension $d$, but only in the limit where $\varepsilon \to 0$. Thus in general we need first to control the singular case $\varepsilon = 0$ by applying first the

<footnote n="58">⁵⁸ Namely, $u(\mathbf{x}) \equiv \tilde{k} r^{2-d}$, satisfying $\Delta u = -4\pi \delta^{(d)}(\mathbf{x})$, where $\tilde{k}$ is defined in (266).</footnote>
Hadamard regularization.<fnref n="59" /> In particular this permits to obtain all the contributions which do not pose a problem. Then we add the specific corrections (typically in the form of poles $1/ε$) due to the use of dimensional regularization. Thus the main step of our strategy consists of computing, in the limit $ε → 0$, the *difference* between the $d$-dimensional result and its Hadamard 3-dimensional counterpart, which will constitute the correction to the Hadamard result. The point is that such difference can be computed solely from the vicinity of the singular points, as the regular parts cancel out in the difference for $ε = 0$.

The Hadamard counterpart of the Poisson integral (274) evaluated at 1 was denoted ($P$)$_1$ and obtained in Eq. (260). Denoting the difference between the two regularizations by means of the script letter $\mathcal{D}$, we pose (for what concerns the point 1)

$$ \mathcal{D}P_1 \equiv P^{(d)}(\mathbf{y}_1) - (P)_1. \qquad (275) $$

That is, $\mathcal{D}P_1$ is what we shall have to add to the Hadamard-regularization result in order to get the $d$-dimensional result. With this method, we only compute the first two terms of the Laurent expansion of $\mathcal{D}P_1$ when $ε → 0$, say $\mathcal{D}P_1 = a_{-1}ε^{-1} + a_0 + \mathcal{O}(ε)$. This is all the information we need to clear up the ambiguities with Hadamard's regularization. The result depends solely on the singular coefficients in the expansion of $F^{(d)}$ around the singularities, Eq. (271). It reads (Blanchet et al. 2004a)

$$ \mathcal{D}P_1 = -\frac{\tilde{k}}{4\pi} \left\{ \sum_{q_0 \le q \le q_1} \left( \frac{1}{q\epsilon} + \ln r'_1 - 1 \right) \langle f_1^{(\epsilon)} \rangle + \sum_{q_0 \le q \le q_1} \left( \frac{1}{(q+1)\epsilon} + \ln s_2 \right) \sum_{\ell=0}^{+\infty} \frac{(-)^{\ell}}{\ell!} \partial_L \left( \frac{1}{r_{12}^{1+\epsilon}} \right) \langle n_2^L f_{-\ell-3,q}^{(\epsilon)} \rangle + \mathcal{O}(\epsilon) \right\}. \qquad (276) $$

where the angular average is defined in $d$ dimensions as

$$ \langle f_1^{(\epsilon)} \rangle = \int d\Omega_1^{(d-1)} f_1^{(\epsilon)}(\mathbf{n}_1), \qquad (277) $$

with $d\Omega_1^{(d-1)}$ the solid angle element around the direction $\mathbf{n}_1$; recall that $\Omega^{(d-1)} = 2\pi^d/\Gamma(d/2)$ is the volume of the unit sphere with $d-1$ dimensions. The first term of (276) corresponds to the contribution of the singularity 1, while the second term concerns the other singularity 2, where we has to perform a multipolar re-expansion of the factor $r_1^{2-d}$ around the point $\mathbf{y}_2$, which reads in STF form<fnref n="60" />

<footnote n="59">⁵⁹ Actually we must employ a crude version of Hadamard’s regularization [in particular avoiding the non-distributivity issue, Eq. (256)] coined the “pure-Hadamard-Schwartz” regularization (Blanchet et al. 2004a).</footnote>

<footnote n="60">⁶⁰ The expansion is STF because $\Delta r^{2-d} = 0$ in $d$ dimensions (in the sense of functions). The Appendix B of Blanchet et al. (2004a) gives a compendium of $d$-dimensional formulae on STF expansions.</footnote>
$$r_1^{2-d} = \sum_{\ell=0}^{+\infty} \frac{(-1)^{\ell}}{\ell!} \partial_L \left( \frac{1}{r_{12}^{1+\epsilon}} \right) r_2^{\ell} n_2^L. \qquad (278)$$

See Eq. (4.24) of Blanchet et al. (2004a) for the similar formula concerning the gradient of the Poisson potential—evidently crucial for the problem of equations of motion.

Note that (276) depends on the two “constants” $\ln r'_1$ and $\ln s_2$ appearing in (260) [remind that $s_1$ disappears from (260)]. One can check that these $\ln r'_1$ and $\ln s_2$ will exactly cancel out the same constants present in the Hadamard calculation, so that the dimensionally regularized equations of motion are finally free of these constants. However the coefficients $f_{p,q}^{(\epsilon)}$ (and $1 \leftrightarrow 2$) in $d$ dimensions depend on the length scale $\ell_0$ introduced in Eq. (262). Taking this dependence into account one can verify that the constants in (276) appear only in the combinations $\ln(r'_1/\ell_0)$ and $\ln(s_2/\ell_0)$.

Next we outline the way we obtain, starting from the computation of the previous difference, the equations of motion in dimensional regularization. As we already alluded, the common “core” of both dimensional and Hadamard regularizations is the pure-Hadamard-Schwartz (pHS) regularization (Blanchet et al. 2004a). This regularization is a specific, minimal Hadamard-type regularization of integrals, based on the partie finie integral (257), together with a minimal treatment of “contact” terms, in which the definition (257) is applied separately to each of the elementary potentials $V$, $V_i$, etc. (and gradients) which enter the PN metric. Furthermore, the regularization of a product of these potentials is assumed to be distributive, i.e., $(FG)_1 = (F)_1(G)_1$ in the case where $F$ and $G$ are given by such elementary potentials; this is thus in contrast with Eq. (256). The pHS regularization also assumes the use of standard Schwartz (1978) distributional derivatives. The interest of the pHS regularization is that the dimensional regularization is equal to it plus the “difference”; see Eq. (279).

Therefore, summing up all contributions, the final result for the acceleration of body 1, regularized with dimensional regularization (DR), reads

$$a_1^{\text{DR}} = a_1^{\text{pHS}} + \mathcal{D}a_1, \qquad (279)$$

where $\mathcal{D}a_1 \sim \sum \mathcal{D}P_1$ is the sum of all the individual differences of Poisson or Poisson-like integrals as computed in Eq. (276) (together with the difference for gradients), and where $a_1^{\text{pHS}}$ is the acceleration obtained with the pHS “core” regularization described above. As obtained in (276), the difference makes a contribution only when a term generates a pole $\sim 1/\epsilon$, in which case the dimensional regularization adds an extra contribution, made of the pole and the finite part associated with the pole—we consistently neglect all terms $\mathcal{O}(\epsilon)$. One must then be especially wary of combinations of terms whose pole parts finally cancel but whose dimensionally regularized finite parts generally do not; these must be evaluated with care. The poles
$\sim 1/\epsilon$ arise at order 3PN in the equations of motion in harmonic coordinates.<fnref n="61" /> Then we have

**Theorem 11** The pole part $\propto 1/\epsilon$ of the dimensional-regularized acceleration (279) can be re-absorbed (i.e. renormalized) into some shifts of the two “bare” world-lines: $\mathbf{y}_1 \to \mathbf{y}_1 + \xi_1$ and $\mathbf{y}_2 \to \mathbf{y}_2 + \xi_2$ (with $\xi_a \sim 1/\epsilon$ say), so that the result, expressed in terms of the “dressed” quantities, is finite:

$$ \mathbf{a}_1 = \lim_{\epsilon \to 0} [\mathbf{a}_1^{\text{DR}} + \delta_{\xi} \mathbf{a}_1], \qquad (280) $$

with $\delta_{\xi} \mathbf{a}_1$ denoting the contribution of the shifts. Furthermore the renormalized acceleration (280) is equivalent to the Hadamard-regularized acceleration if and only if the ambiguous parameters are uniquely fixed.

All the previous discussion concerns the ultra-violet (UV) type divergences appearing in the equations of motion. However starting at the 4PN order the equations of motion also develop infra-red (IR) type divergences. Those too have to be treated with the dimensional regularization. We give more details on the DR of the IR divergences in Sect. 3.1.4.

### 3.1.3 Dimensional regularization of the radiation field

We now address the similar problem concerning the binary’s radiation field—say, to 3PN or 4PN order beyond the Einstein quadrupole formalism. As for the equations of motion, when using the Hadamard regularization for the calculation of the binary’s multipole moments, some ambiguity parameters arise (Blanchet et al. 2002b; Blanchet and Iyer 2005). Again, those are fixed using the dimensional regularization.

Let us first review the multipole expansion in $d$ dimensions. To this end, we consider again the scalar wave equation (264) in $D$-dimensions, but in which the source $S = \Box h$ has a compact spatial support. The multipole expansion $\mathcal{M}(h)$ is the retarded homogeneous solution of the wave equation, i.e. $\Box \mathcal{M}(h) = 0$, valid outside the compact support of the source $S$. It reads

$$ \mathcal{M}(h)(\mathbf{x}, t) = -\frac{1}{4\pi} \sum_{\ell=0}^{+\infty} \frac{(-1)^{\ell}}{\ell!} \partial_L \left[ \frac{\tilde{k}}{r^{\ell-2}} \int_1^{+\infty} dy \gamma_{\frac{1}{2}d}(y) \mathcal{F}_L(t - yr/c) \right], \quad (281) $$

where $\gamma_s(z)$ is part of the Green’s function (265)–(267), and the multipole-moment functions $\mathcal{F}_L$ are STF (in $d$-dimensions) and related to the compact-support source term by

$$ \mathcal{F}_L(u) = \int d^d\mathbf{x} \hat{x}_L S_{[\ell]}(\mathbf{x}, u). \qquad (282) $$

We denote $\hat{x}_L$ the STF product of spatial vectors in $d$ dimensions, and employ the convenient $\ell$-dependent weighted time average

<footnote n="61">⁶¹ The situation in harmonic coordinates is to be contrasted with the calculation in ADM-type coordinates, where it was shown that the pole parts directly cancel out in the total Hamiltonian at 3PN order (Damour et al. 2001a).</footnote>
$$S_{[\ell]}(\mathbf{x}, u) \equiv \int_{-1}^{1} \mathrm{d}z \delta_{\ell}^{(\varepsilon)}(z) S(\mathbf{x}, u + zr/c), \qquad (283a)$$

$$\text{where} \quad \delta_{\ell}^{(\varepsilon)}(z) = \frac{\Gamma(\ell + \frac{3}{2} + \frac{\varepsilon}{2})}{\Gamma(\frac{1}{2})\Gamma(\ell + 1 + \frac{\varepsilon}{2})} (1 - z^2)^{\ell + \frac{3}{2}} \qquad (283b)$$

satisfies the normalization condition $\int_{-1}^{1} \mathrm{d}z \delta_{\ell}^{(\varepsilon)}(z) = 1$, and is a simple $d$-dimensional generalization of Eq. (129). The formal PN expansion of Eq. (283), generalizing the 3 dimensional result (137), reads

$$S_{[\ell]}(\mathbf{x}, u) = \sum_{k=0}^{+\infty} \frac{1}{2^{2k} k!} \frac{\Gamma(\frac{d}{2} + \ell)}{\Gamma(\frac{d}{2} + \ell + k)} \left(\frac{r}{c} \frac{\partial}{\partial t}\right)^{2k} S(\mathbf{x}, u). \qquad (284)$$

In the second step we want to generalize the formalism of Sect. 2.2.3, i.e. solving for the multipole expansion of a source term with non-compact spatial support, given by the stress-energy matter plus gravitation pseudo tensor $\tau^{\alpha\beta}$ of the Einstein field equations in harmonic coordinates (261). We obtain

$$\mathcal{M}(h^{\alpha\beta}) = u^{\alpha\beta} - \frac{4G}{c^4} \sum_{\ell=0}^{+\infty} \frac{(-)^{\ell}}{\ell!} \partial_L \left[ \frac{\tilde{k}}{r^{d-2}} \int_1^{+\infty} dy \gamma_{\frac{1}{2}d}(y) \mathcal{F}_L^{\alpha\beta}(t - yr/c) \right]. \qquad (285)$$

The first term $u^{\alpha\beta}$ takes the same form as the corresponding piece in 3 dimensions, given by Eq. (127a). In the second term the multipole-moment function $\mathcal{F}_L^{\alpha\beta}(u)$ is given by the integral extending over the PN expansion of the pseudo-tensor $\tilde{\tau}^{\alpha\beta}$ (Blanchet et al. 2005a; Henry et al. 2021):

$$\mathcal{F}_L^{\alpha\beta}(u) = \mathrm{FP}_{B=0} \int d^d\mathbf{x} \tilde{r}^B \hat{x}_L \tilde{\tau}_{[\ell]}^{\alpha\beta}(\mathbf{x}, u). \qquad (286)$$

Indeed, we can follow exactly the same procedure as in 3 dimensions (see Sect. 2.2.3), which is appropriate for the non-compact support source term $\tilde{\tau}^{\alpha\beta}$. In particular, we introduced the regularization factor $\tilde{r}^B$ [see Eq. (72)] and the finite part when $B \to 0$ in both terms of Eq. (285). One may wonder why it is necessary to introduce the FP when $B \to 0$ since we have at our disposal the dimensional regularization parameter $d = 3 + \varepsilon$ which should play the role of IR regularization. Actually, applying the FP process “on the top” of dimensional regularization changes the nature of the limit $B \to 0$. While this limit is singular in 3 dimensions, it turns out to be finite in any generic $d$ dimensions. We shall see an example in Sect. 3.1.4 where it is necessary to compute first the finite limit when $B \to 0$, and only then to apply the standard dimensional regularization, of course keeping consistently all poles $\propto 1/\varepsilon$. We call this variant of dimensional regularization the “$B\varepsilon$” regularization.

With the result (286) in hand, we can derive the expressions of the irreducible source moments in $d$ dimensions, generalizing those computed in Theorem 7. Extending the definitions (136) we pose
$$ \bar{\Sigma} = \frac{2}{d-1} \frac{(d-2)\bar{\tau}^{00} + \bar{\tau}^{ii}}{c^2}, \quad \bar{\Sigma}^i = \frac{\bar{\tau}^{i0}}{c}, \quad \bar{\Sigma}^{ij} = \bar{\tau}^{ij}, \qquad (287) $$

where we recall that $\bar{\tau}^{\alpha\beta}$ is the PN expansion of the pseudo-tensor $\tau^{\alpha\beta}$, and obtain the mass-type moment as (Blanchet et al. 2005a)

$$ I_L = \frac{d-1}{2(d-2)} \mathrm{FP}_{B=0} \int d^d\mathbf{x} \tilde{r}^B \left\{ \begin{aligned} & \hat{x}^L \bar{\Sigma}_{[\ell]} - \frac{4(d+2\ell-2)}{c^2(d+\ell-2)(d+2\ell)} \hat{x}^{iL} \bar{\Sigma}^{\ell i}_{[\ell+1]} \\ & + \frac{2(d+2\ell-2)}{c^4(d+\ell-1)(d+\ell-2)(d+2\ell+2)} \hat{x}^{ijL} \bar{\Sigma}^{\ell ij}_{[\ell+2]} \\ & - \frac{4(d-3)(d+2\ell-2)}{c^2(d-1)(d+\ell-2)(d+2\ell)} B \hat{x}^{iL} \frac{\hat{x}^j}{r^2} \bar{\Sigma}^{\ell ij}_{[\ell+1]} \end{aligned} \right\}. \qquad (288) $$

At Newtonian order, this expression reduces to the standard result $I_L = \int d^d\mathbf{x} \rho \hat{x}_L + \mathcal{O}(1/c^2)$ with $\rho = T^{00}/c^2$ denoting the usual Newtonian density. The post-Newtonian corrections can be computed explicitly using Eq. (284). The last term in (288) does not exist in 3 dimensions and plays no role in practical calculations of the mass moment up to 4PN order.

The mass moment is just STF in all its indices so that the symmetry is given by a symmetric Young tableau (with the multi-index $L = i_1 \cdots i_\ell$)

$$ I_L = \boxed{\begin{smallmatrix} i_\ell & \cdots & i_1 \end{smallmatrix}}. \qquad (289) $$

In $d$ dimensions, since we do not dispose of the usual notion of Levi-Civita symbol, it is more tricky to define the proper generalization of the current-type moment. One has (Henry et al. 2021)

$$ J_{i|L} = \mathrm{FP}_{B=0} \mathrm{Sym}_{i|L} \int d^d\mathbf{x} \tilde{r}^B \left\{ -2\hat{x}^L \bar{\Sigma}^i_{[\ell]} + \frac{2(2\ell+d-2)}{c^2(\ell+d-1)(2\ell+d)} \hat{x}^{aL} \bar{\Sigma}^{ia}_{[\ell+1]} \right\}, \qquad (290) $$

where we introduce the specific notation that allows reconstructing the symmetries of the object $J_{i|L} \equiv J_{i|i_1 \cdots i_\ell}$:

$$ \mathrm{Sym}_{i|L} \equiv A_{ii_\ell} \mathrm{TF}_{iL} \mathrm{STF}_L, \qquad (291) $$

where STF is the standard symmetric-trace-free projection with respect to the indices $L$, where TF means removing the traces (but not symmetrizing) of the resulting object with respect to the indices $iL$, and $A_{ii_\ell}$ means the anti-symmetrization with respect to the pair of indices $ii_\ell$ (with the factor $1/2$ included). Thus, $J_{i|L}$ is trace-free, STF with respect to $L-1 = i_1 \cdots i_{\ell-1}$ and anti-symmetric with respect to the pair $i_\ell i$. Finally, besides $I_L$ and $J_{i|L}$, there exists in $d$ dimensions an additional type of irreducible multipole moment, denoted $K_{ij|L} \equiv K_{ij|i_1 \cdots i_\ell}$ and given by Eq. (2.34) of Henry et al. (2021). The symmetries of the moments $J_{i|L}$ and $K_{ij|L}$ are given by the mixed Young tableaux
$$ J_{i|L} = \begin{bmatrix} i_\ell & i_{\ell-1} & \cdots & i_1 \\ i \end{bmatrix}, \qquad K_{ij|L} = \begin{bmatrix} i_\ell & i_{\ell-1} & i_{\ell-2} & \cdots & i_1 \\ j & i \end{bmatrix}, \qquad (292) $$

with the convention that the indices are symmetrized over lines before being antisymmetrized over columns. However, the tensor $K_{ij|L}$ actually vanishes in 3 dimensions. This can be checked by counting the number of its independent components, which turns out to be proportional to $d-3$. In fact, the tensor $K_{ij|L}$ for $\ell=2$ has exactly the same trace-free property, symmetries, and number of independent components as the Weyl tensor. It plays no role with dimensional regularization, which is always performed in the limit $d \to 3$, and shall henceforth be ignored.

The ordinary STF mass-type moment in 3 dimensions is simply recovered as the limit $I_L^{(d=3)} = \lim_{d \to 3} I_L$ (after the UV type renormalization). In particular, we observe that this limit removes the last term in Eq. (288). For the ordinary STF current-type moment in 3 dimensions we have

$$ \lim_{d \to 3} J_{i|L} = \epsilon_{ii'a} J_{aL-1}^{(d=3)} \iff J_L^{(d=3)} = \frac{1}{2} \epsilon_{ab}(i_\ell) \lim_{d \to 3} J_{a|bL-1}, \qquad (293) $$

where $\epsilon_{abi}$ is the usual Levi-Civita symbol in 3 dimensions and underlined indices are excluded from symmetrization denoted by parenthesis. Notice that $J_L^{(d=3)}$, as recovered from (293), not only is symmetric in its indices $L$ but is also automatically trace-free.

As in the problem of equations of motion, the strategy—since the moments (288) and (290) cannot be computed in closed form for any $d$—is to use the expressions of the $d$-dimensional moments to evaluate the difference between the dimensional-regularization result and the previous Hadamard (or pHS) one (Blanchet et al. 2004b, 2005a). The ambiguities arise from the integration regions near the particles, that give rise to poles $\propto 1/\epsilon$, corresponding to logarithmic ultra-violet (UV) divergences in 3 dimensions.<fnref n="62" /> Since the difference is evaluated neglecting $O(\epsilon)$, the compact-support terms in the moments (proportional to the matter source densities $\sigma$, $\sigma_i$ and $\sigma_{ij}$) do not contribute. We are therefore left with evaluating the difference linked with the computation of the non-compact terms in the moments (288) and (290) near the singularities that produce poles in $d$ dimensions.

A generic non-compact support term in 3 dimensions will be given by a function $F(\mathbf{x})$ smooth except at the source points $\mathbf{y}_1$ and $\mathbf{y}_2$ around which it admits a power-like singular expansion of the type (254). The generic contribution in the multipole moments in 3 dimensions is defined by the Hadamard partie finie prescription, given either by (257) or (258), say

$$ H = \mathrm{Pf}_{s_1,s_2} \int d^3\mathbf{x} F(\mathbf{x}), \qquad (294) $$

which is just a function of time. In $d$ dimensions the corresponding function $F^{(d)}(\mathbf{x})$ admits a generalized singular expansion of the type (271); since $\lim_{d \to 3} F^{(d)}(\mathbf{x}) =$

<footnote n="62">⁶² For the moment we neglect IR divergences, which will be discussed in Sect. 3.1.4.</footnote>
$F(\mathbf{x})$ the singular coefficients in the expansion satisfy the constraint (272). In practice, the various coefficients $\int_{p,q} f_{p,q}^{(\varepsilon)}$ in the expansion are computed by specializing the general expressions of the non-linear retarded potentials $V$, $\hat{W}_{ij}$, etc. [see Eqs. (270)], to point particles in $d$ dimensions. Now the analogue of the term (294) in $d$ dimensions reads

$$H^{(d)} = \int d^d \mathbf{x} F^{(d)}(\mathbf{x}). \qquad (295)$$

Given the results of the two regularizations we compute the difference—which will have to be added to the 3-dimensional result—namely

$$\mathcal{D}H \equiv H^{(d)} - H, \qquad (296)$$

and as in Eq. (276) we compute this difference keeping the pole part $\propto \varepsilon^{-1}$ and the finite part, but neglecting terms of order $\mathcal{O}(\varepsilon)$. To this order the difference depends only on the UV behaviour of the integrands, and can therefore be computed “locally”, i.e., in the vicinity of the particles, when $r_1 \to 0$ and $r_2 \to 0$. The outcome is (Blanchet et al. 2005a)

$$\mathcal{D}H = \frac{1}{\varepsilon} \sum_{q_0 \le q \le q_1} \left[ \frac{1}{q+1} + \varepsilon \ln s_1 \right] \langle f_{1-3,q}^{(\varepsilon)} \rangle + 1 \leftrightarrow 2 + \mathcal{O}(\varepsilon), \qquad (297)$$

where the spherical average in $d$ dimensions is defined by (277). This result depends on the two Hadamard regularization scales $s_1$ and $s_2$, and on the scale $\ell_0$ belonging to dimensional regularization and defined by (262); taking into account the dependence of $\int_{p,q} f_{p,q}^{(\varepsilon)}$ upon $\ell_0$, the Hadamard scales enter (297) only through the combinations $\ln(s_1/\ell_0)$ and $\ln(s_2/\ell_0)$.

Having controlled the difference (297), the dimensional regularization of the multipole moments is obtained as the sum of the pHS part and of that difference. Consider the case of the mass quadrupole moment $I_{ij}$ since this is that moments which requests the highest approximation. Denoting $\mathcal{D}I_{ij} \sim \sum \mathcal{D}H$ the sum of all the differences we obtain

$$I_{ij}^{\mathrm{DR}} = I_{ij}^{\mathrm{pHS}} + \mathcal{D}I_{ij}, \qquad (298)$$

which now no longer depends on the Hadamard scales $s_1$ and $s_2$; however it contains poles $\propto 1/\varepsilon$ coming from the differences (297). At this stage all elements are in place to prove (Blanchet et al. 2004b, 2005a):

**Theorem 12** The pole part $\propto 1/\varepsilon$ of the dimensional-regularized quadrupole moment (298) can be renormalized into the same shifts of the world-lines as for the particles accelerations [Eq.(280) in Theorem 11]:

$$I_{ij} = \lim_{\varepsilon \to 0} \left[ I_{ij}^{\mathrm{DR}} + \delta_{\xi} I_{ij} \right], \qquad (299)$$

with $\delta_{\xi} I_{ij}$ denoting the contribution of the shifts.
To 3PN order it has been checked that the finite renormalized quadrupole moment (299) is physically equivalent to the Hadamard-regularized one (Blanchet et al. 2002b; Blanchet and Iyer 2005) if and only if the ambiguity parameters in the Hadamard calculation take uniquely determined values. Of course the procedure is then extended to 4PN order.

### 3.1.4 Infra-red divergences and radiative modes

Besides the UV divergences described in the previous Sects. 3.1.2 and 3.1.3, and which arises at 3PN order, there are also IR divergences, due to the behaviour of spatial integrals at infinity, which start to appear at 4PN order (Damour et al. 2014; Galley et al. 2016; Bernard et al. 2017a; Foffa and Sturani 2012). The IR divergences are present in both the multipole moments of the radiation field and in the local equations of motion (more precisely the Fokker Lagrangian introduced in Sect. 3.2.1).

Consider a non-compact support term $F(\mathbf{x})$ which typically contributes to the binary's multipole moments or the Lagrangian, and is obtained by direct post-Newtonian iteration of the metric in 3 dimensions. The divergences occur from the expansion of $F$ when $r \to +\infty$, which is of the form (for any $\mathcal{N} \in \mathbb{N}$)

$$F(\mathbf{x}) = \sum_{p_0 \le p \le \mathcal{N}} \frac{1}{r^p} \varphi_p(\mathbf{n}) + o\left(\frac{1}{r^\mathcal{N}}\right). \quad (300)$$

The coefficients $\varphi_p$ depend on the unit direction $\mathbf{n} = \mathbf{x}/r$ and on $p$, with minimal value $p_0 \in \mathbb{Z}$ corresponding to growing divergent behaviour with the distance. In $d$ dimensions the same PN iteration yields $F^{(d)}(\mathbf{x})$, which admits a more general expansion of the form

$$F^{(d)}(\mathbf{x}) = \sum_{\substack{p_0 \le p \le \mathcal{N} \\ q_0 \le q \le q_1}} \frac{1}{r^{p+q\varepsilon}} \varphi_{p,q}^{(\varepsilon)}(\mathbf{n}) + o\left(\frac{1}{r^\mathcal{N}}\right). \quad (301)$$

Comparing (300) and (301) in the limit $\varepsilon \to 0$ we have

$$\sum_{q=q_0}^{q_1} \varphi_{p,q}^{(\varepsilon=0)}(\mathbf{n}) = \varphi_p(\mathbf{n}). \quad (302)$$

With such behaviour in hand we want to regularize the three-dimensional divergent integral $K^{(\text{div})} = \int d^3\mathbf{x} F(\mathbf{x})$. Since we are dealing with the IR bound at infinity, we can consider the integration domain $r = |\mathbf{x}| > \mathcal{R}$, where $\mathcal{R}$ is some sufficiently large radius, thus $K^{(\text{div})} = \int_{r > \mathcal{R}} d^3\mathbf{x} F(\mathbf{x})$.

The Hadamard regularization is a good starting point, and we employ the partie finie integral, in the form based on analytic continuation with regularization factor $\tilde{r}^B$, like in Eqs. (7). Thus we define
$$K_R = \underset{B=0}{\mathrm{FP}} \int_{r > R} \mathrm{d}^3 \mathbf{x} \, \tilde{r}^B F(\mathbf{x}). \qquad (303)$$

A straightforward calculation, plugging (300) into (303), yields the Hadamard-regularized version of the integral as

$$K_R = - \sum_{p \neq 3} \frac{\mathcal{R}^{3-p}}{3-p} \langle \varphi_p \rangle_3 - \ln \left( \frac{\mathcal{R}}{r_0} \right) \langle \varphi_3 \rangle_3, \qquad (304)$$

where $\langle \varphi_p \rangle_3 \equiv \int d\Omega \, \varphi_p(\mathbf{n})$ is the angular average in 3 dimensions. As we see the crucial coefficient in the expansion (300) is that for $p = 3$; it generates a logarithmic divergence of the original integral. On the other hand, the dimensional regularization prescription, to be considered in the sense of analytic continuation in $d \in \mathbb{C}$, reads

$$K_R^{(d)} = \int_{r > \mathcal{R}} \mathrm{d}^d \mathbf{x} F^{(d)}(\mathbf{x}). \qquad (305)$$

Inserting the expansion (301), keeping only the pole $\propto 1/\epsilon$ followed by the finite part $\propto \epsilon^0$, and using also Eq. (302), we obtain

$$K_R^{(d)} = - \sum_{p \neq 3} \frac{\mathcal{R}^{3-p}}{3-p} \langle \varphi_p \rangle_3 + \sum_{q_0 \le q \le q_1} \left[ \frac{1}{(q-1)\epsilon} - \ln \mathcal{R} \right] \langle \varphi_{3,q}^{(\epsilon)} \rangle + \mathcal{O}(\epsilon), \qquad (306)$$

where the angular integration in the second term, because of the presence of the pole, is to be performed in $d$ dimensions following (277). A priori the result (306) could contain terms that diverge at infinity. These terms correspond to the value $q = 1$ but do not appear up to the 4PN level. One can check that the result (306) comes out the same with the $B\epsilon$ regularization considered below.

Finally, exactly like in the UV case (Sects. 3.1.2 and 3.1.3), we add the difference between the two prescriptions, i.e.

$$\mathcal{D}K \equiv K_R^{(d)} - K_R. \qquad (307)$$

The first term in (304) and (306) is the same and thus cancels out in the difference, so we obtain

$$\mathcal{D}K = \sum_{q_0 \le q \le q_1} \left[ \frac{1}{(q-1)\epsilon} + \ln r_0 \right] \langle \varphi_{3,q}^{(\epsilon)} \rangle + \mathcal{O}(\epsilon). \qquad (308)$$

As expected, the scale $\mathcal{R}$ has disappeared from the difference. Keeping track of the dependence over $\ell_0$ of the singular coefficients in (301), the remaining logarithm $\ln r_0$ is actually $\ln(r_0/\ell_0)$. The poles $\propto 1/\epsilon$ appearing in (308) can be qualified as “IR
poles”; they appear in both the expression of the (source) multipole moments and in the Fokker Lagrangian for the conservative equations of motion.<fnref n="63" />

We now review that the poles in the “difference” (308) are dropped when we add the “radiative modes” associated with tail (and alike non-linear effects) propagation in the far zone. To see this we start again with the generic d’Alembertian equation in d dimensions, Eq. (264), which we now write as $\Box h = N$ to emphasize that the source $N(\mathbf{x}, t)$ will be a non-linear term in the vacuum Einstein field equations—say a typical term in Eq. (58a). Furthermore since we are in vacuum outside the source, this term will be in the form of a multipole expansion, for instance a quadratic product of two linear homogeneous retarded waves (281). The retarded solution is given by (268) which we recopy here:

$$h = -\frac{\tilde{k}}{4\pi} \int_{1}^{+\infty} dz \gamma_{\frac{1}{2}d}(z) \int d^d \mathbf{x}' \frac{N(\mathbf{x}', t - z|\mathbf{x} - \mathbf{x}'|/c)}{|\mathbf{x} - \mathbf{x}'|^{d-2}}. \quad (309)$$

We split the spatial integration over the source point $\mathbf{x}'$ in (309) into an inner domain corresponding to $r' < r$ (where $r' \equiv |\mathbf{x}'|$ and $r \equiv |\mathbf{x}|$), and the outer domain for which $r' > r$. For the computation of the radiation modes in the conservative Lagrangian and equations of motion, the relevant integration to be considered is that of the outer domain $r' > r$, i.e.

$$h^> = -\frac{\tilde{k}}{4\pi} \int_{1}^{+\infty} dz \gamma_{\frac{1}{2}d}(z) \int_{r' > r} d^d \mathbf{x}' \frac{N(\mathbf{x}', t - z|\mathbf{x} - \mathbf{x}'|/c)}{|\mathbf{x} - \mathbf{x}'|^{d-2}}. \quad (310)$$

Within the integration domain $r' > r$ we are entitled to replace the integrand by its formal expansion when $r \to 0$. Computing the Taylor expansion when $r \to 0$, and transforming it to STF form, we obtain

$$\frac{N(\mathbf{x}', t - z|\mathbf{x} - \mathbf{x}'|/c)}{|\mathbf{x} - \mathbf{x}'|^{d-2}} = \sum_{p=0}^{+\infty} \frac{(-)^p}{p!} \sum_{j=0}^{+\infty} \Delta^{-j} \hat{x}_p \hat{\partial}_p' \Delta^j \left[ \frac{N(\mathbf{y}, t - zr'/c)}{r'^{d-2}} \right]_{y=\mathbf{x}'} \quad (311)$$

where $\hat{\partial}_p'$ is the STF product of $p$ partial derivatives, $\hat{\partial}_1' = \hat{\partial}_1' \cdots \hat{\partial}_p'$ with $\hat{\partial}_i' = \partial/\partial x'^i$, and on the last factor $\mathbf{y} = \mathbf{x}'$ is set after differentiation. We also employ the useful short-hand notation

$$\Delta^{-j} \hat{x}_p = \frac{\Gamma(p + \frac{d}{2})}{\Gamma(p + j + \frac{d}{2})} \frac{r^{2j} \hat{x}_p}{2^{2j} j!} \quad (312)$$

such a notation being motivated by the fact that $\Delta(\Delta^{-j}\hat{x}_p) = \Delta^{-j+1}\hat{x}_p$ in d dimensions. Inserting (311) into (310) we obtain

<footnote n="63">⁶³ Unlike in the EFT approach (Porto 2016; Porto and Rothstein 2017) one does not distinguish between IR poles, say ~1/εIR, and UV poles, ~1/εUV; here ε is simply related to the space dimension by ε = d - 3.</footnote>
$$h^> = -\frac{\tilde{k}}{4\pi} \sum_{p=0}^{+\infty} \frac{(-)^p}{p!} \sum_{j=0}^{+\infty} \frac{\Delta^{-j}\hat{x}_P}{c^{2j}} \int_1^{+\infty} dz \gamma_{\frac{1-d}{2}}(z) \\ \times \int_{r'>r} d^d\mathbf{x}' \hat{c}_P' \left[ \frac{N^{(2j)}(\mathbf{y}, t - zr'/c)}{r'^{d-2}} \right]_{\mathbf{y}=\mathbf{x}'} . \qquad (313)$$

Note that the iterated Laplace operator $\Delta^{ij}$ in (311) has transformed into $2j$ time derivatives of the source as indicated with the superscript (2j).

Next we consider the above general formula in the case where the source term has a definite multipolarity $\ell$, i.e. is of the type $N(\mathbf{x}, t) = \hat{n}_L N(r, t)$. Clearly, after summing up all the multipolar pieces, while assuming the convergence of the multipolar series, there is no restriction on the generality of the solution. We shall denote the corresponding solution by $h_L^>$ since it depends on the multi-index $L = i_1 \cdots i_\ell$. In that case we can explicitly perform the angular integration in $d$ dimensions [see Eq. (B23b) of Blanchet et al. 2004a] and obtain, after simplification by means of a series of integrations by parts over the variable $z$,

$$h_L^> = \frac{-1}{d + 2\ell - 2} \sum_{j=0}^{+\infty} \frac{\Delta^{-j} \hat{x}_L}{c^{2j}} \int_1^{+\infty} dz \gamma_{\frac{1-d}{2}-\ell}(z) \\ \times \int_r^{+\infty} dr' r'^{-\ell+1} N^{(2j)}(r', t - \frac{zr'}{c}). \qquad (314)$$

Note that the index of the $\gamma$-function is $\frac{1-d}{2} - \ell$ instead of $\frac{1-d}{2}$ previously. Now it was shown by Bernard et al. (2017a) that the radiative modes needed for insertion into the Fokker Lagrangian describing the conservative dynamics, are given by the outer domain integral (314), but in which the radial integral from $r$ to $+\infty$ is replaced by the integral from 0 to $+\infty$. Finally the general form of the radiative modes is

$$h_L^{\text{rad}}(\mathbf{x}, t) = \sum_{j=0}^{+\infty} \frac{\Delta^{-j} \hat{x}_L}{c^{2j}} \mathcal{N}^{(2j)}(t), \qquad (315)$$

where the function $\mathcal{N}(t)$ is related to the source term of the (vacuum) Einstein field equations by

$$\mathcal{N}(t) = \frac{-1}{d + 2\ell - 2} \int_1^{+\infty} dz \gamma_{\frac{1-d}{2}-\ell}(z) \frac{\text{FP}}{B=0} \\ \times \int_0^{+\infty} dr' \tilde{r}^{\prime B} r'^{-\ell+1} N(r', t - \frac{zr'}{c}). \qquad (316)$$

From the above form one can directly check that it is an homogeneous solution in $d$ dimensions, i.e. $\Box h_L^{\text{rad}} = 0$. Furthermore, that solution is manifestly regular when $r \to 0$, so it must be identified with an homogeneous anti-symmetric solution of the wave equation in $d$ dimensions, of the type half-retarded minus half-advanced. As such this homogeneous solution can be directly imported into the near zone—via the matching equation, see Eq. (115)—where it contributes to the conservative dynamics described by the Fokker Lagrangian (or the EFT action).
Important in this approach, the “Bε” regularization scheme (Bernard et al. 2017a; Marchand et al. 2018; Larrouturou et al. 2022b) is used in Eq. (316). Thus the regularization factor $\tilde{r}^{\prime B} = (r'/r_0)^B$ is inserted first in order to treat the divergence of the multipole expansion, see Eq. (75). As already mentioned after Eq. (286), this procedure is used even in d dimensions, the reason being that for the radiative modes obtained in Eq. (316), the direct limit $\epsilon \to 0$ is ill defined and must be “protected” by the finite part regularization when $B \to 0$. Following the Bε regularization one considers first the finite part at $B = 0$ in a generic d dimension, and, second, the limit $\epsilon \to 0$ is applied, keeping all the occuring poles as usual. With this method the limit $B \to 0$ is proved to be finite (no poles $\propto 1/B$), so the method actually reduces to the usual dimensional regularization.

We apply the above formulas to a quadratic tail interaction, made of the interaction between the constant mass M and a varying multipole moment $M_L$ or $S_L$. In this case the generic form of the source term is

$$N(r, t) = r^{-k-2\epsilon} \int_{1}^{+\infty} dy y^p \gamma_{-1-\frac{\epsilon}{2}}(y) F(t - yr), \quad (317)$$

where $k, p \in \mathbb{N}$ and the function $F(t)$ stands for some time derivative of either $M_L(t)$ or $S_L(t)$ (and multiplied by M). The non-locality of the source term (317) comes directly from that of the multipolar homogeneous solution (281) in d dimensions. In the quadratic case (317), after a change of integration variable and repeated integrations by parts (Bernard et al. 2017a), the radiative mode (316) becomes

$$\mathcal{N}(t) = \mathrm{FP}_{B=0} \left\{ D_{\ell}^{p,k}(B, \epsilon) \int_{0}^{+\infty} d\tau \, \tilde{\tau}^{B} \tau^{-2\epsilon} F^{(\ell+k-1)}(t-\tau) \right\}, \quad (318)$$

where $\tilde{\tau}^B = (c\tau/r_0)^B$ and one has succeeded in factorizing out a pure numerical coefficient in front given by

$$D_{\ell}^{p,k}(B, \epsilon) = \frac{(-1)^{\ell+k}}{2\ell+1+\epsilon} \frac{\Gamma(2\epsilon-B)}{\Gamma(\ell+k-1+2\epsilon-B)} \\ \times \int_{1}^{+\infty} dy \, y^p \gamma_{-1-\frac{\epsilon}{2}}(y) \int_{1}^{+\infty} dz \, (y+z)^{\ell+k-2+2\epsilon-B} \gamma_{-\ell-1-\frac{\epsilon}{2}}(z). \quad (319)$$

We are ultimately interested in the limit $\epsilon \to 0$, and the integral in (318) is well defined in this limit, i.e. convergent at the bound $\tau = 0$. Since F is a time derivative of a multipole moment, we can assume it to be zero in the remote past so there is no problem with the bound $\tau = +\infty$ of the integral. However the coefficient (319) in front of the integral will develop a pole when $\epsilon \to 0$. In contrast with the IR pole in Eq. (308), such a pole will be a UV pole. The coefficient (319) has been computed in the limit $\epsilon \to 0$ in the Appendix D of Bernard et al. (2017a); as it turns out it is ill-defined for some values of $\ell, p$ or $k$, thus necessitating a “protection” by the B-process which reduces in this case to a mere finite limit $B \to 0$. In conclusion, the UV type pole generated by the coefficient (319) cancels out (without need of a further renormalization) the IR pole coming from the potential modes in Eq. (308).
## 3.2 Newtonian-like equations of motion

### 3.2.1 The Fokker action for point particles

The main methods able to derive the equations of motion of compact binary systems to high post-Newtonian order, namely ADM Hamiltonian, Fokker Lagrangian and EFT, have been reviewed in Sect. 1.2.1. Here we focus on the powerful Fokker Lagrangian method. The starting point is the action $S = S_g + S_m$ given by (46), but in $D = d + 1$ space-time dimensions, and where the gravitational action $S_g$ is written in the Landau and Lifshitz (1971) form and we add the standard gauge-fixing term appropriate to harmonic coordinates,<fnref n="64" />

$$S_g = \frac{c^3}{16\pi G^{(d)}} \int d^{d+1}x \sqrt{-g} \left[ g^{\alpha\beta} \left( \Gamma^{\mu}_{\alpha\nu} \Gamma^{\nu}_{\beta\mu} - \Gamma^{\mu}_{\alpha\beta} \Gamma^{\nu}_{\mu\nu} \right) - \frac{1}{2} g_{\alpha\beta} \Gamma^{\alpha} \Gamma^{\beta} \right]. \quad (320)$$

The matter action $S_m$ takes the form appropriate for $N$ point-particles—we neglect the spins and the body’s internal structure:

$$S_m = - \sum_{a=1}^{N} m_a c \int dt \sqrt{-(g_{\alpha\beta})_a v_a^\alpha v_a^\beta}. \quad (321)$$

Here $m_a$ is the mass of the particles, $v_a^\alpha = dy_a^\alpha/dt = (c, v_a)$ is the coordinate velocity, $y_a^\alpha = (ct, y_a)$ the ordinary trajectory, and $(g_{\alpha\beta})_a$ is the metric evaluated at the location of the particle a. A self-field regularization will be necessary to define the object $(g_{\alpha\beta})_a$ and we shall follow the dimensional regularization scheme. A closed-form expression for the gauge-fixed gravitational action can be written with the help of the gothic metric $g^{\alpha\beta} \equiv \sqrt{-g}g^{\alpha\beta}$ and its inverse $g_{\alpha\beta} = g_{\beta\alpha}/\sqrt{-g}$ as

$$S_g = \frac{c^3}{32\pi G^{(d)}} \int d^{d+1}x \left[ -\frac{1}{2} (g_{\alpha\mu} g_{\beta\nu} - \frac{1}{d-1} g_{\alpha\beta} g_{\mu\nu}) g^{\rho\sigma} \partial_{\rho} g^{\alpha\beta} \partial_{\sigma} g^{\mu\nu} \right. \\ \left. + g_{\alpha\beta} (\partial_{\mu} g^{\alpha\nu} \partial_{\nu} g^{\beta\mu} - \partial_{\mu} g^{\alpha\mu} \partial_{\nu} g^{\beta\nu}) \right]. \quad (322)$$

The relaxed Einstein field equations following from the variation of the gauge-fixed action are of the form

$$\Box h^{\alpha\beta} = \frac{16\pi G^{(d)}}{c^4} |g| T^{\alpha\beta} + \Sigma^{\alpha\beta}, \quad (323)$$

where $h^{\alpha\beta} = g^{\alpha\beta} - \eta^{\alpha\beta}$ denotes the gothic metric deviation, the particle stress-energy tensor takes the standard form

<footnote n="64">⁶⁴ Where $\Gamma^{\alpha}_{\mu\nu} = \frac{1}{2}g^{\alpha\rho}(\partial_{\mu}g_{\nu\rho} + \partial_{\nu}g_{\mu\rho} - \partial_{\rho}g_{\mu\nu})$ is the usual Christoffel symbol and we define $\Gamma^{\alpha} \equiv g^{\mu\nu}\Gamma^{\alpha}_{\mu\nu}$.</footnote>
$$T^{\alpha\beta} = \sum_{a=1}^{N} \frac{m_a v_a^\alpha v_a^\beta}{\sqrt{-(g_{\mu\nu})_a \frac{v_a^\mu v_a^\nu}{c^2}}} \frac{\delta^{(d)}(\mathbf{x} - \mathbf{y}_a)}{\sqrt{-g}}, \quad (324)$$

and the source term differs from the one in harmonic coordinates by terms depending on the “harmonicity” $H^\alpha \equiv \partial_\mu h^{\alpha\mu} = -\sqrt{-g} \Gamma^\alpha$:

$$\Sigma^{\alpha\beta} = \Lambda_{\text{harm}}^{\alpha\beta} - H^{\alpha} H^{\beta} - H^{\mu} \partial_{\mu} h^{\alpha\beta} \\ - \frac{1}{2} g^{\alpha\beta} g_{\mu\nu} H^{\mu} H^{\nu} + 2 g_{\mu\nu} g^{\rho(\alpha} \partial_{\rho} h^{\beta)\mu} H^{\nu}, \quad (325)$$

where $\Lambda_{\text{harm}}^{\alpha\beta}$ is given in d dimensions by Eq. (263).

The Fokker (1929) action is obtained by inserting into the action (321)–(322) an explicit PN solution of the field equations, given as a functional of the particle’s trajectories, i.e., an explicit post-Newtonian metric

$$\bar{g}_{\alpha\beta}(\mathbf{x}, t) \equiv \bar{g}_{\alpha\beta}[\mathbf{x}; \mathbf{y}_b(t), \mathbf{v}_b(t), \dots], \quad (326)$$

where the overbar recalls the post-Newtonian character of this solution. Note that the source term (136) contains all required harmonics $H^\alpha$, which are not assumed to be zero in the PN iteration of the metric (326). Hence the ellipsis in (326) indicate also the dependence on accelerations $a_b(t)$, derivative of accelerations $b_b(t)$, and so on, due to the fact that the equations of motion are “off-shell” at this stage and no replacement of accelerations is performed. Thus, the Fokker generalized action, which is a functional of the particles only, but depending not only on positions and velocities but also on accelerations and their derivatives, is given by

$$S_F[\mathbf{y}_b(t), \mathbf{v}_b(t), \dots] = \int dt \int d^d\mathbf{x} \bar{\mathcal{L}}_g[\mathbf{x}; \mathbf{y}_b(t), \mathbf{v}_b(t), \dots] \\ - \sum_a m_a c \int dt \sqrt{-\bar{g}_{\alpha\beta}[\mathbf{y}_a(t); \mathbf{y}_b(t), \mathbf{v}_b(t), \dots]} v_a^\alpha v_a^\beta, \quad (327)$$

where $\bar{\mathcal{L}}_g$ is the Lagrangian density of the gravitational action (322) computed with the metric (326). As is well known, it is always possible to eliminate from a generalized PN action a contribution that is quadratic in the accelerations by absorbing it into a “double-zero” term which does not contribute to the dynamics (Damour and Schäfer 1985). The argument can be extended to any term polynomial in the accelerations (and their derivatives). The post-Newtonian equations of motion of the particles are obtained as the generalized Lagrange equations

$$\frac{\delta S_F}{\delta \mathbf{y}_b} \equiv \frac{\partial L_F}{\partial \mathbf{y}_b} - \frac{d}{dt} \left( \frac{\partial L_F}{\partial \mathbf{v}_b} \right) + \frac{d^2}{dt^2} \left( \frac{\partial L_F}{\partial \mathbf{a}_b} \right) + \dots = 0, \quad (328)$$

where $L_F$ is the corresponding Lagrangian (such that $S_F = \int dt L_F$). Once they have been constructed, the equations (328) can be order reduced by replacing iteratively
all the higher-order accelerations by their expressions coming from the PN equations of motion themselves.<fnref n="65" />

By definition of the Fokker action, varying it with respect to the post-Newtonian metric (326) yields back the relaxed field equations. Hence (coming back to the gothic variable):

$$ \frac{\delta S_F}{\delta \bar{h}_{\alpha\beta}} = \frac{c^4}{64\pi G^{(d)}} \left[ \Box \bar{h}^{\alpha\beta} - \bar{\Sigma}^{\alpha\beta} - \frac{16\pi G}{c^4} |\bar{g}| \bar{T}^{\alpha\beta} \right] = 0. \quad (329) $$

Now let us suppose that we have constructed the post-Newtonian metric up to some small post-Newtonian remainder, say $\bar{h} = \bar{h}_n + r_n$ (skipping space-time indices for a while), where $\bar{h}_n$ is explicitly known while the small remainder $r_n$ is unknown. Evaluating the Fokker action for this solution we get<fnref n="66" />

$$ S_F[\bar{h}] = S_F[\bar{h}_n] + r_n \frac{\delta S_F}{\delta \bar{h}}[\bar{h}_n] + \mathcal{O}(r_n^2). \quad (330) $$

Actually, because the variation of the action with respect to $\bar{h}$ is zero, Eq. (328), when evaluated for the approximate solution $\bar{h}_n$, the correction term in (330) will acquire an extra small factor $r_n$. More precisely, taking also into account the factor $c^4$ in front of (329),

$$ \frac{\delta S_F}{\delta \bar{h}}[\bar{h}_n] = \mathcal{O}(c^4 r_n) \quad \implies \quad S_F[\bar{h}] = S_F[\bar{h}_n] + \mathcal{O}(c^4 r_n^2). \quad (331) $$

If the unknown remainder is of order $r_n = \mathcal{O}(c^{-n-3})$, which means that we control the post-Newtonian metric up to order $c^{-n-2}$ included, we see that the Fokker action will be controlled up to order

$$ S_F[\bar{h}] = S_F[\bar{h}_n] + \mathcal{O}(c^{-2n-2}), \quad (332) $$

which means that it is controlled up to $n$PN order included. Thus in order to compute the Fokker action to order $n$PN it suffices to insert the metric to order $c^{-n-2}$, which represents for high PN orders an important improvement with respect to the naive expectation $c^{-2n}$. However the previous reasoning neglects the tensorial nature of the metric. The more precise analysis by Bernard et al. (2016) concludes on the same result provided that the components of the metric are decomposed according to

$$ \bar{h}^{\alpha\beta} \rightarrow \begin{cases} \bar{h} = (\bar{h}^{00ii}, \bar{h}^{0i}, \bar{h}^{ij}), \\ \bar{h}^{00ii} = \frac{2}{d-1} [(d-2)\bar{h}^{00} + \bar{h}^{ii}]. \end{cases} \quad (333) $$

That is, in order to obtain the Fokker action accurate to $n$PN order, it suffices to insert (each component of) the previous $\bar{h}$ developed up to order $c^{-n-2}$ included. For instance, in order to obtain the Fokker action (and therefore the full dynamics)

<footnote n="65">⁶⁵ The classical Fokker action is equivalent, in the “tree-level” approximation, to the effective action used by the effective field theory (Foffa and Sturani 2014; Porto 2016).</footnote>

<footnote n="66">⁶⁶ Such writing is a bit sketchy; see Eq. (4.7) of Bernard et al. (2016) for more details.</footnote>
correct at 4PN order one needs only to insert $\bar{h}$ up to the relatively low level to $c^{-6}$, i. e. (since we are interested in the conservative dynamics) $\bar{h}^{00ii}$ and $\bar{h}^{ij}$ to level $c^{-6}$, and $\bar{h}^{0i}$ to level $c^{-5}$. The above technique has been referred to as the “method $n+2$”.

As said in Sect. 3.1.3 the Fokker Lagrangian contains some IR poles coming from the behaviour when $r \to +\infty$ of the gravitational part of the Lagrangian. Thus, one computes the expansion $r \to +\infty$ of the various potentials parametrizing the metric in $d$ dimensions, see e.g. (270), and then apply the formula (308). The potentials are those needed at the 4PN order (say) following the previous method “$n+2$”. In the language of EFT (Porto and Rothstein 2017) the previous calculation corresponds to “potential mode” contribution.

However, there is also to take into account the contribution coming from the conservative part of the 4PN tail effect in $d$ dimensions, which corresponds in the EFT language to the “radiation” contribution and has been reviewed in Sect. 3.1.4. The radiation mode is an imprint of tails in the local conservative dynamics of the source at the 4PN order (Blanchet and Damour 1988). It contains a tail-induced modification of the dissipative radiation reaction force at the relative 1.5PN order, see Eq. (176). The conservative piece modifies the Fokker action of equivalently the EFT at the 4PN order (Foffa and Sturani 2013; Damour et al. 2014; Galley et al. 2016).

Applying Eqs. (318)–(319) to the case of the quadratic interaction between the constant mass M and the varying quadrupole moment $M_{ij}$, one obtains at the leading 4PN order (Bernard et al. 2017a)

$$ \bar{h}^{00ii}_{\text{rad}} = \frac{8G^2M}{15c^{10}} x^i \int_0^{+\infty} d\tau \left[ \ln\left(\frac{c\sqrt{\bar{q}}\tau}{2\ell_0}\right) - \frac{1}{2\epsilon} + \frac{61}{60} \right] M_{ij}^{(7)}(t-\tau), \quad (334a) $$

$$ \bar{h}^{0i}_{\text{rad}} = -\frac{8G^2M}{3c^9} x^i \int_0^{+\infty} d\tau \left[ \ln\left(\frac{c\sqrt{\bar{q}}\tau}{2\ell_0}\right) - \frac{1}{2\epsilon} + \frac{107}{120} \right] M_{ij}^{(6)}(t-\tau), \quad (334b) $$

$$ \bar{h}_{\text{rad}}^{ij} = \frac{8G^2M}{c^8} \int_0^{+\infty} d\tau \left[ \ln\left(\frac{c\sqrt{\bar{q}}\tau}{2\ell_0}\right) - \frac{1}{2\epsilon} + \frac{4}{5} \right] M_{ij}^{(5)}(t-\tau). \quad (334c) $$

We do not indicate the 5PN remainder terms nor the neglected $O(\epsilon)$ terms; G is Newton’s constant in Eq. (262) and we pose $\bar{q} = 4\pi\epsilon\gamma_E$ with $\gamma_E$ denoting the Euler constant.<fnref n="67" /> The tail terms in (334) yield a non-local modification of the equations of motion at 4PN order, and consequently a modification of the Fokker action. One finds

$$ S_{\text{F}}^{\text{rad}} = \frac{2G^2M}{5c^8} \int_{-\infty}^{+\infty} dt M_{ij}^{(3)}(t) \int_0^{+\infty} d\tau \left[ \ln\left(\frac{c\sqrt{\bar{q}}\tau}{2\ell_0}\right) - \frac{1}{2\epsilon} + \frac{41}{60} \right] M_{ij}^{(4)}(t-\tau). \quad (335) $$

An important point to notice is that (335) can be rewritten in a manifestly time-symmetric way; thus the procedure automatically selects the “conservative” part of

<footnote n="67">⁶⁷ Euler’s constant (discovered in 1731) is sometimes called the Euler-Mascheroni constant. Mascheroni’s contribution was to have computed (in 1790) the constant up to 19 decimal places, while Euler had previously computed it to 15 places. Nowadays, one credits the important mathematical contribution, rather than the record of the most numerical digits computed, and it is appropriate to name the constant after Euler alone (Lagarias 2013).</footnote>
the tail at 4PN order—the dissipative part giving no contribution to the action. Hence we have in fact

$$S_F^{\text{rad}} = \frac{G^2 M}{5c^8} \int_{-\infty}^{+\infty} dt M_{ij}^{(3)}(t) \\ \times \int_{0}^{+\infty} d\tau \left[ \ln\left(\frac{c\sqrt{q}\tau}{2\ell_0}\right) - \frac{1}{2\epsilon} + \frac{41}{60} \right] \left( M_{ij}^{(4)}(t-\tau) - M_{ij}^{(4)}(t+\tau) \right), \quad (336)$$

which can elegantly be rewritten by means of the Hadamard partie finie (Pf) integral, defined as the one-dimensional version of Eq. (257), as (Foffa and Sturani 2013; Damour et al. 2014; Galley et al. 2016; Bernard et al. 2017a)

$$S_F^{\text{rad}} = \frac{G^2 M}{5c^8} \text{Pf} \iint \frac{dt \, dt'}{|t-t'|} M_{ij}^{(3)}(t) M_{ij}^{(3)}(t'), \quad (337)$$

with $\tau_0 = \frac{2\ell_0}{c\sqrt{q}} e^{\frac{1-41}{60}}$ playing the role of the Hadamard cut-off scale.<fnref n="68" /> Since the action (337) is non local its variation involves functional derivatives rather than ordinary derivatives. Varying with respect to the particle world-lines $y_a^i(t)$—for a system of point particles with mass $m_a$—one obtains

$$\frac{\delta S_F^{\text{rad}}}{\delta y_a^i(t)} = - \frac{4G^2 M}{5c^8} m_a y_a^i(t) \text{Pf} \int_{-\infty}^{+\infty} \frac{dt'}{|t-t'|} M_{ij}^{(6)}(t'), \quad (338)$$

which coincides with the conservative part of the known 4PN tail contribution in the equations of motion (Blanchet and Damour 1988).

### 3.2.2 The 3.5PN acceleration and 3PN energy

The acceleration of one of the particles, say the particle 1, is presented up to the 3.5PN order, together with the corresponding binary's energy (which is conserved in the absence of radiation reaction). Besides the Hamiltonian and Lagrangian methods reviewed in Sect. 1.2.1, a “direct” post-Newtonian attack of the problem is possible, reducing the metric of an extended regular source, as given in harmonic coordinates to 3.5PN order by Eqs. (161), to the case where the matter tensor is made of delta functions, and then curing the self-field divergences by means of the dimensional regularization (see Sect. 3.1.2). In this approach the equations of motion (in the case of spinless particles) are just the geodesic equations provided to 3.5PN order in Eqs. (165)–(168); the metric therein is the regularized metric generated by the system of particles themselves. The radiation reaction effects at 2.5PN and 3.5PN orders in harmonic coordinates are known from Iyer and Will (1993, 1995), Jaranowski and

$$\text{Pf}_{\tau_0} \int_{-\infty}^{+\infty} dt' \frac{f(t')}{|t-t'|} = \int_{0}^{+\infty} d\tau \ln\left(\frac{\tau}{\tau_0}\right) [f^{(1)}(t-\tau) - f^{(1)}(t+\tau)].$$

<footnote n="68">⁶⁸ For any regular function $f(t)$ tending to zero sufficiently rapidly when $t \to \pm\infty$ we have</footnote>