up at some high post-Newtonian order. Historically these difficulties, even appearing at higher approximations, have cast a doubt on the actual soundness, from a theoretical point of view, of the post-Newtonian expansion. Practically speaking, they posed the question of the reliability of the approximation, when comparing the theory's predictions with very precise experimental results. This was one of the main reason for the famous radiation-reaction controversy raging at the time of the binary pulsar data [182, 418]. In this section we assess the nature of these difficulties – are they purely technical or linked with some fundamental drawback of the approximation scheme? – and eventually resolve them.

1. The first problem we face is that in higher approximations some *divergent* Poisson-type integrals appear. Indeed the post-Newtonian expansion replaces the resolution of a hyperbolic-like d'Alembertian equation by a perturbatively equivalent hierarchy of elliptic-like Poisson equations. Rapidly it is found during the post-Newtonian iteration that the right-hand side of the Poisson equations acquires a non-compact support (it is distributed all over space $\mathbb{R}^3$), and that as a result the standard Poisson integral diverges at the bound of the integral at spatial infinity, i.e., when $r = |x| \to +\infty$, with $t = \text{const}$.

2. The second problem is related with the limitation of the post-Newtonian approximation to the near zone – the region surrounding the source of small extent with respect to the wavelength of the emitted radiation: $r \ll \lambda$. As we have seen, the post-Newtonian expansion assumes from the start that all retardations $r/c$ are small, so it can rightly be viewed as a formal *near-zone* expansion, when $r \to 0$. Note that the fact which makes the Poisson integrals to become typically divergent, namely that the coefficients of the post-Newtonian series blow up at spatial infinity, when $r \to +\infty$, has nothing to do with the actual behaviour of the field at infinity. However, the serious consequence is that it is a *priori* impossible to implement within the post-Newtonian scheme alone the physical information that the matter system is isolated from the rest of the Universe. Most importantly, the no-incoming radiation condition, imposed at past null infinity, cannot be taken directly into account, a *priori*, into the post-Newtonian scheme. In this sense the post-Newtonian approximation is not “self-supporting”, because it necessitates some information taken from outside its own domain of validity.

The divergencies are linked to the fact that the post-Newtonian expansion is actually a singular perturbation, in the sense that the coefficients of the successive powers of $1/c$ are not uniformly valid in space, since they typically blow up at spatial infinity like some powers of $r$. We know for instance that the post-Newtonian expansion cannot be “asymptotically flat” starting at the 2PN or 3PN level, depending on the adopted coordinate system [362]. The result is that the standard Poisson integrals are in general badly-behaving at infinity. Trying to solve the post-Newtonian equations by means of the Poisson integral does not make sense. However, this does not mean that there are no solutions to the problem, but simply that the Poisson integral does not constitute the appropriate solution of the Poisson equation in the context of post-Newtonian expansions.

Here we present, following Refs. [357, 75], a solution of both problems, in the form of a general expression for the near-zone gravitational field, developed to any post-Newtonian order, which has been determined from implementing the matching equation (103). This solution is free of the divergences of Poisson-type integrals we mentioned above, and yields, in particular, some general expression, valid up to any order, of the terms associated with the gravitational radiation reaction force inside the post-Newtonian source.

Though we shall focus our attention on the particular approach advocated in [357, 75], there are other ways to resolve the problems of the post-Newtonian approximation. Notably, an alternative solution to the problem of divergencies, proposed in Refs. [214, 211], is based on an initial-value formulation. In this method the problem of the appearance of divergencies is avoided because
of the finiteness of the causal region of integration, between the initial Cauchy hypersurface and the considered field point. On the other hand, a different approach to the problem of radiation reaction, which does not use a matching procedure, is to work only within a post-Minkowskian iteration scheme without expanding the retardations, see e.g., Ref. [126].

## 5.1 Post-Newtonian iteration in the near zone

We perform the post-Newtonian iteration of the field equations in harmonic coordinates in the near zone of an isolated matter distribution. We deal with a general hydrodynamical fluid, whose stress-energy tensor is smooth, i.e., $T^{\alpha\beta} \in C^{\infty}(\mathbb{R}^4)$. Thus the scheme a priori excludes the presence of singularities and black holes; these will be dealt with in Part B of this article.

We shall now prove [357] that the post-Newtonian expansion can be indefinitely iterated without divergences. Like in Eq. (106) we denote by means of an overline the formal (infinite) post-Newtonian expansion of the field inside the source's near-zone. The general structure of the post-Newtonian expansion is denoted (skipping the space-time indices $\alpha\beta$) as

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

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 3 [see Eq. (53)], 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$.33 For the stress-energy pseudo-tensor appearing in Eq. (106) we have the same type of expansion,

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

Note that the expansion starts with a term of order $c^2$ corresponding to the rest mass-energy ($\bar{\tau}$ has the dimension of an energy density). As usual we shall understand the infinite sums such as (127)–(128) in the sense of formal series, i.e., merely as an ordered collection of 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(\mathbf{x}, t) \in C^{\infty}(\mathbb{R}^4)$.

Inserting the post-Newtonian ansatz (127) into the harmonic-coordinates Einstein field equation (21)–(22) 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 = 16\pi G \bar{\tau}_{m-4} + \partial_t^2 \bar{h}, \qquad (129) $$

where the second term comes from the split of the d'Alembertian operator into a Laplacian and a second time derivative: $\square = \Delta - \frac{1}{c^2}\partial_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}_p$ ($\forall p \le m-1$), and from that show how to infer the next-order coefficient $\bar{h}_m$.

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 a specific process of analytic continuation. For any source term like $\bar{\tau}_m$, we multiply it by the regularization factor $\tilde{\tau}^B$ already extensively used in the construction of the exterior field, thus $B \in \mathbb{C}$ and $\tilde{r} = r/r_0$ is given by Eq. (42). Only then do we apply the usual Poisson integral, which therefore defines a

<footnote n="33">33 For this argument we assume the validity of the matching equation (103) and that the post-Minkowskian series over $n = 1, \dots, \infty$ in Eq. (53) has been formally summed up.</footnote>
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 one can prove [357] that 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 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 → 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:

$$ \Delta^{-1}[\bar{\tau}_m](\mathbf{x}, t) = -\frac{1}{4\pi} \mathcal{F}\mathcal{P}_{B=0} \int \frac{d^3\mathbf{x}'}{|\mathbf{x} - \mathbf{x}'|} \tilde{r}^{\prime B} \bar{\tau}_m(\mathbf{x}', t). \quad (130) $$

The integral extends over all three-dimensional space but with the latter finite-part regularization at infinity denoted $\mathcal{F}\mathcal{P}_{B=0}$ or $\mathcal{F}\mathcal{P}$ for short. The main properties of this generalized Poisson operator is that it solves the Poisson equation,

$$ \Delta(\Delta^{-1}[\bar{\tau}_m]) = \bar{\tau}_m, \quad (131) $$

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

$$ \Delta^{-k-1}[\bar{\tau}_m](\mathbf{x}, t) = -\frac{1}{4\pi} \mathcal{F}\mathcal{P}_{B=0} \int d^3\mathbf{x}' \frac{|\mathbf{x} - \mathbf{x}'|^{2k-1}}{(2k)!} \tilde{r}^{\prime B} \bar{\tau}_m(\mathbf{x}', t). \quad (132) $$

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. (129), and augmented by the most general homogeneous solution of the Poisson equation. Thus, we can write

$$ \bar{h}_m = 16\pi G \Delta^{-1}[\bar{\tau}_{m-4}] + \partial_t^2 \Delta^{-1}[\bar{h}_{m-2}] + \sum_{\ell=0}^{+\infty} \mathcal{B}_L(t) \hat{x}_L. \quad (133) $$

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}_L(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 gravitational field at infinity from the source.

It is now trivial to iterate the process. We substitute for $\bar{h}_{m-2}$ in the right-hand side of Eq. (133) 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 = 0$ or $\bar{h}_1 = 0$. At this point $\bar{h}_m$ is expressed in terms of the previous $\bar{\tau}_p$'s and $\mathcal{B}_L$'s with $p \le m-2$. To finalize the process we introduce what we call the operator of the “instantaneous” potentials and denote $\Box_{\text{inst}}^{-1}$. Our notation is chosen to contrast with the standard operator of the retarded potentials $\Box_{\text{ret}}^{-1}$ defined by Eq. (31). However, beware of the fact that unlike $\Box_{\text{inst}}^{-1}$ the operator $\Box_{\text{inst}}^{-1}$ will be defined only when acting on a post-Newtonian series such as $\bar{\tau}$. Indeed, we pose

$$ \Box_{\text{inst}}^{-1}[\bar{\tau}] = \sum_{k=0}^{+\infty} \left(\frac{\partial}{c\partial t}\right)^{2k} \Delta^{-k-1}[\bar{\tau}], \quad (134) $$

where the k-th iteration of the generalized Poisson operator is defined by Eq. (132). 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 source-free d’Alembertian equation,

$$ \Box (\Box_{\text{inst}}^{-1}[\bar{\tau}]) = \bar{\tau}. \quad (135) $$
On the other hand, the homogeneous solution in Eq. (133) 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(t)$ parametrizing those advanced-minus-retarded free waves. It is very easy to relate if necessary the post-Newtonian expansion of $A_L(t)$ to the functions $mB_L(t)$ previously introduced in Eq. (133). Finally the most general post-Newtonian solution, iterated ad infinitum and without any divergences, is obtained into the form

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

We shall refer to the $A_L(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. (136): There is no means to compute the radiation-reaction functions $A_L(t)$. We are here touching the second problem faced by the standard post-Newtonian approximation.

## 5.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., [6]. 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 [181, 269, 270, 334]. Unfortunately such a 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 the hereditary modification of the radiation reaction force at the 1.5PN relative order [58, 60, 43]. Since the reaction force is itself of order 2.5PN this means that 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 5 and 6 that the matching equation (103) 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$ which have been left unspecified in Eq. (136).

We find [357] that the radiation-reaction functions $A_L$ are composed of the multipole moment functions $F_L$ defined by Eq. (119), which will here characterize “linear-order” radiation reaction effects starting at 2.5PN order, and of an extra piece $R_L$, which will be due to non-linear effects in the radiation reaction and turn out to arise at the 4PN order. Thus,

$$ A_L(t) = F_L(t) + R_L(t). \quad (137) $$

The extra piece $R_L$ is obtained from the multipole expansion of the pseudo-tensor $M(\tau)$.34 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

$$ R_L(t) = FP \int d^3\mathbf{x} \hat{x}_L \int_1^{+\infty} dz \gamma_\ell(z) M(\tau) (\mathbf{x}, t - zr/c). \quad (138) $$

<footnote n="34">34 We mean the fully-fledge $M(\tau)$; i.e., not the formal object $M(\bar{\tau})$.</footnote>
The fact that the multipolar expansion $M(\tau)$ is the source term for the function $R_L$ is the consequence of the matching equation (103). The specific contributions due to $R_L$ in the post-Newtonian metric (136) are associated with tails of waves [58, 43]. Notice that, remarkably, the $FP$ regularization deals with the bound of the integral at $r = 0$, in contrast with Eq. (119) where 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. (119), 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:<fnref n="35" />

$$ \int_{1}^{+\infty} dz \gamma_{\ell}(z) = 1. \qquad (139) $$

The post-Newtonian metric (136) is now fully determined. However, let us now prove a more interesting alternative formulation of it, derived in Ref. [75].

**Theorem 7.** 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

$$ \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{R_L^{\alpha\beta}(t-r/c) - R_L^{\alpha\beta}(t+r/c)}{2r} \right\}. \quad (140) $$

The first term represents a particular solution of the hierarchy of post-Newtonian 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 (138).

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

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

We emphasize that (141) constitutes the definition of a (formal) post-Newtonian 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 (141) should carefully be distinguished from the “global” solution $\Box_{\text{ret}}^{-1}[\tau]$ defined by Eq. (31), with global non-expanded pseudo-tensor $\tau$.

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

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

so constitutes a good prescription for a particular solution of the wave equation – as legitimate as the solution (134). Therefore the two solutions should differ by an homogeneous solution of the

<footnote n="35">³⁵ Though the latter integral is a priori divergent, its value can be determined by invoking complex analytic continuation in $\ell \in \mathbb{C}$.</footnote>
wave equation which is necessarily of the anti-symmetric type (regular inside the source). Detailed investigations [357, 75] 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 (143) $$

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

For computations limited to the 3.5PN order (level of the 1PN correction to the radiation reaction force), the first term in Eq. (140) with the “intuitive” prescription (141) is sufficient. But because of the second term in (140) there is a fundamental breakdown of this scheme at the 4PN order where it becomes necessary to take into account non-linear radiation reaction effects associated with tails. The second term in (140) constitutes a generalization of the tail-transported radiation reaction arising at the 4PN order, i.e., 1.5PN order relative to the dominant radiation reaction order, as determined in Ref. [58]. The tail-transported radiation reaction is required by energy conservation and the presence of tails in the wave zone. The usual radiation reaction terms, up to 3.5PN order, are contained in the first term of Eq. (140), and are parametrized by the same multipole-moment functions $\mathcal{F}_L$ as the exterior multipolar field, as Eq. (143) explicitly shows. In Section 5.4 we shall give an explicit expression of the radiation reaction force showing the usual radiation reaction terms to 3.5PN order, issued from $\mathcal{F}_L$, and exhibiting the above tail-induced 4PN effect, issued from $\mathcal{R}_L$.

Finally note that the post-Newtonian solution, in either form (136) or (140), has been obtained without imposing the condition of harmonic coordinates (21) in an explicit way. We have simply matched together the post-Newtonian and multipolar expansions, satisfying the “relaxed” Einstein field equations (22) in their respective domains, and found that the matching determines uniquely the solution. An important check done in [357, 75], is therefore to verify that the harmonic coordinate condition (21) is indeed satisfied as a consequence of the conservation of the pseudo-tensor (27), so that we really grasp a solution of the full Einstein field equations.

## 5.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}$, ..., 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 basic idea in our post-Newtonian iteration scheme is to use wherever possible a “direct” integration, with the help of some formulas like $\Box_{\text{ret}}^{-1}(\partial_{\mu}V\partial^{\mu}V + V\Box V) = V^2/2$. The 3.5PN harmonic-coordinates metric reads [71]

$$ g_{00} = -1 + \frac{2}{c^2}V - \frac{2}{c^4}V^2 + \frac{8}{c^6}\left(\hat{X} + V_iV_i + \frac{V^3}{6}\right) + \frac{32}{c^8}\left(\hat{T} - \frac{1}{2}V\hat{X} + \hat{R}_iV_i - \frac{1}{2}VV_iV_i - \frac{1}{48}V^4\right) + \mathcal{O}\left(\frac{1}{c^{10}}\right), \quad (144a) $$

$$ 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^2V_i\right) + \mathcal{O}\left(\frac{1}{c^9}\right), \quad (144b) $$
$$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). \quad (144c)$$

All the potentials are generated by the matter stress-energy tensor $T^{\alpha\beta}$ through some convenient definitions recalling Eqs. (124),

$$\sigma = \frac{T^{00} + T^{ii}}{c^2}, \quad (145a)$$

$$\sigma_i = \frac{T^{0i}}{c}, \quad (145b)$$

$$\sigma_{ij} = T^{ij}. \quad (145c)$$

Starting at Newtonian and 1PN orders, $V$ and $V_i$ represent some retarded versions of the usual Newtonian and gravitomagnetic potentials,

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

$$V_i = \Box_{\text{ret}}^{-1} [-4\pi G \sigma_i]. \quad (146b)$$

From the 2PN order we have the potentials

$$\hat{X} = \Box_{\text{ret}}^{-1} \left[ -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 \right], \quad (147a)$$

$$\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 (147b)$$

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

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. (141); thus it involves in principle the finite part regularization $\mathcal{F}\mathcal{P}$ to deal with (IR-type) divergences occurring at high post-Newtonian 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 [68].

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

$$\hat{T} = \Box_{\text{ret}}^{-1} \left[ -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\partial_i V - 2\partial_i V_j \partial_j \hat{R}_i - \partial_i V_j \partial_t \hat{W}_{ij} \right. \\ \left. + V V_i \partial_t \partial_i V + 2V_i \partial_j V_i \partial_j V + \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], \quad (148a)$$

$$\hat{Y}_i = \Box_{\text{ret}}^{-1} \left[ -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) \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 - \frac{3}{2} V_k \partial_i V \partial_k V \\ \left. - \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], \quad (148b)$$

$$\hat{Z}_{ij} = \Box_{\text{ret}}^{-1} \left[ -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_{j} V_{i} - 2\partial_{(i} V_{k} \partial_{k} V_{j)} \right]$$
$$ -\frac{3}{4}\delta_{ij}(\partial_t V)^2 - \delta_{ij}\partial_k V_m(\partial_k V_m - \partial_m V_k)] . \qquad (148c) $$

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. (144).

Note that the 3PN metric (144) 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 (146)–(148); we elaborate more on radiation-reaction effects in the next Section 5.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. (27). These identities read

$$ \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} \left[ \hat{R}_i + V V_i \right] + \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), $$

$$ \partial_t \left\{ V_i + \frac{2}{c^2} \left[ \hat{R}_i + V V_i \right] \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). \quad (149b) $$

For latter applications to systems of compact objects, let us give the geodesic equations of a particle moving in the 3.5PN metric (144).<fnref n="36" /> It is convenient to write these equations as

$$ \frac{dP^i}{dt} = F^i, \qquad (150) $$

where the “linear momentum density” $P^i$ and the “force density” $F^i$ of the particle are given by

$$ P^i = c \frac{g_{i\mu} v^\mu}{\sqrt{-g_{\rho\sigma} v^\rho v^\sigma}}, \qquad (151a) $$

$$ F^i = \frac{c}{2} \frac{\partial_i g_{\mu\nu} v^\mu v^\nu}{\sqrt{-g_{\rho\sigma} v^\rho v^\sigma}}, \qquad (151b) $$

where $v^\mu = (c, v^i)$ with $v^i = dx^i/dt$ being the particle’s ordinary coordinate velocity, and where the metric components are taken at the location of the particle. Notice that we are here viewing the particle as moving in the fixed background metric (144). In Part B of this article, the metric will be generated by the system of particles itself, and we shall have to supplement the computation of the metric at the location of one of these particles by a suitable self-field regularization.

The expressions of both $P^i$ and $F^i$ in terms of the non-linear potentials follow from insertion of the 3.5PN metric coefficients (144). We obtain some complicated-looking (but useful in applications) sums of products of potentials given by

$$ P^i = v^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 + \frac{9}{2} V^2 v^i - 4 V V_i + 4 \hat{W}_{ij} v^j - 8 \hat{R}_i \right) \qquad (152a) $$

<footnote n="36">³⁶ Of course the geodesic equations are appropriate for the motion of particles without spins; for spinning particles one has also to take into account the coupling of the spin to the space-time curvature, see Eq. (377).</footnote>
$$
\begin{align}
&+ \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 + 2 \hat{W}_{ij} v^j v^2 \right. \nonumber \\
&\qquad \left. + 2 \hat{W}_{jk} v^i v^j v^k - 10 V V_i v^2 - 20 V V_j v^i v^j - 4 \hat{R}_i v^2 - 8 \hat{R}_j v^i v^j + \frac{9}{2} V^3 v^i + 12 V_j V_j v^i \right. \nonumber \\
&\qquad \left. + 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), \tag{152b} \\
&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) \nonumber \\
&\qquad + \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. \nonumber \\
&\qquad \left. - 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) \nonumber \\
&\qquad + \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 - 10 V \partial_i V_j v^2 v^j \right. \nonumber \\
&\qquad \left. - 4 \partial_i \hat{R}_j v^2 v^j + \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 \right. \nonumber \\
&\qquad \left. + 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 \right. \nonumber \\
&\qquad \left. - 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 \right. \nonumber \\
&\qquad \left. - 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). \nonumber
\end{align}
$$

Note that it will be supposed that all the accelerations appearing in the potentials and in the
final expression of the equations of motion are order-reduced by means of the equations of motion
themselves. For instance, we see from Eq. (152a) that when computing the time-derivative of $P_i$
we shall meet an acceleration at 1PN order which is 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^{\alpha\mu} = 0$, represent four out of the
ten Einstein field equations, see Section 2.1 for discussion. In the harmonic-coordinate approach the
equations of motion are equivalent to the harmonic gauge conditions $\partial_\mu h^{\alpha\mu} = 0$. Thus, each time
we get an acceleration in some PN expression (including the PN expression of the acceleration
itself), we have also another equation (or the same equation) which tells that the acceleration
is given by another PN expression. The post-Newtonian method assumes that it is legitimate
to replace that acceleration and to re-expand consistently with the PN order. Post-Newtonian
predictions based on such consistent PN order-reduction have been very successful.<fnref n="37" />

5.4 Radiation reaction potentials to 4PN order

We said that the metric (144) contains the correct radiation-reaction terms appropriate for an
isolated system up to the 3.5PN level included. The metric can even be generalized to include the
radiation-reaction terms up to 4PN order. To show this we shall use a particular non-harmonic
coordinate system to describe the radiation reaction terms up to 4PN order, which constitutes a
natural generalization of the Burke & Thorne [114, 113] coordinate system at 2.5PN order. Recall
that at the lowest 2.5PN order the radiation reaction force takes the simple form of Eq. (6), in

<footnote n="37">³⁷ Note, however, that the operation of order-reduction is illicit at the level of the Lagrangian. In fact, it is known that the elimination of acceleration terms 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 [374].</footnote>
which the force $F_i^{\text{reac}} = \rho \partial_i V^{\text{reac}}$ involves only a scalar potential given by

$$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 (153)$$

At such dominant 2.5PN level ("Newtonian" radiation reaction) the source quadrupole moment $Q_{ij}$ is simply given by the usual Newtonian expression (3).

The novel feature when one extends the Newtonian radiation reaction to include the 1PN corrections is that the reaction force 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 [56]. 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 1.5PN order of the lowest-order Burke & Thorne scalar radiation reaction potential (153), in some appropriate coordinate system, following Refs. [43, 47].

At that level (corresponding to 4PN order in the metric), and in this particular coordinate system, it suffices to incorporate some radiation-reaction contributions into the scalar and vectorial potentials $V$ and $V_i$ which parametrize the metric in Eq. (144). We thus pose

$$\mathcal{V} = V^{\text{inst}} + V^{\text{reac}}, \qquad (154a)$$

$$\mathcal{V}_i = V_i^{\text{inst}} + V_i^{\text{reac}}. \qquad (154b)$$

Then the metric, accurate to 4PN order regarding the radiation-reaction contributions – we indicate this by using the symbol $\mathcal{O}^{\text{reac}}$ for the remainders – reads

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

$$g_{0i} = -\frac{4}{c^3}\mathcal{V}_i + \mathcal{O}^{\text{reac}}\left(\frac{1}{c^{10}}\right), \qquad (155b)$$

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

The other contributions, which are conservative (i.e., non radiative), are given up to 3PN order by the metric (144) in which all the potentials take the same form as in Eqs. (146)–(148), 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. (134). This is for instance what we have indicated in Eqs. (154) by writing $V^{\text{inst}}$ and $V_i^{\text{inst}}$. Up to 3.5PN order, in this particular coordinate system, the effect of all these retardations gets replaced by the effect of the radiation-reaction potentials $V^{\text{reac}}$ and $V_i^{\text{reac}}$; furthermore, at the 4PN order there is a modification of the scalar radiation-reaction potential that is imposed by gravitational-wave tails propagating in the wave zone [58]. The explicit form of these potentials is [43, 47]$^{38}$

$$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} \mathbf{x}^2 x^{ij} I_{ij}^{(7)}(t) \right] \\ - \frac{4G^2M}{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), \qquad (156a)$$

<footnote>$^{38}$ Recall the footnote 17 for our notation. In particular $\hat{x}^{ijk}$ in the vector potential denotes the STF combination $\hat{x}^{ijk} = x^{ijk} - \frac{r^2}{5}(x^i \delta^{jk} + x^j \delta^{ki} + x^k \delta^{ij})$ with $x^{ijk} = x^i x^j x^k$.</footnote>
$$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} J_{kl}^{(5)}(t) \right] + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (156b)$$

where the multipole moments $I_L$ and $J_L$ denote the source multipole moments defined in Eqs. (123). Witness the tail integral at 4PN order characterized by a logarithmic kernel; see Section 3.2.

The scalar potential $V^{\text{reac}}$ will obviously reproduce Eq. (153) at the dominant order. However, note that it is crucial to include in Eq. (156a) the 1PN correction in the source quadrupole moment $I_{ij}$. The mass-type moments $I_L$ to 1PN order (and the current-type $J_L$ to Newtonian order), read

$$I_L = \int d^3\mathbf{x} \left\{ \hat{x}_L \sigma + \frac{1}{2c^2(2\ell+3)} \mathbf{x}^2 \hat{x}_L \partial_t^2 \sigma - \frac{4(2\ell+1)}{c^2(\ell+1)(2\ell+3)} \hat{x}_{iL} \partial_t \sigma_i \right\} + \mathcal{O}\left(\frac{1}{c^4}\right), \quad (157a)$$

$$J_L = \int d^3\mathbf{x} \epsilon_{ab<i\ell} \hat{x}_{L-1>a} \sigma_b + \mathcal{O}\left(\frac{1}{c^2}\right). \quad (157b)$$

The matter source densities $\sigma$ and $\sigma_i$ are given in Eqs. (145). Note that the mass multipole moments $I_L$ extend only over the *compact support* of the source even at the 1PN order. Only at the 2PN order will they involve some non-compact supported contributions – i.e., some integrals extending up to infinity [44].

The 3.5PN radiation reaction force in the equations of motion of compact binary systems has been derived by Iyer & Will [258, 259] in an arbitrary gauge, based on the energy and angular momentum balance equations at the relative 1PN order. As demonstrated in Ref. [259] the expressions of the radiation scalar and vector radiation-reaction potentials (156), which are valid in a particular gauge but are here derived from first principles, are fully consistent with the works [258, 259].

With the radiation-reaction potentials (156) in hands, one can *prove* [47] the energy balance equation up to 1.5PN order, namely

$$\frac{dE^{4\text{PN}}}{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). \quad (158)$$

One recognizes in the right-hand side the known positive-definite expression for the energy flux at 1.5PN order. Indeed the effective quadrupole moment which appears in the parenthesis of (158) agrees with the tail-modified *radiative* quadrupole moment $U_{ij}$ parametrizing the field in the far zone; see Eq. (90) where we recall that $M_L$ and $I_L$ are identical up to 2.5PN order.
# Part B: Compact Binary Systems

The problem of the motion and gravitational radiation of compact objects in post-Newtonian approximations is of crucial importance, for at least three reasons listed in the Introduction of this article: Motion of *N* planets in the solar system; gravitational radiation reaction force in binary pulsars; direct detection of gravitational waves from inspiralling compact binaries. As discussed in Section 1.3, the appropriate theoretical description of inspiralling compact binaries is by two structureless point-particles, characterized solely by their masses $m_1$ and $m_2$ (and possibly their spins), and moving on a quasi-circular orbit.

Strategies to detect and analyze the very weak signals from compact binary inspiral involve matched filtering of a set of accurate theoretical template waveforms against the output of the detectors. Many analyses [139, 137, 198, 138, 393, 346, 350, 284, 157, 158, 159, 156, 105, 106, 3, 18, 111] have shown that, in order to get sufficiently accurate theoretical templates, one must include post-Newtonian effects up to the 3PN level or higher. Recall that in practice, the post-Newtonian templates for the inspiral phase have to be matched to numerical-relativity results for the subsequent merger and ringdown phases. The match proceeds essentially through two routes: Either the so-called Hybrid templates obtained by direct matching between the PN expanded waveform and the numerical computations [4, 371], or the Effective-One-Body (EOB) templates [108, 109, 161, 168] that build on post-Newtonian results and extend their realm of validity to facilitate the analytical comparison with numerical relativity [112, 329]. Note also that various post-Newtonian resummation techniques, based on Padé approximants, have been proposed to improve the efficiency of PN templates [157, 158, 161].

## 6 Regularization of the Field of Point Particles

Our aim is to compute the metric (and its gradient needed in the equations of motion) at the 3PN order (say) for a system of two point-like particles. A *priori* one is not allowed to use directly some metric expressions like Eqs. (144) above, which have been derived under the assumption of a continuous (smooth) matter distribution. Applying them to a system of point particles, we find that most of the integrals become divergent at the location of the particles, i.e., when $x \to y_1(t)$ or $y_2(t)$, where $y_1(t)$ and $y_2(t)$ denote the two trajectories. Consequently, we must supplement the calculation by a prescription for how to remove the infinite part of these integrals. At this stage different choices for a “self-field” regularization (which will take care of the infinite self-field of point particles) are possible. In this section we review the:

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

2. Dimensional self-field regularization, an extremely powerful regularization which is free of any ambiguities (at least up to the 3PN level), and therefore permits to uniquely fix the values of the ambiguity parameters coming from Hadamard's regularization. However, dimensional regularization has not yet been implemented to the present problem in the general case (i.e., for an arbitrary space dimension $d \in \mathbb{C}$).

The why and how the final results are unique and independent of the employed self-field regularization (in agreement with the physical expectation) stems from the effacing principle of general relativity [142] – namely that the internal structure of the compact bodies makes a contribution only at the formal 5PN approximation. However, we shall review several alternative computations, independent of the self-field regularization, which confirm the end results.
## 6.1 Hadamard self-field regularization

In most practical computations we employ the Hadamard regularization [236, 381] (see Ref. [382] for an entry to the mathematical literature). Let us present here an account of this regularization, as well as a theory of generalized functions (or pseudo-functions) associated with it, following the detailed investigations in Refs. [70, 72].

Consider the class $\mathcal{F}$ of functions $F(\mathbf{x})$ which are smooth ($C^\infty$) on $\mathbb{R}^3$ except for the two points $\mathbf{y}_1$ and $\mathbf{y}_2$, around which they admit a power-like singular expansion of the type:<fnref n="39" />

$$ \forall \mathcal{N} \in \mathbb{N}, \quad F(\mathbf{x}) = \sum_{a_0 \le a \le \mathcal{N}} r_1^a f_a(\mathbf{n}_1) + o(r_1^\mathcal{N}), \qquad (159) $$

and similarly for the other point 2. Here $r_1 = |\mathbf{x} - \mathbf{y}_1| \to 0$, and the coefficients $f_a$ 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 $a$ of $r_1$ are real, range in discrete steps [i.e., $a \in (a_i)_{i \in \mathbb{N}}$], and are bounded from below ($a_0 \le a$). The coefficients $f_a$ (and $2f_a$) for which $a < 0$ can be referred to as the singular coefficients of $F$. If $F$ and $G$ belong to $\mathcal{F}$ 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 (160) $$

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. (160), the Hadamard partie finie is “non-distributive” in the sense that

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

The non-distributivity of Hadamard’s partie finie is the main source of the appearance of ambiguity parameters at the 3PN order, as discussed in Section 6.2.

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

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

The first term integrates over a domain $S(s)$ defined as $\mathbb{R}^3$ from which the two spherical balls $r_1 \le s$ and $r_2 \le s$ of radius $s$ and centered on the two singularities, denoted $\mathcal{B}(\mathbf{y}_1, s)$ and $\mathcal{B}(\mathbf{y}_2, s)$, are excised: $S(s) = \mathbb{R}^3 \setminus \mathcal{B}(\mathbf{y}_1, s) \cup \mathcal{B}(\mathbf{y}_2, s)$. The other terms, where the value of a function at point 1 takes the meaning (160), 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 the other 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. (162). We shall look for the fate of these constants in the final equations of motion and radiation field. See Ref. [70] for alternative expressions of the partie-finie integral.

We now come to a specific variant of Hadamard’s regularization called the extended Hadamard regularization (EHR) and defined in Refs. [70, 72]. The basic idea is to associate to any $F \in \mathcal{F}$

<footnote n="39">³⁹ The function $F(\mathbf{x})$ depends also on (coordinate) time $t$, through for instance its dependence on the velocities $\mathbf{v}_1(t)$ and $\mathbf{v}_2(t)$, but the time $t$ is purely “spectator” in the regularization process, and thus will not be indicated. See the footnote 20 for the definition of the Landau symbol $o$ for remainders.</footnote>
a pseudo-function, called the partie finie pseudo-function PfF, namely a linear form acting on functions G of F, and which is defined by the duality bracket

$$ \forall G \in \mathcal{F}, \quad \langle \mathrm{Pf}F, G \rangle = \mathrm{Pf} \int d^3\mathbf{x} FG. \qquad (163) $$

When restricted to the set $\mathcal{D}$ of smooth functions, i.e., $C^\infty(\mathbb{R}^4)$, with compact support (obviously we have $\mathcal{D} \subset \mathcal{F}$), the pseudo-function PfF is a distribution in the sense of Schwartz [381]. The product of pseudo-functions coincides, by definition, with the ordinary point-wise product, namely $\mathrm{Pf}F \cdot \mathrm{Pf}G = \mathrm{Pf}(FG)$. In practical computations, we use an interesting pseudo-function, constructed on the basis of the Riesz delta function [365], which plays a role analogous to the Dirac measure in distribution theory, $\delta_1(\mathbf{x}) \equiv \delta(\mathbf{x} - \mathbf{y}_1)$. This is the delta-pseudo-function $\mathrm{Pf}\delta_1$ defined by

$$ \forall F \in \mathcal{F}, \quad \langle \mathrm{Pf}\delta_1, F \rangle = \mathrm{Pf} \int d^3\mathbf{x} \delta_1 F = (F)_1, \qquad (164) $$

where $(F)_1$ is the partie finie of F as given by Eq. (160). From the product of $\mathrm{Pf}\delta_1$ with any PfF we obtain the new pseudo-function $\mathrm{Pf}(F\delta_1)$, that is such that

$$ \forall G \in \mathcal{F}, \quad \langle \mathrm{Pf}(F\delta_1), G \rangle = (FG)_1. \qquad (165) $$

As a general rule, we are not allowed, in consequence of the “non-distributivity” of the Hadamard partie finie, Eq. (161), to replace F within the pseudo-function $\mathrm{Pf}(F\delta_1)$ by its regularized value: $\mathrm{Pf}(F\delta_1) \neq (F)_1 \mathrm{Pf}\delta_1$ in general. It should be noticed that the object $\mathrm{Pf}(F\delta_1)$ has no equivalent in distribution theory.

Next, we treat the spatial derivative of a pseudo-function of the type PfF, namely $\partial_i(\mathrm{Pf}F)$. Essentially, we require [70] that the rule of integration by parts holds. By this we mean that we are allowed to freely operate by parts any duality bracket, with the all-integrated (“surface”) terms always being zero, as in the case of non-singular functions. This requirement is motivated by our will that a computation involving singular functions be as much as possible the same as if we were dealing with regular functions. Thus, by definition,

$$ \forall F, G \in \mathcal{F}, \quad \langle \partial_i(\mathrm{Pf}F), G \rangle = -\langle \partial_i(\mathrm{Pf}G), F \rangle. \qquad (166) $$

Furthermore, we assume that when all the singular coefficients of F vanish, the derivative of PfF reduces to the ordinary derivative, i.e., $\partial_i(\mathrm{Pf}F) = \mathrm{Pf}(\partial_i F)$. Then it is trivial to check that the rule (166) contains as a particular case the standard definition of the distributional derivative [381]. Notably, we see that the integral of a gradient is always zero: $\langle \partial_i(\mathrm{Pf}F), 1 \rangle = 0$. This should certainly be the case if we want to compute a quantity like a Hamiltonian density which is defined only modulo a total divergence. We pose

$$ \partial_i(\mathrm{Pf}F) = \mathrm{Pf}(\partial_i F) + D_i[F], \qquad (167) $$

where $\mathrm{Pf}(\partial_i F)$ represents the “ordinary” derivative and $D_i[F]$ is the distributional term. The following solution of the basic relation (166) was obtained in Ref. [70]:

$$ D_i[F] = 4\pi \mathrm{Pf} \left( n_1^i \left[ \frac{1}{2} r_1 f_{-1} + \sum_{k \ge 0} \frac{1}{r_1^k} f_{-2-k} \right] \delta_1 \right) + 1 \leftrightarrow 2, \qquad (168) $$

where for simplicity we assume that the powers a in the expansion (159) of F are relative integers. The distributional term (168) is of the form $\mathrm{Pf}(G\delta_1)$ plus 1 $\leftrightarrow$ 2; it is generated solely by the
singular coefficients of $F$.<fnref n="40" /> The formula for the distributional term associated with the $ℓ$-th distributional derivative, i.e. $D_L[F] = \partial_L P_f F - P_f \partial_L F$, where $L = i_1 i_2 \cdots i_ℓ$, reads

$$D_L[F] = \sum_{k=1}^{\ell} \partial_{i_1 \dots i_{k-1}} D_{i_k} [\partial_{i_{k+1} \dots i_{\ell}} F]. \qquad (169)$$

We refer to Theorem 4 in Ref. [70] for the definition of another derivative operator, representing the most general derivative satisfying the same properties as the one defined by Eq. (168), and, in addition, the commutation of successive derivatives (or Schwarz lemma).<fnref n="41" />

The distributional derivative defined by (167)–(168) does not satisfy the Leibniz rule for the derivation of a product, in accordance with a general result of Schwartz [380]. Rather, the investigation of Ref. [70] suggests that, in order to construct a consistent theory (using the ordinary point-wise product for pseudo-functions), the Leibniz rule should be weakened, and replaced by the rule of integration by part, Eq. (166), which is in fact nothing but an integrated version of the Leibniz rule. However, the loss of the Leibniz rule *stricto sensu* constitutes one of the reasons for the appearance of the ambiguity parameters at 3PN order, see Section 6.2.

The Hadamard regularization ($F$)$_1$ is defined by Eq. (160) in a preferred spatial hypersurface $t = \text{const}$ of a coordinate system, and consequently is not a *priori* compatible with the Lorentz invariance. Thus we expect that the equations of motion in harmonic coordinates (which manifestly preserve the global Lorentz invariance) should exhibit at some stage a violation of the Lorentz invariance due to the latter regularization. In fact this occurs exactly at the 3PN order. Up to the 2.5PN level, the use of the regularization ($F$)$_1$ is sufficient to get some unambiguous equations of motion which are Lorentz invariant [76]. This problem can be dealt with within Hadamard’s regularization, by introducing a Lorentz-invariant variant of this regularization, denoted $[F]_1$ in Ref. [72]. It consists of performing the Hadamard regularization within the spatial hypersurface that is geometrically orthogonal (in a Minkowskian sense) to the four-velocity of the particle. The regularization $[F]_1$ differs from the simpler regularization ($F$)$_1$ by relativistic corrections of order $1/c^2$ at least. See [72] for the formulas defining this regularization in the form of some infinite power series in $1/c^2$. The regularization $[F]_1$ plays a crucial role in obtaining the equations of motion at the 3PN order in Refs. [69, 71]. In particular, the use of the Lorentz-invariant regularization $[F]_1$ permits to obtain the value of the ambiguity parameter $\omega_{\text{kinetic}}$ in Eq. (170a) below.

## 6.2 Hadamard regularization ambiguities

The standard Hadamard regularization yields some ambiguous results for the computation of certain integrals at the 3PN order, as noticed by Jaranowski & Schäfer [261, 262, 263] in their computation of the equations of motion within the ADM-Hamiltonian formulation of general relativity. By standard Hadamard regularization we mean the regularization based solely on the definitions of the partie finie of a singular function, Eq. (160), and the partie finie of a divergent integral, Eq. (162), and without using a theory of pseudo-functions and generalized distributional derivatives as in Refs. [70, 72]. It was shown in Refs. [261, 262, 263] that there are two and only two types of ambiguous terms in the 3PN Hamiltonian, which were then parametrized by two unknown numerical coefficients called $\omega_{\text{static}}$ and $\omega_{\text{kinetic}}$.

Progressing concurrently, Blanchet & Faye [70, 72] introduced the “extended” Hadamard regularization – the one we outlined in Section 6.1 – and obtained [69, 71] the 3PN equations of motion

<footnote n="40">⁴⁰ The sum over $k$ in Eq. (168) is always finite since there is a maximal order $a_0$ of divergency in Eq. (159).</footnote>

<footnote n="41">⁴¹ It was shown in Ref. [71] that using one or the other of these derivatives results in some equations of motion that differ by a coordinate transformation, and the redefinition of some ambiguity parameter. This indicates that the distributional derivatives introduced in Ref. [70] constitute some technical tools devoid of physical meaning besides precisely the appearance of Hadamard’s ambiguity parameters.</footnote>
complete except for one and only one unknown numerical constant, called $\lambda$. The new extended Hadamard regularization is mathematically well-defined and yields unique results for the computation of any integral in the problem; however, it turned out to be in a sense “incomplete” as it could not determine the value of this constant. The comparison of the result with the work [261, 262], on the basis of the computation of the invariant energy of compact binaries moving on circular orbits, revealed [69] that

$$ \omega_{\text{kinetic}} = \frac{41}{24}, \qquad (170a) $$

$$ \omega_{\text{static}} = -\frac{11}{3}\lambda - \frac{1987}{840}. \qquad (170b) $$

Therefore, the ambiguity $\omega_{\text{kinetic}}$ is fixed, while $\lambda$ is equivalent to the other ambiguity $\omega_{\text{static}}$. Notice that the value (170a) for the kinetic ambiguity parameter $\omega_{\text{kinetic}}$, which is in factor of some velocity dependent terms, is the only one for which the 3PN equations of motion are Lorentz invariant. Fixing up this value was possible because the extended Hadamard regularization [70, 72] was defined in such a way that it keeps the Lorentz invariance.

The value of $\omega_{\text{kinetic}}$ given by Eq. (170a) was recovered in Ref. [162] by directly proving that such value is the unique one for which the global Poincaré invariance of the ADM-Hamiltonian formalism is verified. Since the coordinate conditions associated with the ADM formalism do not manifestly respect the Poincaré symmetry, it was necessary to prove that the 3PN Hamiltonian is compatible with the existence of generators for the Poincaré algebra. By contrast, the harmonic-coordinate conditions preserve the Poincaré invariance, and therefore the associated equations of motion at 3PN order are manifestly Lorentz-invariant, as was found to be the case in Refs. [69, 71].

The appearance of one and only one physical unknown coefficient $\lambda$ in the equations of motion constitutes a quite striking fact, that is related specifically with the use of a Hadamard-type regularization.42 Technically speaking, the presence of the ambiguity parameter $\lambda$ is associated with the non-distributivity of Hadamard's regularization, in the sense of Eq. (161). Mathematically speaking, $\lambda$ is probably related to the fact that it is impossible to construct a distributional derivative operator, such as Eqs. (167)–(168), satisfying the Leibniz rule for the derivation of the product [380]. The Einstein field equations can be written in many different forms, by shifting the derivatives and operating some terms by parts with the help of the Leibniz rule. All these forms are equivalent in the case of regular sources, but since the derivative operator (167)–(168) violates the Leibniz rule they become inequivalent for point particles.

Physically speaking, let us also argue that $\lambda$ has its root in the fact that in a complete computation of the equations of motion valid for two regular extended weakly self-gravitating bodies, many non-linear integrals, when taken individually, start depending, from the 3PN order, on the internal structure of the bodies, even in the “compact-body” limit where the radii tend to zero. However, when considering the full equations of motion, one expects that all the terms depending on the internal structure can be removed, in the compact-body limit, by a coordinate transformation (or by some appropriate shifts of the central world lines of the bodies), and that finally $\lambda$ is given by a pure number, for instance a rational fraction, independent of the details of the internal structure of the compact bodies. From this argument (which could be justified by the effacing principle in general relativity) the value of $\lambda$ is necessarily the one we compute below, Eq. (172), and will be valid for any compact objects, for instance black holes.

The ambiguity parameter $\omega_{\text{static}}$, which is in factor of some static, velocity-independent term, and hence cannot be derived by invoking Lorentz invariance, was computed by Damour, Jaranowski

<footnote n="42">42 Note also that the harmonic-coordinates 3PN equations of motion [69, 71] depend, in addition to $\lambda$, on two arbitrary constants $r'_1$ and $r'_2$ parametrizing some logarithmic terms. These constants are closely related to the constants $s_1$ and $s_2$ in the *partie-finie* integral (162); see Ref. [71] and Eq. (185) below for the precise definition. However, $r'_1$ and $r'_2$ are not “physical” in the sense that they can be removed by a coordinate transformation.</footnote>
& Schäfer [163] by means of dimensional regularization, instead of some Hadamard-type one, within the ADM-Hamiltonian formalism. Their result is

$$ \omega_{\text{static}} = 0. \qquad (171) $$

As argued in [163], clearing up the static ambiguity is made possible by the fact that dimensional regularization, contrary to Hadamard's regularization, respects all the basic properties of the algebraic and differential calculus of ordinary functions: Associativity, commutativity and distributivity of point-wise addition and multiplication, Leibniz's rule, and the Schwarz lemma. In this respect, dimensional regularization is certainly superior to Hadamard's one, which does not respect the distributivity of the product [recall Eq. (161)] and unavoidably violates at some stage the Leibniz rule for the differentiation of a product.

The ambiguity parameter $\lambda$ is fixed from the result (171) and the necessary link (170b) provided by the equivalence between the harmonic-coordinates and ADM-Hamiltonian formalisms [69, 164]. However, $\lambda$ has also been computed directly by Blanchet, Damour & Esposito-Farèse [61] applying dimensional regularization to the 3PN equations of motion in harmonic coordinates (in the line of Refs. [69, 71]). The end result,

$$ \lambda = -\frac{1987}{3080}, \qquad (172) $$

is in full agreement with Eq. (171).43 Besides the independent confirmation of the value of $\omega_{\text{static}}$ or $\lambda$, the work [61] provides also a confirmation of the consistency of dimensional regularization, since the explicit calculations are entirely different from the ones of Ref. [163]: Harmonic coordinates instead of ADM-type ones, work at the level of the equations of motion instead of the Hamiltonian, a different form of Einstein's field equations which is solved by a different iteration scheme.

Let us comment that 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 in Refs. [280, 234] which 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.44 At the 3PN order we expect that the extended-body program should give the value of the regularization parameter $\lambda$ – probably after a coordinate transformation to remove the terms depending on the internal structure. Ideally, its value should also be confirmed by independent and more physical methods like those of Refs. [407, 281, 172].

An important work, in several aspects more physical than the formal use of regularizations, is the one of Itoh & Futamase [255, 253, 254], following previous investigations in Refs. [256, 257]. These authors derived the 3PN equations of motion in harmonic coordinates by means of a particular variant of the famous “surface-integral” method $\partial la$ Einstein, Infeld & Hoffmann [184]. The aim is to describe extended relativistic compact binary systems in the so-called strong-field point particle limit which has been defined in Ref. [212]. This approach is interesting because it is based on the physical notion of extended compact bodies in general relativity, and is free of the problems of ambiguities. The end result of Refs. [255, 253] is in agreement with the 3PN harmonic

<footnote n="43">43 One may wonder why the value of $\lambda$ is a complicated rational fraction while $\omega_{\text{static}}$ is so simple. This is because $\omega_{\text{static}}$ was introduced precisely to measure the amount of ambiguities of certain integrals, while by contrast, $\lambda$ was introduced [see Eq. (185)] as an unknown constant entering the relation between the arbitrary scales $r'_1, r'_2$ on the one hand, and $s_1, s_2$ on the other hand, which has a priori nothing to do with the ambiguous part of integrals.</footnote>

<footnote n="44">44 See however some comments on the latter work in Ref. [145], pp. 168–169.</footnote>
coordinates equations of motion [69, 71] and is unambiguous, as it does directly determine the ambiguity parameter $\lambda$ to exactly the value (172).

The 3PN equations of motion in harmonic coordinates or, more precisely, the associated 3PN Lagrangian, were also derived by Foffa & Sturani [203] using another important approach, coined the effective field theory (EFT) [223]. Their result is fully compatible with the value (172) for the ambiguity parameter $\lambda$; however, in contrast with the surface-integral method of Refs. [255, 253], this does not check the method of regularization because the EFT approach is also based on dimensional self-field regularization.

We next consider the problem of the binary's radiation field, where the same phenomenon occurs, with the appearance of some Hadamard regularization ambiguity parameters at 3PN order. More precisely, Blanchet, Iyer & Joguet [81], computing the 3PN compact binary's mass quadrupole moment $I_{ij}$, found it necessary to introduce three Hadamard regularization constants $\xi$, $\kappa$, and $\zeta$, which are independent of the equation-of-motion related constant $\lambda$. The total gravitational-wave flux at 3PN order, in the case of circular orbits, was found to depend on a single combination of the latter constants, $\theta = \xi + 2\kappa + \zeta$, and the binary's orbital phase, for circular orbits, involved only the linear combination of $\theta$ and $\lambda$ given by $\hat{\theta} = \theta - \frac{7}{3}\lambda$, as shown in [73].

Dimensional regularization (instead of Hadamard's) has next been applied in Refs. [62, 63] to the computation of the 3PN radiation field of compact binaries, leading to the following unique determination of the ambiguity parameters:<fnref n="45" />

$$ \xi = -\frac{9871}{9240}, \qquad (173a) $$

$$ \kappa = 0, \qquad (173b) $$

$$ \zeta = -\frac{7}{33}. \qquad (173c) $$

These values represent the end result of dimensional regularization. However, several alternative calculations provide a check, independent of dimensional regularization, for all the parameters (173). One computes [80] the 3PN binary's mass dipole moment $I_i$ using Hadamard's regularization, and identifies $I_i$ with the 3PN center of mass vector position $G_i$, already known as a conserved integral associated with the Poincaré invariance of the 3PN equations of motion in harmonic coordinates [174]. This yields $\xi + \kappa = -9871/9240$ in agreement with Eqs. (173). Next, one considers [65] the limiting physical situation where the mass of one of the particles is exactly zero (say, $m_2 = 0$), and the other particle moves with uniform velocity. Technically, the 3PN quadrupole moment of a boosted Schwarzschild black hole is computed and compared with the result for $I_{ij}$ in the limit $m_2 = 0$. The result is $\zeta = -7/33$, and represents a direct verification of the global Poincaré invariance of the wave generation formalism (the parameter $\zeta$ representing the analogue for the radiation field of the parameter $\omega_{\text{kinetic}}$). Finally, one proves [63] that $\kappa = 0$ by showing that there are no dangerously divergent diagrams corresponding to non-zero $\kappa$-values, where a diagram is meant here in the sense of Ref. [151].

The determination of the parameters (173) completes the problem of the general relativistic prediction for the templates of inspiralling compact binaries up to 3.5PN order. The numerical values of these parameters indicate, following measurement-accuracy analyses [105, 106, 159, 156], that the 3.5PN order should provide an excellent approximation for both the on-line search and the subsequent off-line analysis of gravitational wave signals from inspiralling compact binaries in the LIGO and VIRGO detectors.

<footnote n="45">⁴⁵ The result for $\xi$ happens to be amazingly related to the one for $\lambda$ by a cyclic permutation of digits; compare $3\xi = -9871/3080$ with $\lambda = -1987/3080$.</footnote>
## 6.3 Dimensional regularization of the equations of motion

As reviewed in Section 6.2, work at 3PN order using Hadamard’s self-field regularization showed the appearance of ambiguity parameters, due to an incompleteness of the Hadamard regularization employed for curing the infinite self field of point particles. We give here more details on the determination using dimensional regularization of the ambiguity parameter $\lambda$ [or equivalently $\omega_{\text{static}}$, see Eq. (170b)] which appeared in the 3PN equations of motion.

Dimensional regularization was invented as a means to preserve the gauge symmetry of perturbative quantum field theories [391, 91, 100, 131]. 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, we shall prove that dimensional regularization determines the kinetic ambiguity parameter $\omega_{\text{kinetic}}$ (and its radiation-field analogue $\zeta$), and is therefore able to correctly keep track of the global Lorentz-Poincaré invariance of the gravitational field of isolated systems. The dimensional regularization is also an important ingredient of the EFT approach to equations of motion and gravitational radiation [223].

The Einstein field equations in $d+1$ space-time dimensions, relaxed by the condition of harmonic coordinates $\partial_{\mu}h^{\alpha\mu} = 0$, take exactly the same form as given in Eqs. (18)–(23). In particular the box operator $\square$ now denotes the flat space-time d’Alembertian operator in $d+1$ dimensions with signature $(-1, 1, 1, \dots)$. The gravitational constant $G$ is related to the usual three-dimensional Newton’s constant $G_N$ by

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

where $\ell_0$ denotes an arbitrary length scale. The explicit expression of the gravitational source term $\Lambda^{\alpha\beta}$ involves some $d$-dependent coefficients, and is given by

$$\begin{aligned} \Lambda^{\alpha\beta} = & -h^{\mu\nu}\partial_{\mu}^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})\left(g_{\lambda\tau}g_{\epsilon\pi} - \frac{1}{d-1}g_{\tau\epsilon}g_{\lambda\pi}\right)\partial_{\mu}h^{\lambda\pi}\partial_{\nu}h^{\tau\epsilon}. \end{aligned} \qquad (175)$$

When $d=3$ we recover Eq. (24). In the following we assume, as usual 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 shall often pose $\varepsilon \equiv d-3$.

We parametrize the 3PN metric in $d$ dimensions 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 Section 5.3. Those are obtained by post-Newtonian iteration of the $d$-dimensional field equations, starting from appropriate definitions of matter source densities generalizing Eqs. (145):

$$\sigma = \frac{2}{d-1} \frac{(d-2)T^{00} + T^{ii}}{c^2}, \qquad (176a)$$

$$\sigma_i = \frac{T^{0i}}{c}, \qquad (176b)$$

$$\sigma_{ij} = T^{ij}. \qquad (176c)$$

As a result all the expressions of Section 5.3 acquire some explicit $d$-dependent coefficients. For instance we find [61]

$$V = \square_{\text{ret}}^{-1} [-4\pi G\sigma], \qquad (177a)$$
$$ \hat{W}_{ij} = \Box_{\text{ret}}^{-1} \left[ -4\pi G \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 (177b) $$

Here $\Box_{\text{ret}}^{-1}$ means the retarded integral in $d+1$ space-time dimensions, which admits, though, no simple expression generalizing Eq. (31) in physical $(t, \mathbf{x})$ space.<fnref n="46" />

As reviewed in Section 6.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 in powers and inverse powers of $r_1 \equiv |\mathbf{x} - \mathbf{y}_1|$ and $r_2 \equiv |\mathbf{x} - \mathbf{y}_2|$, given by Eq. (178). In $d$ spatial dimensions, there is an analogue of the function $F$, which results from the post-Newtonian 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, as it reads

$$ 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}}). \qquad (178) $$

The coefficients $f_{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$ and $q_1$ 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 at 3PN 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_{p,q}^{(\varepsilon=0)}(\mathbf{n}_1) = f_p(\mathbf{n}_1). \qquad (179) $$

For the problem at hand, we essentially have to deal with the regularization of Poisson integrals, or iterated Poisson integrals (and their gradients needed in the equations of motion), 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,

$$ P^{(d)}(\mathbf{x}') = \Delta^{-1} [F^{(d)}(\mathbf{x})] \equiv -\frac{\tilde{k}}{4\pi} \int \frac{d^d\mathbf{x}}{|\mathbf{x} - \mathbf{x}'|^{d-2}} F^{(d)}(\mathbf{x}), \qquad (180) $$

where $\tilde{k}$ is a constant related to the usual Eulerian $\Gamma$-function by<fnref n="47" />

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

We need 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 allow simply to get $[P^{(d)}(\mathbf{x}')]_{\mathbf{x}'=\mathbf{y}_1}$ by replacing $\mathbf{x}'$ by $\mathbf{y}_1$ inside the explicit integral (180). 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}). \qquad (182) $$

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

<footnote n="46">⁴⁶ In higher approximations there will be also IR divergences and one should really employ the $d$-dimensional version of Eq. (141).</footnote>

<footnote n="47">⁴⁷ 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>
It is not possible at present to compute the equations of motion in the general d-dimensional case, but only in the limit where $\epsilon \to 0$ [163, 61]. The main technical step of our strategy consists of computing, in the limit $\epsilon \to 0$, the difference between the d-dimensional Poisson potential (182), and its Hadamard 3-dimensional counterpart given by $(P)_1$, where the Hadamard partie finie is defined by Eq. (160). But we must be precise when defining the Hadamard partie finie of a Poisson integral. Indeed, the definition (160) stricto sensu is applicable when the expansion of the function $F$, for $r_1 \to 0$, does not involve logarithms of $r_1$; see Eq. (160). However, the Poisson integral $P(\mathbf{x}')$ of $F(\mathbf{x})$ will typically involve such logarithms at the 3PN order, namely some $\ln r_1'$ where $r_1' \equiv |\mathbf{x}' - \mathbf{y}_1|$ formally tends to zero (hence $\ln r_1'$ is formally infinite). The proper way to define the Hadamard partie finie in this case is to include the $\ln r_1'$ into its definition; we arrive at [70]

$$ (P)_1 = -\frac{1}{4\pi} \mathrm{Pf}_{r_1', s_2'} \int \frac{d^3\mathbf{x}}{r_1'} F(\mathbf{x}) - (r_1^2 F)_1. \qquad (183) $$

The first term follows from Hadamard's partie finie integral (162); the second one is given by Eq. (160). Notice that in this result the constant $s_1$ entering the partie finie integral (162) has been “replaced” by $r_1'$, which plays the role of a new regularization constant (together with $r_2'$ for the other particle), and which ultimately parametrizes the final Hadamard regularized 3PN equations of motion. It was shown that $r_1'$ and $r_2'$ are unphysical, in the sense that they can be removed by a coordinate transformation [69, 71]. On the other hand, the constant $s_2$ remaining in the result (183) is the source for the appearance of the physical ambiguity parameter $\lambda$. Denoting the difference between the dimensional and Hadamard 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 (184) $$

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. However, we shall only compute the first two terms of the Laurent expansion of $\mathcal{D}P_1$ when $\epsilon \to 0$, say $\mathcal{D}P_1 = a_{-1} \epsilon^{-1} + a_0 + \mathcal{O}(\epsilon)$. This is the information we need to clear up the ambiguity parameter. We insist that the difference $\mathcal{D}P_1$ comes exclusively from the contribution of terms developing some poles $\propto 1/\epsilon$ in the d-dimensional calculation.

Next we outline the way we obtain, starting from the computation of the “difference”, the 3PN equations of motion in dimensional regularization, and show how the ambiguity parameter $\lambda$ is determined. By contrast to $r_1'$ and $r_2'$ which are pure gauge, $\lambda$ is a genuine physical ambiguity, introduced in Refs. [70, 71] as the single unknown numerical constant parametrizing the ratio between $s_2$ and $r_2'$ [where $s_2$ is the constant left in Eq. (183)] as

$$ \ln\left(\frac{r_2'}{s_2}\right) = \frac{159}{308} + \lambda \frac{m_1 + m_2}{m_2} \quad (\text{and} 1 \leftrightarrow 2), \qquad (185) $$

where $m_1$ and $m_2$ are the two masses. The terms corresponding to the $\lambda$-ambiguity in the acceleration $\mathbf{a}_1 = dv_1/dt$ of particle 1 read simply

$$ \Delta \mathbf{a}_1[\lambda] = -\frac{44\lambda}{3} \frac{G_N^4 m_1 m_2^2 (m_1 + m_2)}{r_{12}^5 c^6} \mathbf{n}_{12}, \qquad (186) $$

where the relative distance between particles is denoted $\mathbf{y}_1 - \mathbf{y}_2 \equiv r_{12} \mathbf{n}_{12}$ (with $\mathbf{n}_{12}$ being the unit vector pointing from particle 2 to particle 1). We start from the end result of Ref. [71] for the 3PN harmonic coordinates acceleration $\mathbf{a}_1$ in Hadamard's regularization, abbreviated as HR. Since the result was obtained by means of the specific extended variant of Hadamard's regularization (in short EHR, see Section 6.1) we write it as

$$ \mathbf{a}_1^{(\mathrm{HR})} = \mathbf{a}_1^{(\mathrm{EHR})} + \Delta \mathbf{a}_1[\lambda], \qquad (187) $$
where $a_1^{(\text{EHR})}$ is a fully determined functional of the masses $m_1$ and $m_2$, the relative distance $r_{12} n_{12}$, the coordinate velocities $v_1$ and $v_2$, and also the gauge constants $r'_1$ and $r'_2$. The only ambiguous term is the second one and is given by Eq. (186).

Our strategy is to extract from both the dimensional and Hadamard regularizations their common core part, obtained by applying the so-called “pure-Hadamard-Schwartz” (pHS) regularization. Following the definition in Ref. [61], the pHS regularization is a specific, minimal Hadamard-type regularization of integrals, based on the partie finie integral (162), together with a minimal treatment of “contact” terms, in which the definition (162) is applied separately to each of the elementary potentials $V$, $V_i$, etc. (and gradients) that enter the post-Newtonian 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. (161). The pHS regularization also assumes the use of standard Schwartz distributional derivatives [381]. The interest of the pHS regularization is that the dimensional regularization is equal to it plus the “difference”; see Eq. (190).

To obtain the pHS-regularized acceleration we need to subtract from the EHR result a series of contributions, which are specific consequences of the use of EHR [70, 72]. For instance, one of these contributions corresponds to the fact that in the EHR the distributional derivative is given by Eqs. (167)–(168) which differs from the Schwartz distributional derivative in the pHS regularization. Hence we define

$$a_1^{(\text{pHS})} = a_1^{(\text{EHR})} - \sum \delta a_1, \qquad (188)$$

where the $\delta a_1$’s denote the extra terms following from the EHR prescriptions. The pHS-regularized acceleration (188) constitutes essentially the result of the first stage of the calculation of $a_1$, as reported in Ref. [193].

The next step consists of evaluating the Laurent expansion, in powers of $\epsilon = d - 3$, of the difference between the dimensional regularization and the pHS (3-dimensional) computation. As we reviewed above, this 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 $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, and must be evaluated with care. We denote the above defined difference by

$$\mathcal{D}a_1 = \sum \mathcal{D}P_1. \qquad (189)$$

It is made of the sum of all the individual differences of Poisson or Poisson-like integrals as computed in Eq. (184). The total difference (189) depends on the Hadamard regularization scales $r'_1$ and $s_2$ (or equivalently on $\lambda$ and $r'_1, r'_2$), and on the parameters associated with dimensional regularization, namely $\epsilon$ and the characteristic length scale $\ell_0$ introduced in Eq. (174). Finally, the result is the explicit computation of the $\epsilon$-expansion of the dimensional regularization (DR) acceleration as

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

With this result we can prove two theorems [61].

**Theorem 8.** The pole part $\propto 1/\epsilon$ of the DR acceleration (190) can be re-absorbed (i.e. renormalized) into some shifts of the two “bare” world-lines: $y_1 \to y_1 + \xi_1$ and $y_2 \to y_2 + \xi_2$, with $\xi_{1,2} \propto 1/\epsilon$ say, so that the result, expressed in terms of the “dressed” quantities, is finite when $\epsilon \to 0$.

The situation in harmonic coordinates is to be contrasted with the calculation in ADM-type coordinates within the Hamiltonian formalism, where it was shown that all pole parts directly cancel
out in the total 3PN Hamiltonian: No renormalization of the world-lines is needed [163]. The central result is then:

**Theorem 9.** The renormalized (finite) DR acceleration is physically equivalent to the Hadamard-regularized (HR) acceleration (end result of Ref. [71]), in the sense that

$$a_1^{(\text{HR})} = \lim_{\epsilon \to 0} \left[ a_1^{(\text{DR})} + \delta_{\xi} a_1 \right], \qquad (191)$$

where $\delta_{\xi} a_1$ denotes the effect of the shifts on the acceleration, if and only if the HR ambiguity parameter $\lambda$ entering the harmonic-coordinates equations of motion takes the unique value (172).

The precise shifts $\xi_1$ and $\xi_2$ needed in Theorem 9 involve not only a pole contribution $\propto 1/\epsilon$, but also a finite contribution when $\epsilon \to 0$. Their explicit expressions read:48

$$\xi_1 = \frac{11}{3} \frac{G_N^2 m_1^2}{c^6} \left[ \frac{1}{\epsilon} - 2 \ln \left( \frac{r_1 \bar{q}^{1/2}}{\ell_0} \right) - \frac{327}{1540} \right] a_1^{\text{N}} \quad (\text{together with } 1 \leftrightarrow 2), \qquad (192)$$

where $G_N$ is Newton's constant, $\ell_0$ is the characteristic length scale of dimensional regularization, cf. Eq. (174), where $a_1^N$ is the Newtonian acceleration of the particle 1 in $d$ dimensions, and $\bar{q} \equiv 4\pi e^{\gamma_E}$ depends on Euler's constant $\gamma_E \simeq 0.577$.

## 6.4 Dimensional regularization of the radiation field

We now address the similar problem concerning the binary's radiation field – to 3PN order beyond Einstein's quadrupole formalism (2)–(3). As reviewed in Section 6.2, three ambiguity parameters: $\xi$, $\kappa$ and $\zeta$, have been shown to appear in the 3PN expression of the quadrupole moment [81, 80].

To apply dimensional regularization, we must use as in Section 6.3 the $d$-dimensional post-Newtonian iteration leading to potentials such as those in Eqs. (177); and we have to generalize to $d$ dimensions some key results of the wave generation formalism of Part A. Essentially, we need the $d$-dimensional analogues of the multipole moments of an isolated source $I_L$ and $J_L$ in Eqs. (123). Here we report the result we find in the case of the mass-type moment:

$$I_L^{(d)}(t) = \frac{d-1}{2(d-2)} \mathcal{FP} \int d^d\mathbf{x} \left\{ \hat{x}_L \sum_{[\ell]} \frac{4(d+2\ell-2)}{c^2(d+\ell-2)(d+2\ell)} \hat{x}_{aL} \sum_{[\ell+1]}^{(1)} \mathbf{a}(\mathbf{x}, t) + \frac{2(d+2\ell-2)}{c^4(d+\ell-1)(d+\ell-2)(d+2\ell+2)} \hat{x}_{abL} \sum_{[\ell+2]}^{(2)} \mathbf{a}(\mathbf{x}, t) \right\}, \qquad (193)$$

in which we denote, generalizing Eqs. (124),

$$\Sigma = \frac{2}{d-1} \frac{(d-2)\bar{\tau}^{00} + \bar{\tau}^{ii}}{c^2}, \qquad (194a)$$

$$\Sigma_i = \frac{\bar{\tau}^{0i}}{c}, \qquad (194b)$$

$$\Sigma_{ij} = \bar{\tau}^{ij}, \qquad (194c)$$

and where for any source densities the underscript $[\ell]$ means the infinite series

$$\Sigma(\mathbf{x}, t) = \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} \Sigma(\mathbf{x}, t). \qquad (195)$$

<footnote n="48">48 When working at the level of the equations of motion (not considering the metric outside the world-lines), the effect of shifts can be seen as being induced by a coordinate transformation of the bulk metric as in Ref. [71].</footnote>
The latter definition represents the d-dimensional version of the post-Newtonian expansion series (126). At Newtonian order, the expression (193) reduces to the standard result $I_L^{(d)} = \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 ambiguity parameters $\xi$, $\kappa$ and $\zeta$ come from the Hadamard regularization of the mass quadrupole moment $I_{ij}$ at the 3PN order. The terms corresponding to these ambiguities were found to be

$$ \Delta I_{ij}[\xi, \kappa, \zeta] = \frac{44 G_N^2 m_1^3}{3 c^6} \left[ \left( \xi + \kappa \frac{m_1 + m_2}{m_1} \right) y_1^{\langle i} a_1^{j\rangle} + \zeta v_1^{\langle i} v_1^{j\rangle} \right] + 1 \leftrightarrow 2, \quad (196) $$

where $y_1$, $v_1$ and $a_1$ denote the first particle's position, velocity and acceleration (and the brackets $\langle\rangle$ surrounding indices refer to the STF projection). Like in Section 6.3, we express both the Hadamard and dimensional results in terms of the more basic pHS regularization. The first step of the calculation [80] is therefore to relate the Hadamard-regularized quadrupole moment $I_{ij}^{(\text{HR})}$, for general orbits, to its pHS part:

$$ I_{ij}^{(\text{HR})} = I_{ij}^{(\text{pHS})} + \Delta I_{ij} \left[ \xi + \frac{1}{22}, \kappa, \zeta + \frac{9}{110} \right]. \quad (197) $$

In the right-hand side we find both the pHS part, and the effect of adding the ambiguities, with some numerical shifts of the ambiguity parameters ($\xi \to \xi + 1/22$, $\zeta \to \zeta + 9/110$) due to the difference between the specific Hadamard-type regularization scheme used in Ref. [81] and the pHS one. The pHS part is free of ambiguities but depends on the gauge constants $r'_1$ and $r'_2$ introduced in the harmonic-coordinates equations of motion [69, 71].

We next use the d-dimensional moment (193) to compute the difference between the dimensional regularization (DR) result and the pHS one [62, 63]. As in the work on equations of motion, we find that the ambiguities arise solely from the terms in the integration regions near the particles, that give rise to poles $\propto 1/\epsilon$, corresponding to logarithmic ultra-violet (UV) divergences in 3 dimensions. The infra-red (IR) region at infinity, i.e., $|x| \to +\infty$, does not contribute to the difference between DR and pHS. The compact-support terms in the integrand of Eq. (193), proportional to the matter source densities $\sigma$, $\sigma_a$, and $\sigma_{ab}$, are also found not to contribute to the difference. We are therefore left with evaluating the difference linked with the computation of the non-compact terms in the expansion of the integrand of (193) near the singularities that produce poles in d dimensions.

Let $F^{(d)}(\mathbf{x})$ be the non-compact part of the integrand of the quadrupole moment (193) (with indices $L = ij$), where $F^{(d)}$ includes the appropriate multipolar factors such as $\hat{x}_{ij}$, so that

$$ I_{ij}^{(d)} = \int d^d\mathbf{x} F^{(d)}(\mathbf{x}). \quad (198) $$

We do not indicate that we are considering here only the non-compact part of the moments. Near the singularities the function $F^{(d)}(\mathbf{x})$ admits a singular expansion of the type (178). In practice, the various coefficients $f_{p,q}^{(\epsilon)}$ are computed by specializing the general expressions of the non-linear retarded potentials $V$, $V_a$, $\hat{W}_{ab}$, etc. (valid for general extended sources) to point particles in d dimensions. On the other hand, the analogue of Eq. (198) in 3 dimensions is

$$ I_{ij} = \mathrm{Pf} \int d^3\mathbf{x} F(\mathbf{x}), \quad (199) $$

where Pf refers to the Hadamard partie finie defined in Eq. (162). The difference $\mathcal{D}I_{ij}$ between the DR evaluation of the d-dimensional integral (198) and its corresponding three-dimensional evaluation (199), reads then

$$ \mathcal{D}I_{ij} = I_{ij}^{(d)} - I_{ij}. \quad (200) $$
Such 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$. We find that Eq. (200) depends on two constant scales $s_1$ and $s_2$ coming from Hadamard’s partie finie (162), and on the constants belonging to dimensional regularization, i.e., $\epsilon = d-3$ and $\ell_0$ defined by Eq. (174). The dimensional regularization of the 3PN quadrupole moment is then obtained as the sum of the pHS part, and of the difference computed according to Eq. (200), namely

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

An important fact, hidden in our too-compact notation (201), is that the sum of the two terms in the right-hand side of Eq. (201) does not depend on the Hadamard regularization scales $s_1$ and $s_2$. Therefore it is possible without changing the sum to re-express these two terms (separately) by means of the constants $r'_1$ and $r'_2$ instead of $s_1$ and $s_2$, where $r'_1, r'_2$ are the two fiducial scales entering the Hadamard-regularization result (197). This replacement being made the pHS term in Eq. (201) is exactly the same as the one in Eq. (197). At this stage all elements are in place to prove the following theorem [62, 63].

**Theorem 10.** The DR quadrupole moment (201) is physically equivalent to the Hadamard-regularized one (end result of Refs. [81, 80]), in the sense that

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

where $\delta_{\xi}I_{ij}$ denotes the effect of the same shifts as determined in Theorems 8 and 9, if and only if the HR ambiguity parameters $\xi, \kappa$ and $\zeta$ take the unique values reported in Eqs. (173). Moreover, the poles $1/\epsilon$ separately present in the two terms in the brackets of Eq. (202) cancel out, so that the physical (“dressed”) DR quadrupole moment is finite and given by the limit when $\epsilon \to 0$ as shown in Eq. (202).

This theorem finally provides an unambiguous determination of the 3PN radiation field by dimensional regularization. Furthermore, as reviewed in Section 6.2, several checks of this calculation could be done, which provide independent confirmations for the ambiguity parameters. Such checks also show the powerfulness of dimensional regularization and its validity for describing the classical general-relativistic dynamics of compact bodies.

# 7 Newtonian-like Equations of Motion

## 7.1 The 3PN acceleration and energy for particles

We present the acceleration of one of the particles, say the particle 1, at the 3PN order, as well as the 3PN energy of the binary, which is conserved in the absence of radiation reaction. To get this result we used essentially a “direct” post-Newtonian method (issued from Ref. [76]), which consists of reducing the 3PN metric of an extended regular source, worked out in Eqs. (144), to the case where the matter tensor is made of delta functions, and then curing the self-field divergences by means of the Hadamard regularization technique. The equations of motion are simply the 3PN geodesic equations explicitly provided in Eqs. (150)–(152); the metric therein is the regularized metric generated by the system of particles itself. Hadamard’s regularization permits to compute all the terms but one, and the Hadamard ambiguity parameter $\lambda$ is obtained from dimensional regularization; see Section 6.3. We also add the 3.5PN terms in harmonic coordinates which are known from Refs. [258, 259, 260, 336, 278, 322, 254]. These correspond to radiation reaction effects at relative 1PN order (see Section 5.4 for discussion on radiation reaction up to 1.5PN order).
Though the successive post-Newtonian approximations are really a consequence of general relativity, the final equations of motion must be interpreted in a Newtonian-like fashion. That is, once a convenient general-relativistic (Cartesian) coordinate system is chosen, we should express the results in terms of the coordinate positions, velocities, and accelerations of the bodies, and view the trajectories of the particles as taking place in the absolute Euclidean space of Newton. But because the equations of motion are actually relativistic, they must:

1. Stay manifestly invariant – at least in harmonic coordinates – when we perform a global post-Newtonian-expanded Lorentz transformation;

2. Possess the correct “perturbative” limit, given by the geodesics of the (post-Newtonian- expanded) Schwarzschild metric, when one of the masses tends to zero;

3. Be conservative, i.e., to admit a Lagrangian or Hamiltonian formulation, when the gravita- tional radiation reaction is turned off.

We denote by $r_{12} = |\mathbf{y}_1(t) - \mathbf{y}_2(t)|$ the harmonic-coordinate distance between the two particles, with $\mathbf{y}_1 = (y_1^i)$ and $\mathbf{y}_2 = (y_2^i)$, by $\mathbf{n}_{12} = (\mathbf{y}_1 - \mathbf{y}_2)/r_{12}$ the corresponding unit direction, and by $\mathbf{v}_1 = d\mathbf{y}_1/dt$ and $\mathbf{a}_1 = d\mathbf{v}_1/dt$ the coordinate velocity and acceleration of the particle 1 (and idem for 2). Sometimes we pose $\mathbf{v}_{12} = \mathbf{v}_1 - \mathbf{v}_2$ for the relative velocity. The usual Euclidean scalar product of vectors is denoted with parentheses, e.g., $(n_{12}v_1) = \mathbf{n}_{12} \cdot \mathbf{v}_1$ and $(v_1v_2) = \mathbf{v}_1 \cdot \mathbf{v}_2$. The equations of the body 2 are obtained by exchanging all the particle labels 1 ↔ 2 (remembering that $\mathbf{n}_{12}$ and $\mathbf{v}_{12}$ change sign in this operation):

$$
\begin{equation}
\begin{aligned}
\mathbf{a}_1 = & -\frac{Gm_2}{r_{12}^2} \mathbf{n}_{12} \\
& + \frac{1}{c^2} \left\{ \left[ \frac{5G^2 m_1 m_2}{r_{12}^3} + \frac{4G^2 m_2^2}{r_{12}^3} + \frac{Gm_2}{r_{12}^2} \left( \frac{3}{2}(n_{12}v_2)^2 - v_1^2 + 4(v_1v_2) - 2v_2^2 \right) \right] \mathbf{n}_{12} \right. \\
& \qquad \left. + \frac{Gm_2}{r_{12}^2} (4(n_{12}v_1) - 3(n_{12}v_2)) \mathbf{v}_{12} \right\} \\
& + \frac{1}{c^4} \left\{ \left[ -\frac{57G^3 m_1^2 m_2}{4r_{12}^4} - \frac{69G^3 m_1 m_2^2}{2r_{12}^4} - \frac{9G^3 m_2^3}{r_{12}^4} \right. \right. \\
& \qquad \left. + \frac{Gm_2}{r_{12}^2} \left( -\frac{15}{8}(n_{12}v_2)^4 + \frac{3}{2}(n_{12}v_2)^2 v_1^2 - 6(n_{12}v_2)^2 (v_1v_2) - 2(v_1v_2)^2 + \frac{9}{2}(n_{12}v_2)^2 v_2^2 \right. \right. \\
& \qquad \left. + 4(v_1v_2)v_2^2 - 2v_2^4 \right) \right. \\
& \qquad \left. + \frac{G^2 m_1 m_2}{r_{12}^3} \left( \frac{39}{2}(n_{12}v_1)^2 - 39(n_{12}v_1)(n_{12}v_2) + \frac{17}{2}(n_{12}v_2)^2 - \frac{15}{4}v_1^2 - \frac{5}{2}(v_1v_2) + \frac{5}{4}v_2^2 \right) \right. \\
& \qquad \left. + \frac{G^2 m_2^2}{r_{12}^3} (2(n_{12}v_1)^2 - 4(n_{12}v_1)(n_{12}v_2) - 6(n_{12}v_2)^2 - 8(v_1v_2) + 4v_2^2) \right] \mathbf{n}_{12} \\
& + \left[ \frac{G^2 m_2^2}{r_{12}^3} (-2(n_{12}v_1) - 2(n_{12}v_2)) + \frac{G^2 m_1 m_2}{r_{12}^3} \left( -\frac{63}{4}(n_{12}v_1) + \frac{55}{4}(n_{12}v_2) \right) \right. \\
& \qquad \left. + \frac{Gm_2}{r_{12}^2} \left( -6(n_{12}v_1)(n_{12}v_2)^2 + \frac{9}{2}(n_{12}v_2)^3 + (n_{12}v_2)v_1^2 - 4(n_{12}v_1)(v_1v_2) \right. \right. \\
& \qquad \left. + 4(n_{12}v_2)(v_1v_2) + 4(n_{12}v_1)v_2^2 - 5(n_{12}v_2)v_2^2 \right) \right] \mathbf{v}_{12} \right\}
\end{aligned}
\end{equation}
$$
$$
\begin{aligned}
&+ \frac{1}{c^5} \left\{ \left[ \frac{208 G^3 m_1 m_2^2}{15 r_{12}^4} (n_{12} v_{12}) - \frac{24 G^3 m_1^2 m_2}{5 r_{12}^4} (n_{12} v_{12}) + \frac{12 G^2 m_1 m_2}{5 r_{12}^3} (n_{12} v_{12}) v_{12}^2 \right] \mathbf{n}_{12} \right. \\
&\qquad \left. + \left[ \frac{8 G^3 m_1^2 m_2}{5 r_{12}^4} - \frac{32 G^3 m_1 m_2^2}{5 r_{12}^4} - \frac{4 G^2 m_1 m_2}{5 r_{12}^3} v_{12}^2 \right] \mathbf{v}_{12} \right\} \\
&+ \frac{1}{c^6} \left\{ \left[ \frac{G m_2}{r_{12}^2} \left( \frac{35}{16} (n_{12} v_2)^6 - \frac{15}{8} (n_{12} v_2)^4 v_1^2 + \frac{15}{2} (n_{12} v_2)^4 (v_1 v_2) + 3 (n_{12} v_2)^2 (v_1 v_2)^2 \right. \right. \\
&\qquad \left. - \frac{15}{2} (n_{12} v_2)^4 v_2^2 + \frac{3}{2} (n_{12} v_2)^2 v_1^2 v_2^2 - 12 (n_{12} v_2)^2 (v_1 v_2) v_2^2 - 2 (v_1 v_2)^2 v_2^2 \right. \\
&\qquad \left. + \frac{15}{2} (n_{12} v_2)^2 v_2^4 + 4 (v_1 v_2) v_2^4 - 2 v_2^6 \right) \right. \\
&\qquad \left. + \frac{G^2 m_1 m_2}{r_{12}^3} \left( - \frac{171}{8} (n_{12} v_1)^4 + \frac{171}{2} (n_{12} v_1)^3 (n_{12} v_2) - \frac{723}{4} (n_{12} v_1)^2 (n_{12} v_2)^2 \right. \right. \\
&\qquad \left. + \frac{383}{2} (n_{12} v_1) (n_{12} v_2)^3 - \frac{455}{8} (n_{12} v_2)^4 + \frac{229}{4} (n_{12} v_1)^2 v_1^2 \right. \\
&\qquad \left. - \frac{205}{2} (n_{12} v_1) (n_{12} v_2) v_1^2 + \frac{191}{4} (n_{12} v_2)^2 v_1^2 - \frac{91}{8} v_1^4 - \frac{229}{2} (n_{12} v_1)^2 (v_1 v_2) \right. \\
&\qquad \left. + 244 (n_{12} v_1) (n_{12} v_2) (v_1 v_2) - \frac{225}{2} (n_{12} v_2)^2 (v_1 v_2) + \frac{91}{2} v_1^2 (v_1 v_2) \right. \\
&\qquad \left. - \frac{177}{4} (v_1 v_2)^2 + \frac{229}{4} (n_{12} v_1)^2 v_2^2 - \frac{283}{2} (n_{12} v_1) (n_{12} v_2) v_2^2 \right. \\
&\qquad \left. + \frac{259}{4} (n_{12} v_2)^2 v_2^2 - \frac{91}{4} v_1^2 v_2^2 + 43 (v_1 v_2) v_2^2 - \frac{81}{8} v_2^4 \right) \right. \\
&\qquad \left. + \frac{G^2 m_2^2}{r_{12}^3} \left( - 6 (n_{12} v_1)^2 (n_{12} v_2)^2 + 12 (n_{12} v_1) (n_{12} v_2)^3 + 6 (n_{12} v_2)^4 \right. \right. \\
&\qquad \left. + 4 (n_{12} v_1) (n_{12} v_2) (v_1 v_2) + 12 (n_{12} v_2)^2 (v_1 v_2) + 4 (v_1 v_2)^2 \right. \\
&\qquad \left. - 4 (n_{12} v_1) (n_{12} v_2) v_2^2 - 12 (n_{12} v_2)^2 v_2^2 - 8 (v_1 v_2) v_2^2 + 4 v_2^4 \right) \right. \\
&\qquad \left. + \frac{G^3 m_2^3}{r_{12}^4} \left( -(n_{12} v_1)^2 + 2 (n_{12} v_1) (n_{12} v_2) + \frac{43}{2} (n_{12} v_2)^2 + 18 (v_1 v_2) - 9 v_2^2 \right) \right. \\
&\qquad \left. + \frac{G^3 m_1 m_2^2}{r_{12}^4} \left( \frac{415}{8} (n_{12} v_1)^2 - \frac{375}{4} (n_{12} v_1) (n_{12} v_2) + \frac{1113}{8} (n_{12} v_2)^2 - \frac{615}{64} (n_{12} v_{12})^2 \pi^2 \right. \right. \\
&\qquad \left. + 18 v_1^2 + \frac{123}{64} \pi^2 v_{12}^2 + 33 (v_1 v_2) - \frac{33}{2} v_2^2 \right) \right. \\
&\qquad \left. + \frac{G^3 m_1^2 m_2}{r_{12}^4} \left( - \frac{45887}{168} (n_{12} v_1)^2 + \frac{24025}{42} (n_{12} v_1) (n_{12} v_2) - \frac{10469}{42} (n_{12} v_2)^2 + \frac{48197}{840} v_1^2 \right. \right. \\
&\qquad \left. - \frac{36227}{420} (v_1 v_2) + \frac{36227}{840} v_2^2 + 110 (n_{12} v_{12})^2 \ln \left( \frac{r_{12}}{r_1'} \right) - 22 v_{12}^2 \ln \left( \frac{r_{12}}{r_1'} \right) \right) \right. \\
&\qquad \left. + \frac{16 G^4 m_2^4}{r_{12}^5} + \frac{G^4 m_1^2 m_2^2}{r_{12}^5} \left( 175 - \frac{41}{16} \pi^2 \right) + \frac{G^4 m_1^3 m_2}{r_{12}^5} \left( - \frac{3187}{1260} + \frac{44}{3} \ln \left( \frac{r_{12}}{r_1'} \right) \right) \right. \\
&\qquad \left. + \frac{G^4 m_1 m_2^3}{r_{12}^5} \left( \frac{110741}{630} - \frac{41}{16} \pi^2 - \frac{44}{3} \ln \left( \frac{r_{12}}{r_2'} \right) \right) \right] \mathbf{n}_{12} \right. \\
&\qquad \left. + \left[ \frac{G m_2}{r_{12}^2} \left( \frac{15}{2} (n_{12} v_1) (n_{12} v_2)^4 - \frac{45}{8} (n_{12} v_2)^5 - \frac{3}{2} (n_{12} v_2)^3 v_1^2 + 6 (n_{12} v_1) (n_{12} v_2)^2 (v_1 v_2) \right] \right\}
\end{aligned}
$$
$$
\begin{aligned}
& -6(n_{12}v_2)^3(v_1v_2) - 2(n_{12}v_2)(v_1v_2)^2 - 12(n_{12}v_1)(n_{12}v_2)^2v_2^2 + 12(n_{12}v_2)^3v_2^2 \\
& \quad + (n_{12}v_2)v_1^2v_2^2 - 4(n_{12}v_1)(v_1v_2)v_2^2 + 8(n_{12}v_2)(v_1v_2)v_2^2 + 4(n_{12}v_1)v_2^4 \\
& \quad - 7(n_{12}v_2)v_2^4) \\
& + \frac{G^2m_2^2}{r_{12}^3} \left( -2(n_{12}v_1)^2(n_{12}v_2) + 8(n_{12}v_1)(n_{12}v_2)^2 + 2(n_{12}v_2)^3 + 2(n_{12}v_1)(v_1v_2) \right. \\
& \quad \left. \quad + 4(n_{12}v_2)(v_1v_2) - 2(n_{12}v_1)v_2^2 - 4(n_{12}v_2)v_2^2 \right) \\
& + \frac{G^2m_1m_2}{r_{12}^3} \left( -\frac{243}{4}(n_{12}v_1)^3 + \frac{565}{4}(n_{12}v_1)^2(n_{12}v_2) - \frac{269}{4}(n_{12}v_1)(n_{12}v_2)^2 \right. \\
& \quad \left. \quad - \frac{95}{12}(n_{12}v_2)^3 + \frac{207}{8}(n_{12}v_1)v_1^2 - \frac{137}{8}(n_{12}v_2)v_1^2 - 36(n_{12}v_1)(v_1v_2) \right. \\
& \quad \left. \quad + \frac{27}{4}(n_{12}v_2)(v_1v_2) + \frac{81}{8}(n_{12}v_1)v_2^2 + \frac{83}{8}(n_{12}v_2)v_2^2 \right) \\
& + \frac{G^3m_2^3}{r_{12}^4} (4(n_{12}v_1) + 5(n_{12}v_2)) \\
& + \frac{G^3m_1m_2^2}{r_{12}^4} \left( -\frac{307}{8}(n_{12}v_1) + \frac{479}{8}(n_{12}v_2) + \frac{123}{32}(n_{12}v_{12})\pi^2 \right) \\
& + \frac{G^3m_1^2m_2}{r_{12}^4} \left( \frac{31397}{420}(n_{12}v_1) - \frac{36227}{420}(n_{12}v_2) - 44(n_{12}v_{12})\ln\left(\frac{r_{12}}{r_1'}\right) \right) v_{12} \\
& + \frac{1}{c^7} \left\{ \left[ \frac{G^4m_3^3m_2}{r_{12}^5} \left( \frac{3992}{105}(n_{12}v_1) - \frac{4328}{105}(n_{12}v_2) \right) \right. \right. \\
& \quad \left. \quad + \frac{G^4m_1^2m_2^2}{r_{12}^6} \left( -\frac{13576}{105}(n_{12}v_1) + \frac{2872}{21}(n_{12}v_2) \right) - \frac{3172}{21} \frac{G^4m_1m_2^3}{r_{12}^6} (n_{12}v_{12}) \right. \\
& \quad \left. \quad + \frac{G^3m_1^2m_2}{r_{12}^4} \left( 48(n_{12}v_1)^3 - \frac{696}{5}(n_{12}v_1)^2(n_{12}v_2) + \frac{744}{5}(n_{12}v_1)(n_{12}v_2)^2 - \frac{288}{5}(n_{12}v_2)^3 \right. \right. \\
& \quad \left. \quad - \frac{4888}{105}(n_{12}v_1)v_1^2 + \frac{5056}{105}(n_{12}v_2)v_1^2 + \frac{2056}{21}(n_{12}v_1)(v_1v_2) \right. \\
& \quad \left. \quad - \frac{2224}{21}(n_{12}v_2)(v_1v_2) - \frac{1028}{21}(n_{12}v_1)v_2^2 + \frac{5812}{105}(n_{12}v_2)v_2^2 \right) \\
& \quad + \frac{G^3m_1m_2^2}{r_{12}^4} \left( -\frac{582}{5}(n_{12}v_1)^3 + \frac{1746}{5}(n_{12}v_1)^2(n_{12}v_2) - \frac{1954}{5}(n_{12}v_1)(n_{12}v_2)^2 \right. \\
& \quad \left. \quad + 158(n_{12}v_2)^3 + \frac{3568}{105}(n_{12}v_{12})v_1^2 - \frac{2864}{35}(n_{12}v_1)(v_1v_2) \right. \\
& \quad \left. \quad + \frac{10048}{105}(n_{12}v_2)(v_1v_2) + \frac{1432}{35}(n_{12}v_1)v_2^2 - \frac{5752}{105}(n_{12}v_2)v_2^2 \right) \\
& \quad + \frac{G^2m_1m_2}{r_{12}^3} \left( -56(n_{12}v_{12})^5 + 60(n_{12}v_1)^3v_{12}^2 - 180(n_{12}v_1)^2(n_{12}v_2)v_{12}^2 \right. \\
& \quad \left. \quad + 174(n_{12}v_1)(n_{12}v_2)^2v_{12}^2 - 54(n_{12}v_2)^3v_{12}^2 - \frac{246}{35}(n_{12}v_{12})v_1^4 \right. \\
& \quad \left. \quad + \frac{1068}{35}(n_{12}v_1)v_1^2(v_1v_2) - \frac{984}{35}(n_{12}v_2)v_1^2(v_1v_2) - \frac{1068}{35}(n_{12}v_1)(v_1v_2)^2 \right)
\end{aligned}
$$
$$
\begin{equation}
\begin{aligned}
& + \frac{180}{7}(n_{12}v_2)(v_1v_2)^2 - \frac{534}{35}(n_{12}v_1)v_1^2v_2^2 + \frac{90}{7}(n_{12}v_2)v_1^2v_2^2 \\
& + \frac{984}{35}(n_{12}v_1)(v_1v_2)v_2^2 - \frac{732}{35}(n_{12}v_2)(v_1v_2)v_2^2 - \frac{204}{35}(n_{12}v_1)v_2^4 \\
& + \frac{24}{7}(n_{12}v_2)v_2^4 \biggr] \mathbf{n}_{12} \\
& + \left[ -\frac{184}{21}\frac{G^4m_1^3m_2}{r_{12}^5} + \frac{6224}{105}\frac{G^4m_1^2m_2^2}{r_{12}^6} + \frac{6388}{105}\frac{G^4m_1m_2^3}{r_{12}^6} \right. \\
& \qquad \left. + \frac{G^3m_1^2m_2}{r_{12}^4}\left(\frac{52}{15}(n_{12}v_1)^2 - \frac{56}{15}(n_{12}v_1)(n_{12}v_2) - \frac{44}{15}(n_{12}v_2)^2 - \frac{132}{35}v_1^2 + \frac{152}{35}(v_1v_2)\right. \right. \\
& \qquad \left. - \frac{48}{35}v_2^2 \right) \\
& \qquad + \frac{G^3m_1m_2^2}{r_{12}^4}\left(\frac{454}{15}(n_{12}v_1)^2 - \frac{372}{5}(n_{12}v_1)(n_{12}v_2) + \frac{854}{15}(n_{12}v_2)^2 - \frac{152}{21}v_1^2\right. \\
& \qquad \left. + \frac{2864}{105}(v_1v_2) - \frac{1768}{105}v_2^2\right) \\
& \qquad + \frac{G^2m_1m_2}{r_{12}^3}\left(60(n_{12}v_1)^4 - \frac{348}{5}(n_{12}v_1)^2v_{12}^2 + \frac{684}{5}(n_{12}v_1)(n_{12}v_2)v_{12}^2\right. \\
& \qquad \left. - 66(n_{12}v_2)^2v_{12}^2 + \frac{334}{35}v_1^4 - \frac{1336}{35}v_1^2(v_1v_2) + \frac{1308}{35}(v_1v_2)^2 + \frac{654}{35}v_1^2v_2^2\right. \\
& \qquad \left. - \frac{1252}{35}(v_1v_2)v_2^2 + \frac{292}{35}v_2^4\right) \biggr] v_{12} \Biggr\} + \mathcal{O}\left(\frac{1}{c^8}\right). \tag{203}
\end{aligned}
\end{equation}
$$

The 2.5PN and 3.5PN terms are associated with gravitational radiation reaction.<fnref n="49" /> The 3PN harmonic-coordinates equations of motion depend on two arbitrary length scales $r'_1$ and $r'_2$ associated with the logarithms present at the 3PN order. It has been proved in Ref. [71] that $r'_1$ and $r'_2$ are merely linked with the choice of coordinates – we can refer to $r'_1$ and $r'_2$ as “gauge constants”. In our approach [69, 71], the harmonic coordinate system is not uniquely fixed by the coordinate condition $\partial_\mu h^{\alpha\mu} = 0$. In fact there are infinitely many “locally-defined” harmonic coordinate systems. For general smooth matter sources, as in the general formalism of Part A, we expect the existence and uniqueness of a global harmonic coordinate system. But here we have some point-particles, with delta-function singularities, and in this case we do not have the notion of a global coordinate system. We can always change the harmonic coordinates by means of the gauge vector $\eta^\alpha = \delta x^\alpha$, satisfying $\Delta\eta^\alpha = 0$ except at the location of the two particles (we assume that the transformation is at the 3PN level, so we can consider simply a flat-space Laplace equation). More precisely, we can show that the logarithms appearing in Eq. (203), together with the constants $r'_1$ and $r'_2$ therein, can be removed by the coordinate transformation associated with the 3PN gauge vector (with $r_1 = |\mathbf{x} - \mathbf{y}_1(t)|$ and $r_2 = |\mathbf{x} - \mathbf{y}_2(t)|$; and $\partial^\alpha = \eta^{\alpha\mu}\partial_\mu$):

$$ \eta^{\alpha} = -\frac{22}{3} \frac{G^2 m_1 m_2}{c^6} \partial^{\alpha} \left[ \frac{G m_1}{r_2} \ln \left( \frac{r_{12}}{r'_1} \right) + \frac{G m_2}{r_1} \ln \left( \frac{r_{12}}{r'_2} \right) \right]. \quad (204) $$

Therefore, the arbitrariness in the choice of the constants $r'_1$ and $r'_2$ is innocuous on the physical

$$ \frac{1}{\pi} \int \frac{d^3 \mathbf{x}}{r_1^2 r_2^2} = \frac{\pi^2}{r_{12}}. $$

<footnote n="49">⁴⁹ Notice the dependence upon the irrational number $\pi^2$. Technically, the $\pi^2$ terms arise from non-linear interactions involving some integrals such as</footnote>
point of view, because the physical results must be gauge invariant. Indeed we shall verify that $r'_1$ and $r'_2$ cancel out in our final results.

When retaining the “even” relativistic corrections at the 1PN, 2PN and 3PN orders, and neglecting the “odd” radiation reaction terms at the 2.5PN and 3.5PN orders, we find that the equations of motion admit a conserved energy (and a Lagrangian, as we shall see); that energy can be straightforwardly obtained by guess-work starting from Eq. (203), with the result

$$
\begin{aligned}
E = & \frac{m_1 v_1^2}{2} - \frac{G m_1 m_2}{2 r_{12}} \\
& + \frac{1}{c^2} \left\{ \frac{G^2 m_1^2 m_2}{2 r_{12}^2} + \frac{3 m_1 v_1^4}{8} + \frac{G m_1 m_2}{r_{12}} \left( -\frac{1}{4}(n_{12} v_1)(n_{12} v_2) + \frac{3}{2} v_1^2 - \frac{7}{4}(v_1 v_2) \right) \right\} \\
& + \frac{1}{c^4} \left\{ -\frac{G^3 m_1^3 m_2}{2 r_{12}^3} - \frac{19 G^3 m_1^2 m_2^2}{8 r_{12}^3} + \frac{5 m_1 v_1^6}{16} \right. \\
& \qquad \left. + \frac{G m_1 m_2}{r_{12}} \left( \frac{3}{8}(n_{12} v_1)^3 (n_{12} v_2) + \frac{3}{16}(n_{12} v_1)^2 (n_{12} v_2)^2 - \frac{9}{8}(n_{12} v_1)(n_{12} v_2) v_1^2 \right. \right. \\
& \qquad \left. - \frac{13}{8}(n_{12} v_2)^2 v_1^2 + \frac{21}{8} v_1^4 + \frac{13}{8}(n_{12} v_1)^2 (v_1 v_2) + \frac{3}{4}(n_{12} v_1)(n_{12} v_2)(v_1 v_2) \right. \\
& \qquad \left. - \frac{55}{8} v_1^2 (v_1 v_2) + \frac{17}{8}(v_1 v_2)^2 + \frac{31}{16} v_1^2 v_2^2 \right) \\
& + \frac{G^2 m_1^2 m_2}{r_{12}^2} \left( \frac{29}{4}(n_{12} v_1)^2 - \frac{13}{4}(n_{12} v_1)(n_{12} v_2) + \frac{1}{2}(n_{12} v_2)^2 - \frac{3}{2} v_1^2 + \frac{7}{4} v_2^2 \right) \\
& \left. + \frac{1}{c^6} \left\{ \frac{35 m_1 v_1^8}{128} \right. \right. \\
& \qquad \left. + \frac{G m_1 m_2}{r_{12}} \left( -\frac{5}{16}(n_{12} v_1)^5 (n_{12} v_2) - \frac{5}{16}(n_{12} v_1)^4 (n_{12} v_2)^2 - \frac{5}{32}(n_{12} v_1)^3 (n_{12} v_2)^3 \right. \right. \\
& \qquad \left. + \frac{19}{16}(n_{12} v_1)^3 (n_{12} v_2) v_1^2 + \frac{15}{16}(n_{12} v_1)^2 (n_{12} v_2)^2 v_1^2 + \frac{3}{4}(n_{12} v_1)(n_{12} v_2)^3 v_1^2 \right. \\
& \qquad \left. + \frac{19}{16}(n_{12} v_2)^4 v_1^2 - \frac{21}{16}(n_{12} v_1)(n_{12} v_2) v_1^4 - 2(n_{12} v_2)^2 v_1^4 \right. \\
& \qquad \left. + \frac{55}{16} v_1^6 - \frac{19}{16}(n_{12} v_1)^4 (v_1 v_2) - (n_{12} v_1)^3 (n_{12} v_2)(v_1 v_2) \right. \\
& \qquad \left. - \frac{15}{32}(n_{12} v_1)^2 (n_{12} v_2)^2 (v_1 v_2) + \frac{45}{16}(n_{12} v_1)^2 v_1^2 (v_1 v_2) \right. \\
& \qquad \left. + \frac{5}{4}(n_{12} v_1)(n_{12} v_2) v_1^2 (v_1 v_2) + \frac{11}{4}(n_{12} v_2)^2 v_1^2 (v_1 v_2) - \frac{139}{16} v_1^4 (v_1 v_2) \right. \\
& \qquad \left. - \frac{3}{4}(n_{12} v_1)^2 (v_1 v_2)^2 + \frac{5}{16}(n_{12} v_1)(n_{12} v_2)(v_1 v_2)^2 + \frac{41}{8} v_1^2 (v_1 v_2)^2 + \frac{1}{16}(v_1 v_2)^3 \right. \\
& \qquad \left. - \frac{45}{16}(n_{12} v_1)^2 v_1^2 v_2^2 - \frac{23}{32}(n_{12} v_1)(n_{12} v_2) v_1^2 v_2^2 + \frac{79}{16} v_1^4 v_2^2 - \frac{161}{32} v_1^2 (v_1 v_2) v_2^2 \right) \\
& + \frac{G^2 m_1^2 m_2}{r_{12}^2} \left( -\frac{49}{8}(n_{12} v_1)^4 + \frac{75}{8}(n_{12} v_1)^3 (n_{12} v_2) - \frac{187}{8}(n_{12} v_1)^2 (n_{12} v_2)^2 \right. \\
& \qquad \left. + \frac{247}{24}(n_{12} v_1)(n_{12} v_2)^3 + \frac{49}{8}(n_{12} v_1)^2 v_1^2 + \frac{81}{8}(n_{12} v_1)(n_{12} v_2) v_1^2 \right. \\
& \qquad \left. - \frac{21}{4}(n_{12} v_2)^2 v_1^2 + \frac{11}{2} v_1^4 - \frac{15}{2}(n_{12} v_1)^2 (v_1 v_2) - \frac{3}{2}(n_{12} v_1)(n_{12} v_2)(v_1 v_2) \right)
\end{aligned}
$$
$$
\begin{equation}
\begin{aligned}
& + \frac{21}{4}(n_{12}v_2)^2(v_1v_2) - 27v_1^2(v_1v_2) + \frac{55}{2}(v_1v_2)^2 + \frac{49}{4}(n_{12}v_1)^2v_2^2 \\
& - \frac{27}{2}(n_{12}v_1)(n_{12}v_2)v_2^2 + \frac{3}{4}(n_{12}v_2)^2v_2^2 + \frac{55}{4}v_1^2v_2^2 - 28(v_1v_2)v_2^2 + \frac{135}{16}v_2^4) \\
& + \frac{3G^4m_1^4m_2}{8r_{12}^4} + \frac{G^4m_1^3m_2^2}{r_{12}^4}\left(\frac{9707}{420} - \frac{22}{3}\ln\left(\frac{r_{12}}{r_1}\right)\right) \\
& + \frac{G^3m_1^2m_2^2}{r_{12}^3}\left(\frac{547}{12}(n_{12}v_1)^2 - \frac{3115}{48}(n_{12}v_1)(n_{12}v_2) - \frac{123}{64}(n_{12}v_1)(n_{12}v_{12})\pi^2 - \frac{575}{18}v_1^2\right. \\
& \qquad \left. + \frac{41}{64}\pi^2(v_1v_{12}) + \frac{4429}{144}(v_1v_2)\right) \\
& + \frac{G^3m_1^3m_2}{r_{12}^3}\left(-\frac{44627}{840}(n_{12}v_1)^2 + \frac{32027}{840}(n_{12}v_1)(n_{12}v_2) + \frac{3}{2}(n_{12}v_2)^2 + \frac{24187}{2520}v_1^2\right. \\
& \qquad \left. - \frac{27967}{2520}(v_1v_2) + \frac{5}{4}v_2^2 + 22(n_{12}v_1)(n_{12}v_{12})\ln\left(\frac{r_{12}}{r_1}\right) - \frac{22}{3}(v_1v_{12})\ln\left(\frac{r_{12}}{r_1}\right)\right)\right\} \\
& + 1 \leftrightarrow 2 + \mathcal{O}\left(\frac{1}{c^7}\right).
\end{aligned}
\tag{205}
\end{equation}
$$

To the terms given above, we must add the same terms but corresponding to the relabelling $1 \leftrightarrow 2$. Actually, this energy is not conserved because of the radiation reaction. Thus its time derivative, as computed by means of the 3PN equations of motion themselves (i.e., by order-reducing all the accelerations), is purely equal to the 2.5PN effect,

$$
\begin{equation}
\begin{aligned}
\frac{dE}{dt} ={}& \frac{4}{5} \frac{G^2 m_1^2 m_2}{c^5 r_{12}^3} \left[ (v_1 v_{12}) \left( -v_{12}^2 + 2 \frac{G m_1}{r_{12}} - 8 \frac{G m_2}{r_{12}} \right) + (n_{12} v_1) (n_{12} v_{12}) \left( 3 v_{12}^2 - 6 \frac{G m_1}{r_{12}} + \frac{52}{3} \frac{G m_2}{r_{12}} \right) \right] \\
& + 1 \leftrightarrow 2 + \mathcal{O}\left(\frac{1}{c^7}\right).
\end{aligned}
\tag{206}
\end{equation}
$$

The resulting energy balance equation can be better expressed by transfering to the left-hand side certain 2.5PN terms so that we recognize in the right-hand side the familiar form of a total energy flux. Posing

$$
E^{2.5\text{PN}} = E + \frac{4G^2 m_1^2 m_2}{5c^5 r_{12}^2} (n_{12}v_1) \left[ v_{12}^2 - \frac{2G(m_1 - m_2)}{r_{12}} \right] + 1 \leftrightarrow 2, \quad (207)
$$

we find agreement with the standard Einstein quadrupole formula (4):

$$
\frac{dE^{2.5\text{PN}}}{dt} = - \frac{G}{5c^5} \frac{d^3 Q_{ij}}{dt^3} \frac{d^3 Q_{ij}}{dt^3} + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (208)
$$

where the Newtonian trace-free quadrupole moment reads $Q_{ij} = m_1(y_1^i y_1^j - \frac{1}{3}\delta^{ij}y_1^i) + 1 \leftrightarrow 2$. We refer to [258, 259] for the discussion of the energy balance equation up to the next 3.5PN order. See also Eq. (158) for the energy balance equation at relative 1.5PN order for general fluid systems.

7.2 Lagrangian and Hamiltonian formulations

The conservative part of the equations of motion in harmonic coordinates (203) is derivable from
a *generalized* Lagrangian, depending not only on the positions and velocities of the bodies, but
also on their accelerations: **a**₁ = d**v**₁/dt and **a**₂ = d**v**₂/dt. As shown in Ref. [147], the accelera-
tions in the harmonic-coordinates Lagrangian occur already from the 2PN order. This fact is in
accordance with a general result [308] that N-body equations of motion cannot be derived from an ordinary Lagrangian beyond the 1PN level, provided that the gauge conditions preserve the manifest Lorentz invariance. Note that we can always arrange for the dependence of the Lagrangian upon the accelerations to be linear, at the price of adding some so-called “multi-zero” terms to the Lagrangian, which do not modify the equations of motion (see, e.g., Ref. [169]). At the 3PN level, we find that the Lagrangian also depends on accelerations. It is notable that these accelerations are sufficient – there is no need to include derivatives of accelerations. Note also that the Lagrangian is not unique because we can always add to it a total time derivative $dF/dt$, where F is any function depending on the positions and velocities, without changing the dynamics. We find [174]

$$ L^{\text{harm}} = \frac{Gm_1m_2}{2r_{12}} + \frac{m_1v_1^2}{2} + \frac{1}{c^2} \left\{ -\frac{G^2m_1^2m_2}{2r_{12}^2} + \frac{m_1v_1^4}{8} + \frac{Gm_1m_2}{r_{12}} \left( -\frac{1}{4}(n_{12}v_1)(n_{12}v_2) + \frac{3}{2}v_1^2 - \frac{7}{4}(v_1v_2) \right) \right\} \\ + \frac{1}{c^4} \left\{ \frac{G^3m_1^3m_2}{2r_{12}^3} + \frac{19G^3m_1^2m_2^2}{8r_{12}^3} + \frac{G^2m_1^2m_2}{r_{12}^2} \left( \frac{7}{2}(n_{12}v_1)^2 - \frac{7}{2}(n_{12}v_1)(n_{12}v_2) + \frac{1}{2}(n_{12}v_2)^2 + \frac{1}{4}v_1^2 - \frac{7}{4}(v_1v_2) + \frac{7}{4}v_2^2 \right) \right. \\ + \frac{Gm_1m_2}{r_{12}} \left( \frac{3}{16}(n_{12}v_1)^2(n_{12}v_2)^2 - \frac{7}{8}(n_{12}v_2)^2v_1^2 + \frac{7}{8}v_1^4 + \frac{3}{4}(n_{12}v_1)(n_{12}v_2)(v_1v_2) \right. \\ \left. - 2v_1^2(v_1v_2) + \frac{1}{8}(v_1v_2)^2 + \frac{15}{16}v_1^2v_2^2 \right) + \frac{m_1v_1^6}{16} \\ + Gm_1m_2 \left( -\frac{7}{4}(a_1v_2)(n_{12}v_2) - \frac{1}{8}(n_{12}a_1)(n_{12}v_2)^2 + \frac{7}{8}(n_{12}a_1)v_2^2 \right) \\ + \frac{1}{c^6} \left\{ \frac{G^2m_1^2m_2}{r_{12}^2} \left( \frac{13}{18}(n_{12}v_1)^4 + \frac{83}{18}(n_{12}v_1)^3(n_{12}v_2) - \frac{35}{6}(n_{12}v_1)^2(n_{12}v_2)^2 - \frac{245}{24}(n_{12}v_1)^2v_1^2 \right. \right. \\ + \frac{179}{12}(n_{12}v_1)(n_{12}v_2)v_1^2 - \frac{235}{24}(n_{12}v_2)^2v_1^2 + \frac{373}{48}v_1^4 + \frac{529}{24}(n_{12}v_1)^2(v_1v_2) \\ - \frac{97}{6}(n_{12}v_1)(n_{12}v_2)(v_1v_2) - \frac{719}{24}v_1^2(v_1v_2) + \frac{463}{24}(v_1v_2)^2 - \frac{7}{24}(n_{12}v_1)^2v_2^2 \\ \left. - \frac{1}{2}(n_{12}v_1)(n_{12}v_2)v_2^2 + \frac{1}{4}(n_{12}v_2)^2v_2^2 + \frac{463}{48}v_1^2v_2^2 - \frac{19}{2}(v_1v_2)v_2^2 + \frac{45}{16}v_2^4 \right) \\ + \frac{5m_1v_1^8}{128} \\ + Gm_1m_2 \left( \frac{3}{8}(a_1v_2)(n_{12}v_1)(n_{12}v_2)^2 + \frac{5}{12}(a_1v_2)(n_{12}v_2)^3 + \frac{1}{8}(n_{12}a_1)(n_{12}v_1)(n_{12}v_2)^3 \right. \\ + \frac{1}{16}(n_{12}a_1)(n_{12}v_2)^4 + \frac{11}{4}(a_1v_1)(n_{12}v_2)v_1^2 - (a_1v_2)(n_{12}v_2)v_1^2 \\ - 2(a_1v_1)(n_{12}v_2)(v_1v_2) + \frac{1}{4}(a_1v_2)(n_{12}v_2)(v_1v_2) \\ + \frac{3}{8}(n_{12}a_1)(n_{12}v_2)^2(v_1v_2) - \frac{5}{8}(n_{12}a_1)(n_{12}v_1)^2v_2^2 + \frac{15}{8}(a_1v_1)(n_{12}v_2)v_2^2 \\ - \frac{15}{8}(a_1v_2)(n_{12}v_2)v_2^2 - \frac{1}{2}(n_{12}a_1)(n_{12}v_1)(n_{12}v_2)v_2^2 \\ \left. - \frac{5}{16}(n_{12}a_1)(n_{12}v_2)^2v_2^2 \right) \right\} $$
$$
\begin{equation}
\begin{aligned}
& + \frac{G^2 m_1^2 m_2}{r_{12}} \left( -\frac{235}{24} (a_2 v_1) (n_{12} v_1) - \frac{29}{24} (n_{12} a_2) (n_{12} v_1)^2 - \frac{235}{24} (a_1 v_2) (n_{12} v_2) \right. \\
& \qquad \left. - \frac{17}{6} (n_{12} a_1) (n_{12} v_2)^2 + \frac{185}{16} (n_{12} a_1) v_1^2 - \frac{235}{48} (n_{12} a_2) v_1^2 \right. \\
& \qquad \left. - \frac{185}{8} (n_{12} a_1) (v_1 v_2) + \frac{20}{3} (n_{12} a_1) v_2^2 \right) \\
& + \frac{G m_1 m_2}{r_{12}} \left( -\frac{5}{32} (n_{12} v_1)^3 (n_{12} v_2)^3 + \frac{1}{8} (n_{12} v_1) (n_{12} v_2)^3 v_1^2 + \frac{5}{8} (n_{12} v_2)^4 v_1^2 \right. \\
& \qquad \left. - \frac{11}{16} (n_{12} v_1) (n_{12} v_2) v_1^4 + \frac{1}{4} (n_{12} v_2)^2 v_1^4 + \frac{11}{16} v_1^6 \right. \\
& \qquad \left. - \frac{15}{32} (n_{12} v_1)^2 (n_{12} v_2)^2 (v_1 v_2) + (n_{12} v_1) (n_{12} v_2) v_1^2 (v_1 v_2) \right. \\
& \qquad \left. + \frac{3}{8} (n_{12} v_2)^2 v_1^2 (v_1 v_2) - \frac{13}{16} v_1^4 (v_1 v_2) + \frac{5}{16} (n_{12} v_1) (n_{12} v_2) (v_1 v_2)^2 \right. \\
& \qquad \left. + \frac{1}{16} (v_1 v_2)^3 - \frac{5}{8} (n_{12} v_1)^2 v_1^2 v_2^2 - \frac{23}{32} (n_{12} v_1) (n_{12} v_2) v_1^2 v_2^2 + \frac{1}{16} v_1^4 v_2^2 \right. \\
& \qquad \left. - \frac{1}{32} v_1^2 (v_1 v_2) v_2^2 \right) \\
& - \frac{3 G^4 m_1^4 m_2}{8 r_{12}^4} + \frac{G^4 m_1^3 m_2^2}{r_{12}^4} \left( -\frac{9707}{420} + \frac{22}{3} \ln \left( \frac{r_{12}}{r_1'} \right) \right) \\
& + \frac{G^3 m_1^2 m_2^2}{r_{12}^3} \left( \frac{383}{24} (n_{12} v_1)^2 - \frac{889}{48} (n_{12} v_1) (n_{12} v_2) - \frac{123}{64} (n_{12} v_1) (n_{12} v_2) \pi^2 - \frac{305}{72} v_1^2 \right. \\
& \qquad \left. + \frac{41}{64} \pi^2 (v_1 v_{12}) + \frac{439}{144} (v_1 v_2) \right) \\
& + \frac{G^3 m_1^3 m_2}{r_{12}^3} \left( -\frac{8243}{210} (n_{12} v_1)^2 + \frac{15541}{420} (n_{12} v_1) (n_{12} v_2) + \frac{3}{2} (n_{12} v_2)^2 + \frac{15611}{1260} v_1^2 \right. \\
& \qquad \left. - \frac{17501}{1260} (v_1 v_2) + \frac{5}{4} v_2^2 + 22 (n_{12} v_1) (n_{12} v_2) \ln \left( \frac{r_{12}}{r_1'} \right) \right. \\
& \qquad \left. - \frac{22}{3} (v_1 v_{12}) \ln \left( \frac{r_{12}}{r_1'} \right) \right) \Biggr\} + 1 \leftrightarrow 2 + \mathcal{O} \left( \frac{1}{c^7} \right) .
\end{aligned}
\tag{209}
\end{equation}
$$

Witness the accelerations occurring at the 2PN and 3PN orders; see also the gauge-dependent logarithms of $r_{12}/r_1'$ and $r_{12}/r_2'$. We refer to [174] for the explicit expressions of the ten conserved quantities corresponding to the integrals of energy [also given in Eq. (205)], linear and angular momenta, and center-of-mass position. Notice that while it is strictly forbidden to replace the accelerations by the equations of motion in the Lagrangian, this can and *should* be done in the final expressions of the conserved integrals derived from that Lagrangian.

Now we want to exhibit a transformation of the particles' dynamical variables – or *contact* transformation, as it is called in the jargon – which transforms the 3PN harmonic-coordinates Lagrangian (209) into a new Lagrangian, valid in some ADM or ADM-like coordinate system, and such that the associated Hamiltonian coincides with the 3PN Hamiltonian that has been obtained by Jaranowski & Schäfer [261, 262]. In ADM coordinates the Lagrangian will be ordinary, depending only on the positions and velocities of the bodies. Let this contact transformation be $\mathbf{Y}_1(t) = \mathbf{y}_1(t) + \delta\mathbf{y}_1(t)$ and $1 \leftrightarrow 2$, where $\mathbf{Y}_1$ and $\mathbf{y}_1$ denote the trajectories in ADM and harmonic coordinates, respectively. For this transformation to be able to remove all the accelerations in the
initial Lagrangian $L^{\text{harm}}$ up to the 3PN order, we determine [174] it to be necessarily of the form

$$ \delta \mathbf{y}_1 = \frac{1}{m_1} \left[ \frac{\partial L^{\text{harm}}}{\partial \boldsymbol{a}_1} + \frac{\partial F}{\partial \boldsymbol{v}_1} + \frac{1}{c^6} \boldsymbol{X}_1 \right] + \mathcal{O} \left( \frac{1}{c^8} \right) \quad (\text{and idem } 1 \leftrightarrow 2), \qquad (210) $$

where $F$ is a freely adjustable function of the positions and velocities, made of 2PN and 3PN terms, and where $\mathbf{X}_1$ represents a special correction term, that is purely of order 3PN. The point is that once the function $F$ is specified there is a unique determination of the correction term $\mathbf{X}_1$ for the contact transformation to work (see Ref. [174] for the details). Thus, the freedom we have is entirely encoded into the function $F$, and the work then consists in showing that there exists a unique choice of $F$ for which our Lagrangian $L^{\text{harm}}$ is physically equivalent, via the contact transformation (210), to the ADM Hamiltonian of Refs. [261, 262]. An interesting point is that not only the transformation must remove all the accelerations in $L^{\text{harm}}$, but it should also cancel out all the logarithms $\ln(r_{12}/r'_1)$ and $\ln(r_{12}/r'_2)$, because there are no logarithms in ADM coordinates. The result we find, which can be checked to be in full agreement with the expression of the gauge vector in Eq. (204), is that $F$ involves the logarithmic terms

$$ F = \frac{22}{3} \frac{G^3 m_1 m_2}{c^6 r_{12}^2} \left[ m_1^2 (n_{12} v_1) \ln \left( \frac{r_{12}}{r'_1} \right) - m_2^2 (n_{12} v_2) \ln \left( \frac{r_{12}}{r'_2} \right) \right] + \dots, \qquad (211) $$

together with many other non-logarithmic terms (indicated by dots) that are entirely specified by the isometry of the harmonic and ADM descriptions of the motion. For this particular choice of $F$ the ADM Lagrangian reads

$$ L^{\text{ADM}} = L^{\text{harm}} + \frac{\delta L^{\text{harm}}}{\delta y_1^i} \delta y_1^i + \frac{\delta L^{\text{harm}}}{\delta y_2^i} \delta y_2^i + \frac{dF}{dt} + \mathcal{O} \left( \frac{1}{c^8} \right). \qquad (212) $$

Inserting into this equation all our explicit expressions we find

$$ \begin{aligned} L^{\text{ADM}} = & \frac{Gm_1 m_2}{2R_{12}} + \frac{1}{2} m_1 V_1^2 \\ & + \frac{1}{c^2} \left\{ -\frac{G^2 m_1^2 m_2}{2R_{12}^2} + \frac{1}{8} m_1 V_1^4 + \frac{Gm_1 m_2}{R_{12}} \left( -\frac{1}{4}(N_{12}V_1)(N_{12}V_2) + \frac{3}{2}V_1^2 - \frac{7}{4}(V_1V_2) \right) \right\} \\ & + \frac{1}{c^4} \left\{ \frac{G^3 m_1^3 m_2}{4R_{12}^3} + \frac{5G^3 m_1^2 m_2^2}{8R_{12}^3} + \frac{m_1 V_1^6}{16} \right. \\ & \qquad \left. + \frac{G^2 m_1^2 m_2}{R_{12}^2} \left( \frac{15}{8}(N_{12}V_1)^2 + \frac{11}{8}V_1^2 - \frac{15}{4}(V_1V_2) + 2V_2^2 \right) \right. \\ & \qquad \left. + \frac{Gm_1 m_2}{R_{12}} \left( \frac{3}{16}(N_{12}V_1)^2 (N_{12}V_2)^2 - \frac{1}{4}(N_{12}V_1)(N_{12}V_2)V_1^2 - \frac{5}{8}(N_{12}V_2)^2 V_1^2 + \frac{7}{8}V_1^4 \right. \right. \\ & \qquad \left. \left. + \frac{3}{4}(N_{12}V_1)(N_{12}V_2)(V_1V_2) - \frac{7}{4}V_1^2(V_1V_2) + \frac{1}{8}(V_1V_2)^2 + \frac{11}{16}V_1^2 V_2^2 \right) \right\} \\ & + \frac{1}{c^6} \left\{ \frac{5m_1 V_1^8}{128} - \frac{G^4 m_1^4 m_2}{8R_{12}^4} + \frac{G^4 m_1^3 m_2^2}{R_{12}^4} \left( -\frac{227}{24} + \frac{21}{32}\pi^2 \right) \right. \\ & \qquad \left. + \frac{Gm_1 m_2}{R_{12}} \left( -\frac{5}{32}(N_{12}V_1)^3 (N_{12}V_2)^3 + \frac{3}{16}(N_{12}V_1)^2 (N_{12}V_2)^2 V_1^2 \right. \right. \\ & \qquad \left. \left. + \frac{9}{16}(N_{12}V_1)(N_{12}V_2)^3 V_1^2 - \frac{3}{16}(N_{12}V_1)(N_{12}V_2)V_1^4 - \frac{5}{16}(N_{12}V_2)^2 V_1^4 \right. \right. \\ & \qquad \left. \left. + \frac{11}{16}V_1^6 - \frac{15}{32}(N_{12}V_1)^2 (N_{12}V_2)^2 (V_1V_2) + \frac{3}{4}(N_{12}V_1)(N_{12}V_2)V_1^2 (V_1V_2) \right) \right\} \end{aligned} $$
$$
\begin{equation}
\begin{aligned}
& -\frac{1}{16}(N_{12}V_2)^2 V_1^2 (V_1 V_2) - \frac{21}{16}V_1^4 (V_1 V_2) + \frac{5}{16}(N_{12}V_1)(N_{12}V_2)(V_1 V_2)^2 \\
& \quad + \frac{1}{8}V_1^2 (V_1 V_2)^2 + \frac{1}{16}(V_1 V_2)^3 - \frac{5}{16}(N_{12}V_1)^2 V_1^2 V_2^2 \\
& \quad - \frac{9}{32}(N_{12}V_1)(N_{12}V_2)V_1^2 V_2^2 + \frac{7}{8}V_1^4 V_2^2 - \frac{15}{32}V_1^2 (V_1 V_2) V_2^2) \\
& + \frac{G^2 m_1^2 m_2}{R_{12}^2} \left( -\frac{5}{12}(N_{12}V_1)^4 - \frac{13}{8}(N_{12}V_1)^3(N_{12}V_2) - \frac{23}{24}(N_{12}V_1)^2 (N_{12}V_2)^2 \right. \\
& \quad \left. + \frac{13}{16}(N_{12}V_1)^2 V_1^2 + \frac{1}{4}(N_{12}V_1)(N_{12}V_2)V_1^2 + \frac{5}{6}(N_{12}V_2)^2 V_1^2 + \frac{21}{16}V_1^4 \right. \\
& \quad \left. - \frac{1}{2}(N_{12}V_1)^2 (V_1 V_2) + \frac{1}{3}(N_{12}V_1)(N_{12}V_2)(V_1 V_2) - \frac{97}{16}V_1^2 (V_1 V_2) \right. \\
& \quad \left. + \frac{341}{48}(V_1 V_2)^2 + \frac{29}{24}(N_{12}V_1)^2 V_2^2 - (N_{12}V_1)(N_{12}V_2)V_2^2 + \frac{43}{12}V_1^2 V_2^2 \right. \\
& \quad \left. - \frac{71}{8}(V_1 V_2) V_2^2 + \frac{47}{16}V_2^4 \right) \\
& + \frac{G^3 m_1^2 m_2^2}{R_{12}^3} \left( \frac{73}{16}(N_{12}V_1)^2 - 11(N_{12}V_1)(N_{12}V_2) + \frac{3}{64}\pi^2(N_{12}V_1)(N_{12}V_2) \right. \\
& \quad \left. - \frac{265}{48}V_1^2 - \frac{1}{64}\pi^2(V_1 V_2) + \frac{59}{8}(V_1 V_2) \right) \\
& + \frac{G^3 m_1^3 m_2}{R_{12}^3} \left( -5(N_{12}V_1)^2 - \frac{1}{8}(N_{12}V_1)(N_{12}V_2) + \frac{173}{48}V_1^2 - \frac{27}{8}(V_1 V_2) + \frac{13}{8}V_2^2 \right) \Biggr\} \\
& + 1 \leftrightarrow 2 + \mathcal{O}\left(\frac{1}{c^7}\right).
\end{aligned}
\tag{213}
\end{equation}
$$

The notation is the same as in Eq. (209), except that we use upper-case letters to denote the ADM-coordinates positions and velocities; thus, for instance $N_{12} = (\mathbf{Y}_1 - \mathbf{Y}_2)/R_{12}$ and $(N_{12}V_1) = N_{12} \cdot \mathbf{V}_1$. The Hamiltonian is simply deduced from the latter Lagrangian by applying the usual Legendre transformation. Posing $\mathbf{P}_1 = \partial L^{\text{ADM}}/\partial \mathbf{V}_1$ and $1 \leftrightarrow 2$, we get [261, 262, 263, 162, 174]

$$
\begin{equation}
\begin{aligned}
H^{\text{ADM}} = & -\frac{Gm_1 m_2}{2R_{12}} + \frac{P_1^2}{2m_1} \\
& + \frac{1}{c^2} \left\{ -\frac{P_1^4}{8m_1^3} + \frac{G^2 m_1^2 m_2}{2R_{12}^2} + \frac{Gm_1 m_2}{R_{12}} \left( \frac{1}{4} \frac{(N_{12}P_1)(N_{12}P_2)}{m_1 m_2} - \frac{3}{2} \frac{P_1^2}{m_1^2} + \frac{7}{4} \frac{(P_1 P_2)}{m_1 m_2} \right) \right\} \\
& + \frac{1}{c^4} \left\{ \frac{P_1^6}{16m_1^5} - \frac{G^3 m_1^3 m_2}{4R_{12}^3} - \frac{5G^3 m_1^2 m_2^2}{8R_{12}^3} \right. \\
& \qquad \left. + \frac{G^2 m_1^2 m_2}{R_{12}^2} \left( -\frac{3}{2} \frac{(N_{12}P_1)(N_{12}P_2)}{m_1 m_2} + \frac{19}{4} \frac{P_1^2}{m_1^2} - \frac{27}{4} \frac{(P_1 P_2)}{m_1 m_2} + \frac{5P_2^2}{2m_1^2} \right) \right. \\
& \qquad \left. + \frac{Gm_1 m_2}{R_{12}} \left( -\frac{3}{16} \frac{(N_{12}P_1)^2 (N_{12}P_2)^2}{m_1^2 m_2^2} + \frac{5}{8} \frac{(N_{12}P_2)^2 P_1^2}{m_1^2 m_2^2} \right. \right. \\
& \qquad \left. \left. + \frac{5}{8} \frac{P_1^4}{m_1^4} - \frac{3}{4} \frac{(N_{12}P_1)(N_{12}P_2)(P_1 P_2)}{m_1^2 m_2^2} - \frac{1}{8} \frac{(P_1 P_2)^2}{m_1^2 m_2^2} - \frac{11}{16} \frac{P_1^2 P_2^2}{m_1^2 m_2^2} \right) \right\} \\
& + \frac{1}{c^6} \left\{ -\frac{5P_1^8}{128m_1^7} + \frac{G^4 m_1^4 m_2}{8R_{12}^4} + \frac{G^4 m_1^3 m_2^2}{R_{12}^4} \left( \frac{227}{24} - \frac{21}{32}\pi^2 \right) \right\}
\end{aligned}
\tag{214}
\end{equation}
$$
$$
\begin{align}
&+ \frac{G^3 m_1^2 m_2^2}{R_{12}^3} \left( -\frac{43 (N_{12} P_1)^2}{16 m_1^2} + \frac{119 (N_{12} P_1) (N_{12} P_2)}{16 m_1 m_2} - \frac{3 \pi^2 (N_{12} P_1)^2}{64 m_1^2} \right. \nonumber \\
&\qquad \left. + \frac{3}{64} \pi^2 \frac{(N_{12} P_1) (N_{12} P_2)}{m_1 m_2} - \frac{473 P_1^2}{48 m_1^2} + \frac{1}{64} \pi^2 \frac{P_1^2}{m_1^2} + \frac{143 (P_1 P_2)}{16 m_1 m_2} \right. \nonumber \\
&\qquad \left. - \frac{1}{64} \pi^2 \frac{(P_1 P_2)}{m_1 m_2} \right) \nonumber \\
&+ \frac{G^3 m_1^3 m_2}{R_{12}^3} \left( \frac{5 (N_{12} P_1)^2}{4 m_1^2} + \frac{21 (N_{12} P_1) (N_{12} P_2)}{8 m_1 m_2} - \frac{425 P_1^2}{48 m_1^2} + \frac{77 (P_1 P_2)}{8 m_1 m_2} - \frac{25 P_2^2}{8 m_2^2} \right) \nonumber \\
&+ \frac{G^2 m_1^2 m_2}{R_{12}^2} \left( \frac{5 (N_{12} P_1)^4}{12 m_1^4} - \frac{3 (N_{12} P_1)^3 (N_{12} P_2)}{2 m_1^3 m_2} + \frac{10 (N_{12} P_1)^2 (N_{12} P_2)^2}{3 m_1^2 m_2^2} \right. \nonumber \\
&\qquad \left. + \frac{17 (N_{12} P_1)^2 P_1^2}{16 m_1^4} - \frac{15 (N_{12} P_1) (N_{12} P_2) P_1^2}{8 m_1^3 m_2} - \frac{55 (N_{12} P_2)^2 P_1^2}{12 m_1^2 m_2^2} \right. \nonumber \\
&\qquad \left. + \frac{P_1^4}{16 m_1^4} - \frac{11 (N_{12} P_1)^2 (P_1 P_2)}{8 m_1^3 m_2} + \frac{125 (N_{12} P_1) (N_{12} P_2) (P_1 P_2)}{12 m_1^2 m_2^2} \right. \nonumber \\
&\qquad \left. - \frac{115 P_1^2 (P_1 P_2)}{16 m_1^3 m_2} + \frac{25 (P_1 P_2)^2}{48 m_1^2 m_2^2} - \frac{193 (N_{12} P_1)^2 P_2^2}{48 m_1^2 m_2^2} + \frac{371 P_1^2 P_2^2}{48 m_1^2 m_2^2} \right. \nonumber \\
&\qquad \left. - \frac{27 P_2^4}{16 m_2^4} \right) \nonumber \\
&+ \frac{G m_1 m_2}{R_{12}} \left( \frac{5 (N_{12} P_1)^3 (N_{12} P_2)^3}{32 m_1^3 m_2^3} + \frac{3 (N_{12} P_1)^2 (N_{12} P_2)^2 P_1^2}{16 m_1^4 m_2^2} \right. \nonumber \\
&\qquad \left. - \frac{9 (N_{12} P_1) (N_{12} P_2)^3 P_1^2}{16 m_1^3 m_2^3} - \frac{5 (N_{12} P_2)^2 P_1^4}{16 m_1^4 m_2^2} - \frac{7 P_1^6}{16 m_1^6} \right. \nonumber \\
&\qquad \left. + \frac{15 (N_{12} P_1)^2 (N_{12} P_2)^2 (P_1 P_2)}{32 m_1^3 m_2^3} + \frac{3 (N_{12} P_1) (N_{12} P_2) P_1^2 (P_1 P_2)}{4 m_1^4 m_2^2} \right. \nonumber \\
&\qquad \left. + \frac{1 (N_{12} P_2)^2 P_1^2 (P_1 P_2)}{16 m_1^3 m_2^3} - \frac{5 (N_{12} P_1) (N_{12} P_2) (P_1 P_2)^2}{16 m_1^3 m_2^3} \right. \nonumber \\
&\qquad \left. + \frac{P_1^2 (P_1 P_2)^2}{8 m_1^4 m_2^2} - \frac{1 (P_1 P_2)^3}{16 m_1^3 m_2^3} - \frac{5 (N_{12} P_1)^2 P_1^2 P_2^2}{16 m_1^4 m_2^2} \right. \nonumber \\
&\qquad \left. + \frac{7 (N_{12} P_1) (N_{12} P_2) P_1^2 P_2^2}{32 m_1^3 m_2^3} + \frac{1 P_1^4 P_2^2}{2 m_1^4 m_2^2} + \frac{1 P_1^2 (P_1 P_2) P_2^2}{32 m_1^3 m_2^3} \right) \right\} \tag{214}
\end{align}
$$

Arguably, the results given by the ADM-Hamiltonian formalism (for the problem at hand) look simpler than their harmonic-coordinate counterparts. Indeed, the ADM Lagrangian is ordinary – no accelerations – and there are no logarithms nor associated gauge constants $r'_1$ and $r'_2$.<fnref n="50" /> Of course, one is free to describe the binary motion in whatever coordinates one likes, and the two formalisms, harmonic (209) and ADM (213)–(214), describe rigorously the same physics. On the other hand, the higher complexity of the harmonic-coordinates Lagrangian (209) enables one to perform more tests of the computations, notably by inquiring about the future of the constants $r'_1$ and $r'_2$, that we know must disappear from physical quantities such as the center-of-mass energy and the total gravitational-wave flux.

<footnote n="50">⁵⁰ On the other hand, the ADM-Hamiltonian formalism provides a limited description of the gravitational radiation field, compared to what will be done using harmonic coordinates in Section 9.</footnote>
## 7.3 Equations of motion in the center-of-mass frame

In this section we translate the origin of coordinates to the binary's center-of-mass by imposing the vanishing of the binary's mass dipole moment: $I_i = 0$ in the notation of Part A. Actually the dipole moment is computed as the center-of-mass conserved integral associated with the boost symmetry of the 3PN equations of motion [174, 79]. This condition results in the 3PN-accurate relationship between the individual positions in the center-of-mass frame $\mathbf{y}_1$ and $\mathbf{y}_2$, and the relative position $\mathbf{x} \equiv \mathbf{y}_1 - \mathbf{y}_2$ and velocity $\mathbf{v} \equiv \mathbf{v}_1 - \mathbf{v}_2 = d\mathbf{x}/dt$ (formerly denoted $\mathbf{y}_{12}$ and $\mathbf{v}_{12}$). We shall also use the orbital separation $r \equiv |\mathbf{x}|$, together with $\mathbf{n} = \mathbf{x}/r$ and $\dot{\mathbf{r}} \equiv \mathbf{n} \cdot \mathbf{v}$. Mass parameters are: The total mass $m = m_1 + m_2$ (to be distinguished from the ADM mass denoted by M in Part A); the relative mass difference $\Delta = (m_1 - m_2)/m$; the reduced mass $\mu = m_1 m_2/m$; and the very useful symmetric mass ratio

$$ \nu \equiv \frac{\mu}{m} \equiv \frac{m_1 m_2}{(m_1 + m_2)^2}. \qquad (215) $$

The usefulness of this ratio lies in its interesting range of variation: $0 < \nu \le 1/4$, with $\nu = 1/4$ in the case of equal masses, and $\nu \to 0$ in the test-mass limit for one of the bodies. Thus $\nu$ is numerically rather small and may be viewed as a small expansion parameter. We also pose $X_1 = m_1/m$ and $X_2 = m_2/m$ so that $\Delta = X_1 - X_2$ and $\nu = X_1 X_2$.

For reference we give the 3PN-accurate expressions of the individual positions in the center-of-mass frame in terms of relative variables. They are in the form

$$ \mathbf{y}_1 = [X_2 + \nu \Delta \mathcal{P}] \mathbf{x} + \nu \Delta \mathcal{Q} \mathbf{v} + \mathcal{O}\left(\frac{1}{c^7}\right), \qquad (216a) $$

$$ \mathbf{y}_2 = [-X_1 + \nu \Delta \mathcal{P}] \mathbf{x} + \nu \Delta \mathcal{Q} \mathbf{v} + \mathcal{O}\left(\frac{1}{c^7}\right), \qquad (216b) $$

where all post-Newtonian corrections, beyond Newtonian order, are proportional to the mass ratio $\nu$ and the mass difference $\Delta$. The two dimensionless coefficients $\mathcal{P}$ and $\mathcal{Q}$ read

$$ \begin{aligned} \mathcal{P} = & \frac{1}{c^2} \left\{ \frac{v^2}{2} - \frac{Gm}{2r} \right\} \\ & + \frac{1}{c^4} \left\{ \frac{3v^4}{8} - \frac{3\nu v^4}{2} + \frac{Gm}{r} \left( -\frac{\dot{r}^2}{8} + \frac{3\dot{r}^2\nu}{4} + \frac{19v^2}{8} + \frac{3\nu v^2}{2} \right) + \frac{G^2m^2}{r^2} \left( \frac{7}{4} - \frac{\nu}{2} \right) \right\} \\ & + \frac{1}{c^6} \left\{ \frac{5v^6}{16} - \frac{11\nu v^6}{4} + 6\nu^2 v^6 \right. \\ & \qquad \left. + \frac{Gm}{r} \left( \frac{\dot{r}^4}{16} - \frac{5\dot{r}^4\nu}{8} + \frac{21\dot{r}^4\nu^2}{16} - \frac{5\dot{r}^2v^2}{16} + \frac{21\dot{r}^2\nu v^2}{16} \right. \right. \\ & \qquad \left. \left. - \frac{11\dot{r}^2\nu^2 v^2}{2} + \frac{53v^4}{16} - 7\nu v^4 - \frac{15\nu^2 v^4}{2} \right) \right. \\ & \qquad \left. + \frac{G^2m^2}{r^2} \left( -\frac{7\dot{r}^2}{3} + \frac{73\dot{r}^2\nu}{8} + 4\dot{r}^2\nu^2 + \frac{101v^2}{12} - \frac{33\nu v^2}{8} + 3\nu^2 v^2 \right) \right. \\ & \qquad \left. + \frac{G^3m^3}{r^3} \left( -\frac{14351}{1260} + \frac{\nu}{8} - \frac{\nu^2}{2} + \frac{22}{3} \ln\left(\frac{r}{r_0''}\right) \right) \right\}, \end{aligned} \qquad (217a) $$

$$ \begin{aligned} \mathcal{Q} = & \frac{1}{c^4} \left\{ -\frac{7Gm\dot{r}}{4} \right\} + \frac{1}{c^5} \left\{ \frac{4Gmv^2}{5} - \frac{8G^2m^2}{5r} \right\} \\ & + \frac{1}{c^6} \left\{ Gm\dot{r} \left( \frac{5\dot{r}^2}{12} - \frac{19\dot{r}^2\nu}{24} - \frac{15v^2}{8} + \frac{21\nu v^2}{4} \right) + \frac{G^2m^2\dot{r}}{r} \left( -\frac{235}{24} - \frac{21\nu}{4} \right) \right\}. \end{aligned} \qquad (217b) $$
Up to 2.5PN order there is agreement with the circular-orbit limit of Eqs. (6.4) in Ref. [45]. Notice the 2.5PN radiation-reaction term entering the coefficient $Q$; such 2.5PN term is explicitly displayed for circular orbits in Eqs. (224) below. In Eqs. (217) the logarithms at the 3PN order appear only in the coefficient $P$. They contain a particular combination $r_0''$ of the two gauge-constants $r_1'$ and $r_2'$ defined by

$$ \Delta \ln r_0'' = X_1^2 \ln r_1' - X_2^2 \ln r_2', \quad (218) $$

and which happens to be different from a similar combination $r_0'$ we shall find in the equations of relative motion, see Eq. (221).

The 3PN and even 3.5PN center-of-mass equations of motion are obtained by replacing in the 3.5PN equations of motion (203) in a general frame, the positions and velocities by their center-of-mass expressions (216)–(217), applying as usual the order-reduction of all accelerations where necessary. We write the relative acceleration in the center-of-mass frame in the form

$$ \frac{d\mathbf{v}}{dt} = -\frac{Gm}{r^2} \left[ (1 + \mathcal{A}) \mathbf{n} + \mathcal{B} \mathbf{v} \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (219) $$

and find that the coefficients $\mathcal{A}$ and $\mathcal{B}$ are [79]

$$
\begin{aligned}
\mathcal{A} ={}& \frac{1}{c^2} \left\{ -\frac{3\dot{r}^2\nu}{2} + v^2 + 3\nu v^2 - \frac{Gm}{r}(4+2\nu) \right\} \\
& + \frac{1}{c^4} \left\{ \frac{15\dot{r}^4\nu}{8} - \frac{45\dot{r}^4\nu^2}{8} - \frac{9\dot{r}^2\nu v^2}{2} + 6\dot{r}^2\nu^2 v^2 + 3\nu v^4 - 4\nu^2 v^4 \right. \\
& \qquad \left. + \frac{Gm}{r} \left( -2\dot{r}^2 - 25\dot{r}^2\nu - 2\dot{r}^2\nu^2 - \frac{13\nu v^2}{2} + 2\nu^2 v^2 \right) + \frac{G^2 m^2}{r^2} \left( 9 + \frac{87\nu}{4} \right) \right\} \\
& + \frac{1}{c^5} \left\{ -\frac{24\dot{r}^2 v^2}{5} \frac{Gm}{r} - \frac{136\dot{r}\nu}{15} \frac{G^2 m^2}{r^2} \right\} \\
& + \frac{1}{c^6} \left\{ -\frac{35\dot{r}^6\nu}{16} + \frac{175\dot{r}^6\nu^2}{16} - \frac{175\dot{r}^6\nu^3}{16} + \frac{15\dot{r}^4\nu v^2}{2} - \frac{135\dot{r}^4\nu^2 v^2}{4} + \frac{255\dot{r}^4\nu^3 v^2}{8} \right. \\
& \qquad \left. - \frac{15\dot{r}^2\nu v^4}{2} + \frac{237\dot{r}^2\nu^2 v^4}{8} - \frac{45\dot{r}^2\nu^3 v^4}{2} + \frac{11\nu v^6}{4} - \frac{49\nu^2 v^6}{4} + 13\nu^3 v^6 \right. \\
& \qquad \left. + \frac{Gm}{r} \left( 79\dot{r}^4\nu - \frac{69\dot{r}^4\nu^2}{2} - 30\dot{r}^4\nu^3 - 121\dot{r}^2\nu v^2 + 16\dot{r}^2\nu^2 v^2 + 20\dot{r}^2\nu^3 v^2 + \frac{75\nu v^4}{4} \right. \right. \\
& \qquad \left. \left. + 8\nu^2 v^4 - 10\nu^3 v^4 \right) \right. \\
& \qquad \left. + \frac{G^2 m^2}{r^2} \left( \dot{r}^2 + \frac{32573\dot{r}^2\nu}{168} + \frac{11\dot{r}^2\nu^2}{8} - 7\dot{r}^2\nu^3 + \frac{615\dot{r}^2\nu\pi^2}{64} - \frac{26987\nu v^2}{840} + \nu^3 v^2 \right. \right. \\
& \qquad \left. \left. - \frac{123\nu\pi^2 v^2}{64} - 110\dot{r}^2\nu \ln\left(\frac{r}{r_0'}\right) + 22\nu v^2 \ln\left(\frac{r}{r_0'}\right) \right) \right. \\
& \qquad \left. + \frac{G^3 m^3}{r^3} \left( -16 - \frac{437\nu}{4} - \frac{71\nu^2}{2} + \frac{41\nu\pi^2}{16} \right) \right\} \\
& + \frac{1}{c^7} \left\{ \frac{Gm}{r} \dot{r} \left( \frac{366}{35}\nu v^4 + 12\nu^2 v^4 - 114\nu^2\nu\dot{r}^2 - 12\nu^2 v^2\dot{r}^2 + 112\nu\dot{r}^4 \right) \right. \\
& \qquad \left. + \frac{G^2 m^2}{r^2} \dot{r} \left( \frac{692}{35}\nu v^2 - \frac{724}{15}\nu v^2 + \frac{294}{5}\nu\dot{r}^2 + \frac{376}{5}\nu^2\dot{r}^2 \right) \right. \\
& \qquad \left. + \frac{G^3 m^3}{r^3} \dot{r} \left( \frac{3956}{35}\nu + \frac{184}{5}\nu^2 \right) \right\},
\end{aligned}
\quad (220a)
$$
$$
\begin{equation}
\begin{aligned}
\mathcal{B} ={}& \frac{1}{c^2} \{-4\dot{r} + 2\dot{r}\nu\} \\
& + \frac{1}{c^4} \left\{ \frac{9\dot{r}^3\nu}{2} + 3\dot{r}^3\nu^2 - \frac{15\dot{r}\nu v^2}{2} - 2\dot{r}\nu^2 v^2 + \frac{Gm}{r} \left( 2\dot{r} + \frac{41\dot{r}\nu}{2} + 4\dot{r}\nu^2 \right) \right\} \\
& + \frac{1}{c^5} \left\{ \frac{8\nu v^2}{5} \frac{Gm}{r} + \frac{24\nu}{5} \frac{G^2 m^2}{r^2} \right\} \\
& + \frac{1}{c^6} \left\{ -\frac{45\dot{r}^5\nu}{8} + 15\dot{r}^5\nu^2 + \frac{15\dot{r}^5\nu^3}{4} + 12\dot{r}^3\nu v^2 - \frac{111\dot{r}^3\nu^2 v^2}{4} - 12\dot{r}^3\nu^3 v^2 - \frac{65\dot{r}\nu v^4}{8} \right. \\
& \qquad \left. + 19\dot{r}\nu^2 v^4 + 6\dot{r}\nu^3 v^4 \right. \\
& \qquad \left. + \frac{Gm}{r} \left( \frac{329\dot{r}^3\nu}{6} + \frac{59\dot{r}^3\nu^2}{2} + 18\dot{r}^3\nu^3 - 15\dot{r}\nu v^2 - 27\dot{r}\nu^2 v^2 - 10\dot{r}\nu^3 v^2 \right) \right. \\
& \qquad \left. + \frac{G^2 m^2}{r^2} \left( -4\dot{r} - \frac{18169\dot{r}\nu}{840} + 25\dot{r}\nu^2 + 8\dot{r}\nu^3 - \frac{123\dot{r}\nu\pi^2}{32} + 44\dot{r}\nu \ln\left(\frac{r}{r_0'}\right) \right) \right\} \\
& + \frac{1}{c^7} \left\{ \frac{Gm}{r} \left( -\frac{626}{35}\nu v^4 - \frac{12}{5}\nu^2 v^4 + \frac{678}{5}\nu v^2 \dot{r}^2 + \frac{12}{5}\nu^2 v^2 \dot{r}^2 - 120\nu \dot{r}^4 \right) \right. \\
& \qquad \left. + \frac{G^2 m^2}{r^2} \left( \frac{164}{21}\nu v^2 + \frac{148}{5}\nu^2 v^2 - \frac{82}{3}\nu \dot{r}^2 - \frac{848}{15}\nu^2 \dot{r}^2 \right) \right. \\
& \qquad \left. + \frac{G^3 m^3}{r^3} \left( -\frac{1060}{21}\nu - \frac{104}{5}\nu^2 \right) \right\}.
\end{aligned}
\tag{220b}
\end{equation}
$$

Up to the 2.5PN order the result agrees with Ref. [302]. The 3.5PN term is issued from Refs. [258, 259, 260, 336, 278, 322, 254]. At the 3PN order we have some gauge-dependent logarithms containing a constant $r'_0$ which is the “logarithmic barycenter” of the two constants $r'_1$ and $r'_2$:

$$ \ln r'_{0} = X_{1} \ln r'_{1} + X_{2} \ln r'_{2}. \tag{221} $$

The logarithms in Eqs. (220), together with the constant $r'_0$ therein, can be removed by applying the gauge transformation (204), while still staying within the class of harmonic coordinates. The resulting modification of the equations of motion will affect only the coefficients of the 3PN order in Eqs. (220); let us denote them by $\mathcal{A}_{3PN}$ and $\mathcal{B}_{3PN}$. The new values of these coefficients, obtained after removal of the logarithms by the latter harmonic gauge transformation, will be denoted $\mathcal{A}_{3PN}^{\text{MH}}$ and $\mathcal{B}_{3PN}^{\text{MH}}$. Here MH stands for the *modified harmonic* coordinate system, differing from the SH (*standard harmonic*) coordinate system containing logarithms at the 3PN order in the coefficients $\mathcal{A}_{3PN}$ and $\mathcal{B}_{3PN}$. See Ref. [9] for a full description of the coordinate transformation between SH and MH coordinates for various quantities. We have [320, 9]

$$
\begin{equation}
\begin{aligned}
\mathcal{A}_{3\text{PN}}^{\text{MH}} ={}& \frac{1}{c^6} \left\{ -\frac{35\dot{r}^6\nu}{16} + \frac{175\dot{r}^6\nu^2}{16} - \frac{175\dot{r}^6\nu^3}{16} + \frac{15\dot{r}^4\nu v^2}{2} - \frac{135\dot{r}^4\nu^2 v^2}{4} + \frac{255\dot{r}^4\nu^3 v^2}{8} - \frac{15\dot{r}^2\nu v^4}{2} \right. \\
& \left. + \frac{237\dot{r}^2\nu^2 v^4}{8} - \frac{45\dot{r}^2\nu^3 v^4}{2} + \frac{11\nu v^6}{4} - \frac{49\nu^2 v^6}{4} + 13\nu^3 v^6 \right. \\
& \left. + \frac{Gm}{r} \left( 79\dot{r}^4\nu - \frac{69\dot{r}^4\nu^2}{2} - 30\dot{r}^4\nu^3 - 121\dot{r}^2\nu v^2 + 16\dot{r}^2\nu^2 v^2 + 20\dot{r}^2\nu^3 v^2 + \frac{75\nu v^4}{4} \right. \right. \\
& \left. \left. + 8\nu^2 v^4 - 10\nu^3 v^4 \right) \right. \\
& \left. + \frac{G^2 m^2}{r^2} \left( \dot{r}^2 + \frac{22717\dot{r}^2\nu}{168} + \frac{11\dot{r}^2\nu^2}{8} - 7\dot{r}^2\nu^3 + \frac{615\dot{r}^2\nu\pi^2}{64} - \frac{20827\nu v^2}{840} + \nu^3 v^2 \right. \right. \\
& \left. \left. - \frac{123\nu\pi^2 v^2}{64} \right) \right\}
\end{aligned}
\tag{222}
\end{equation}
$$
$$ + \frac{G^3 m^3}{r^3} \left( -16 - \frac{1399 \nu}{12} - \frac{71 \nu^2}{2} + \frac{41 \nu \pi^2}{16} \right) \right\} , \qquad (222a) $$

$$ \begin{aligned} B_{3\text{PN}}^{\text{MH}} = & \frac{1}{c^6} \left\{ -\frac{45 \dot{r}^5 \nu}{8} + 15 \dot{r}^5 \nu^2 + \frac{15 \dot{r}^5 \nu^3}{4} + 12 \dot{r}^3 \nu \nu^2 - \frac{111 \dot{r}^3 \nu^2 \nu^2}{4} - 12 \dot{r}^3 \nu^3 \nu^2 - \frac{65 \dot{r} \nu \nu^4}{8} \right. \\ & \left. \qquad + 19 \dot{r} \nu^2 \nu^4 + 6 \dot{r} \nu^3 \nu^4 \right. \\ & \left. \qquad + \frac{Gm}{r} \left( \frac{329 \dot{r}^3 \nu}{6} + \frac{59 \dot{r}^3 \nu^2}{2} + 18 \dot{r}^3 \nu^3 - 15 \dot{r} \nu \nu^2 - 27 \dot{r} \nu^2 \nu^2 - 10 \dot{r} \nu^3 \nu^2 \right) \right. \\ & \left. \qquad + \frac{G^2 m^2}{r^2} \left( -4 \dot{r} - \frac{5849 \dot{r} \nu}{840} + 25 \dot{r} \nu^2 + 8 \dot{r} \nu^3 - \frac{123 \dot{r} \nu \pi^2}{32} \right) \right\} . \end{aligned} \qquad (222b) $$

Again, the other terms in the equations of motion (219)–(220) are unchanged. These gauge-transformed coefficients in MH coordinates are useful because they do not yield the usual complications associated with logarithms. However, they must be handled with care in applications such as in Ref. [320], because one must ensure that all other quantities in the problem (energy, angular momentum, gravitational-wave fluxes, etc.) are defined in the same specific MH gauge avoiding logarithms. In the following we shall no longer use the MH coordinate system leading to Eqs. (222), except when constructing the generalized quasi-Keplerian representation of the 3PN motion in Section 10.2. Therefore all expressions we shall derive below, notably all those concerning the radiation field, are valid in the SH coordinate system in which the equations of motion are fully given by Eq. (203) or, in the center-of-mass frame, by Eqs. (219)–(220).

For future reference let also give the 3PN center-of-mass Hamiltonian in ADM coordinates derived in Refs. [261, 262, 162]. In the center-of-mass frame the conjugate variables are the relative separation $X = Y_1 - Y_2$ and the conjugate momentum (per unit reduced mass) $P$ such that $\mu P = P_1 = -P_2$ where $P_1$ and $P_2$ are defined in Section 7.2). Posing $N \equiv X/R$ with $R \equiv |X|$, together with $P^2 \equiv P^2$ and $P_R \equiv N \cdot P$, we have

$$ \begin{aligned} \frac{H^{\text{ADM}}}{\mu} = & \frac{P^2}{2} - \frac{Gm}{R} \\ & + \frac{1}{c^2} \left\{ -\frac{P^4}{8} + \frac{3 \nu P^4}{8} + \frac{Gm}{R} \left( -\frac{P_R^2 \nu}{2} - \frac{3 P^2}{2} - \frac{\nu P^2}{2} \right) + \frac{G^2 m^2}{2 R^2} \right\} \\ & + \frac{1}{c^4} \left\{ \frac{P^6}{16} - \frac{5 \nu P^6}{16} + \frac{5 \nu^2 P^6}{16} \right. \\ & \qquad \left. + \frac{Gm}{R} \left( -\frac{3 P_R^4 \nu^2}{8} - \frac{P_R^2 P^2 \nu^2}{4} + \frac{5 P^4}{8} - \frac{5 \nu P^4}{2} - \frac{3 \nu^2 P^4}{8} \right) \right. \\ & \qquad \left. + \frac{G^2 m^2}{R^2} \left( \frac{3 P_R^2 \nu}{2} + \frac{5 P^2}{2} + 4 \nu P^2 \right) \right. \\ & \qquad \left. + \frac{G^3 m^3}{R^3} \left( -\frac{1}{4} - \frac{3 \nu}{4} \right) \right\} \\ & + \frac{1}{c^6} \left\{ -\frac{5 P^8}{128} + \frac{35 \nu P^8}{128} - \frac{35 \nu^2 P^8}{64} + \frac{35 \nu^3 P^8}{128} \right. \\ & \qquad \left. + \frac{Gm}{R} \left( -\frac{5 P_R^6 \nu^3}{16} + \frac{3 P_R^4 P^2 \nu^2}{16} - \frac{3 P_R^4 P^2 \nu^3}{16} + \frac{P_R^2 P^4 \nu^2}{8} \right) \right. \\ & \qquad \left. - \frac{3 P_R^2 P^4 \nu^3}{16} - \frac{7 P^6}{16} + \frac{21 \nu P^6}{8} - \frac{53 \nu^2 P^6}{16} - \frac{5 \nu^3 P^6}{16} \right) \right\} \\ & + \frac{G^2 m^2}{R^2} \left( \frac{5 P_R^4 \nu}{12} + \frac{43 P_R^4 \nu^2}{12} + \frac{17 P_R^2 P^2 \nu}{16} \right) \end{aligned} $$
$$ 
\begin{aligned}
& + \left( \frac{15 P_R^2 P^2 \nu^2}{8} - \frac{27 P^4}{16} + \frac{17 \nu P^4}{2} + \frac{109 \nu^2 P^4}{16} \right) \\
& + \frac{G^3 m^3}{R^3} \left( -\frac{85 P_R^2 \nu}{16} - \frac{7 P_R^2 \nu^2}{4} - \frac{25 P^2}{8} - \frac{335 \nu P^2}{48} \right. \\
& \qquad \left. - \frac{23 \nu^2 P^2}{8} - \frac{3 P_R^2 \nu \pi^2}{64} + \frac{\nu P^2 \pi^2}{64} \right) \\
& + \frac{G^4 m^4}{R^4} \left( \frac{1}{8} + \frac{109 \nu}{12} - \frac{21 \nu \pi^2}{32} \right) \right\} + \mathcal{O} \left( \frac{1}{c^8} \right).
\end{aligned} 
\quad (223) $$

## 7.4 Equations of motion and energy for quasi-circular orbits

Most inspiralling compact binaries will have been circularized by the time they become visible by the detectors LIGO and VIRGO; see Section 1.2. In the case of orbits that are circular – apart from the gradual radiation-reaction inspiral – the complicated equations of motion simplify drastically, since we have $\dot{r} = (nv) = \mathcal{O}(1/c^5)$. For circular orbits, up to the 2.5PN order, the relation between center-of-mass variables and the relative ones reads

$$ \mathbf{y}_1 = \mathbf{x} \left[ X_2 + 3\gamma^2\nu\Delta \right] - \frac{4}{5} \frac{G^2\nu m^2 \Delta}{rc^5} \mathbf{v} + \mathcal{O} \left( \frac{1}{c^6} \right), \quad (224a) $$

$$ \mathbf{y}_2 = \mathbf{x} \left[ -X_1 + 3\gamma^2\nu\Delta \right] - \frac{4}{5} \frac{G^2\nu m^2 \Delta}{rc^5} \mathbf{v} + \mathcal{O} \left( \frac{1}{c^6} \right), \quad (224b) $$

where we recall $X_1 = m_1/m$, $X_2 = m_2/m$ and $\Delta = X_1 - X_2$. See Eqs. (216)–(217) for more general formulas. To conveniently display the successive post-Newtonian corrections, we employ the post-Newtonian parameter

$$ \gamma \equiv \frac{Gm}{rc^2} = \mathcal{O} \left( \frac{1}{c^2} \right). \quad (225) $$

Notice that there are no corrections of order 1PN in Eqs. (224) for circular orbits; the dominant term is of order 2PN, i.e., is proportional to $\gamma^2 = \mathcal{O}(1/c^4)$. See Ref. [79] for a systematic calculation of Eqs. (224) to higher order.

The relative acceleration $\mathbf{a} \equiv \mathbf{a}_1 - \mathbf{a}_2$ of two bodies moving on a circular orbit at the 3.5PN order is then given by

$$ \mathbf{a} = -\Omega^2 \mathbf{x} - \frac{32}{5} \frac{G^3 m^3 \nu}{c^5 r^4} \left[ 1 + \gamma \left( -\frac{743}{336} - \frac{11}{4} \nu \right) \right] \mathbf{v} + \mathcal{O} \left( \frac{1}{c^8} \right), \quad (226) $$

where $\mathbf{x} \equiv \mathbf{y}_1 - \mathbf{y}_2$ is the relative separation (in harmonic coordinates) and $\Omega$ denotes the angular frequency of the circular motion. The second term in Eq. (226), opposite to the velocity $\mathbf{v} \equiv \mathbf{v}_1 - \mathbf{v}_2$, represents the radiation reaction force up to 3.5PN order, which comes from the reduction of the coefficients of $1/c^5$ and $1/c^7$ in Eqs. (220). The radiation-reaction force is responsible for the secular decrease of the separation $r$ and increase of the orbital frequency $\Omega$:

$$ \dot{r} = -\frac{64}{5} \frac{G^3 m^3 \nu}{r^3 c^5} \left[ 1 + \gamma \left( -\frac{1751}{336} - \frac{7}{4} \nu \right) \right], \quad (227a) $$

$$ \dot{\Omega} = \frac{96}{5} \frac{Gm\nu}{r^3} \gamma^{5/2} \left[ 1 + \gamma \left( -\frac{2591}{336} - \frac{11}{12}\nu \right) \right]. \quad (227b) $$

Concerning conservative effects, the main content of the 3PN equations (226) is the relation between the frequency $\Omega$ and the orbital separation $r$, which is given by the following generalized
version of Kepler's third law:

$$ \Omega^2 = \frac{Gm}{r^3} \left\{ 1 + (-3 + \nu)\gamma + \left(6 + \frac{41}{4}\nu + \nu^2\right)\gamma^2 + \left(-10 + \left[-\frac{75707}{840} + \frac{41}{64}\pi^2 + 22 \ln\left(\frac{r}{r'_0}\right)\right]\nu + \frac{19}{2}\nu^2 + \nu^3\right)\gamma^3 \right\} + \mathcal{O}\left(\frac{1}{c^8}\right). \quad (228) $$

The length scale $r'_0$ is given in terms of the two gauge-constants $r'_1$ and $r'_2$ by Eq. (221). As for the energy, it is inferred from the circular-orbit reduction of the general result (205). We have

$$ E = -\frac{\mu c^2 \gamma}{2} \left\{ 1 + \left(-\frac{7}{4} + \frac{1}{4}\nu\right)\gamma + \left(-\frac{7}{8} + \frac{49}{8}\nu + \frac{1}{8}\nu^2\right)\gamma^2 + \left(-\frac{235}{64} + \left[\frac{46031}{2240} - \frac{123}{64}\pi^2 + \frac{22}{3}\ln\left(\frac{r}{r'_0}\right)\right]\nu + \frac{27}{32}\nu^2 + \frac{5}{64}\nu^3\right)\gamma^3 \right\} + \mathcal{O}\left(\frac{1}{c^8}\right). \quad (229) $$

This expression is that of a physical observable $E$; however, it depends on the choice of a coordinate system, as it involves the post-Newtonian parameter $\gamma$ defined from the harmonic-coordinate separation $r$. But the numerical value of $E$ should not depend on the choice of a coordinate system, so $E$ must admit a frame-invariant expression, the same in all coordinate systems. To find it we re-express $E$ with the help of the following frequency-related parameter $x$, instead of the post-Newtonian parameter $\gamma$:<fnref n="51" />

$$ x \equiv \left( \frac{Gm\Omega}{c^3} \right)^{2/3} = \mathcal{O} \left( \frac{1}{c^2} \right). \quad (230) $$

We readily obtain from Eq. (228) the expression of $\gamma$ in terms of $x$ at 3PN order,

$$ \gamma = x \left\{ 1 + \left(1 - \frac{\nu}{3}\right)x + \left(1 - \frac{65}{12}\nu\right)x^2 + \left(1 + \left[-\frac{2203}{2520} - \frac{41}{192}\pi^2 - \frac{22}{3}\ln\left(\frac{r}{r'_0}\right)\right]\nu + \frac{229}{36}\nu^2 + \frac{1}{81}\nu^3\right)x^3 + \mathcal{O}\left(\frac{1}{c^8}\right) \right\}, \quad (231) $$

that we substitute back into Eq. (229), making all appropriate post-Newtonian re-expansions. As a result, we gladly discover that the logarithms together with their associated gauge constant $r'_0$ have cancelled out. Therefore, our result is [160, 69]

$$ E = -\frac{\mu c^2 x}{2} \left\{ 1 + \left(-\frac{3}{4} - \frac{1}{12}\nu\right)x + \left(-\frac{27}{8} + \frac{19}{8}\nu - \frac{1}{24}\nu^2\right)x^2 + \left(-\frac{675}{64} + \left[\frac{34445}{576} - \frac{205}{96}\pi^2\right]\nu - \frac{155}{96}\nu^2 - \frac{35}{5184}\nu^3\right)x^3 + \mathcal{O}\left(\frac{1}{c^8}\right) \right\}. \quad (232) $$

For circular orbits one can check that there are no terms of order $x^{7/2}$ in Eq. (232), so this result is actually valid up to the 3.5PN order. We shall discuss in Section 11 how the effects of the spins of the two black holes affect the latter formula.

The formula (232) has been extended to include the logarithmic terms $\propto \ln x$ at the 4PN and 5PN orders [67, 289], that are due to tail effects occurring in the near zone, see Sections 5.2 and

<footnote n="51">⁵¹ This parameter is an invariant in a large class of coordinate systems – those for which the metric becomes asymptotically Minkowskian far from the system: $g_{\alpha\beta} \to \text{diag}(-1, 1, 1, 1)$.</footnote>
5.4. Adding also the Schwarzschild test-mass limit<fnref n="52" /> up to 5PN order, we get:

$$E = -\frac{\mu c^2 x}{2} \left\{ 1 + \left(-\frac{3}{4} - \frac{\nu}{12}\right)x + \left(-\frac{27}{8} + \frac{19}{8}\nu - \frac{\nu^2}{24}\right)x^2 \right. \\ \left. + \left(-\frac{675}{64} + \left[\frac{34445}{576} - \frac{205}{96}\pi^2\right]\nu - \frac{155}{96}\nu^2 - \frac{35}{5184}\nu^3\right)x^3 \right. \\ \left. + \left(-\frac{3969}{128} + \nu e_4(\nu) + \frac{448}{15}\nu \ln x\right)x^4 \right. \\ \left. + \left(-\frac{45927}{512} + \nu e_5(\nu) + \left[-\frac{4988}{35} - \frac{656}{5}\nu\right]\nu \ln x\right)x^5 \right\} + \mathcal{O}\left(\frac{1}{c^{12}}\right). \quad (233)$$

We can write also a similar expression for the angular momentum,

$$J = \frac{G \mu m}{c x^{1/2}} \left\{ 1 + \left(\frac{3}{2} + \frac{\nu}{6}\right)x + \left(\frac{27}{8} - \frac{19}{8}\nu + \frac{\nu^2}{24}\right)x^2 \right. \\ \left. + \left(\frac{135}{16} + \left[-\frac{6889}{144} + \frac{41}{24}\pi^2\right]\nu + \frac{31}{24}\nu^2 + \frac{7}{1296}\nu^3\right)x^3 \right. \\ \left. + \left(\frac{2835}{128} + \nu j_4(\nu) - \frac{64}{3}\nu \ln x\right)x^4 \right. \\ \left. + \left(\frac{15309}{256} + \nu j_5(\nu) + \left[\frac{9976}{105} + \frac{1312}{15}\nu\right]\nu \ln x\right)x^5 \right\} + \mathcal{O}\left(\frac{1}{c^{12}}\right). \quad (234)$$

For circular orbits the energy $E$ and angular momentum $J$ are known to be linked together by the so-called “thermodynamic” relation

$$\frac{\partial E}{\partial \Omega} = \Omega \frac{\partial J}{\partial \Omega}, \quad (235)$$

which is actually just one aspect of the “first law of binary black hole mechanics” that we shall discuss in more details in Section 8.3.

We have introduced in Eqs. (233)–(234) some non-logarithmic 4PN and 5PN coefficients $e_4(\nu)$, $j_4(\nu)$ and $e_5(\nu)$, $j_5(\nu)$, which can however be proved to be *polynomials* in the symmetric mass ratio $\nu$.<fnref n="53" /> Recent works on the 4PN approximation to the equations of motion by means of both EFT methods [204] and the traditional ADM-Hamiltonian approach [264, 265], and complemented by an analytic computation of the gravitational self-force in the small mass ratio $\nu$ limit [36], have yielded the next-order 4PN coefficient as ($\gamma_E$ being Euler’s constant)

$$e_4(\nu) = -\frac{123671}{5760} + \frac{9037}{1536}\pi^2 + \frac{1792}{15}\ln 2 + \frac{896}{15}\gamma_E \\ + \left[-\frac{498449}{3456} + \frac{3157}{576}\pi^2\right]\nu + \frac{301}{1728}\nu^2 + \frac{77}{31104}\nu^3. \quad (236)$$

The numerical value $e_4(0) \simeq 153.88$ was predicted before thanks to a comparison with numerical self-force calculations [289, 287].

$$E^{\text{Schw}} = \mu c^2 \left[ \frac{1 - 2x}{\sqrt{1 - 3x}} - 1 \right].$$

$$j_4(\nu) = -\frac{5}{7}e_4(\nu) + \frac{64}{35},$$

$$j_5(\nu) = -\frac{2}{3}e_5(\nu) - \frac{4988}{945} - \frac{656}{135}\nu.$$

<footnote n="52">⁵² Namely,</footnote>

<footnote n="53">⁵³ From the thermodynamic relation (235) we necessarily have the relations</footnote>
## 7.5 The 2.5PN metric in the near zone

The near-zone metric is given by Eqs. (144) for general post-Newtonian matter sources. For point-particles binaries all the potentials $V$, $V_i$, ... parametrizing the metric must be computed and iterated for delta-function sources. Up to the 2.5PN order it is sufficient to cure the divergences due to singular sources by means of the Hadamard self-field regularization. Let us point out that the computation is greatly helped – and indeed is made possible at all – by the existence of the following solution $g$ of the elementary Poisson equation

$$ \Delta g = \frac{1}{r_1 r_2}, \qquad (237) $$

which takes the very nice closed analytic form [202]

$$ g = \ln S, \qquad (238a) $$

$$ S \equiv r_1 + r_2 + r_{12}, \qquad (238b) $$

where $r_a = |\mathbf{x} - \mathbf{y}_a|$ and $r_{12} = |\mathbf{y}_1 - \mathbf{y}_2|$. Furthermore, to obtain the metric at the 2.5PN order, the solutions of even more difficult elementary Poisson equations are required. Namely we meet

$$ \Delta K_1 = 2 \partial_i \partial_j \left( \frac{1}{r_2} \right) \partial_i \partial_j \ln r_1, \qquad (239a) $$

$$ \Delta H_1 = 2 \partial_i \partial_j \left( \frac{1}{r_1} \right) \partial_i \partial_j g, \qquad (239b) $$

with $a_i \partial_i$ denoting the partial derivatives with respect to the source points $y_a^i$ (and as usual $\partial_i$ being the partial derivative with respect to the field point $x^i$). It is quite remarkable that the solutions of the latter equations are known in closed analytic form. By combining several earlier results from Refs. [120, 324, 377], one can write these solutions into the form [64, 76]

$$ K_1 = \left( \frac{1}{2} \Delta - \Delta_1 \right) \left[ \frac{\ln r_1}{r_2} \right] + \frac{1}{2} \Delta_2 \left[ \frac{\ln r_{12}}{r_2} \right] + \frac{r_2}{2r_{12}^2 r_1^2} + \frac{1}{r_{12}^2 r_2}, \qquad (240a) $$

$$ H_1 = \frac{1}{2} \Delta_1 \left[ \frac{g}{r_1} + \frac{\ln r_1}{r_{12}} - \Delta_1 \left( \frac{r_1 + r_{12}}{2} g \right) \right] - \frac{r_2}{2r_1^2 r_{12}^2} + \partial_i \partial_j \left[ \frac{\ln r_{12}}{r_1} + \frac{\ln r_1}{2r_{12}} \right] + \frac{1}{2} \Delta_2 \left[ \frac{\ln r_{12}}{r_1} \right], \qquad (240b) $$

where $\Delta_a$ are the Laplacians with respect to the two source points.

We report here the complete expression of the 2.5PN metric in harmonic coordinates valid at any field point in the near zone. Posing $g_{\alpha\beta} = \eta_{\alpha\beta} + k_{\alpha\beta}$ we have [76]

$$ k_{00} = \frac{2Gm_1}{c^2 r_1} + \frac{1}{c^4} \left[ \frac{Gm_1}{r_1} \left( -(n_1 v_1)^2 + 4v_1^2 \right) - 2 \frac{G^2 m_1^2}{r_1^2} \right. \\ + G^2 m_1 m_2 \left( -\frac{2}{r_1 r_2} - \frac{r_1}{2r_{12}^3} + \frac{r_1^2}{2r_2 r_{12}^3} - \frac{5}{2r_2 r_{12}} \right) + \frac{4G^2 m_1 m_2}{3c^5 r_{12}^2} (n_{12} v_{12}) \\ + \frac{1}{c^6} \left[ \frac{Gm_1}{r_1} \left( \frac{3}{4}(n_1 v_1)^4 - 3(n_1 v_1)^2 v_1^2 + 4v_1^4 \right) + \frac{G^2 m_1^2}{r_1^2} (3(n_1 v_1)^2 - v_1^2) + 2 \frac{G^3 m_1^3}{r_1^3} \right. \\ + G^2 m_1 m_2 \left( v_1^2 \left( \frac{3r_1^3}{8r_{12}^5} - \frac{3r_1^2 r_2}{8r_{12}^5} - \frac{3r_1 r_2^2}{8r_{12}^5} + \frac{3r_2^3}{8r_{12}^5} - \frac{37r_1}{8r_{12}^3} + \frac{r_1^2}{r_2 r_{12}^3} + \frac{3r_2}{8r_{12}^3} \right. \right. \\ \left. + \frac{2r_2^2}{r_1 r_{12}^3} + \frac{6}{r_1 r_{12}} - \frac{5}{r_2 r_{12}} - \frac{8r_{12}}{r_1 r_2 S} + \frac{16}{r_{12} S} \right) \right] $$
$$
\begin{align*}
&+ (v_1 v_2) \left( \frac{8}{r_1 r_2} - \frac{3r_1^3}{4r_{12}^5} + \frac{3r_1^2 r_2}{4r_{12}^5} + \frac{13r_1}{4r_{12}^3} - \frac{2r_1^2}{r_2 r_{12}^3} - \frac{6}{r_1 r_{12}} - \frac{16}{r_1 S} - \frac{12}{r_{12} S} \right) \\
&+ (n_{12} v_1)^2 \left( -\frac{15r_1^3}{8r_{12}^5} + \frac{15r_1^2 r_2}{8r_{12}^5} + \frac{15r_1 r_2^2}{8r_{12}^5} - \frac{15r_2^3}{8r_{12}^5} + \frac{57r_1}{8r_{12}^3} - \frac{3r_1^2}{4r_2 r_{12}^3} - \frac{33r_2}{8r_{12}^3} \right. \\
&\qquad \left. + \frac{7}{4r_2 r_{12}} - \frac{16}{S^2} - \frac{16}{r_{12} S} \right) \\
&+ (n_{12} v_1)(n_{12} v_2) \left( \frac{15r_1^3}{4r_{12}^5} - \frac{15r_1^2 r_2}{4r_{12}^5} - \frac{9r_1}{4r_{12}^3} + \frac{12}{S^2} + \frac{12}{r_{12} S} \right) \\
&+ (n_1 v_1)^2 \left( \frac{2}{r_1 r_2} - \frac{r_1}{4r_{12}^3} - \frac{3r_2^2}{4r_1 r_{12}^3} + \frac{7}{4r_1 r_{12}} - \frac{8}{S^2} - \frac{8}{r_1 S} \right) \\
&+ (n_1 v_1)(n_1 v_2) \left( \frac{r_1}{r_{12}^3} + \frac{16}{S^2} + \frac{16}{r_1 S} \right) - (n_1 v_2)^2 \left( \frac{8}{S^2} + \frac{8}{r_1 S} \right) \\
&+ (n_{12} v_1)(n_1 v_1) \left( -\frac{3r_1^2}{r_{12}^4} + \frac{3r_2^2}{2r_{12}^4} + \frac{3}{2r_{12}^2} + \frac{16}{S^2} \right) + \frac{16(n_1 v_2)(n_2 v_1)}{S^2} \\
&+ (n_{12} v_2)(n_1 v_1) \left( \frac{3r_1^2}{r_{12}^4} - \frac{3r_2^2}{2r_{12}^4} + \frac{13}{2r_{12}^2} - \frac{40}{S^2} \right) - \frac{12(n_1 v_1)(n_2 v_2)}{S^2} \\
&+ (n_{12} v_1)(n_1 v_2) \left( \frac{3r_1^2}{2r_{12}^4} + \frac{4}{r_{12}^2} + \frac{16}{S^2} \right) + (n_{12} v_2)(n_1 v_2) \left( \frac{-3r_1^2}{2r_{12}^4} - \frac{3}{r_{12}^2} + \frac{16}{S^2} \right) \\
&+ G^3 m_1^2 m_2 \left( \frac{4}{r_1^3} + \frac{1}{2r_2^3} + \frac{9}{2r_1^2 r_2} - \frac{r_1^3}{4r_{12}^6} + \frac{3r_1^4}{16r_2 r_{12}^6} - \frac{r_1^2 r_2}{8r_{12}^6} + \frac{r_1 r_2^2}{4r_{12}^6} - \frac{r_2^3}{16r_{12}^6} + \frac{5r_1}{4r_{12}^4} \right. \\
&\qquad \left. - \frac{23r_1^2}{8r_2 r_{12}^4} + \frac{43r_2}{8r_{12}^4} - \frac{5r_2^2}{2r_1 r_{12}^4} - \frac{3}{r_{12}^3} + \frac{3r_1}{2r_2 r_{12}^3} + \frac{r_2}{r_1 r_{12}^3} - \frac{5r_2^2}{r_1^2 r_{12}^3} + \frac{4r_2^3}{r_1^3 r_{12}^3} + \frac{3}{2r_1 r_{12}^2} \right. \\
&\qquad \left. - \frac{r_1^2}{4r_2^3 r_{12}^2} + \frac{3}{16r_2 r_{12}^2} + \frac{15r_2}{4r_1^2 r_{12}^2} - \frac{4r_2^2}{r_1^3 r_{12}^2} + \frac{5}{r_1^2 r_{12}} + \frac{5}{r_1 r_2 r_{12}} - \frac{4r_2}{r_1^3 r_{12}} - \frac{r_{12}^2}{4r_1^2 r_2^3} \right) \\
&+ \frac{1}{c^7} \left[ G^2 m_1 m_2 \left( (n_{12} v_{12})^2 (n_1 v_1) \left( -\frac{8r_1^3}{r_{12}^5} - \frac{16r_1}{r_{12}^3} \right) + (n_{12} v_{12})^2 (n_1 v_2) \left( \frac{8r_1^3}{r_{12}^5} + \frac{5r_1}{r_{12}^3} \right) \right. \right. \\
&\qquad \left. + (n_{12} v_{12})^3 \left( -\frac{7r_1^4}{2r_{12}^6} + \frac{7r_1^2 r_2^2}{2r_{12}^6} - \frac{11r_1^2}{r_{12}^4} - \frac{37}{4r_{12}^2} \right) \right. \\
&\qquad \left. + (n_{12} v_1)(n_{12} v_{12})^2 \left( \frac{20r_1^2}{r_{12}^4} - \frac{11}{2r_{12}^2} \right) - 4(n_{12} v_{12})(n_1 v_1)^2 \frac{r_1^2}{r_{12}^4} \right. \\
&\qquad \left. + 4(n_{12} v_{12})(n_1 v_1)(n_1 v_2) \frac{r_1^2}{r_{12}^4} + (n_{12} v_1)^2 (n_1 v_1) \frac{r_1}{r_{12}^3} \right. \\
&\qquad \left. + 22(n_{12} v_1)(n_{12} v_{12})(n_1 v_1) \frac{r_1}{r_{12}^3} - (n_{12} v_1)^2 (n_1 v_2) \frac{r_1}{r_{12}^3} \right. \\
&\qquad \left. + 4(n_{12} v_1)(n_{12} v_{12})(n_1 v_2) \frac{r_1}{r_{12}^3} + 11(n_{12} v_1)^2 (n_{12} v_{12}) \frac{1}{2r_{12}^2} \right. \\
&\qquad \left. + (n_1 v_2)^2 v_{12}^2 \left( -\frac{8r_1^3}{5r_{12}^5} - \frac{2r_1}{3r_{12}^3} \right) + (n_1 v_1)^2 v_{12}^2 \left( \frac{8r_1^3}{5r_{12}^5} + \frac{11r_1}{3r_{12}^3} \right) \right. \\
&\qquad \left. - (n_{12} v_1)^2 v_{12}^2 \left( \frac{4r_1^2}{r_{12}^4} + \frac{5}{2r_{12}^2} \right) - (n_{12} v_{12}) v_1^2 \left( \frac{12r_1^2}{r_{12}^4} + \frac{5}{2r_{12}^2} \right) \right. \\
&\qquad \left. + (n_{12} v_{12}) v_{12}^2 \left( \frac{3r_1^4}{2r_{12}^6} - \frac{3r_1^2 r_2^2}{2r_{12}^6} + \frac{7r_1^2}{r_{12}^4} + \frac{27}{4r_{12}^2} \right) \right]
\end{align*}
$$
$$
\begin{equation}
\begin{aligned}
& -29(n_1 v_1) v_1^2 \frac{r_1}{3r_{12}^3} + (n_1 v_2) v_1^2 \frac{r_1}{r_{12}^3} + 5(n_{12} v_1) v_1^2 \frac{1}{r_{12}^2} \\
& \quad + (n_{12} v_{12}) (v_1 v_2) \left( \frac{12r_1^2}{r_{12}^4} + \frac{3}{r_{12}^2} \right) + 8(n_1 v_1) (v_1 v_2) \frac{r_1}{r_{12}^3} \\
& \quad + 2(n_1 v_2) (v_1 v_2) \frac{r_1}{3r_{12}^3} - 5(n_{12} v_1) (v_1 v_2) \frac{1}{r_{12}^2} \\
& + G^3 m_1^2 m_2 \left( (n_1 v_{12}) \left( -\frac{8r_1^3}{15r_{12}^6} + \frac{8r_1 r_2^2}{15r_{12}^6} - \frac{16r_1}{3r_{12}^4} + \frac{8}{r_{12}^3} - \frac{8r_2^2}{r_1^2 r_{12}^3} + \frac{8}{r_1^2 r_{12}^2} \right) \right. \\
& \quad \left. + (n_{12} v_1) \left( -\frac{4r_1^2}{3r_{12}^5} + \frac{4r_2^2}{3r_{12}^5} + \frac{20}{3r_{12}^3} \right) + 8(n_1 v_1) \frac{r_1}{3r_{12}^4} \right. \\
& \quad \left. + (n_{12} v_{12}) \left( \frac{3}{r_1^3} - \frac{4r_1^2}{3r_{12}^5} + \frac{68r_2^2}{15r_{12}^5} + \frac{3r_1}{r_{12}^4} - \frac{6r_2^2}{r_1 r_{12}^4} + \frac{3r_2^4}{r_1^3 r_{12}^4} - \frac{76}{3r_{12}^3} \right. \right. \\
& \quad \left. \left. + \frac{2}{r_1 r_{12}^2} - \frac{6r_2^2}{r_1^3 r_{12}^2} \right) \right] + 1 \leftrightarrow 2 + \mathcal{O}\left(\frac{1}{c^8}\right),
\end{aligned}
\tag{241a}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
k_{0i} = & -\frac{4Gm_1}{c^3 r_1} v_1^i + \frac{1}{c^5} \left[ n_1^i \left( -\frac{G^2 m_1^2}{r_1^2} (n_1 v_1) + \frac{G^2 m_1 m_2}{S^2} \left( -16(n_{12} v_1) + 12(n_{12} v_2) \right. \right. \right. \\
& \left. \left. -16(n_2 v_1) + 12(n_2 v_2) \right) \right) \\
& + n_{12}^i G^2 m_1 m_2 \left( -6(n_{12} v_{12}) \frac{r_1}{r_{12}^3} - 4(n_1 v_1) \frac{1}{r_{12}^2} + 12(n_1 v_1) \frac{1}{S^2} \right. \\
& \left. \left. -16(n_1 v_2) \frac{1}{S^2} + 4(n_{12} v_1) \frac{1}{S} \left( \frac{1}{S} + \frac{1}{r_{12}} \right) \right) \right. \\
& + v_1^i \left( \frac{Gm_1}{r_1} (2(n_1 v_1)^2 - 4v_1^2) + \frac{G^2 m_1^2}{r_1^2} + G^2 m_1 m_2 \left( \frac{3r_1}{r_{12}^3} - \frac{2r_2}{r_{12}^3} \right) \right. \\
& \left. \left. + G^2 m_1 m_2 \left( -\frac{r_2^2}{r_1 r_{12}^3} - \frac{3}{r_1 r_{12}} + \frac{8}{r_2 r_{12}} - \frac{4}{r_{12} S} \right) \right) \right] \\
& + \frac{1}{c^6} \left[ n_{12}^i \left( G^2 m_1 m_2 \left( -10(n_{12} v_{12})^2 \frac{r_1^2}{r_{12}^4} - 12(n_{12} v_{12}) (n_1 v_1) \frac{r_1}{r_{12}^3} + 2v_{12}^2 \frac{r_1^2}{r_{12}^4} - 4 \frac{v_1^2}{r_{12}^2} \right) \right. \right. \\
& \left. \left. + G^3 m_1^2 m_2 \left( \frac{2r_1^2}{3r_{12}^5} - \frac{2r_2^2}{3r_{12}^5} - \frac{2}{r_{12}^3} \right) \right) \right. \\
& \left. + v_1^i \frac{G^2 m_1 m_2}{r_{12}^3} \left( \frac{16(n_1 v_{12}) r_1}{3} - 4(n_{12} v_2) r_{12} \right) \right. \\
& \left. + v_{12}^i \frac{G^2 m_1 m_2}{r_{12}^2} \left( -2(n_{12} v_1) + 6(n_{12} v_2) \frac{r_1^2}{r_{12}^2} \right) \right] + 1 \leftrightarrow 2 + \mathcal{O}\left(\frac{1}{c^7}\right),
\end{aligned}
\tag{241b}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
k_{ij} ={}& \frac{2Gm_1}{c^2 r_1} \delta^{ij} + \frac{1}{c^4} \left[ \delta^{ij} \left( -\frac{Gm_1}{r_1} (n_1 v_1)^2 + \frac{G^2 m_1^2}{r_1^2} \right. \right. \\
& \left. \left. + G^2 m_1 m_2 \left( \frac{2}{r_1 r_2} - \frac{r_1}{2r_{12}^3} + \frac{r_1^2}{2r_2 r_{12}^3} - \frac{5}{2r_1 r_{12}} + \frac{4}{r_{12} S} \right) \right) \right. \\
& \left. + 4 \frac{Gm_1}{r_1} v_1^i v_1^j + \frac{G^2 m_1^2}{r_1^2} n_1^i n_1^j - 4 G^2 m_1 m_2 n_{12}^i n_{12}^j \left( \frac{1}{S^2} + \frac{1}{r_{12} S} \right) \right. \\
& \left. + \frac{4G^2 m_1 m_2}{S^2} \left( n_1^{(i} n_2^{j)} + 2 n_1^{(i} n_{12}^{j)} \right) \right]
\end{aligned}
\tag{241c}
\end{equation}
$$
$$ + \frac{G^2 m_1 m_2}{c^5 r_{12}^2} \left( -\frac{2}{3} (n_{12} v_{12}) \delta^{ij} - 6 (n_{12} v_{12}) n_{12}^i n_{12}^j + 8 n_{12}^{(i} v_{12}^{j)} \right) + 1 \leftrightarrow 2 + \mathcal{O}\left(\frac{1}{c^6}\right). \quad (241c) $$

Here we pose $S = r_1 + r_2 + r_{12}$ and $1 \leftrightarrow 2$ refers to the same quantity but with all particle labels exchanged. To higher order one needs the solution of elementary equations still more intricate than (239) and the 3PN metric valid in closed form all over the near zone is not currently known.

Let us also display the latter 2.5PN metric computed at the location of the particle 1 for instance, thanks to the Hadamard self-field regularization, i.e., in the sense of the Hadamard partie finie defined by Eq. (160). We get

$$ (k_{00})_1 = \frac{2Gm_2}{c^2 r_{12}} + \frac{Gm_2}{c^4 r_{12}} \left( 4v_2^2 - (n_{12}v_2)^2 - 3\frac{Gm_1}{r_{12}} - 2\frac{Gm_2}{r_{12}} \right) \\ + \frac{8G^2 m_1 m_2}{3c^5 r_{12}^2} (n_{12}v_{12}) + \frac{Gm_2}{c^6 r_{12}} \left( \frac{3}{4}(n_{12}v_2)^4 - 3(n_{12}v_2)^2 v_2^2 + 4v_2^4 \right) \\ + \frac{G^2 m_1 m_2}{c^6 r_{12}^2} \left( -\frac{87}{4}(n_{12}v_1)^2 + \frac{47}{2}(n_{12}v_1)(n_{12}v_2) - \frac{55}{4}(n_{12}v_2)^2 + \frac{23}{4}v_1^2 - \frac{39}{2}(v_1v_2) \right) \\ + \frac{47}{4} \frac{G^2 m_1 m_2}{c^6 r_{12}^2} v_2^2 + \frac{Gm_2}{c^6 r_{12}} \left\{ \frac{Gm_2}{r_{12}} [3(n_{12}v_2)^2 - v_2^2] - \frac{G^2 m_1^2}{r_{12}^2} + \frac{17}{2} \frac{G^2 m_1 m_2}{r_{12}^2} + 2 \frac{G^2 m_2^2}{r_{12}^2} \right\} \\ + \frac{G^2 m_1 m_2}{c^7 r_{12}^2} \left\{ -20(n_{12}v_1)^3 + 40(n_{12}v_1)^2(n_{12}v_2) - 36(n_{12}v_1)(n_{12}v_2)^2 + 16(n_{12}v_2)^3 \right. \\ + \frac{296}{15}(n_{12}v_1)v_1^2 - \frac{116}{15}(n_{12}v_2)v_1^2 - \frac{104}{5}(n_{12}v_1)(v_1v_2) + \frac{232}{15}(n_{12}v_2)(v_1v_2) \\ + \frac{56}{15}(n_{12}v_1)v_2^2 - \frac{52}{5}(n_{12}v_2)v_2^2 + \frac{Gm_1}{r_{12}} \left( -\frac{64}{5}(n_{12}v_1) + \frac{104}{5}(n_{12}v_2) \right) \\ \left. + \frac{Gm_2}{r_{12}} \left( -\frac{144}{5}(n_{12}v_1) + \frac{392}{15}(n_{12}v_2) \right) \right\} + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (242a) $$

$$ (k_{0i})_1 = -\frac{4Gm_2}{c^3 r_{12}} v_2^i + \frac{Gm_2}{c^5 r_{12}} \left\{ n_{12}^i \left[ \frac{Gm_1}{r_{12}} (10(n_{12}v_1) + 2(n_{12}v_2)) - \frac{Gm_2}{r_{12}} (n_{12}v_2) \right] \right\} \\ + \frac{Gm_2}{c^5 r_{12}} \left\{ 4 \frac{Gm_1}{r_{12}} v_1^i + v_2^i \left( 2(n_{12}v_2)^2 - 4v_2^2 - 2 \frac{Gm_1}{r_{12}} + \frac{Gm_2}{r_{12}} \right) \right\} \\ + \frac{G^2 m_1 m_2}{c^6 r_{12}^2} \left\{ n_{12}^i (10(n_{12}v_1)^2 - 8(n_{12}v_1)(n_{12}v_2) - 2(n_{12}v_2)^2 - 6v_1^2 \right. \\ \left. + 4(v_1v_2) + 2v_2^2 - \frac{8}{3} \frac{Gm_1}{r_{12}} + \frac{4}{3} \frac{Gm_2}{r_{12}}) \right. \\ \left. - 8(n_{12}v_1)v_1^i + v_2^i \left( \frac{20}{3}(n_{12}v_1) + \frac{4}{3}(n_{12}v_2) \right) \right\} + \mathcal{O}\left(\frac{1}{c^7}\right), \quad (242b) $$

$$ (k_{ij})_1 = \frac{2Gm_2}{c^2 r_{12}} \delta^{ij} + \frac{Gm_2}{c^4 r_{12}} \delta^{ij} \left( -(n_{12}v_2)^2 + \frac{Gm_1}{r_{12}} + \frac{Gm_2}{r_{12}} \right) \\ + \frac{Gm_2}{c^4 r_{12}} \left\{ n_{12}^{ij} \left( -8 \frac{Gm_1}{r_{12}} + \frac{Gm_2}{r_{12}} \right) + 4v_2^{ij} \right\} - \frac{4G^2 m_1 m_2}{3c^5 r_{12}^2} \delta^{ij} (n_{12}v_{12}) \\ + \frac{G^2 m_1 m_2}{c^5 r_{12}^2} \left\{ -12n_{12}^{(i} (n_{12}v_{12}) + 16n_{12}^{(i} v_{12}^{j)} \right\} + \mathcal{O}\left(\frac{1}{c^6}\right). \quad (242c) $$

When regularized at the location of the particles, the metric can be computed to higher order, for instance 3PN. We shall need it when we compute the so-called redshift observable in Sections (8.3) and (8.4); indeed, see Eq. (276).
# 8 Conservative Dynamics of Compact Binaries

## 8.1 Concept of innermost circular orbit

Having in hand the conserved energy $E(x)$ for circular orbits given by Eq. (232), or even more accurate by (233), we define the innermost circular orbit (ICO) as the minimum, when it exists, of the energy function $E(x)$ – see e.g., Ref. [51]. Notice that the ICO is not defined as a point of dynamical general-relativistic instability. Hence, we prefer to call this point the ICO rather than, strictly speaking, an innermost stable circular orbit or ISCO. A study of the dynamical stability of circular binary orbits in the post-Newtonian approximation is reported in Section 8.2.

The previous definition of the ICO is motivated by the comparison with some results of numerical relativity. Indeed we shall confront the prediction of the standard (Taylor-based) post-Newtonian approximation with numerical computations of the energy of binary black holes under the assumptions of conformal flatness for the spatial metric and of exactly circular orbits [228, 232, 133, 121]. The latter restriction is implemented by requiring the existence of an “helical” Killing vector (HKV), which is time-like inside the light cylinder associated with the circular motion, and space-like outside. The HKV will be defined in Eq. (273) below. In the numerical approaches of Refs. [228, 232, 133, 121] there are no gravitational waves, the field is periodic in time, and the gravitational potentials tend to zero at spatial infinity within a restricted model equivalent to solving five out of the ten Einstein field equations (the so-called Isenberg–Wilson–Mathews approximation; see Ref. [228] for a discussion). Considering an evolutionary sequence of equilibrium configurations the circular-orbit energy $E(\Omega)$ and the ICO of binary black holes are obtained numerically (see also Refs. [92, 229, 301] for related calculations of binary neutron stars and strange quark stars).

Since the numerical calculations [232, 133] have been performed in the case of two corotating black holes, which are spinning essentially with the orbital angular velocity, we must for the comparison include within our post-Newtonian formalism the effects of spins appropriate to two Kerr black holes rotating at the orbital rate. The total relativistic masses of the two Kerr black holes (with $a = 1, 2$ labelling the black holes) are given by$^{54}$

$$m_a^2 = \mu_a^2 + \frac{S_a^2}{4\mu_a^2}. \qquad (243)$$

We assume the validity of the Christodoulou mass formula for Kerr black holes [127, 129]; i.e., we neglect the influence of the companion. Here $S_a$ is the spin, related to the usual Kerr parameter by $S_a = m_a a_a$, and $\mu_a \equiv m_a^{\text{irr}}$ is the irreducible mass, not to be confused with the reduced mass of the binary system, and given by $4\pi\mu_a = \sqrt{A_a}$ ($A_a$ is the hole’s surface area). The angular velocity of the black hole, defined by the angular velocity of the outgoing photons that remain for ever at the location of the light-like horizon, is

$$\omega_a = \left. \frac{\partial m_a}{\partial S_a} \right|_{\mu_a} = \frac{S_a}{4m_a \mu_a^2}. \qquad (244)$$

We shall give in Eq. (284) below a more general formulation of the “internal structure” of the black holes. Combining Eqs. (243)–(244) we obtain $m_a$ and $S_a$ as functions of $\mu_a$ and $\omega_a$,

$$m_a = \frac{\mu_a}{\sqrt{1 - 4\mu_a^2 \omega_a^2}}, \qquad (245a)$$

$$S_a = \frac{4\mu_a^3 \omega_a}{\sqrt{1 - 4\mu_a^2 \omega_a^2}}. \qquad (245b)$$

<footnote>$^{54}$ In all of Section 8 we pose $G = 1 = c$.</footnote>
In the limit of slow rotation we get

$$S_a = I_a \omega_a + \mathcal{O}(\omega_a^3), \qquad (246a)$$

$$m_a = \mu_a + \frac{1}{2} I_a \omega_a^2 + \mathcal{O}(\omega_a^4), \qquad (246b)$$

where $I_a = 4\mu_a^3$ is the moment of inertia of the black hole. We see that the total mass-energy $m_a$ involves the irreducible mass augmented by the usual kinetic energy of the spin.

We now need the relation between the rotation frequency $\omega_a$ of each of the corotating black holes and the orbital frequency $\Omega$ of the binary system. Indeed $\Omega$ is the basic variable describing each equilibrium configuration calculated numerically in Refs. [232, 133], with the irreducible masses held constant along the numerical evolutionary sequences. Here we report the result of an investigation of the condition for corotation based on the first law of mechanics for spinning black holes [55], which concluded that the corotation condition at 2PN order reads

$$\omega_a = \Omega \left\{ 1 - \nu x + \nu \left( -\frac{3}{2} + \frac{\nu}{3} \right) x^2 + \mathcal{O}(x^3) \right\}, \qquad (247)$$

where $x$ denotes the post-Newtonian parameter (230) and $\nu$ the symmetric mass ratio (215). The condition (247) is issued from the general relation which will be given in Eq. (285). Interestingly, notice that $\omega_1 = \omega_2$ up to the rather high 2PN order. In the Newtonian limit $x \to 0$ or the test-particle limit $\nu \to 0$ we simply have $\omega_a = \Omega$, in agreement with physical intuition.

To take into account the spin effects our first task is to replace all the masses entering the energy function (232) by their equivalent expressions in terms of $\omega_a$ and the irreducible masses $\mu_a$, and then to replace $\omega_a$ in terms of $\Omega$ according to the corotation prescription (247).$^{55}$ It is clear that the leading contribution is that of the spin kinetic energy given in Eq. (246b), and it comes from the replacement of the rest mass-energy $m = m_1 + m_2$. From Eq. (246b) this effect is of order $\Omega^2$ in the case of corotating binaries, which means by comparison with Eq. (232) that it is equivalent to an “orbital” effect at the 2PN order (i.e., $\propto x^2$). Higher-order corrections in Eq. (246b), which behave at least like $\Omega^4$, will correspond to the orbital 5PN order at least and are negligible for the present purpose. In addition there will be a subdominant contribution, of the order of $\Omega^{8/3}$ equivalent to 3PN order, which comes from the replacement of the masses into the Newtonian part, proportional to $x \propto \Omega^{2/3}$, of the energy $E$; see Eq. (232). With the 3PN accuracy we do not need to replace the masses that enter into the post-Newtonian corrections in $E$, so in these terms the masses can be considered to be the irreducible ones.

Our second task is to include the specific relativistic effects due to the spins, namely the spin-orbit (SO) interaction and the spin-spin (SS) one. In the case of spins $S_1$ and $S_2$ aligned parallel to the orbital angular momentum (and right-handed with respect to the sense of motion) the SO energy reads

$$E_{\text{SO}} = -m \nu (m\Omega)^{5/3} \left[ \left( \frac{4 m_1^2}{3 m^2} + \nu \right) \frac{S_1}{m_1^2} + \left( \frac{4 m_2^2}{3 m^2} + \nu \right) \frac{S_2}{m_2^2} \right]. \qquad (248)$$

We shall review in Section 11 the most up-to-date results for the spin-orbit energy and related quantities; here we are simply employing the leading-order formula obtained in Refs. [27, 28, 275, 271] and given by the first term in Eq. (415). We immediately infer from this formula that in the case of corotating black holes the SO effect is equivalent to a 3PN orbital effect and thus must be retained with the present accuracy. With this approximation, the masses in Eq. (248) can be replaced by the irreducible ones. As for the SS interaction (still in the case of spins aligned with

<footnote>$^{55}$ Note that this is an iterative process because the masses in Eq. (247) are themselves to be replaced by the irreducible masses.</footnote>
the orbital angular momentum) it is given by

$$E_{SS} = \mu \nu (m\Omega)^2 \frac{S_1 S_2}{m_1^2 m_2^2}. \qquad (249)$$

The SS effect can be neglected here because it is of order 5PN for corotating systems. Summing up all the spin contributions to 3PN order we find that the supplementary energy due to the corotating spins is [51, 55]$^{56}$

$$\Delta E^{\text{corot}} = \mu x_{\mu} \left[ (2 - 6\eta) x_{\mu}^2 + \eta(-10 + 25\eta) x_{\mu}^3 + \mathcal{O}(x_{\mu}^4) \right]. \qquad (250)$$

The total mass $\mu = \mu_1 + \mu_2$, the symmetric mass ratio $\eta = \mu_1\mu_2/\mu^2$, and the dimensionless invariant post-Newtonian parameter $x_\mu = (\mu\Omega)^{2/3}$ are now expressed in terms of the irreducible masses $\mu_a$, rather than the masses $m_a$. The complete 3PN energy of the corotating binary is finally given by the sum of Eqs. (232) and (250), in which all the masses are now understood as being the irreducible ones, which must be assumed to stay constant when the binary evolves for the comparison with the numerical calculation.

The left panel of Figure 1 shows the results for $E_{\text{ICO}}$ in the case of irrotational and corotational binaries. Since $\Delta E^{\text{corot}}$, given by Eq. (250), is at least of order 2PN, the result for $1\text{PN}^{\text{corot}}$ is the same as for 1PN in the irrotational case; then, obviously, $2\text{PN}^{\text{corot}}$ takes into account only the leading 2PN corotation effect, i.e., the spin kinetic energy given by Eq. (246b), while $3\text{PN}^{\text{corot}}$ involves also, in particular, the corotational SO coupling at the 3PN order. In addition we present the numerical point obtained by numerical relativity under the assumptions of conformal flatness and of helical symmetry [228, 232]. As we can see the 3PN points, and even the 2PN ones, are in good agreement with the numerical value. The fact that the 2PN and 3PN values are so close to each other is a good sign of the convergence of the expansion. In fact one might say that the role of the 3PN approximation is merely to “confirm” the value already given by the 2PN one (but of course, had we not computed the 3PN term, we would not be able to trust very much the 2PN value). As expected, the best agreement we obtain is for the 3PN approximation and in the case of corotation, i.e., the point $3\text{PN}^{\text{corot}}$. However, the 1PN approximation is clearly not precise enough, but this is not surprising in the highly relativistic regime of the ICO. The right panel of Figure 1 shows other very interesting comparisons with numerical relativity computations [133, 121], done not only for the case of corotational binaries but also in the irrotational (non-spinning) case. Witness in particular the almost perfect agreement between the standard 3PN point (PN standard, shown with a green triangle) and the numerical quasi-equilibrium point (QE, red triangle) in the case of irrotational non-spinning (NS) binaries.

However, we recall that the numerical works [228, 232, 133, 121] assume that the spatial metric is conformally flat, which is incompatible with the post-Newtonian approximation starting from the 2PN order (see [196] for a discussion). Nevertheless, the agreement found in Figure 1 constitutes an appreciable improvement of the previous situation, because the first estimations of the ICO in post-Newtonian theory [274] and numerical relativity [132, 342, 29] disagreed with each other, and do not match with the present 3PN results.

## 8.2 Dynamical stability of circular orbits

In this section, following Ref. [79], we shall investigate the problem of the stability, against dynamical perturbations, of circular orbits at the 3PN order. We propose to use two different methods, one based on a linear perturbation at the level of the center-of-mass equations of motion (219)–(220) in (standard) harmonic coordinates, the other one consisting of perturbing the Hamiltonian

<footnote>$^{56}$ In Ref. [51] it was assumed that the corotation condition was given by the leading-order result $\omega_a = \Omega$. The 1PN correction in Eq. (247) modifies the 3PN terms in Eq. (250) with respect to the result of Ref. [51].</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1:** The binding energy $E_{ICO}$ versus $\Omega_{ICO}$ in the equal-mass case ($\nu = 1/4$). Left panel: Comparison with the numerical relativity result of Gourgoulhon, Grandclément et al. [228, 232] valid in the corotating case (marked by a star). Points indicated by $nPN$ are computed from the minimum of Eq. (232), and correspond to irrotational binaries. Points denoted by $nPN^{corot}$ come from the minimum of the sum of Eqs. (232) and (250), and describe corotational binaries. Note the very good convergence of the standard (Taylor-expanded) PN series. Right panel: Numerical relativity results of Cook, Pfeiffer et al. [133, 121] for quasi-equilibrium (QE) configurations and various boundary conditions for the lapse function, in the non-spinning (NS), leading-order non spinning (LN) and corotating (CO) cases. The point from [228, 232] (HKV-GGB) is also reported as in the left panel, together with IVP, the initial value approach with effective potential [132, 342], as well as standard PN predictions from the left panel and non-standard (EOB) ones. The agreement between the QE computation and the standard non-resummed 3PN point is excellent especially in the irrotational NS case.</figcaption></figure>

equations in ADM coordinates for the center-of-mass Hamiltonian (223). We shall find a criterion for the stability of circular orbits and shall present it in an invariant way – the same in different coordinate systems. We shall check that our two methods agree on the result.

We deal first with the perturbation of the equations of motion, following Kidder, Will & Wise- man [275] (see their Section III.A). We introduce polar coordinates ($r, \varphi$) in the orbital plane and pose $u \equiv \dot{r}$ and $\Omega \equiv \dot{\varphi}$. Then Eq. (219) yields the system of equations

$$ \dot{r} = u, \qquad (251a) $$

$$ \dot{u} = -\frac{Gm}{r^2} [1 + A + Bu] + r\Omega^2, \qquad (251b) $$

$$ \dot{\Omega} = -\Omega \left[ \frac{Gm}{r^2} B + \frac{2u}{r} \right], \qquad (251c) $$

where $A$ and $B$ are given by Eqs. (220) as functions of $r, u$ and $\Omega$ (through $v^2 = u^2 + r^2\Omega^2$).

In the case of an orbit that is circular apart from the adiabatic inspiral at the 2.5PN order (we neglect the radiation-reaction damping effects), we have $\dot{r}_0 = \dot{u}_0 = \dot{\Omega}_0 = 0$ hence $u_0 = 0$. In this section we shall indicate quantities associated with the circular orbit, which constitutes the zero-th approximation in our perturbation scheme, using the subscript 0. Hence Eq. (251b) gives
the angular velocity $\Omega_0$ of the circular orbit as

$$ \Omega_0^2 = \frac{Gm}{r_0^3} (1 + A_0). \qquad (252) $$

Solving iteratively this relation at the 3PN order using the equations of motion (219)–(220), we obtain $\Omega_0$ as a function of the circular-orbit radius $r_0$ in standard harmonic coordinates; the result agrees with Eq. (228).$^{57}$

We now investigate the linear perturbation around the circular orbit defined by the constants $r_0$, $u_0 = 0$ and $\Omega_0$. We pose

$$ r = r_0 + \delta r, \qquad (253a) $$

$$ u = \delta u, \qquad (253b) $$

$$ \Omega = \Omega_0 + \delta\Omega, \qquad (253c) $$

where $\delta r$, $\delta u$ and $\delta \Omega$ denote the linear perturbations of the circular orbit. Then a system of linear equations readily follows:

$$ \dot{\delta}r = \delta u, \qquad (254a) $$

$$ \dot{\delta}u = \alpha_0 \delta r + \beta_0 \delta \Omega, \qquad (254b) $$

$$ \dot{\delta}\Omega = \gamma_0 \delta u, \qquad (254c) $$

where the coefficients, which solely depend on the unperturbed circular orbit (hence the added subscript 0), read as [275]

$$ \alpha_0 = 3\Omega_0^2 - \frac{Gm}{r_0^2} \left(\frac{\partial A}{\partial r}\right)_0, \qquad (255a) $$

$$ \beta_0 = 2r_0\Omega_0 - \frac{Gm}{r_0^2} \left(\frac{\partial A}{\partial \Omega}\right)_0, \qquad (255b) $$

$$ \gamma_0 = -\Omega_0 \left[ \frac{2}{r_0} + \frac{Gm}{r_0^2} \left( \frac{\partial B}{\partial u} \right)_0 \right]. \qquad (255c) $$

In obtaining these equations we use the fact that $A$ is a function of the square $u^2$ through $v^2 = u^2 + r^2\Omega^2$, so that $\partial A/\partial u$ is proportional to $u$ and thus vanishes in the unperturbed configuration (because $u = \delta u$). On the other hand, since the radiation reaction is neglected, $B$ is also proportional to $u$ [see Eq. (220b)], so only $\partial B/\partial u$ can contribute at the zero-th perturbative order. Now by examining the fate of perturbations that are proportional to some $e^{i\sigma t}$, we arrive at the condition for the frequency $\sigma$ of the perturbation to be real, and hence for stable circular orbits to exist, as being [275]

$$ \hat{C}_0 \equiv -\alpha_0 - \beta_0 \gamma_0 > 0. \qquad (256) $$

Substituting into this $A$ and $B$ at the 3PN order we then arrive at the orbital-stability criterion

$$ \hat{C}_0 = \frac{Gm}{r_0^3} \left\{ 1 + \frac{Gm}{r_0 c^2} (-9 + \nu) + \frac{G^2 m^2}{r_0^2 c^4} \left( 30 + \frac{65}{4}\nu + \nu^2 \right) + \frac{G^3 m^3}{r_0^3 c^6} \left( -70 + \left[ -\frac{29927}{840} - \frac{451}{64}\pi^2 + 22 \ln\left(\frac{r_0}{r_0'}\right) \right]\nu + \frac{19}{2}\nu^2 + \nu^3 \right) + \mathcal{O}\left(\frac{1}{c^8}\right) \right\}, \qquad (257) $$

<footnote>$^{57}$ One should not confuse the circular-orbit radius $r_0$ with the constant $r'_0$ entering the logarithm at the 3PN order in Eq. (228) and which is defined by Eq. (221).</footnote>
where we recall that $r_0$ is the radius of the orbit in harmonic coordinates.

Our second method is to use the Hamiltonian equations associated with the 3PN center-of-mass Hamiltonian in ADM coordinates $H^{ADM}$ given by Eq. (223). We introduce the polar coordinates $(R, \Psi)$ in the orbital plane – we assume that the orbital plane is equatorial, given by $\Theta = \frac{\pi}{2}$ in the spherical coordinate system $(R, \Theta, \Psi)$ – and make the substitution

$$P^2 = P_R^2 + \frac{P_\Psi^2}{R^2}. \qquad (258)$$

This yields a reduced Hamiltonian that is a function of $R$, $P_R$ and $P_\Psi$, and describes the motion in polar coordinates in the orbital plane; henceforth we denote it by $\mathcal{H} = \mathcal{H}[R, P_R, P_\Psi] \equiv H^{ADM}/\mu$. The Hamiltonian equations then read

$$\frac{dR}{dt} = \frac{\partial \mathcal{H}}{\partial P_R}, \qquad (259a)$$

$$\frac{d\Psi}{dt} = \frac{\partial \mathcal{H}}{\partial P_\Psi}, \qquad (259b)$$

$$\frac{dP_R}{dt} = -\frac{\partial \mathcal{H}}{\partial R}, \qquad (259c)$$

$$\frac{dP_\Psi}{dt} = 0. \qquad (259d)$$

Evidently the constant $P_\Psi$ is nothing but the conserved angular-momentum integral. For circular orbits we have $R = R_0$ (a constant) and $P_R = 0$, so

$$\frac{\partial \mathcal{H}}{\partial R} [R_0, 0, P_\Psi^0] = 0, \qquad (260)$$

which gives the angular momentum $P_\Psi^0$ of the circular orbit as a function of $R_0$, and

$$\Omega_0 = \left(\frac{d\Psi}{dt}\right)_0 = \frac{\partial \mathcal{H}}{\partial P_\Psi} [R_0, 0, P_\Psi^0], \qquad (261)$$

which yields the angular frequency of the circular orbit $\Omega_0$, which is evidently the same numerical quantity as in Eq. (252), but is here expressed in terms of the separation $R_0$ in ADM coordinates. The last equation, which is equivalent to $R = \text{const} = R_0$, is

$$\frac{\partial \mathcal{H}}{\partial P_R} [R_0, 0, P_\Psi^0] = 0. \qquad (262)$$

It is automatically verified because $\mathcal{H}$ is a quadratic function of $P_R$ and hence $\partial \mathcal{H}/\partial P_R$ is zero for circular orbits.

We consider now a perturbation of the circular orbit defined by

$$P_R = \delta P_R, \qquad (263a)$$

$$P_\Psi = P_\Psi^0 + \delta P_\Psi, \qquad (263b)$$

$$R = R_0 + \delta R, \qquad (263c)$$

$$\Omega = \Omega_0 + \delta\Omega. \qquad (263d)$$

The Hamiltonian equations (259), worked out at the linearized order, read as

$$\delta \dot{P}_R = -\pi_0 \delta R - \rho_0 \delta P_\Psi, \qquad (264a)$$