<figure><img src="image_1.png" /></figure>

HAL
open science

# Equations of motion of self-gravitating *N*-body systems in the first post-Minkowskian approximation

Luc Blanchet, Athanassios S. Fokas

► To cite this version:

Luc Blanchet, Athanassios S. Fokas. Equations of motion of self-gravitating *N*-body systems in the first post-Minkowskian approximation. Physical Review D, 2018, 98 (8), pp.084005. 10.1103/Phys-RevD.98.084005. hal-01903005

HAL Id: hal-01903005

https://hal.science/hal-01903005v1

Submitted on 2 Oct 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Equations of motion of self-gravitating *N*-body systems
in the first post-Minkowskian approximation

Luc Blanchet<fnref n="1" />,\* and Athanassios S. Fokas²,³,†

¹GReCO, Institut d'Astrophysique de Paris, UMR 7095, CNRS, Sorbonne Universités
& UPMC Université Paris 6, 98bis Boulevard Arago, 75014 Paris, France

²Department of Applied Mathematics and Theoretical Physics, University of Cambridge,
Cambridge, CB3 0WA, United Kingdom

³Viterbi School of Engineering, University of Southern California,
Los Angeles, California, 90089-2560, USA

<figure><img src="image_3.png" /></figure>

(Received 21 June 2018; published 8 October 2018)

We revisit the problem of the equations of motion of a system of *N* self-interacting massive particles (without spins) in the first post-Minkowskian (1PM) approximation of general relativity. We write the equations of motion, gravitational field, and associated conserved integrals of the motion in a form suitable for comparison with recently published post-Newtonian (PN) results at the 4PN order. We show that the Lagrangian associated with the equations of motion in harmonic coordinates is a generalized one, and we compute all the terms linear in *G* up to 5PN order. We discuss the Hamiltonian in the frame of the center of mass and exhibit a canonical transformation connecting it to previous results directly obtained with the Hamiltonian formalism of general relativity. Finally we recover the known result for the gravitational scattering angle of two particles at the 1PM order.

DOI: 10.1103/PhysRevD.98.084005

# I. INTRODUCTION

The post-Minkowskian (PM) expansion is one of the most important approximation schemes in general relativity. Assuming that the gravitational field is weak, it expands in nonlinearities around the Minkowski background, with small dimensionless expansion parameter $r_s/r$, where $r_s \sim Gm/c^2$ is the Schwarzschild radius of the source, and no restriction on the slowness parameter $v/c$.<fnref n="1" /> It has been developed starting from the 1916 Einstein paper (linearized gravity) [1], by many pioneering works from the 1960s on [2–9]. Notably the approximation scheme has been systematically investigated [4,10,11] and the equations of motion of systems of particles have been obtained up to 2PM order (quadratic in $G$) [5,8,12,13] (in particular, the work [13] introduced a regularization procedure for treating ultraviolet divergences due to point particles at 2PM order).

On the other hand, the post-Newtonian (PN) approximation, or expansion in powers of $v/c$ (the other very important approximation scheme in general relativity), has played a crucial role in the data analysis of the recent gravitational wave events. While the gravitational wave-form generated by compact (black hole or neutron star)

binary systems is known to 3.5PN order beyond the quadrupole formula [14], the equations of motion have been derived to the 4PN order by means of two independent approaches: the canonical Hamiltonian formalism in Arnowitt-Deser-Misner (ADM) coordinates [15–20] and the Fokker Lagrangian in harmonic coordinates [21–25]. In addition, the effective field theory approach has obtained partial results at 4PN order (with missing $G^4$ terms in Newton’s constant) [26–29].

The most well-known physical application of the PM approximation is for unbound orbits, i.e., the scattering problem with high relative velocities and small deviation angle. The gravitational scattering of two particles has been solved at the 2PM order [5,7,8,12]. More recently, Ledvinka et al. [30] published a closed-form expression for the Hamiltonian of *N* particles in the 1PM approximation. Very recently new works appeared on the gravitational scattering angle, the Hamiltonian of two point-particles at the 2PM order, and the link with the effective-one-body formalism [31,32].

The aim of the present paper is to revisit the PM approximation in order to check some terms in the recent PN derivations of the equations of motion of compact binary systems in harmonic coordinates. For this purpose we restrict attention to the 1PM approximation, keeping only linear terms in $G$. Of course, the 1PM approximation is well known, but we find that some work is necessary in order to make possible the comparison with PN results.

<footnote>\*luc.blanchet@iap.fr</footnote>

<footnote>†tf227@cam.ac.uk</footnote>

<footnote n="1">¹Here, *r* denotes the typical size and *v* the typical velocity of the source. Henceforth, small *n*PM terms ~(*r*<sub>s</sub>/*r*)<sup>*n*</sup> will be said to be of order *G*<sup>*n*</sup> in Newton’s constant.</footnote>
In particular, we need to explicitly perform the expansion of retardations consistently with the approximation, so as to recast the 1PM equations of motion into a PN-like form (valid up to any PN order but order linear in G only). In addition we want to provide closed-form expressions for the associated conserved quantities like energy.

Concerning the Lagrangian, the problem is not so straightforward because in harmonic coordinates the Lagrangian is a generalized one, depending not only on positions and velocities but also on the accelerations of the bodies. We shall not find a closed-form expression for the harmonic-coordinates Lagrangian at the 1PM order in the general case, but propose a method to obtain it at any PN order (and linear order in G). In the two-body case with equal masses we give explicit expressions for the Lagrangian and Hamiltonian, and we shall connect, via a canonical transformation, the Hamiltonian to the one obtained in [30]. Finally we have been able to check all the terms linear in G in many of the recently obtained 4PN results [21–25].

The plan of this paper is as follows. In Sec. II we present the derivation of the 1PM gravitational field and equations of motion of N particles in a PN-like form. In Sec. III we obtain the conserved energy and linear momentum in closed analytic form. The special case of two equal masses is investigated in Sec. IV. We investigate the (acceleration-dependent) Lagrangian in harmonic coordinates in Sec. V and present new formulas for computing it at any PN order. We discuss the corresponding Hamiltonian formalism in Sec. VI. Further material is relegated to some Appendices: the scattering of two particles at 1PM order in Appendix A (we recover the known result), the rather academic case of bounded circular orbits and their stability in Appendix B,<fnref n="2" /> and the terms of order G in the harmonic-coordinates Lagrangian up to 5PN order in Appendix C.

## II. DERIVATION OF THE EQUATIONS OF MOTION

In this paper we restrict attention to the linearized Einstein field equations in harmonic coordinates. Thus, our theory is based on the linearized Einstein-Hilbert action of general relativity and on the standard linearized harmonic gauge-fixing term. We denote by $h^{\mu\nu} = \sqrt{-gg^{\mu\nu}} - \eta^{\mu\nu}$ the “gothic” metric deviation from Minkowski’s metric, where $g^{\mu\nu}$ is the inverse of the usual covariant metric $g_{\rho\sigma}$ and $g = \det(g_{\rho\sigma})$. Furthermore, we raise and lower indices using the Minkowski metric, for instance, $h = \eta^{\mu\nu}h_{\mu\nu}$. Thus, we have the action

$$S = \int \frac{d^4 x}{c} \left[ \frac{c^4}{64\pi G} h^{\mu\nu} \left( \eta_{\mu(\rho}\eta_{\sigma)\nu} - \frac{1}{2}\eta_{\mu\nu}\eta_{\rho\sigma} \right) \Box h^{\rho\sigma} \right. \\ \left. - \frac{1}{2} \left( h_{\mu\nu} - \frac{1}{2}h \right) T^{\mu\nu} \right]. \qquad (2.1)$$

Here $\Box$ is the flat space-time d’Alembertian operator, and $T^{\mu\nu}$ is the stress-energy tensor of matter fields in special relativity, corresponding to a system of N point particles without spins, namely,

$$T^{\mu\nu}(x) = \sum_{a=1}^{N} m_a c^2 \int_{-\infty}^{+\infty} d\tau_a u_a^{\mu} u_a^{\nu} \delta^{(4)}(x - y_a). \quad (2.2)$$

The sum runs over the N particles; $\delta^{(4)}$ is the four-dimensional Dirac distribution; $m_a$’s are the constant masses of particles; the worldlines of the particles are denoted by $y_a^\mu$; and the particles’ special-relativistic proper time is $\tau_a$, so that $d\tau_a = \sqrt{-\eta_{\mu\nu}dy_a^\mu dy_a^\nu/c^2}$. The dimensionless four-velocities are defined by $cu_a^\mu = dy_a^\mu/d\tau_a$ and normalized to $\eta_{\mu\nu}u_a^\mu u_a^\nu = -1$. The linearized harmonic gauge-fixed Einstein field equations deduced from the action (2.1) are given by

$$\Box h^{\mu\nu} = \frac{16\pi G}{c^4} T^{\mu\nu}. \qquad (2.3)$$

We solve the latter equations for the case of the explicit stress-energy tensor of particles given by (2.2) by means of the standard Lienard-Wiechert integration procedure [33]: adopting a parametrization of the trajectories by means of the coordinate time $t$, i.e., such that $y_a^\mu(t) = (ct, y_a(t))$, the retarded time $t_a^{\text{ret}}$ on the trajectory $a$, which is associated with the propagation from the particle $a$ to the field point $x^\mu = (ct, \mathbf{x})$, is given by the implicit equation

$$t_a^{\text{ret}}(\mathbf{x}, t) = t - \frac{1}{c} |\mathbf{x} - \mathbf{y}_a(t_a^{\text{ret}}(\mathbf{x}, t))|, \quad (2.4)$$

where $|\cdot|$ denotes the usual Euclidean norm. Denoting by $r_a^{\text{ret}} = |\mathbf{x} - \mathbf{y}_a(t_a^{\text{ret}})|$ the spatial interval, and by $\mathbf{n}_a^{\text{ret}} = (\mathbf{x} - \mathbf{y}_a(t_a^{\text{ret}}))/r_a^{\text{ret}}$ the unit spatial direction of propagation from the source point to the field point, the retarded solution of Eq. (2.3) is given by

$$h^{\mu\nu}(\mathbf{x}, t) = -\frac{4}{c^2} \sum_a \frac{G m_a u_a^\mu u_a^\nu}{r_a^{\text{ret}} \gamma_a (1 - \mathbf{n}_a^{\text{ret}} \cdot \mathbf{v}_a / c)} \cdot \quad (2.5)$$

The factor

$$k_{\mu}^{a} u_{a}^{\mu} = \gamma_{a} \left( -1 + \frac{\mathbf{n}_{a}^{\text{ret}} \cdot \mathbf{v}_{a}}{c} \right), \qquad (2.6)$$

which occurs in the denominator of (2.5), is the standard redshift factor, where $k_a^\mu$ is the Minkowski null vector

<footnote n="2">²We keep in mind that in the case of bound orbits, depending on the scale, the 1PM approximation may miss crucial physical terms present in GR, of order $v/c$ to the power 4 or higher.</footnote>
between source and field points, given by $k_a^\mu = [x^\mu - y_a^\mu(t_a^{\text{ret}})]/r_a^{\text{ret}} = (1, n_a^{\text{ret}})$, where $\gamma_a = u_a^0$ is the usual Lorentz factor, and the ordinary coordinate velocity is defined by $v_a^\mu = cu_a^\mu/\gamma_a = (c, v_a)$, with $n_a^{\text{ret}} \cdot v_a$ denoting the ordinary Euclidean scalar product.

The velocities in Eq. (2.5) should also be computed at the retarded time $t_a^{\text{ret}}$. However, here we shall repeatedly use the fact that the accelerations are of order $G$. Therefore their contributions in Eq. (2.5) will be of order $G^2$ and hence can be neglected within the present approximation, which is strictly confined to the linear order in $G$. Thus, in all our computations we can always assume that the four-velocity $u_a^\mu$ and Lorentz factor $\gamma_a$ are constant. Furthermore, neglecting terms of order $G^2$, we can solve the retardation equation (2.4), to reexpress the retarded time $r_a^{\text{ret}}$, distance $r_a^{\text{ret}}$, and direction $n_a^{\text{ret}}$, in terms of their values at the current time $t$, i.e., the “instantaneous” distance $r_a = |\mathbf{x} - \mathbf{y}_a(t)|$ and direction $\mathbf{n}_a = [\mathbf{x} - \mathbf{y}_a(t)]/r_a$. In this way, we find

$$r_a^{\text{ret}} = \gamma_a r_a \left[ \sqrt{1 + (n_a u_a)^2} + (n_a u_a) \right], \quad (2.7a)$$

$$\mathbf{n}_a^{\text{ret}} = \frac{\mathbf{v}_a}{c} + \frac{\mathbf{n}_a}{\gamma_a} \left[ \sqrt{1 + (n_a u_a)^2} - (n_a u_a) \right]. \quad (2.7b)$$

In what follows, for convenience of notation we use parentheses to denote the ordinary scalar product, i.e., $(n_a u_a) = \mathbf{n}_a \cdot \mathbf{u}_a = \gamma_a(n_a v_a)/c$. Substituting these formulas into Eq. (2.5) we find that the equivalent instantaneous expression of the field takes the form

$$h^{\mu\nu} = -\frac{4}{c^2} \sum_a \frac{Gm_a u_a^\mu u_a^\nu}{r_a \sqrt{1 + (n_a u_a)^2}}. \quad (2.8)$$

Next, in preparation for the derivation of the geodesic equation, we compute the space-time derivatives of the field, which can be easily achieved using the explicit expression (2.8). Recalling that within the present approximation the four-velocities are constant, we find

$$\partial_{\rho} h^{\mu\nu} = \frac{4}{c^2} \sum_a \frac{Gm_a u_a^\mu u_a^\nu}{r_a^2 [1 + (n_a u_a)^2]^{3/2}} (n_\rho^a + u_\rho^a (n_a u_a)). \quad (2.9)$$

Here we have extended the definition of the unit direction $\mathbf{n}_a$ by posing $n_a^\mu = (0, \mathbf{n}_a)$; also, we use $u_a^\mu = \gamma_a(1, v_a/c)$. An immediate check of the result (2.9) is that the divergence of $h^{\mu\nu}$ is seen to be zero: $H^\mu = \partial_\nu h^{\mu\nu} = 0$. This is consistent with the harmonic-coordinates field equations (2.3) together with the conservation law of the matter tensor, which is identical with the conservation law in special relativity, namely, $\partial_\nu T^{\mu\nu} = 0$.

The expressions (2.8) and (2.9) are valid everywhere except at the singular location of the particles. Nevertheless, we can easily extend the validity of these equations by using a self-field regularization, which deals

with the infinite self-field of particles. For the purpose, it is sufficient to simply discard the self-field contribution from the sum of the particles.<fnref n="3" /> Therefore, at the particle $a$, the expression (2.9) becomes

$$(\partial_{\rho} h^{\mu\nu})_a = \frac{4}{c^2} \sum_{b \neq a} \frac{Gm_b u_b^{\mu} u_b^{\nu}}{r_{ab}^2 [1 + (n_{ab} u_b)^2]^{3/2}} (n_{\rho}^{ab} + u_{\rho}^{ab} (n_{ab} u_b)), \quad (2.10)$$

where now the sum runs over all particles except $a$, and where $r_{ab} = |\mathbf{y}_a - \mathbf{y}_b|$ and $n_{ab}^0 = 0$, $n_{ab}^i = [\mathbf{y}_a - \mathbf{y}_b]/r_{ab}$. The equation of motion is obtained by just inserting the latter expression into the geodesic equation. At linearized order, which is consistent with our approximation, it is given by

$$\frac{d u_a^\mu}{d \tau_a} = \frac{c}{2} u_a^\rho u_a^\sigma [(\partial^\mu k_{\rho\sigma})_a - 2(\partial_\rho k_\sigma^\mu)_a - u_a^\mu u_a^\lambda (\partial_\lambda k_{\rho\sigma})_a], \quad (2.11)$$

where we have used the trace-reversed variable $k^{\mu\nu} = \frac{1}{2}\eta^{\mu\nu}h - h^{\mu\nu}$ that represents the linear perturbation in the ordinary covariant metric $g_{\mu\nu}$. Hence, we find that the complete equation of motion at the 1PM approximation is given by

$$\frac{d u_a^\mu}{d \tau_a} = -\frac{1}{c} \sum_{b \neq a} \frac{Gm_b}{r_{ab}^2 [1 + (n_{ab} u_b)^2]^{3/2}} \times \left[ (2\epsilon_{ab}^2 - 1) n_{ab}^\mu \right. \\ \left. \qquad + (2\epsilon_{ab}^2 + 1)(-(n_{ab} u_a) + \epsilon_{ab}(n_{ab} u_b)) u_a^\mu \right. \\ \left. \qquad + (4\epsilon_{ab}(n_{ab} u_a) - (2\epsilon_{ab}^2 + 1)(n_{ab} u_b)) u_b^\mu \right]. \quad (2.12)$$

We have used $\epsilon_{ab} = -(u_a u_b)$ as a shorthand notation<fnref n="4" />; parentheses represent four-dimensional scalar products, for instance, $(n_{ab} u_a) = n_{ab}^\mu u_a^\mu = \gamma_a \mathbf{n}_{ab} \cdot \mathbf{v}_a/c$ and $\epsilon_{ab} = \gamma_a \gamma_b (1 - \mathbf{v}_a \cdot \mathbf{v}_b/c^2)$. We often use parentheses to represent both three-dimensional and four-dimensional scalar products, depending on the context. Using (2.12) we can verify that $(u_a \frac{d u_a}{d \tau_a}) = 0$, which is, of course, the consequence of the normalization $(u_a u_a) = -1$.

Finally, we present the noncovariant form of the equations of motion, introducing the ordinary velocity and

$$\epsilon_{ab} = \frac{s_{ab} - m_a^2 - m_b^2}{2m_a m_b} = - \frac{(p_a p_b)}{m_a m_b c^2},$$

<footnote n="3">³See [13] for an explicit self-field regularization procedure applied in the PM approximation.</footnote>

<footnote n="4">⁴This notation agrees with the energy function used for two-body problems in particle physics, in particular,</footnote>

<footnote>where $s_{ab} = -(p_a + p_b)^2/c^2$ is the Mandelstam variable; see, e.g., [9,31,32].</footnote>
acceleration fields and the relevant Lorentz factors. It turns out that this form is relatively simple,

$$ \frac{dv_a}{dt} = -\gamma_a^{-2} \sum_{b \neq a} \frac{Gm_b}{r_{ab}^2 y_{ab}^{3/2}} \left[ (2\epsilon_{ab}^2 - 1)\mathbf{n}_{ab} + \gamma_b \left( -4\epsilon_{ab}\gamma_a(n_{ab}v_a) + (2\epsilon_{ab}^2 + 1)\gamma_b(n_{ab}v_b) \right) \frac{\mathbf{v}_{ab}}{c^2} \right], \quad (2.13) $$

where for convenience of notation we use $y_{ab} = 1 + \gamma_b^2(n_{ab}v_b)^2/c^2$ (notice that $y_{ab} \neq y_{ba}$). We note that the second term is proportional to the relative velocity $\mathbf{v}_{ab} = \mathbf{v}_a - \mathbf{v}_b$ and that there are no terms proportional to individual velocities $\mathbf{v}_a$ or $\mathbf{v}_b$.

We have verified that the formula (2.13) exactly reproduces the post-Newtonian results in the case of two particles ($N=2$) at linear order in $G$ and up to the 4PN order [23,25].

### III. CONSERVED ENERGY AND LINEAR MOMENTUM

The equations of motion at the 1PM approximation, i.e., Eqs. (2.12) and (2.13) are conservative and thus admit conserved integrals of energy, angular momentum, and linear momentum. Indeed, the gravitational radiation reaction dissipative effects in the equations of motion are at least of second order in $G$; see, e.g., [14]. To find the conserved energy we proceed in a standard way; namely, we form the combination $W = \sum_a m_a \gamma_a^3 (v_a a_a)$, where $a_a$ denotes the acceleration (2.13). By definition, $W$ is the time derivative of the special-relativistic total energy, $W = \frac{d}{dt}[\sum_a m_a \gamma_a c^2]$. Also, after replacing $a_a$ via Eq. (2.13), we need to rewrite the resulting expression in the form of the total time derivative of some potential up to order $G$, say, $W = -\frac{dV}{dt}$. Then, the conserved energy $E$, satisfying $\frac{dE}{dt} = 0$, is given by

$$ E = \sum_a m_a \gamma_a c^2 + V. \quad (3.1) $$

We find this energy in two steps. First, a part of the terms in $W$ can be readily integrated thanks to the easily checked identity

$$ \frac{d}{dt} \left( \frac{1}{r_{ab} y_{ab}^{1/2}} \right) = \frac{\frac{\gamma_b}{\gamma_a} \epsilon_{ab} (n_{ab} v_b) - (n_{ab} v_a)}{r_{ab}^2 y_{ab}^{3/2}}, \quad (3.2) $$

where $y_{ab}$ was defined via the equation $y_{ab} = 1 + \gamma_b^2(n_{ab}v_b)^2/c^2$ and we recall that $\epsilon_{ab} = \gamma_a \gamma_b [1 - (v_a v_b)/c^2]$. In this way, several terms occurring in the potential $V$ are obtained,

$$ V = \sum_a \sum_{b \neq a} G m_a m_b \left[ \gamma_a \frac{2\epsilon_{ab}^2 + 1 - 4\frac{\gamma_b}{\gamma_a} \epsilon_{ab}}{r_{ab} y_{ab}^{1/2}} + \frac{\gamma_b^2}{\gamma_a} (2\epsilon_{ab}^2 - 1)(D_{ab} v_b) \right]. \quad (3.3) $$

However, we find that the last term corresponds to a more complicated structure of the energy and is given by the scalar product $(D_{ab}v_b) = \mathbf{D}_{ab} \cdot \mathbf{v}_b = D_{ab}^i v_b^i$, between the velocity $v_b^i$ (which is constant) and some elementary solution of the equation

$$ \frac{dD_{ab}^i}{dt} = \frac{n_{ab}^i}{r_{ab}^2 y_{ab}^{3/2}}. \quad (3.4) $$

We have obtained an analytic closed-form solution of Eq. (3.4) using the following approach: recalling that within our approximation the velocities are constant and trajectories are straight lines, we can always obtain the unit direction vector $n_{ab}^i = (y_a^i - y_b^i)/r_{ab}$ as a simple function of the relative distance $r_{ab}$, namely,

$$ n_{ab}^i = \frac{v_{ab}^i}{v_{ab}^2} \sqrt{v_{ab}^2 - \frac{C_{ab}^2}{r_{ab}^2}} + \frac{\mu_{ab}^i}{r_{ab}}, \quad (3.5) $$

where $C_{ab}$ and the three-vector $\mu_{ab} = (\mu_{ab}^i)$ are constants of the motion. The first constant, $C_{ab}$, is nothing but the specific angular momentum (i.e., the angular momentum per unit mass) of the relative motion of the two particles $a$ and $b$; the second constant, $\mu_{ab}$, is related to the initial positions at time $t=0$ of the two particles on their straight-line trajectories.<fnref n="5" /> We note that $C_{ab} = C_{ba}$ whereas $\mu_{ab} = -\mu_{ba}$. From Eq. (3.5) we obtain the scalar products $(n_{ab}v_a)$ and $(n_{ab}v_b)$ as

$$ (n_{ab} v_a) = \frac{(v_{ab} v_a)}{v_{ab}^2} \sqrt{v_{ab}^2 - \frac{C_{ab}^2}{r_{ab}^2}} + \frac{(\mu_{ab} v_a)}{r_{ab}}, \quad (3.6) $$

together with $a \leftrightarrow b$. We recall our notations $(v_{ab}v_a) = v_a^2 - (v_a v_b)$, $(v_{ab}v_b) = -v_b^2 + (v_a v_b)$, and $v_{ab}^2 = v_a^2 - 2(v_a v_b) + v_b^2$. Next we obtain, using Eq. (3.6), the time derivative of the relative distance $r_{ab}$ as

$$ \frac{dr_{ab}}{dt} = (n_{ab} v_a) - (n_{ab} v_b) = \sqrt{v_{ab}^2 - \frac{C_{ab}^2}{r_{ab}^2}}. \quad (3.7) $$

$$ \mu_{ab} = \mathbf{y}_{ab}^0 - \frac{(\mathbf{y}_{ab}^0 v_{ab}) \mathbf{y}_{ab}}{v_{ab}^2}. $$

<footnote n="5">⁵More precisely, let the trajectories be $\mathbf{y}_a = \mathbf{v}_a t + \mathbf{y}_a^0$, and pose $\mathbf{y}_{ab}^0 = \mathbf{y}_a^0 - \mathbf{y}_b^0$ for the initial relative position at $t=0$. Then we have</footnote>

<footnote>Note that this vector is perpendicular to the relative velocity: $(\mu_{ab} v_{ab}) = 0$.</footnote>
The above formulas are valid in that portion of the trajectories for which $r_{ab}$ increases. Equation (3.7) gives us a one-to-one relationship between time $t$ and the distance $r_{ab}$, and it allows one to perform a change of variable in Eq. (3.4) from the variable $t$ to the distance $r_{ab}$ itself, and hence we obtain the ordinary differential equation

$$ \frac{dD_{ab}^i}{dr_{ab}} = \frac{\frac{v_{ab}^i}{v_{ab}^i} + \frac{\mu_{ab}^i}{\sqrt{v_{ab}^2 r_{ab}^2 - C_{ab}^2}}}{r_{ab}^2 \left( 1 + \frac{\gamma_b^2}{c^2} \left[ \frac{(v_{ab}v_b)}{v_{ab}^2} \sqrt{v_{ab}^2 - \frac{C_{ab}^2}{r_{ab}^2} + \frac{(\mu_{ab}v_b)}{r_{ab}}} \right]^2 \right)^{3/2}} \quad (3.8) $$

This equation can be integrated in closed analytic form:

$$ D_{ab}^i = \frac{1}{C_{ab}^2 + \frac{\gamma_b^2}{c^2} \left[ (\mu_{ab} v_b)^2 v_{ab}^2 + C_{ab}^2 \frac{(v_{ab} v_b)^2}{v_{ab}^2} \right]} \times \frac{\mu_{ab}^i \sqrt{v_{ab}^2 r_{ab}^2 - C_{ab}^2 - C_{ab}^2 \frac{v_{ab}^i}{v_{ab}^2}}}{r_{ab} \left( 1 + \frac{\gamma_b^2}{c^2} \left[ \frac{(v_{ab} v_b)}{v_{ab}^2} \sqrt{v_{ab}^2 - \frac{C_{ab}^2}{r_{ab}^2} + \frac{(v_{ab} v_b)}{r_{ab}}} \right]^2 \right)^{1/2}} \quad (3.9) $$

We note that the first factor in the above equation, within our approximation, is a pure constant. Finally, this form as it stands is not very useful, and we express the constants $C_{ab}$ and $\mu_{ab}^i$ in terms of their earlier expressions as functions of the original variables $n_{ab}^i$, $v_{ab}^i$; the scalar products $(n_{ab}v_a)$ and $(n_{ab}v_b)$; and $(n_{ab}v_{ab}) = \dot{r}_{ab} \equiv dr_{ab}/dt$. This gives the following interesting structure for the most “complicated” term in the energy:

$$ D_{ab}^i = \frac{[\dot{r}_{ab} + \frac{\gamma_b^2}{c^2}(n_{ab}v_b)(v_{ab}v_b)]n_{ab}^i - y_{ab}v_{ab}^i}{r_{ab}y_{ab}^{1/2}[(v_{ab}^2 - \dot{r}_{ab}^2)y_{ab} + \frac{\gamma_b^2}{c^2}(\dot{r}_{ab}(n_{ab}v_b) - (v_{ab}v_b))^2]'} \quad (3.10) $$

plus an integration constant which can be absorbed into the definition of $E$. Hence, combining Eqs. (3.3) and (3.10), we have a closed-form expression for the conserved energy (3.1), thoroughly given by the expression for the potential $V$ at order $G$:

$$ V = \sum_a \sum_{b \neq a} \frac{G m_a m_b}{r_{ab} y_{ab}^{1/2}} \left\{ \begin{aligned} & \gamma_a \left( 2 \epsilon_{ab}^2 + 1 - 4 \frac{\gamma_b}{\gamma_a} \epsilon_{ab} \right) \\ & + \frac{\gamma_b^2}{\gamma_a} (2 \epsilon_{ab}^2 - 1) \\ & \times \frac{\dot{r}_{ab} (n_{ab} v_b) - (v_{ab} v_b)}{(v_{ab}^2 - \dot{r}_{ab}^2) y_{ab} + \frac{\gamma_b^2}{c^2} (\dot{r}_{ab} (n_{ab} v_b) - (v_{ab} v_b))^2} \end{aligned} \right\} \quad (3.11) $$

We have verified that this energy reproduces in the post-Newtonian limit for two particles the 4PN results up to order $G$ [23,25].

We apply next the same method to the integral of linear momentum. Namely, we start by forming the combination $\frac{d}{dt}[\sum_a m_a \gamma_a \mathbf{v}_a] \equiv \sum_a m_a \gamma_a [\mathbf{a}_a + \gamma_a^2 (v_a a_a) \mathbf{v}_a / c^2]$, since $\sum_a m_a \gamma_a \mathbf{v}_a$ represents the total linear momentum in special relativity. Then, we replace the accelerations by the equations of motion and are able to transform the result into a total time derivative thanks to the integration formulas (3.2) and, most importantly, (3.10). As a result we obtain the conserved total linear momentum, such that $d\mathbf{P}/dt = \mathbf{0}$, as

$$ P^i = \sum_a m_a \gamma_a v_a^i + \Pi^i, \quad (3.12) $$

in which the terms of order $G$ are given by

$$ \Pi^i = \sum_a \sum_{b \neq a} \frac{G m_a m_b}{r_{ab} y_{ab}^{1/2}} \left\{ \begin{aligned} & (2 \epsilon_{ab}^2 + 1) \gamma_a \frac{v_a^i}{c^2} - 4 \epsilon_{ab} \gamma_b \frac{v_b^i}{c^2} \\ & + \frac{2 \epsilon_{ab}^2 - 1}{\gamma_a} \left( \delta^{ij} + \frac{\gamma_b^2}{c^2} v_b^i v_b^j \right) \\ & \times \frac{[\dot{r}_{ab} + \frac{\gamma_b^2}{c^2} (n_{ab} v_b) (v_{ab} v_b)] n_{ab}^j - y_{ab} v_{ab}^j}{(v_{ab}^2 - \dot{r}_{ab}^2) y_{ab} + \frac{\gamma_b^2}{c^2} (\dot{r}_{ab} (n_{ab} v_b) - (v_{ab} v_b))^2} \end{aligned} \right\} \quad (3.13) $$

As is well known, the main physical situation meaningfully described by the PM approximation is the scattering of particles. The scattering of two particles has been worked out long ago up to the 2PM order [5,7,8]; see Eq. (4.78) in [8]. In Appendix A we shall check that the previous formalism, in particular the integration formulas (3.9)–(3.10), correctly recovers the known result for the scattering of two particles at the 1PM order.

## IV. THE TWO-BODY CASE WITH EQUAL MASSES

In this section we discuss the equations of motion in the case of two bodies ($N=2$) with equal masses. For equal masses we necessarily have $v_1^i = -v_2^i$ in the frame of the center of mass, which is an important simplification of the general case. We denote

$$ m_1 = m_2 = \frac{m}{2}, \quad \mathbf{v}_1^i = -\mathbf{v}_2^i = \frac{\mathbf{v}}{2}, \quad \text{and} \quad \gamma_1 = \gamma_2 = \gamma, \quad (4.1) $$

where $m = m_1 + m_2$ is the total mass, $\mathbf{v} = \mathbf{v}_1 - \mathbf{v}_2$ is the relative velocity, and $\gamma = (1 - \frac{v^2}{4c^2})^{-1/2}$ is the common Lorentz factor of the two bodies. Notice that $\epsilon_{12} = 2\gamma^2 - 1$. We also pose $r = r_{12} = |\mathbf{y}_1 - \mathbf{y}_2|$ and $\mathbf{n} = \mathbf{n}_{12} = (\mathbf{y}_1 - \mathbf{y}_2)/r$, such that $\dot{r} = (nv)$. We also denote,
consistently with our previous notation in Sec. II, $y = 1 + \frac{\gamma^2(nv)^2}{4c^2}$ (note that $y = y_{12} = y_{21}$ for equal masses).

Inserting (4.1) in the basic equations of motion (2.12)–(2.13), it follows that the motion of two equal masses in the 1PM approximation is characterized by the single equation<fnref n="6" />

$$ \frac{d}{d\tau}(\gamma\mathbf{v}) = -\frac{Gm}{r^2 y^{3/2}} \left[ (8\gamma^4 - 8\gamma^2 + 1)\mathbf{n} - \frac{\gamma^2}{4}(16\gamma^6 - 8\gamma^4 + 6\gamma^2 - 1)\frac{\dot{\mathbf{r}}\mathbf{v}}{c^2} \right]. \quad (4.2) $$

Recalling the definition of the force, we obtain

$$ f = m \frac{d}{dt}(\gamma \mathbf{v}) = -\frac{Gm^2}{r^2 y^{3/2}} \left[ (8\gamma^3 - 8\gamma + \gamma^{-1})\mathbf{n} - \frac{1}{4}(16\gamma^7 - 8\gamma^5 + 6\gamma^3 - \gamma)\frac{\dot{\mathbf{r}}\mathbf{v}}{c^2} \right]. \quad (4.3) $$

Of course, in the small velocity limit, we recover the form of the usual Newtonian gravitational force $f_N = -\frac{Gm^2}{r^2} \mathbf{n}$.

Next we derive the evolution equation for $\gamma^2$. We use the identity $\frac{d}{dt}(\gamma\mathbf{v} \cdot \gamma\mathbf{v}) = 2\gamma\mathbf{v} \cdot \frac{d}{dt}(\gamma\mathbf{v})$, together with $\gamma\mathbf{v} \cdot \gamma\mathbf{v} = \gamma^2 v^2 = 4c^2(\gamma^2 - 1)$, to arrive at

$$ \frac{d\gamma^2}{dt} = \frac{8Gm\gamma^2}{r^2 y^{3/2} c^2} \dot{\mathbf{r}} \left( \gamma^6 - \frac{3}{2}\gamma^4 + \frac{3}{8}\gamma^2 + \frac{1}{16} \right), \quad (4.4) $$

or equivalently, changing the variable $t \to r$,

$$ \frac{d\gamma^2}{dr} = \frac{8Gm}{r^2 y^{3/2} c^2} \left( \gamma^8 - \frac{3}{2}\gamma^6 + \frac{3}{8}\gamma^4 + \frac{\gamma^2}{16} \right). \quad (4.5) $$

Similarly we derive also the evolution equation for $y$. We consider the identity $\frac{d}{dt}(\gamma\dot{\mathbf{r}}) = \frac{dn}{dt} \cdot \gamma\mathbf{v} + \mathbf{n} \cdot \frac{d}{dt}(\gamma\mathbf{v})$ and the fact that $v^2 = \dot{\mathbf{r}}^2 + r\mathbf{v} \cdot \frac{dn}{dt}$, which gives, after insertion of the equation of motion (4.2),

$$ \frac{dy}{dt} = \dot{\mathbf{r}} \left\{ \frac{2}{r}(\gamma^2 - y) + \frac{8Gm}{r^2 c^2 y^{1/2}} \left( \gamma^6 - \frac{\gamma^4}{2} + \frac{3\gamma^2}{8} - \frac{1}{16} \right) - \frac{8Gm}{r^2 c^2 y^{3/2}} \left( \gamma^6 - \frac{\gamma^2}{8} \right) \right\}. \quad (4.6) $$

Hence, again with changing variable from $t$ to $r$, we obtain

$$ \frac{dv^i}{dt} = - \frac{Gm}{r^2 y^{3/2}} \left[ \left( 8\gamma^2 - 8 + \frac{1}{\gamma^2} \right) n^i - \left( 4\gamma^4 - \frac{1}{2} \right) (nv) \frac{v^i}{c^2} \right]. $$

$$ \begin{aligned} \frac{dy}{dr} + \frac{2}{r}y - \frac{2}{r}\gamma^2 &= \frac{8Gm}{r^2 c^2 y^{3/2}} \left[ y \left( \gamma^6 - \frac{\gamma^4}{2} + \frac{3}{8}\gamma^2 - \frac{1}{16} \right) - \left( \gamma^6 - \frac{\gamma^2}{8} \right) \right]. \end{aligned} \quad (4.7) $$

Equations (4.5) and (4.7) form a system of two coupled ordinary differential equations which determine the two unknown functions $\gamma$ and $y$ in terms of $r$. Actually, it will be shown below that it is possible to express explicitly $y$ in terms of $\gamma$ and $r$. Before deriving this expression we will derive the expressions for the conserved energy and conserved angular momentum.

The energy $E$ is defined by (3.1); hence $E = mc^2\gamma + V$ for the two-equal-masses case. Hence, conservation of energy implies $\frac{dV}{dr} = -mc^2\frac{d\gamma}{dr}$. Thus, using Eq. (4.5) we find

$$ \frac{dV}{dr} = -\frac{8Gm^2}{r^2 y^{3/2}} \left( \gamma^7 - \frac{3}{2}\gamma^5 + \frac{3}{8}\gamma^3 + \frac{\gamma}{16} \right). \quad (4.8) $$

In order to integrate the above equation we first observe that since $d\gamma/dt = O(G)$, $\gamma$ can be treated as constant. Furthermore, we will employ the identity

$$ \frac{d}{dr} \left( \frac{1}{ry^{1/2}} \right) = - \frac{\gamma^2}{r^2 y^{3/2}} + O(G), \quad (4.9) $$

which is a consequence of Eq. (4.7) when neglecting terms of order $G$. The integration of Eq. (4.8) is then straightforward and we obtain

$$ V = \frac{Gm^2}{ry^{1/2}} (2\gamma^2 - 1) \left( 2\gamma^3 - 2\gamma - \frac{1}{4\gamma} \right). \quad (4.10) $$

This result is perfectly consistent with the two-equal-mass case of the general analysis done in the previous section; see Eq. (3.11).

The angular momentum $J$ is given by the usual special-relativistic expression plus a 1PM correction $K$ of order $G$, i.e., $J = m_1\gamma_1 y_1 \times v_1 + m_2\gamma_2 y_2 \times v_2 + K$, which gives for the case at hand (denoting $C = x \times v$)

$$ J = \frac{m\gamma}{4} C + K. \quad (4.11) $$

Thus $\frac{dK}{dt} = -\frac{m}{4}\frac{d}{dt}(\gamma C) = -\frac{1}{4}\mathbf{x} \times \mathbf{f}$, and then replacing the force by the right-hand side (rhs) of Eq. (4.3) we obtain

$$ \frac{dK}{dt} = -C \frac{Gm^2}{c^2 r^2 y^{3/2}} \left( 2\gamma^7 - \gamma^5 + 2\gamma^3 - \frac{\gamma}{8} \right). $$

We can integrate this equation treating both $\gamma$ and $C$ as constants to order $G$ and noticing that $C = 2cr\sqrt{1 - y/\gamma^2}$ in our notation. Thus, using also the identity (4.9), we find

<footnote n="6">⁶Or, equivalently,</footnote>
$$ \mathbf{K} = \mathbf{C} \frac{Gm^2}{c^2 r y^{1/2}} \left( \gamma^5 - \frac{\gamma^3}{2} + \frac{3}{8}\gamma - \frac{1}{16\gamma} \right). \quad (4.12) $$

Let $x = \gamma^2$. It is remarkable that the polynomial in $x$ appearing in the rhs of (4.5) can be written in terms of the two polynomials appearing in the rhs of (4.7):

$$ x^4 - \frac{3}{2}x^3 + \frac{3}{8}x^2 + \frac{x}{16} = x \left( x^3 - \frac{x^2}{2} + \frac{3}{8}x - \frac{1}{16} \right) - \left( x^3 - \frac{x}{8} \right). \quad (4.13) $$

Thus, Eq. (4.5) can be rewritten in the form

$$ \frac{dx}{dr} = \frac{8Gm}{r^2 y^{3/2} c^2} \left[ x \left( x^3 - \frac{x^2}{2} + \frac{3}{8}x - \frac{1}{16} \right) - \left( x^3 - \frac{x}{8} \right) \right]. \quad (4.14) $$

Subtracting Eqs. (4.7) and (4.14) we find

$$ \frac{d}{dr}(x-y) + \frac{2}{r}(x-y) = \frac{8Gm}{r^2 y^{3/2} c^2}(x-y) \left( x^3 - \frac{x^2}{2} + \frac{3}{8}x - \frac{1}{16} \right); \quad (4.15) $$

hence

$$ \frac{\frac{d}{dr}[r^2(x-y)]}{r^2(x-y)} = \frac{8Gm}{r^2 y^{3/2} c^2} \left( x^3 - \frac{x^2}{2} + \frac{3}{8}x - \frac{1}{16} \right). \quad (4.16) $$

Dividing the above equation and (4.5) we find

$$ \frac{d[r^2(x-y)]}{r^2(x-y)} = \frac{x^3 - \frac{x^2}{2} + \frac{3}{8}x - \frac{1}{16}}{x^4 - \frac{3x^3}{2} + \frac{3}{8}x^2 + \frac{x}{16}} dx. \quad (4.17) $$

Using the fact that the denominator of the rhs of (4.17) factorizes,<fnref n="7" /> we can finally integrate (4.17) in the nice form

$$ r^2(x-y) = k \frac{(x^2 - x - \frac{1}{8})^{4/3}}{x(x - \frac{1}{2})^{2/3}}, \quad (4.18) $$

where $k$ is a constant. In summary, the basic equation for $\gamma$ is Eq. (4.5), where $y$ is given by

$$ y = \gamma^2 - k \frac{(\gamma^4 - \gamma^2 - \frac{1}{8})^{4/3}}{r^2 \gamma^2 (\gamma^2 - \frac{1}{2})^{2/3}}. \quad (4.19) $$

Of course, Eq. (4.19) is physically meaningful only at the 1PM order, i.e., at leading order in $G$, and the constant $k$ is

$$ x^4 - \frac{3x^3}{2} + \frac{3}{8}x^2 + \frac{x}{16} = x \left( x - \frac{1}{2} \right) \left( x^2 - x - \frac{1}{8} \right). $$

not independent from the two fundamental integrals of motion $E$ and $J$. An easy calculation shows that

$$ k = \frac{4J^2 c^2}{E^2} \frac{\bar{\gamma}^4 (\bar{\gamma}^2 - \frac{1}{2})^{2/3}}{(\bar{\gamma}^4 - \bar{\gamma}^2 - \frac{1}{8})^{4/3}}, \quad (4.20) $$

where $\bar{\gamma} = \frac{E}{mc^2}$ denotes the constant Lorentz factor of special relativity. The integral of motion (4.19) expanded at order $G$ is equivalent to saying that $\gamma = \bar{\gamma} + \delta\gamma$, where

$$ \delta\gamma = - \frac{Gm}{r\bar{\gamma}^{1/2}} (2\bar{\gamma}^2 - 1) \left( 2\bar{\gamma}^3 - 2\bar{\gamma} - \frac{1}{4\bar{\gamma}} \right), \quad (4.21) $$

in agreement with the earlier result (4.10). Note that, here, $\bar{\gamma}^{1/2} = \bar{\gamma} \sqrt{1 - \frac{4J^2c^2}{r^2E^2}}$ with this approximation.

## V. LAGRANGIAN FORMALISM IN HARMONIC COORDINATES

### A. The general case

In this section we look for a Lagrangian associated with our general equations of motion, for any $N$, in harmonic coordinates; see Eqs. (2.12)–(2.13). The Lagrangian will be given by the special-relativistic expression plus terms of order $G$, and we neglect higher-order terms in $G$. It is known that the Lagrangian in harmonic coordinates is a generalized one, depending not only on positions and velocities $\mathbf{y}_a$, $\mathbf{v}_a$ but also on the accelerations $\mathbf{a}_a = d\mathbf{v}_a/dt$ of the particles. In a PN expansion, the accelerations appear at order 2PN [34,35] and they are contained in terms linear in $G$; see Eq. (209) in [14]. Of course, replacing the accelerations by the equations of motion would yield negligible terms of order $G^2$. However, it is not allowed to replace accelerations in a Lagrangian while remaining in the same coordinate system. Such replacement is equivalent to a shift in the particles’ trajectories (or “contact” transformation); i.e., the new Lagrangian is physically equivalent to the original one but written in a different coordinate system [36]. We shall confirm this result within our 1PM framework. Furthermore, by employing the technique of double-zero (or multiple-zero) terms [35], it is sufficient to consider a Lagrangian that is linear in accelerations. Indeed, the procedure can work for any PN order and is thus formally valid at the 1PM order. Therefore, we look for the Lagrangian in the form

$$ L[y, v, a] = -\sum_a \frac{m_a c^2}{\gamma_a} + \lambda + \sum_a q_a^i a_a^i, \quad (5.1) $$

where we symbolize the functional dependence of $L$ by $L[y, v, a] \equiv L[\{\mathbf{y}_a, \mathbf{v}_a, \mathbf{a}_a\}]$, where $\lambda$ and $q_a^i$ are of order $G$ and both depend only on positions and velocities, i.e., $\lambda[y, v]$ and $q_a^i[y, v]$. The linear dependence in accelerations is made explicit in (5.1). Denoting by $p_a^i$ and $q_a^i$ the

<footnote n="7">⁷Namely,</footnote>
conjugate momenta associated with the positions $y_a^i$ and velocities $v_a^i$, i.e.,

$$p_a^i = \frac{\delta L}{\delta v_a^i} = \frac{\partial L}{\partial v_a^i} - \frac{d}{dt} \left( \frac{\partial L}{\partial a_a^i} \right), \quad (5.2a)$$

$$q_a^i = \frac{\delta L}{\delta a_a^i} = \frac{\partial L}{\partial a_a^i}, \quad (5.2b)$$

the equations of motion take the form

$$\frac{\mathrm{d}p_a^i}{\mathrm{d}t} = \frac{\partial L}{\partial y_a^i}. \quad (5.3)$$

The conserved energy is given by the generalized Legendre transformation (see, e.g., [37,38])

$$E = \sum_a (p_a^i v_a^i + q_a^i a_a^i) - L. \quad (5.4)$$

In both Eqs. (5.3) and (5.4) the accelerations have to be replaced by the equations of motion. For instance, the term $q_a^i a_a^i$ in $E$ is second order in $G$ and can be neglected at 1PM order (but note that it cancels anyway with the same term in $L$). Using the form of the Lagrangian (5.1) we obtain the equations of motion as

$$f_a^i = \frac{\delta \lambda}{\delta y_a^i} + \ddot{q}_a^i. \quad (5.5)$$

Here we define the force $f_a^i = m_a \frac{d}{dt}(\gamma_a v_a^i)$ and $\frac{\delta \lambda}{\delta y_a^i} = \frac{\partial \lambda}{\partial y_a^i} - \frac{d}{dt} \left( \frac{\partial \lambda}{\partial v_a^i} \right)$, and the dots refer to time derivatives.<fnref n="8" />

On the other hand, the integral of the energy is given by (3.1) with

$$V = \sum_a v_a^i \frac{\partial \lambda}{\partial v_a^i} - \lambda - \sum_a v_a^i \dot{q}_a^i. \quad (5.6)$$

The left-hand sides (LHS) of Eqs. (5.5) and (5.6) are known, since we have determined $f_a^i$ and $V$ in (2.12)–(2.13) and (3.11). However, the two equations are not independent, since $f_a^i$ and $V$ satisfy the constraint

$$\frac{\mathrm{d}V}{\mathrm{d}t} + \sum_a v_a^i f_a^i = 0. \quad (5.7)$$

Our strategy is to first find a particular solution, say $(\hat{\lambda}, \hat{q}_a^i)$, such that the set of conjugate momenta $\hat{q}_a^i$ obey the equation

$$\sum_a v_a^i \hat{q}_a^i = 0. \quad (5.8)$$

To order $G$, the same constraint is also satisfied by the time derivative $\dot{q}_a^i$; therefore, for the solution $(\hat{\lambda}, \hat{q}_a^i)$, Eq. (5.6) reduces to an ordinary Legendre transformation,

$$V = \sum_a v_a^i \frac{\partial \hat{\lambda}}{\partial v_a^i} - \hat{\lambda}. \quad (5.9)$$

In order to determine $\hat{\lambda}$, we note that the potential $V$ given in (3.11) reduces in the limit $c \to +\infty$ to the Newtonian approximation, namely, $V = U + \mathcal{O}(1/c^2)$, where the Newtonian potential is

$$U = -\sum_{a<b} \frac{Gm_a m_b}{r_{ab}}. \quad (5.10)$$

If we subtract the Newtonian limit $U$, we get a quantity which tends to zero when $c \to +\infty$ like $\mathcal{O}(1/c^2)$. Then, it is straightforward to show that a well-behaved solution of Eq. (5.9) is

$$\hat{\lambda} = -U + \frac{1}{c} \int_c^{+\infty} ds \left[ V \left( y_b, \frac{v_b}{s} \right) - U(y_b) \right]. \quad (5.11)$$

Indeed, we insert into Eq. (3.11) all the relevant factors $c$ and make the replacement $c \to s$, and then integrate over the “speed of light” $s$ from the physical value $c$ up to infinity. The bound $s \to \infty$ corresponds to the Newtonian limit and we see from the definition of the Newtonian potential (5.10) that the integral is convergent. The first term in Eq. (5.11) represents the Newtonian approximation with the correct minus sign for a Lagrangian, and the integral represents formally the complete series of PN corrections, but resummed in the PM approximation. The result (5.11) can be rewritten in a simpler way as the “partie finie” (PF) of the integral in the sense of Hadamard [39], taking care of the divergence at infinity:

$$\hat{\lambda} = \mathrm{PF} \frac{1}{c} \int_c^{+\infty} ds V \left( y_b, \frac{v_b}{s} \right). \quad (5.12)$$

Note that for this very simple type of divergence $\sim s^0 + \mathcal{O}(s^{-2})$, the PF does not depend on an arbitrary constant. Unfortunately, given the complicated structure of $V$ in Eq. (3.11), we have not been able to perform explicitly the integration and obtain a closed-form expression in the general case. But we shall discuss in Sec. IV the equal-mass case for which an analytical form exists.

Equations (5.11) and (5.12) give a solution of Eq. (5.6), but we still have to adjust $\hat{q}_a^i$ in order to satisfy the equations of motion (5.5). Thus, we look for $\hat{q}_a^i$ satisfying

$$\ddot{\hat{q}}_a^i = f_a^i - \frac{\delta \hat{\lambda}}{\delta y_a^i}, \quad (5.13)$$

where the rhs is now known. Since it is of order $G$, we can integrate following the same method as used in Sec. III. That is, we transform (5.13) into an ordinary differential

<footnote n="8">⁸We have $f^i = 4f_1^i = -4f_2^i$ in the notation of Eq. (4.3).</footnote>
equation using the change of variable $t \to r_{ab}$ and then integrate twice to determine $\dot{q}_a^i$, followed by $\dot{q}_a^i$. The solution automatically satisfies the constraint (5.8) by virtue of (5.7).

Finally we have found a particular Lagrangian associated with our equations of motion, in the form of the particular solution $(\hat{\lambda}, \dot{q}_a^i)$. Now the general solution $(\lambda, q_a^i)$ can be obtained by adding an arbitrary total time-derivative $dF/dt$, where $F$ is a function of the positions $y_a^i$ and velocities $v_a^i$; indeed, by definition of the functional derivative, $\delta(dF/dt)/\delta y_a^i = 0$. Hence the general solution (for the class of harmonic-coordinates Lagrangians that are linear in accelerations) reads

$$ \lambda = \hat{\lambda} + \sum_a v_a^i \frac{\partial F}{\partial y_a^i}, \qquad (5.14a) $$

$$ q_a^i = \hat{q}_a^i + \frac{\partial F}{\partial v_a^i}. \qquad (5.14b) $$

At the 1PM order we find that the Lagrangian in harmonic coordinates irreducibly depends on accelerations; i.e., it is impossible to determine $F$ such that $q_a^i = 0$. However, we know that the accelerations in a Lagrangian can always be eliminated by appropriate shifts of the trajectories. In the present case, the shifts are given by $Y_a^i = y_a^i + \xi_a^i$, where $\xi_a^i$ is of order $G$ and is determined by the conjugate momenta $q_a^i$. We find (adapting techniques from Refs. [35–37,40])

$$ \xi_a^i = -\frac{1}{m_a \gamma_a} \left[ q_a^i - \frac{v_a^i (v_a q_a)}{c^2} \right]. \qquad (5.15) $$

Then, an ordinary Lagrangian, dynamically equivalent to the harmonic-coordinates one (5.1), but valid for the shifted trajectories $Y_a^i$ instead, is given to order $G$ by

$$ L'[Y, V] = -\sum_a \frac{m_a c^2}{\gamma_a} + \lambda, \qquad (5.16) $$

where $\lambda$ and $\gamma_a$ denote the same functionals as before, but expressed in terms of the shifted variables, i.e., $\lambda[Y, V]$ and $\gamma_a[V]$. The conserved energy is given by $E = \sum_a m_a c^2 \gamma_a + V'$ but with the modified potential function

$$ V' = V + \sum_a V_a^i q_a^i. \qquad (5.17) $$

Hence, in this construction, the particular solution $\hat{\lambda}$ found in Eqs. (5.11)–(5.12) represents the equivalent ordinary Lagrangian in shifted variables for which the energy functional remains the same ($V' = V$); compare (5.8) and (5.17).

In spite of the fact that it has not been possible to find a closed-form expression for the integral (5.11) in the general case, the relevant integral can easily be worked out in the PN limit $c \to \infty$. We start from the known 4PN expansion of the potential $V$ and explicitly perform the integration

(5.11) term by term, to obtain the corresponding 4PN expansion of $\hat{\lambda}$. Then, we derive the 4PN result for $\dot{q}_a^i$ by the method indicated previously. Finally, we find a unique total time derivative, with some function $F_{PN}$ given in the form of a PN expansion, so that the Lagrangian agrees up to order $G$ with the 4PN Lagrangian published in Refs. [23,25]. In conclusion, the 4PN limit of our Lagrangian is correct. The function $F_{PN}$ is given by Eqs. (C1)–(C2) in Appendix C. Furthermore we have extended the analysis to higher order and give the terms of order $G$ in the harmonic-coordinates Lagrangian up to 5PN order in Appendix C; see (C5)–(C6) there.

## B. The two-body equal-mass case

We apply the previous procedure for two bodies ($N = 2$) with equal masses, following the notation (4.1). The harmonic-coordinates Lagrangian becomes in this case

$$ L = -\frac{mc^2}{\gamma} + \lambda + q^i a^i, \qquad (5.18) $$

where $q^i = \frac{1}{2}(q_1^i - q_2^i)$ is the conjugate momentum of the relative velocity $v^i$, as defined by $q^i = \partial L/\partial a^i$. We follow the method of Sec. VA to determine $\lambda$ and $q^i$. We first look for a particular solution $(\hat{\lambda}, \dot{q}^i)$ such that $v^i \dot{q}^i = 0$. Then, we determine $\hat{\lambda}$ from the potential function $V$ determined for the case at hand in Sec. IV as

$$ V = \frac{Gm^2}{4ry^{1/2}} (2\gamma^2 - 1) \left( 8\gamma^3 - 8\gamma - \frac{1}{\gamma} \right). \qquad (5.19) $$

Substituting (5.19) into Eq. (5.12), we can integrate and obtain a rather complicated result:

$$ \begin{aligned} \hat{\lambda} = & \frac{Gm^2}{4r\gamma} \left(1 + \frac{8(\gamma^2 - 1)^2}{y-1}\right) \sqrt{y} + \frac{Gm^2}{4r} \frac{\gamma^2 - 1}{\gamma\sqrt{\gamma^2 - y}} \\ & \times \left(\arctan \sqrt{\frac{y}{\gamma^2 - y}} - \frac{\pi}{2}\right) \\ & - \frac{2Gm^2}{r} \frac{(\gamma^2 - 1)(\gamma^2 - y)}{\gamma(y-1)^{3/2}} \ln(\sqrt{y} + \sqrt{y-1}). \end{aligned} \qquad (5.20) $$

Note that the circular orbit limit of this expression, $(nv) \to 0$ or $y \to 1$, is well defined, as can be verified directly in Eq. (5.20).

As discussed in Sec. VA, the Lagrangian $\hat{L} = -\frac{mc^2}{\gamma} + \hat{\lambda}$ represents an equivalent Lagrangian which is ordinary (no accelerations) and valid in a coordinate system shifted with respect to the harmonic-coordinates system, with shift $\xi^i = -\frac{4}{m\gamma}\dot{q}^i$; see Eq. (5.15). We shall investigate in Sec. VI A the connection of this Lagrangian to the 1PM Hamiltonian of [30].

Next, we have to determine $\dot{q}^i$ using the method of Sec. VA. The computation is straightforward and we obtain an expression which contains transcendental functions like
in (5.20). We do not give the complete expression because we can now adjust the arbitrary function F in Eqs. (5.14) in order to simplify the result as much as possible. In fact we find that all the transcendental functions (arctan x and ln x) in both $\hat{\lambda}$ and $\hat{q}^i$ can be removed by the following choice:

$$F = -\frac{Gm^2}{8} \sqrt{\frac{y-1}{\gamma^2-y}} \left( \arctan \sqrt{\frac{y}{\gamma^2-y}} - \frac{\pi}{2} \right) - \frac{Gm^2 \gamma^2 - y}{2(y-1)} \left[ \ln \left( \sqrt{y} + \sqrt{y-1} \right) - \sqrt{y-1} \right]. \quad (5.21)$$

In the second term, we have subtracted from the logarithm the term $\sqrt{y-1}$ in order to ensure the well-defined circular orbit limit $y \to 1$. Finally, adding up the total time derivative $dF/dt$ [see Eqs. (5.14)], we obtain the expressions

$$\lambda = \frac{Gm^2 4y^2 - 8y + 1 + \sqrt{y}(8\gamma^4 - 8\gamma^2 + 1) + 4\gamma^4}{4r \gamma(y + \sqrt{y})}, \quad (5.22a)$$

$$q^i = \frac{Gm^2}{16} \left[ \frac{4\gamma^4 - 4y\gamma^2 + 1 + \sqrt{y}(8\gamma^4 - 8\gamma^2 + 1)}{\gamma(y + \sqrt{y})} n^i + \frac{4\gamma^4 - 4y\gamma^2 + 1 - \sqrt{y}(8\gamma^2 - 1)}{2(y + \sqrt{y})} \sqrt{y-1} v^i \right], \quad (5.22b)$$

which completely specify our acceleration-dependent Lagrangian (5.1). By varying this Lagrangian we can recover the equations of motion (4.2)–(4.3).

# VI. HAMILTONIAN FORMALISM

## A. Equal-mass case

Let us rewrite in full form the ordinary Lagrangian we obtained in the previous section, $\hat{L} = -\frac{mc^2}{\gamma} + \hat{\lambda}$, and which we recall is not valid in harmonic coordinates but in a shifted coordinate system $X^i = x^i + \xi^i$; see Eq. (5.15). Using the more common notations $V^2$ and $(NV) = \dot{R}$ (with $c=1$ and removing the “hat”), we find

$$L[X, V] = -m\sqrt{1 - \frac{V^2}{4}} + \frac{Gm^2}{8R} \left( 1 + \frac{8V^4}{(NV)^2(4 - V^2)} \right) \sqrt{4 + (NV)^2 - V^2} \\ + \frac{Gm^2}{8R} \frac{V^2}{\sqrt{V^2 - (NV)^2}} \left( \arctan \sqrt{\frac{4 + (NV)^2 - V^2}{V^2 - (NV)^2}} - \frac{\pi}{2} \right) \\ - \frac{Gm^2}{R} \frac{V^2(V^2 - (NV)^2)}{(NV)^3} \ln \left( \frac{(NV) + \sqrt{4 + (NV)^2 - V^2}}{\sqrt{4 - V^2}} \right). \quad (6.1)$$

The Hamiltonian is easily constructed by an ordinary Legendre transformation, and we obtain, up to order G (posing $\bar{m} = \sqrt{m^2 + 4P^2}$)

$$H[X, P] = \bar{m} - \frac{G}{4R\bar{m}} \left( m^2 + \frac{32P^4}{(NP)^2} \right) \sqrt{m^2 + 4(NP)^2} \\ - \frac{Gm^2}{2R\bar{m}} \frac{P^2}{\sqrt{P^2 - (NP)^2}} \left( \arctan \sqrt{\frac{m^2 + 4(NP)^2}{4(P^2 - (NP)^2)}} - \frac{\pi}{2} \right) \\ + \frac{4Gm^2}{R\bar{m}} \frac{P^2(P^2 - (NP)^2)}{(NP)^3} \ln \left( \frac{2(NP) + \sqrt{m^2 + 4(NP)^2}}{m} \right). \quad (6.2)$$

It is very interesting to compare this Hamiltonian with the one obtained by Ledvinka, Schäfer, and Bičák (LSB) [30]. These authors obtained a nice, closed-form expression for the Hamiltonian in the general case; see Eq. (11) in Ref. [30]. On the other hand, the specific method explained in Sec. V did not lead to a closed-form expression for the Lagrangian except in the equal-mass case, yielding the results (6.1)–(6.2).

We want to find a canonical transformation between our Hamiltonian (6.2) and the LSB one [30]. In particular, this canonical transformation should be able to remove the transcendental functions (arctan x and ln x) from the Hamiltonian (6.2), since there are no such functions in the LSB Hamiltonian. This is easy to do, as we have already seen in Eq. (5.21), where the transcendental functions at the
level of the Lagrangian can be removed by adding an appropriate total time derivative to the Lagrangian.

To find the canonical transformation we adopt reduced canonical variables $(R, \Phi, P_R, P_\Phi)$ by posing $P^2 = P_R^2 + P_\Phi^2/R^2$ and $P_R = (NP)$. Then, we look for a canonical transformation

$$ (R, \Phi, P_R, P_\Phi) \rightarrow (R', \Phi', P'_R, P'_\Phi), \quad (6.3) $$

associated with some generating function $M^{\text{gen}}(R, \Phi, P'_R, P'_\Phi)$. We assume no explicit time dependence of the generating function so that

$$ H(R, \Phi, P_R, P_\Phi) = H'(R', \Phi', P'_R, P'_\Phi). \quad (6.4) $$

Furthermore, since the special-relativistic limits of both Hamiltonians agree, the generating function should be of the type

$$ M^{\text{gen}} = RP'_{R} + \Phi P'_{\Phi} + F^{\text{gen}}(R, P'_{R}, P'_{\Phi}), \quad (6.5) $$

where the function $F^{\text{gen}}(R, P'_R, P'_\Phi)$ is of order $G$. Finally, both Hamiltonians are rotationally invariant, so the function $F^{\text{gen}}$ should not depend on $\Phi$. To order $G$ the transformation laws between the two sets of variables are as follows:

$$ R' = R + \frac{\partial F^{\text{gen}}}{\partial P_R}(R, P_R, P_\Phi), \quad (6.6a) $$

$$ \Phi' = \Phi + \frac{\partial F^{\text{gen}}}{\partial P_{\Phi}}(R, P_{R}, P_{\Phi}), \quad (6.6b) $$

$$ P'_{R} = P_{R} - \frac{\partial F^{\text{gen}}}{\partial R}(R, P_{R}, P_{\Phi}), \quad (6.6c) $$

and, in addition, $P'_\Phi = P_\Phi$, which implies that the conserved angular momentum of both Hamiltonians is the same to order $G$. Inserting Eqs. (6.6) into (6.4) we obtain (to order $G$)

$$ H'(R, \Phi, P_R, P_\Phi) = H(R, \Phi, P_R, P_\Phi) - \frac{4}{\bar{m}} \left( P_R \frac{\partial F^{\text{gen}}}{\partial R} + \frac{P_\Phi^2}{R^3} \frac{\partial F^{\text{gen}}}{\partial P_R} \right). \quad (6.7) $$

We proceed in two steps. We first apply a generating function in order to remove the transcendental functions from Eq. (6.2). As the effect of a canonical transformation with generating function $F^{\text{gen}}$ is the same as the effect of adding a total time derivative in the Lagrangian formalism, we may use for the generating function the one already computed in Eq. (5.21); hence $F^{\text{gen}} = F$, with

$$ F = -\frac{Gm^2 RP_R}{8P_\Phi} \left( \arctan \left( \frac{R\sqrt{m^2 + 4P_R^2}}{2P_\Phi} \right) - \frac{\pi}{2} \right) - \frac{Gm^2 P_\Phi^2}{2R^2 P_R^2} \left[ \ln \left( \frac{2P_R + \sqrt{m^2 + 4P_R^2}}{m} \right) - \frac{2P_R}{m} \right]. \quad (6.8) $$

Under this canonical transformation our Hamiltonian (6.2) is very much simplified. Posing $y = 1 + 4(NP)^2/m^2$ (consistently with our notation in previous sections) we get

$$ H' = \bar{m} + \frac{4Gm}{R\bar{m}} \left( \frac{P^4 - (NP)^4}{(NP)^2} + \frac{1}{\sqrt{y}} \left[ -\frac{m^2}{16} + (NP)^2 - 2P^2 - 8\frac{P^4}{m^2} - \frac{P^4}{(NP)^2} \right] \right). \quad (6.9) $$

On the other hand, the LSB Hamiltonian, Eq. (11) in [30] for two equal masses, is

$$ H_{\text{LSB}} = \bar{m} + \frac{Gm}{R} \left( -7P^2 - (NP)^2 + \frac{m^3}{4\bar{m}\sqrt{y}} \right), \quad (6.10) $$

so it must be a simple task to find the canonical transformation linking the two Hamiltonians. Indeed, we find that there exists such a canonical transformation, but which is not so simple as it does not admit a closed-form expression, involving the two-variable hypergeometric Appell function $F_1(\alpha; \beta, \beta'; \gamma; x, y)$ [41]. We finally get

$$ \begin{aligned} \delta F_{\text{LSB}} = & -\frac{G\bar{m}P_R}{4} - \frac{GmP_\Phi^2}{R^2P_R} \\ & + \frac{Gm}{\sqrt{y}} \left( P_R + \frac{4P_R^3}{m^2} + \frac{P_\Phi^2}{R^2P_R} + \frac{4P_R P_\Phi^2}{m^2R^2} \right) \\ & - GP \left( 2\bar{m} + \frac{m^2}{\bar{m}} \right) \ln(R[\bar{m}P_R + mP\sqrt{y}]) \\ & + 2G\bar{m}P \ln(RP[P + P_R]) \\ & + \frac{GmPP_R^3R^3}{3P_\Phi^3} F_1 \left( \frac{3}{2}; \frac{1}{2}, 1; \frac{5}{2}; -\frac{\bar{m}^2P_R^2R^2}{m^2P_\Phi^2}, -\frac{P_R^2R^2}{P_\Phi^2} \right). \end{aligned} \quad (6.11) $$

Hence, the complete canonical transformation linking our Hamiltonian to the LSB one in the equal-mass case is generated by $F + \delta F_{\text{LSB}}$. However, as we said, we have not been able to connect the two Hamiltonians in the more general case.

## B. Frame of the center of mass

In special relativity the linear momentum of noninteracting particles is $m_a\gamma_a v_a^t$, and we have obtained in Sec. III the total linear momentum at the 1PM order (in the case of two particles with arbitrary masses):
$$P^i = m_1\gamma_1 v_1^i + m_2\gamma_2 v_2^i + \Pi^i, \qquad (6.12)$$

where $\Pi^i$ is the correction of order G given by Eq. (3.13).

The frame of the center of mass (CM) is defined by $\mathbf{P} = 0$. In this frame the dynamics will be described by dynamical variables $(\mathbf{x}, \mathbf{v})$, where $\mathbf{x} = \mathbf{y}_1 - \mathbf{y}_2$ and $\mathbf{v} = \mathbf{v}_1 - \mathbf{v}_2$ are the relative variables. The CM Lagrangian in harmonic coordinates will be a functional of $\mathbf{x}, \mathbf{v}$ and also the acceleration $\mathbf{a}$. In the Hamiltonian formalism the CM dynamics is described by canonical variables $(\mathbf{x}, \mathbf{p})$, where $\mathbf{p}$ is the CM linear momentum $\mathbf{p} = \mathbf{p}_1 = -\mathbf{p}_2$. In special relativity we have $\mathbf{v}_1 = \mathbf{p}/\bar{m}_1$ and $\mathbf{v}_2 = -\mathbf{p}/\bar{m}_2$, where we recall that $\bar{m}_a = \sqrt{m_a^2 + p^2}$ (and pose $c = 1$). At the 1PM order, let us pose with full generality, for the individual velocities,

$$v_1^i = \frac{p^i + Y^i}{\bar{m}_1} + W^i, \qquad (6.13a)$$

$$v_2^i = -\frac{p^i + Y^i}{\bar{m}_2} + W^i, \qquad (6.13b)$$

where $Y^i$ and $W^i$ represent the 1PM corrections of order G, with $Y^i$ (like $p^i$) changing sign when we exchange the two particles' labels, and with $W^i$ staying invariant under such exchange $1 \leftrightarrow 2$. The relative velocity is given by

$$v^i = \frac{p^i + Y^i}{\bar{\mu}}, \qquad (6.14)$$

where the associated reduced mass is $\bar{\mu} = \bar{m}_1\bar{m}_2/\bar{m}$, where $\bar{m} = \bar{m}_1 + \bar{m}_2$ is the total mass. From Eqs. (6.13) one readily derives (to order G)

$$m_1\gamma_1 v_1^i = p^i + X^i + \bar{m}_1 \left( W^i + \frac{p^i}{m_1^2} (pW) \right), \qquad (6.15a)$$

$$m_2\gamma_2 v_2^i = -p^i - X^i + \bar{m}_2 \left( W^i + \frac{p^i}{m_2^2} (pW) \right), \qquad (6.15b)$$

where $(pW)$ denotes the scalar product as usual, and where we introduce the intermediate notation $X^i$ defined by the (easily invertible) relation<fnref n="9" />

$$Y^i = X^i - (pX) \frac{\partial \ln \bar{\mu}}{\partial p^i}. \qquad (6.16)$$

The expressions (6.15) are ready for insertion into the total linear momentum (6.12) and the CM prescription $P^i = 0$.

$$\frac{\partial \ln \bar{\mu}}{\partial p^i} = \frac{1 - 3\bar{\nu} p^i}{\bar{\nu}^2 \bar{m}}.$$

Although $X^i$ is clearly left undetermined at this stage, the CM prescription specifies $W^i$ in terms of the known 1PM correction in the linear momentum, given by (3.13), and we find

$$W^i = - \frac{1}{\bar{m}_1 + \bar{m}_2} \left[ \Pi^i - \frac{\bar{m}_1 m_2^2 + \bar{m}_2 m_1^2}{\bar{m}_1^3 m_2^2 + \bar{m}_2^3 m_1^2} (p\Pi) p^i \right]. \qquad (6.17)$$

The CM Hamiltonian is then obtained starting from the expression of the energy $E = m_1\gamma_1 + m_2\gamma_2 + V$, where $V$ is given by (3.11), by replacing the velocities $v_1^i, v_2^i$ in terms of the CM canonical momentum $p^i$ using (6.13). In this way, we obtain the CM Hamiltonian to order G in the form

$$H = \bar{m}_1 + \bar{m}_2 + \left( \frac{1}{\bar{m}_1} + \frac{1}{\bar{m}_2} \right) (pX) + K, \qquad (6.18)$$

where $(pX)$ is still undetermined, but $K$ is known at this stage:

$$K = - \frac{\bar{m}_1^2 m_2^2 - \bar{m}_2^2 m_1^2}{\bar{m}_1^3 m_2^2 + \bar{m}_2^3 m_1^2} (p\Pi) + V, \qquad (6.19)$$

or, after an explicit calculation,

$$K = - \frac{\bar{m}_1 \bar{m}_2^2 + p^2}{\bar{m}_1^3 m_2^2 + \bar{m}_2^3 m_1^2} \times (\bar{m}_1^2 \bar{m}_2^2 - (3\bar{m}_1^2 + 4\bar{m}_1 \bar{m}_2 + 3\bar{m}_2^2) p^2 + p^4) \times \left[ \frac{G \bar{m}_1^2}{r\sqrt{m_1^2 + (np)^2}} + \frac{G \bar{m}_2^2}{r\sqrt{m_2^2 + (np)^2}} \right]. \qquad (6.20)$$

Now with the previous form (6.18) of the Hamiltonian, we can determine $X^i$ (or equivalently $Y^i$), recalling the definition (6.14) of the velocity, by using the Hamiltonian equation $v^i = \partial H/\partial p^i$, which is easily seen to be equivalent to the equation

$$p^i \frac{\partial X^j}{\partial p^i} + \bar{\mu} \frac{\partial K}{\partial p^i} = 0. \qquad (6.21)$$

The latter equation, after contraction with $p^i$, implies in turn the partial differential equation

$$p^i \frac{\partial (pX)}{\partial p^i} - (pX) = -\bar{\mu} p^i \frac{\partial K}{\partial p^i}, \qquad (6.22)$$

which can in principle be solved to obtain $(pX)$ and the Hamiltonian (6.18). However, the source term in the rhs involves the partial derivatives of the expression (6.20), which are rather cumbersome to compute in practice, so we prefer solving the equivalent equation obeyed by the auxiliary quantity $\sigma = (pX) + \bar{\mu}K$, which is

<footnote n="9">⁹Denoting $\bar{\nu} = \bar{\mu}/\bar{m}$ we have</footnote>
$$p^i \frac{\partial \sigma}{\partial p^i} - \sigma = -\bar{\mu} K \left( 1 - p^i \frac{\partial \ln \bar{\mu}}{\partial p^i} \right), \quad (6.23)$$

and whose source term $S = -\bar{\mu}K(1 - p^i \frac{\partial \ln \bar{\mu}}{\partial p^i})$ is simpler. Not surprisingly, the latter equation has the same form as the one found in the Lagrangian formalism, Eq. (5.9), and we solve it by means of the same technique, i.e., (with PF denoting the Hadamard *partie finie*),

$$\sigma = \text{PF} \frac{1}{c} \int_{c}^{+\infty} ds S \left( \mathbf{x}, \frac{\mathbf{p}}{s} \right). \quad (6.24)$$

Like for the Lagrangian formalism we could not obtain a closed-form expression for the integral (6.24) in the general case, but only in the equal-mass case $m_1 = m_2$, in which case we simply recover the Hamiltonian (6.2).

## ACKNOWLEDGMENTS

One of us (L. B.) would like to thank Gilles Esposito-Farèse for interesting discussions and for comments on a preliminary version of the manuscript.

## APPENDIX A: GRAVITATIONAL SCATTERING OF TWO PARTICLES

We compute the scattering of two particles with arbitrary masses using the formalism developed in Secs. II and III. The total change of linear momentum of the particle 1 between the “in” and “out” states at $t = \pm \infty$ is given by

$$\Delta p_1^i = m_1 c \Delta u_1^i = m_1 c \int_{-\infty}^{+\infty} d\tau_1 \frac{du_1^i}{d\tau_1}. \quad (A1)$$

The equations of motion of the particle 1 in the field of particle 2 are given by Eq. (2.12); therefore the total change in the spatial components of the linear momentum is given by the equation

$$\Delta p_1^i = -Gm_1 m_2 \gamma_1^{-1} \int_{-\infty}^{+\infty} \frac{dt}{r_{12}^2 y_{12}^{3/2}} \left[ (2\epsilon_{12}^2 - 1)n_{12}^i \right. \\ \left. \qquad + (2\epsilon_{12}^2 + 1) \left( -(n_{12}u_1) + \epsilon_{12}(n_{12}u_2) \right) u_1^i \right. \\ \left. \qquad + \left( 4\epsilon_{12}(n_{12}u_1) - (2\epsilon_{12}^2 + 1)(n_{12}u_2) \right) u_2^i \right], \quad (A2)$$

whereas the total change in the zeroth component of the momentum is $\Delta p_1^0 = \frac{v_1^i}{c} \Delta p_1^i$.

Recall that to order $G$ we can integrate the motion on straight lines; therefore, all velocity-dependent factors in the integrand of (A2) are in fact constant, except those involving the separation $r_{12}$ and unit direction $n_{12}^i$. Therefore, we see immediately that the integral (A2) can be expressed by means of the solution $D_{12}^i$ of the

elementary equation (3.4). Denoting the total change in this quantity by $\Delta D_{12}^i = D_{12}^i(+\infty) - D_{12}^i(-\infty)$, we obtain

$$\Delta p_1^i = -Gm_1 m_2 \gamma_1^{-1} [(2\epsilon_{12}^2 - 1)\Delta D_{12}^i \\ + (2\epsilon_{12}^2 + 1)(-(\Delta D_{12}u_1) + \epsilon_{12}(\Delta D_{12}u_2))u_1^i \\ + (4\epsilon_{12}(\Delta D_{12}u_1) - (2\epsilon_{12}^2 + 1)(\Delta D_{12}u_2))u_2^i], \quad (A3)$$

where we denote, for instance, $(u_1\Delta D_{12}) = \Delta D_{12}^i u_1^i = \gamma_1 \Delta D_{12} \cdot \mathbf{v}_1/c$.

The explicit expression of the solution $D_{12}^i$ has been obtained in Eq. (3.9); again, since the velocities are constant, the only time-dependent variable in this expression is the separation $r_{12}$. At $t = +\infty$ in the future we have $r_{12} \to +\infty$ and we readily find the limit from Eq. (3.9). When $t = -\infty$ we must choose the part of the trajectory for which $r_{12}$ decreases rather than increases, and the limit just changes sign, adding finally a factor 2. Thus, the total change in this quantity is given by

$$\Delta D_{12}^i = \frac{2 \mu_{12}^i |v_{12}|}{C_{12}^2 + \frac{\gamma_2^2}{c^2} \left[ (\mu_{12} v_2)^2 v_{12}^2 + C_{12}^2 \frac{(v_{12} v_2)^2}{v_{12}^4} \right]} \sqrt{1 + \frac{\gamma_2^2 (v_{12} v_2)^2}{c^2 v_{12}^4}} \quad (A4)$$

We then return to the original variables, as we did in (3.10). Since the result (A4) is constant, it can be evaluated at any point along the trajectory. We choose to evaluate it at the minimal distance of approach of the two particles, characterized by $\dot{r}_{12} = (n_{12}v_{12}) = 0$. At this point, $r_{12} = b_{12}$ represents the impact parameter of the incoming particles, and $n_{12}^i = b_{12}^i/b_{12}$ is the associated unit direction between the particles at that point. Furthermore, we adopt a center-of-mass frame. Since we are already in a small term of order $G$, we can use the special-relativistic notion of center of mass and denote the CM linear momentum by $p^i = m_1\gamma_1 v_1^i = -m_2\gamma_2 v_2^i$. In this frame, not only do we have $(b_{12}v_{12}) = 0$ at the point of minimal approach, but also $(b_{12}v_1) = (b_{12}v_2) = 0$. The expression (A4) then drastically simplifies, and we obtain

$$\Delta D_{12}^i = \frac{2}{\sqrt{v_{12}^2 + \frac{\gamma_2^2}{c^2} (v_{12} v_2)^2}} \frac{b_{12}^i}{b_{12}^2}. \quad (A5)$$

Note in particular that the scalar products $(u_1\Delta D_{12})$ and $(u_2\Delta D_{12})$ vanish, so only the first term in Eq. (A3) will contribute. On the other hand, we can simplify the result by noting that $v_{12}^2 + \frac{\gamma_2^2}{c^2} (v_{12}v_2)^2 = c^2\gamma_1^{-2}(\epsilon_{12}^2 - 1)$. Finally, we obtain the total change in linear momentum during the scattering encounter as
$$ \Delta \mathbf{p}_1 = - \frac{2Gm_1m_2}{c} \frac{2\epsilon_{12}^2 - 1}{\sqrt{\epsilon_{12}^2 - 1}} \frac{\mathbf{b}_{12}}{b_{12}^2}. \quad (A6) $$

This is a manifestly covariant form [recall that $\epsilon_{12} = -(u_1 u_2)$], which agrees with published results in the literature, like the recent derivation in Fourier space, Eq. (58) of Ref. [31]. In terms of the center-of-mass linear momentum $\mathbf{p} = \bar{m}_1 \mathbf{v}_1 = -\bar{m}_2 \mathbf{v}_2$, where we have posed $\bar{m}_1 = \sqrt{m_1^2 + \mathbf{p}^2/c^2}$ and $\bar{m}_2 = \sqrt{m_2^2 + \mathbf{p}^2/c^2}$, we have

$$ \Delta \mathbf{p}_1 = - \frac{2G}{|\mathbf{p}|} \frac{\bar{m}_1^2 \bar{m}_2^2}{\bar{m}_1 + \bar{m}_2} \left[ 1 + \left( \frac{1}{\bar{m}_1^2} + \frac{1}{\bar{m}_2^2} + \frac{4}{\bar{m}_1 \bar{m}_2} \right) \frac{\mathbf{p}^2}{c^2} + \frac{\mathbf{p}^4}{\bar{m}_1^2 \bar{m}_2^2 c^4} \right] \frac{\mathbf{b}_{12}}{b_{12}^2}. \quad (A7) $$

The latter form agrees with results from Refs. [8,30]. Actually the scattering of two particles has even been computed up to 2PM order in Refs. [5,7,8,12].

## APPENDIX B: CIRCULAR ORBITS AND THEIR STABILITY

Although the most well-known physical application of the PM approximation is for unbound orbits, typically in scattering situations (see Appendix A), it may be mathematically interesting to investigate the case of bound circular orbits. Note that in the case of circular orbits, the terms of order $G^2$ corresponding to 2PM corrections turn out to be also of the order 2PN in the post-Newtonian expansion; these terms are thus neglected with the present 1PM approximation.

We consider two particles with equal masses $m_1 = m_2 = m/2$; see the notation in Sec. IV. The motion takes place in the fixed orbital plane described by polar coordinates $(r, \varphi)$. Posing $C = r^2 \dot{\varphi} = |\mathbf{x} \times \mathbf{v}|$, the equations of motion, equivalent to those presented in Sec. IV, take the form

$$ \ddot{r} - \frac{C^2}{r^3} = - \frac{Gm}{r^2 y^{3/2}} \left[ 8\gamma^2 - 8 + \frac{1}{\gamma^2} - \left(4\gamma^4 - \frac{1}{2}\right)\dot{r}^2 \right], \quad (B1a) $$

$$ \frac{\dot{C}}{C} = \frac{Gm}{r^2 y^{3/2}} \left( 4\gamma^4 - \frac{1}{2} \right) \dot{r}. \quad (B1b) $$

By integration of these equations we recover the conserved energy $E$ and angular momentum $J$ already obtained in Sec. IV, namely,

$$ E = mc^2\gamma + \frac{Gm^2}{4ry^{1/2}}(2\gamma^2 - 1)\left(8\gamma^3 - 8\gamma - \frac{1}{\gamma}\right), \quad (B2a) $$

$$ J = \frac{C}{4} \left[ m\gamma + \frac{Gm^2}{rc^2y^{1/2}} \left( 4\gamma^5 - 2\gamma^3 + \frac{3}{2}\gamma - \frac{1}{4\gamma} \right) \right]. \quad (B2b) $$

The circular orbit is defined by $\dot{r} = \ddot{r} = 0$; hence $r = r_0$ and $C = C_0$ are constant. We denote the orbital frequency of the circular orbit by $\omega_0 = \dot{\varphi}_0 = C_0/r_0^2$, given by

$$ \omega_0^2 = \frac{Gm}{r_0^3} \left( 8\gamma_0^2 - 8 + \frac{1}{\gamma_0^2} \right). \quad (B3) $$

Since we have by definition $v_0^2 = r_0^2 \omega_0^2 = 4c^2(1 - 1/\gamma_0^2)$, Eq. (B3) gives the relation linking the orbital frequency $\omega_0$ to the radius $r_0$ for circular orbits. It can also be written as a relation between $u_0 = \frac{Gm}{r_0 c^2}$ (which represents in fact the small dimensionless PM parameter) and the relativistic factor $\gamma_0$:

$$ u_0 = \frac{4(\gamma_0^2 - 1)}{8\gamma_0^4 - 8\gamma_0^2 + 1}. \quad (B4) $$

We observe that the limit $u_0 \to 0$ (or large separation $r_0 \to \infty$) corresponds to two opposite regimes. One is the usual PN limit which is such that $\gamma_0 \to 1$. However, there is also the ultrarelativistic (UR) regime for which $\gamma_0 \to \infty$. In other words, the large separation regime is composed of two branches,

$$ \gamma_0^2 = \frac{1 + 2u_0 - \epsilon\sqrt{1 - 4u_0 + 2u_0^2}}{4u_0}, \quad (B5) $$

with $\epsilon = +1$ in the PN regime and $\epsilon = -1$ in the UR regime; see the left panel of Fig. 1. Note that in between these two regimes, there is a minimal value for the radius of a circular orbit, which is given by $r_{\min} = (2 + \sqrt{2}) \frac{Gm}{c^2}$. Along the UR branch $u_0^{-1}$ tends to infinity like $u_0^{-1} = \frac{c^2 r_0}{Gm} \sim 2\gamma_0^2$.

Next we investigate the stability of the circular orbits, under a small linear perturbation [still assuming the linear theory (2.1) which is not GR]. We insert the perturbation ansatz $r = r_0 + \delta r$ and $C = C_0 + \delta C$ into Eqs. (B1) and readily obtain, to first order,

$$ \frac{\delta C}{C_0} = \frac{Gm}{r_0^2} \left( 4\gamma_0^4 - \frac{1}{2} \right) \delta r. \quad (B6) $$

On the other hand, the radial equation arises from a combination of (B1a) with (B6), and we find the harmonic oscillator,

$$ \ddot{r} + \Omega_0^2 \delta r = 0, \quad (B7) $$

with characteristic frequency

$$ \Omega_0^2 = \omega_0^2 \left[ 1 + \frac{2(\gamma_0^2 - 1)(8\gamma_0^2 - 1)(16\gamma_0^6 - 40\gamma_0^4 + 22\gamma_0^2 - 1)}{(8\gamma_0^4 - 8\gamma_0^2 + 1)^2} \right]. \quad (B8) $$

To study the stability of the circular orbits, we just have to look for the sign of the characteristic frequency squared
<figure><img src="image_4.png" /><figcaption>FIG. 1. (Left panel) The circular orbits in the two branches corresponding to $\epsilon = 1$ (PN regime) and $\epsilon = -1$ (UR regime) in Eq. (B5). The point with a minimal circular radius corresponds to $r_{\min} = (2 + \sqrt{2}) \frac{Gm}{c^2}$ and $\gamma_{\min} = \frac{1}{2}(4 + \sqrt{2})^{1/2}$. (Right panel) Stability criterion of circular orbits. The unstable region ($\Omega_0^2 < 0$) is for $1.02 \lesssim \gamma_0 \lesssim 1.21$ and surrounds the minimal radius $r_{\min}$ in the figure on the left (indeed $\gamma_{\min} \approx 1.16$).</figcaption></figure>

$\Omega_0^2$. From the right panel of Fig. 1 we see that $\Omega_0^2$ is negative, and thus the orbit is unstable, when $\gamma_- \le \gamma_0 \le \gamma_+$, where $\gamma_- \approx 1.01678$ and $\gamma_+ \approx 1.21395$. This corresponds via Eq. (B5) to the point $r_- \sim 9.45 \frac{Gm}{c^2}$ on the PN branch $\epsilon = 1$ and to the point $r_+ \sim 3.48 \frac{Gm}{c^2}$ on the UR branch $\epsilon = -1$. In particular, we note that the point with minimal circular radius $r_{\min} \approx 3.14 \frac{Gm}{c^2}$ is within the unstable region.

Therefore, we conclude that the two PN and UR regimes are in fact separated by the unstable region $\gamma_- \le \gamma_0 \le \gamma_+$. This means, for instance, that if we start with a circular orbit in the PN regime with $\gamma_0 \approx 1$ and increase the velocity, we shall never be able to reach the UR regime $\gamma_0 \to \infty$, since the orbit will first become unstable when reaching the point $\gamma_0 = \gamma_-$. Similarly, starting in the UR regime, it would be

impossible to reach the PN regime since the orbit will be unstable when $\gamma_0 = \gamma_+$.

## APPENDIX C: THE 5PN HARMONIC-COORDINATES LAGRANGIAN TO ORDER G

The total time derivative which is to be added to the construction of the 1PM Lagrangian in Sec. VA in order to ensure the equivalence with the 4PN result [25] is of the form

$$F_{\text{PN}} = \frac{1}{c^4} F_{\text{2PN}} + \frac{1}{c^6} F_{\text{3PN}} + \frac{1}{c^6} F_{\text{4PN}}, \quad (\text{C1})$$

where the coefficients are explicitly given by

$$F_{\text{2PN}} = \frac{Gm_1 m_2}{24} \left[ (n_{12}v_1)^2 (n_{12}v_2) + 14(n_{12}v_1)(v_1v_2) - 7(n_{12}v_2)v_1^2 \right] + 1 \leftrightarrow 2, \quad (\text{C2a})$$

$$F_{\text{3PN}} = -\frac{Gm_1 m_2}{240} \left[ 3(n_{12}v_1)^4(n_{12}v_2) + 6(n_{12}v_1)^3(n_{12}v_2)^2 + 20(n_{12}v_1)^3(v_1v_2) \right. \\ \left. - 15(n_{12}v_1)^2(n_{12}v_2)v_1^2 + 36(n_{12}v_1)^2(n_{12}v_2)(v_1v_2) - 24(n_{12}v_1)(n_{12}v_2)^2v_1^2 \right. \\ \left. - 90(n_{12}v_1)v_1^2(v_1v_2) + 90(n_{12}v_1)v_1^2v_2^2 + 12(n_{12}v_1)(v_1v_2)^2 \right. \\ \left. - 30(n_{12}v_2)^3v_1^2 - 132(n_{12}v_2)v_1^4 + 144(n_{12}v_2)v_1^2(v_1v_2) \right] + 1 \leftrightarrow 2, \quad (\text{C2b})$$

$$F_{\text{4PN}} = \frac{Gm_1 m_2}{2688} \left[ 15(n_{12}v_1)^6(n_{12}v_2) + 30(n_{12}v_1)^5(n_{12}v_2)^2 + 78(n_{12}v_1)^5(v_1v_2) \right. \\ - 60(n_{12}v_1)^4(n_{12}v_2)^3 - 81(n_{12}v_1)^4(n_{12}v_2)v_1^2 + 132(n_{12}v_1)^4(n_{12}v_2)(v_1v_2) \\ - 132(n_{12}v_1)^3(n_{12}v_2)^2v_1^2 + 180(n_{12}v_1)^3(n_{12}v_2)^2(v_1v_2) - 308(n_{12}v_1)^3v_1^2(v_1v_2) \\ + 308(n_{12}v_1)^3v_1^2v_2^2 - 24(n_{12}v_1)^3(v_1v_2)^2 - 162(n_{12}v_1)^2(n_{12}v_2)^3v_1^2 \\ + 159(n_{12}v_1)^2(n_{12}v_2)v_1^4 - 312(n_{12}v_1)^2(n_{12}v_2)v_1^2(v_1v_2) - 552(n_{12}v_1)^2(n_{12}v_2)v_1^2v_2^2 \\ \left. + 240(n_{12}v_1)^2(n_{12}v_2)(v_1v_2)^2 - 180(n_{12}v_1)(n_{12}v_2)^4v_1^2 + 186(n_{12}v_1)(n_{12}v_2)^2v_1^4 \right]$$
$$ 
\begin{aligned}
& -252(n_{12}v_1)(n_{12}v_2)^2 v_1^2(v_1v_2) + 738(n_{12}v_1)v_1^4(v_1v_2) - 738(n_{12}v_1)v_1^4v_2^2 \\
& -72(n_{12}v_1)v_1^2(v_1v_2)^2 + 1188(n_{12}v_1)v_1^2(v_1v_2)v_2^2 - 72(n_{12}v_1)(v_1v_2)^3 \\
& + 78(n_{12}v_2)^5v_1^2 - 228(n_{12}v_2)^3v_1^4 - 272(n_{12}v_2)^3v_1^2(v_1v_2) - 225(n_{12}v_2)v_1^6 \\
& + 588(n_{12}v_2)v_1^4(v_1v_2) + 1116(n_{12}v_2)v_1^4v_2^2 + 48(n_{12}v_2)v_1^2(v_1v_2)^2] + 1 \leftrightarrow 2.
\end{aligned} 
\quad (C2c) $$

The harmonic-coordinates Lagrangian is specified by the couple $(\lambda, q_a^i)$, where $\lambda$ represents the ordinary part of the Lagrangian and $q_a^i$ denotes the coefficients of the accelerations; see Eq. (5.1). Up to the 5PN level and linear in $G$, we have

$$ \lambda = \lambda_N + \frac{1}{c^2}\lambda_{1PN} + \frac{1}{c^4}\lambda_{2PN} + \frac{1}{c^6}\lambda_{3PN} + \frac{1}{c^8}\lambda_{4PN} + \frac{1}{c^{10}}\lambda_{5PN} + \mathcal{O}\left(\frac{1}{c^{12}}\right), \quad (C3) $$

$$ q_1^i = \frac{1}{c^4} q_{1,2PN}^i + \frac{1}{c^6} q_{1,3PN}^i + \frac{1}{c^8} q_{1,4PN}^i + \frac{1}{c^{10}} q_{1,5PN}^i + \mathcal{O}\left(\frac{1}{c^{12}}\right), \quad (C4) $$

together with $q_2^i$ obtained by label exchange. We explicitly find, after adjustment with the function $F_{PN}$ obtained previously,

$$ \lambda_N = \frac{Gm_1m_2}{r_{12}}, \quad (C5a) $$

$$ \lambda_{1PN} = -\frac{Gm_1m_2}{4r_{12}}[(n_{12}v_1)(n_{12}v_2) - 6v_1^2 + 7(v_1v_2)] + 1 \leftrightarrow 2, \quad (C5b) $$

$$ \lambda_{2PN} = \frac{Gm_1m_2}{16r_{12}} [3(n_{12}v_1)^2(n_{12}v_2)^2 + 12(n_{12}v_1)(n_{12}v_2)(v_1v_2) - 14(n_{12}v_2)^2 v_1^2 \\ 
+ 14v_1^4 - 32v_1^2(v_1v_2) + 15v_1^2v_2^2 + 2(v_1v_2)^2] + 1 \leftrightarrow 2, \quad (C5c) $$

$$ \lambda_{3PN} = -\frac{Gm_1m_2}{32r_{12}} [5(n_{12}v_1)^3(n_{12}v_2)^3 + 15(n_{12}v_1)^2(n_{12}v_2)^2(v_1v_2) + 20(n_{12}v_1)^2v_1^2v_2^2 \\ 
- 4(n_{12}v_1)(n_{12}v_2)^3v_1^2 + 22(n_{12}v_1)(n_{12}v_2)v_1^4 - 32(n_{12}v_1)(n_{12}v_2)v_1^2(v_1v_2) \\ 
+ 23(n_{12}v_1)(n_{12}v_2)v_1^2v_2^2 - 10(n_{12}v_1)(n_{12}v_2)(v_1v_2)^2 - 20(n_{12}v_2)^4v_1^2 \\ 
- 8(n_{12}v_2)^2v_1^4 - 12(n_{12}v_2)^2v_1^2(v_1v_2) - 22v_1^6 + 26v_1^4(v_1v_2) - 2v_1^4v_2^2 \\ 
+ v_1^2(v_1v_2)v_2^2 - 2(v_1v_2)^3] + 1 \leftrightarrow 2, \quad (C5d) $$

$$ \lambda_{4PN} = \frac{Gm_1m_2}{256r_{12}} [70(n_{12}v_1)^5(n_{12}v_2)^3 - 35(n_{12}v_1)^4(n_{12}v_2)^4 - 30(n_{12}v_1)^4(n_{12}v_2)^2(v_1v_2) \\ 
- 40(n_{12}v_1)^3(n_{12}v_2)^3v_1^2 + 120(n_{12}v_1)^3(n_{12}v_2)^3(v_1v_2) + 396(n_{12}v_1)^3(n_{12}v_2)v_1^2v_2^2 \\ 
- 120(n_{12}v_1)^3(n_{12}v_2)(v_1v_2)^2 - 84(n_{12}v_1)^2(n_{12}v_2)^2v_1^2v_2^2 + 44(n_{12}v_1)^2v_1^2(v_1v_2)v_2^2 \\ 
+ 40(n_{12}v_1)^2(v_1v_2)^3 - 150(n_{12}v_1)(n_{12}v_2)^5v_1^2 + 114(n_{12}v_1)(n_{12}v_2)^3v_1^4 \\ 
- 80(n_{12}v_1)(n_{12}v_2)^3v_1^2(v_1v_2) - 370(n_{12}v_1)(n_{12}v_2)v_1^4v_2^2 + 72(n_{12}v_1)(n_{12}v_2)v_1^2(v_1v_2)^2 \\ 
+ 280(n_{12}v_1)(n_{12}v_2)v_1^2(v_1v_2)v_2^2 - 48(n_{12}v_1)(n_{12}v_2)(v_1v_2)^3 - 30(n_{12}v_2)^4v_1^4 \\ 
- 106(n_{12}v_2)^4v_1^2(v_1v_2) - 88(n_{12}v_2)^2v_1^6 + 78(n_{12}v_2)^2v_1^4(v_1v_2) + 92(n_{12}v_2)^2v_1^4v_2^2 \\ 
- 120(n_{12}v_2)^2v_1^2(v_1v_2)^2 + 150v_1^8 - 320v_1^6(v_1v_2) + 160v_1^6v_2^2 + 192v_1^4(v_1v_2)^2 \\ 
- 198v_1^4(v_1v_2)v_2^2 + 3v_1^4v_2^4 + 8v_1^2(v_1v_2)^3 - 8v_1^2(v_1v_2)^2v_2^2 + 8(v_1v_2)^4] + 1 \leftrightarrow 2, \quad (C5e) $$
$$
\begin{aligned}
\lambda_{5\text{PN}} = & -\frac{Gm_1m_2}{4608r_{12}} [126(n_{12}v_2)(n_{12}v_1)^9 + 126(n_{12}v_2)^2(n_{12}v_1)^8 + 434(v_1v_2)(n_{12}v_1)^8 \\
& + 126(n_{12}v_2)^3(n_{12}v_1)^7 - 714(n_{12}v_2)v_1^2(n_{12}v_1)^7 + 336(n_{12}v_2)(v_1v_2)(n_{12}v_1)^7 \\
& + 126(n_{12}v_2)^4(n_{12}v_1)^6 + 192(v_1v_2)^2(n_{12}v_1)^6 - 602(n_{12}v_2)^2v_1^2(n_{12}v_1)^6 \\
& + 266(n_{12}v_2)^2(v_1v_2)(n_{12}v_1)^6 - 1846v_1^2(v_1v_2)(n_{12}v_1)^6 + 1846v_1^2v_2^2(n_{12}v_1)^6 \\
& + 63(n_{12}v_2)^5(n_{12}v_1)^5 + 1526(n_{12}v_2)v_1^4(n_{12}v_1)^5 - 388(n_{12}v_2)(v_1v_2)^2(n_{12}v_1)^5 \\
& - 504(n_{12}v_2)^3v_1^2(n_{12}v_1)^5 + 224(n_{12}v_2)^3(v_1v_2)(n_{12}v_1)^5 - 968(n_{12}v_2)v_1^2(v_1v_2)(n_{12}v_1)^5 \\
& + 886(n_{12}v_2)v_1^2v_2^2(n_{12}v_1)^5 + 76(v_1v_2)^3(n_{12}v_1)^4 + 1010(n_{12}v_2)^2v_1^4(n_{12}v_1)^4 \\
& - 520(n_{12}v_2)^2(v_1v_2)^2(n_{12}v_1)^4 - 820v_1^2(v_1v_2)^2(n_{12}v_1)^4 - 420(n_{12}v_2)^4v_1^2(n_{12}v_1)^4 \\
& + 105(n_{12}v_2)^4(v_1v_2)(n_{12}v_1)^4 + 3082v_1^4(v_1v_2)(n_{12}v_1)^4 - 440(n_{12}v_2)^2v_1^2(v_1v_2)(n_{12}v_1)^4 \\
& - 3082v_1^4v_2^2(n_{12}v_1)^4 + 850(n_{12}v_2)^2v_1^2v_2^2(n_{12}v_1)^4 + 2874v_1^2(v_1v_2)v_2^2(n_{12}v_1)^4 \\
& - 1546(n_{12}v_2)v_1^6(n_{12}v_1)^3 + 256(n_{12}v_2)(v_1v_2)^3(n_{12}v_1)^3 + 650(n_{12}v_2)^3v_1^4(n_{12}v_1)^3 \\
& - 280(n_{12}v_2)^3(v_1v_2)^2(n_{12}v_1)^3 + 736(n_{12}v_2)v_1^2(v_1v_2)^2(n_{12}v_1)^3 - 350(n_{12}v_2)^5v_1^2(n_{12}v_1)^3 \\
& + 992(n_{12}v_2)v_1^4(v_1v_2)(n_{12}v_1)^3 - 160(n_{12}v_2)^3v_1^2(v_1v_2)(n_{12}v_1)^3 - 1160(n_{12}v_2)v_1^4v_2^2(n_{12}v_1)^3 \\
& + 420(n_{12}v_2)^3v_1^2v_2^2(n_{12}v_1)^3 + 144(n_{12}v_2)v_1^2(v_1v_2)v_2^2(n_{12}v_1)^3 - 16(v_1v_2)^4(n_{12}v_1)^2 \\
& - 738(n_{12}v_2)^2v_1^6(n_{12}v_1)^2 + 168(n_{12}v_2)^2(v_1v_2)^3(n_{12}v_1)^2 - 96v_1^2(v_1v_2)^3(n_{12}v_1)^2 \\
& + 410(n_{12}v_2)^4v_1^4(n_{12}v_1)^2 + 1272v_1^4(v_1v_2)^2(n_{12}v_1)^2 + 576(n_{12}v_2)^2v_1^2(v_1v_2)^2(n_{12}v_1)^2 \\
& + 2640v_1^4v_2^4(n_{12}v_1)^2 - 294(n_{12}v_2)^6v_1^2(n_{12}v_1)^2 - 2550v_1^6(v_1v_2)(n_{12}v_1)^2 \\
& + 222(n_{12}v_2)^2v_1^4(v_1v_2)(n_{12}v_1)^2 - 50(n_{12}v_2)^4v_1^2(v_1v_2)(n_{12}v_1)^2 + 2550v_1^6v_2^2(n_{12}v_1)^2 \\
& - 870(n_{12}v_2)^2v_1^4v_2^2(n_{12}v_1)^2 + 1416v_1^2(v_1v_2)^2v_2^2(n_{12}v_1)^2 - 3912v_1^4(v_1v_2)v_2^2(n_{12}v_1)^2 \\
& + 12(n_{12}v_2)^2v_1^2(v_1v_2)v_2^2(n_{12}v_1)^2 + 758(n_{12}v_2)v_1^8(n_{12}v_1) - 24(n_{12}v_2)(v_1v_2)^4(n_{12}v_1) \\
& - 348(n_{12}v_2)^3v_1^6(n_{12}v_1) - 192(n_{12}v_2)v_1^2(v_1v_2)^3(n_{12}v_1) + 260(n_{12}v_2)^5v_1^4(n_{12}v_1) \\
& - 372(n_{12}v_2)v_1^4(v_1v_2)^2(n_{12}v_1) + 296(n_{12}v_2)^3v_1^2(v_1v_2)^2(n_{12}v_1) + 351(n_{12}v_2)v_1^4v_2^4(n_{12}v_1) \\
& - 252(n_{12}v_2)^7v_1^2(n_{12}v_1) - 504(n_{12}v_2)v_1^6(v_1v_2)(n_{12}v_1) + 64(n_{12}v_2)^3v_1^4(v_1v_2)(n_{12}v_1) \\
& - 56(n_{12}v_2)^5v_1^2(v_1v_2)(n_{12}v_1) + 730(n_{12}v_2)v_1^6v_2^2(n_{12}v_1) - 694(n_{12}v_2)^3v_1^4v_2^2(n_{12}v_1) \\
& - 168(n_{12}v_2)v_1^2(v_1v_2)^2v_2^2(n_{12}v_1) - 264(n_{12}v_2)v_1^4(v_1v_2)v_2^2(n_{12}v_1) - 2394v_1^{10} \\
& - 8(v_1v_2)^5 + 1316(n_{12}v_2)^2v_1^8 - 32v_1^2(v_1v_2)^4 - 844(n_{12}v_2)^4v_1^6 - 36v_1^4(v_1v_2)^3 \\
& - 56(n_{12}v_2)^2v_1^2(v_1v_2)^3 + 596(n_{12}v_2)^6v_1^4 - 3292v_1^6(v_1v_2)^2 + 888(n_{12}v_2)^2v_1^4(v_1v_2)^2 \\
& - 524(n_{12}v_2)^4v_1^2(v_1v_2)^2 - 3450v_1^6v_2^4 + 3465v_1^4(v_1v_2)v_2^4 - 434(n_{12}v_2)^8v_1^2 \\
& + 5770v_1^8(v_1v_2) - 2132(n_{12}v_2)^2v_1^6(v_1v_2) + 1292(n_{12}v_2)^4v_1^4(v_1v_2) - 788(n_{12}v_2)^6v_1^2(v_1v_2) \\
& - 3250v_1^8v_2^2 + 2178(n_{12}v_2)^2v_1^2v_2^2 - 8v_1^2(v_1v_2)^3v_2^2 - 2054(n_{12}v_2)^4v_1^4v_2^2 \\
& - 3444v_1^4(v_1v_2)^2v_2^2 + 6742v_1^6(v_1v_2)v_2^2 - 3498(n_{12}v_2)^2v_1^4(v_1v_2)v_2^2] + 1 \leftrightarrow 2,
\end{aligned}
\tag{C5f}
$$

and

$$
q_{1,2\text{PN}}^i = -\frac{Gm_1m_2}{8} [((n_{12}v_2)^2 - 7v_1^2)n_{12}^i + 14(n_{12}v_2)v_2^i], \qquad (C6a)
$$
$$ q_{1,3\text{PN}}^i = -\frac{Gm_1m_2}{48} \left[ \begin{aligned} & (30(n_{12}v_1)^2v_2^2 - 6(n_{12}v_1)(n_{12}v_2)^3 + 24(n_{12}v_1)(n_{12}v_2)v_2^2 \\ & \quad - 3(n_{12}v_2)^4 - 18(n_{12}v_2)^2(v_1v_2) + 15(n_{12}v_2)^2v_2^2)n_{12}^i \\ & + (-132(n_{12}v_2)v_1^2 + 96(n_{12}v_2)(v_1v_2) - 90(n_{12}v_2)v_2^2)v_1^i \\ & + (-18(n_{12}v_1)(n_{12}v_2)^2 - 20(n_{12}v_2)^3 + 48(n_{12}v_2)v_1^2 - 12(n_{12}v_2)(v_1v_2) + 90(n_{12}v_2)v_2^2)v_2^i \end{aligned} \right], \quad \text{(C6b)} $$

$$ q_{1,4\text{PN}}^i = -\frac{Gm_1m_2}{384} \left[ \begin{aligned} & (60(n_{12}v_1)^3(n_{12}v_2)^3 - 180(n_{12}v_1)^3(n_{12}v_2)v_2^2 - 162(n_{12}v_1)^2(n_{12}v_2)^2v_2^2 \\ & \quad - 204(n_{12}v_1)^2(v_1v_2)v_2^2 + 30(n_{12}v_1)(n_{12}v_2)^5 + 120(n_{12}v_1)(n_{12}v_2)^3(v_1v_2) \\ & \quad - 132(n_{12}v_1)(n_{12}v_2)^3v_2^2 + 276(n_{12}v_1)(n_{12}v_2)v_1^2v_2^2 - 120(n_{12}v_1)(n_{12}v_2)(v_1v_2)^2 \\ & \quad - 168(n_{12}v_1)(n_{12}v_2)(v_1v_2)v_2^2 + 186(n_{12}v_1)(n_{12}v_2)v_2^4 + 15(n_{12}v_2)^6 \\ & \quad + 66(n_{12}v_2)^4(v_1v_2) - 81(n_{12}v_2)^4v_2^2 - 156(n_{12}v_2)^2(v_1v_2)v_2^2 \\ & \quad + 159(n_{12}v_2)^2v_2^4 + 294(v_1v_2)v_2^4 - 225v_2^6)n_{12}^i \\ & \quad + (276(n_{12}v_1)^2(n_{12}v_2)v_2^2 - 78(n_{12}v_2)^5 + 228(n_{12}v_2)^3v_1^2 \\ & \quad + 308(n_{12}v_2)^3v_2^2 - 1116(n_{12}v_2)v_1^2v_2^2 - 24(n_{12}v_2)(v_1v_2)^2 \\ & \quad + 792(n_{12}v_2)(v_1v_2)v_2^2 - 738(n_{12}v_2)v_2^4)v_1^i \\ & \quad + (-68(n_{12}v_1)^3v_2^2 + 60(n_{12}v_1)^2(n_{12}v_2)^3 - 120(n_{12}v_1)^2(n_{12}v_2)(v_1v_2) \\ & \quad - 84(n_{12}v_1)^2(n_{12}v_2)v_2^2 + 66(n_{12}v_1)(n_{12}v_2)^4 - 156(n_{12}v_1)(n_{12}v_2)^2v_2^2 \\ & \quad + 294(n_{12}v_1)v_2^4 + 78(n_{12}v_2)^5 - 24(n_{12}v_2)^3(v_1v_2) \\ & \quad - 308(n_{12}v_2)^3v_2^2 - 24(n_{12}v_2)v_1^2(v_1v_2) + 396(n_{12}v_2)v_1^2v_2^2 \\ & \quad - 72(n_{12}v_2)(v_1v_2)^2 - 72(n_{12}v_2)(v_1v_2)v_2^2 + 738(n_{12}v_2)v_2^4)v_2^i \end{aligned} \right], \quad \text{(C6c)} $$

$$ q_{1,5\text{PN}}^i = \frac{Gm_1m_2}{5760} \left[ \begin{aligned} & (140(n_{12}v_1)^7(n_{12}v_2) + 245(n_{12}v_1)^6(n_{12}v_2)^2 + 560(n_{12}v_1)^6(v_1v_2) \\ & \quad - 560(n_{12}v_1)^6v_2^2 + 315(n_{12}v_1)^5(n_{12}v_2)^3 - 645(n_{12}v_1)^5(n_{12}v_2)v_1^2 \\ & \quad + 780(n_{12}v_1)^5(n_{12}v_2)(v_1v_2) - 810(n_{12}v_1)^5(n_{12}v_2)v_2^2 + 350(n_{12}v_1)^4(n_{12}v_2)^4 \\ & \quad - 900(n_{12}v_1)^4(n_{12}v_2)^2v_1^2 + 825(n_{12}v_1)^4(n_{12}v_2)^2(v_1v_2) - 1050(n_{12}v_1)^4(n_{12}v_2)^2v_2^2 \\ & \quad - 1855(n_{12}v_1)^4v_1^2(v_1v_2) + 1855(n_{12}v_1)^4v_1^2v_2^2 - 190(n_{12}v_1)^4(v_1v_2)^2 \\ & \quad - 690(n_{12}v_1)^4(v_1v_2)v_2^2 + 880(n_{12}v_1)^4v_2^4 + 350(n_{12}v_1)^3(n_{12}v_2)^5 \\ & \quad - 900(n_{12}v_1)^3(n_{12}v_2)^3v_1^2 + 800(n_{12}v_1)^3(n_{12}v_2)^3(v_1v_2) - 1250(n_{12}v_1)^3(n_{12}v_2)^3v_2^2 \\ & \quad + 1010(n_{12}v_1)^3(n_{12}v_2)v_1^4 - 1600(n_{12}v_1)^3(n_{12}v_2)v_1^2(v_1v_2) + 2010(n_{12}v_1)^3(n_{12}v_2)v_1^2v_2^2 \\ & \quad - 900(n_{12}v_1)^3(n_{12}v_2)(v_1v_2)^2 - 560(n_{12}v_1)^3(n_{12}v_2)(v_1v_2)v_2^2 + 1300(n_{12}v_1)^3(n_{12}v_2)v_2^4 \\ & \quad + 315(n_{12}v_1)^2(n_{12}v_2)^6 - 750(n_{12}v_1)^2(n_{12}v_2)^4v_1^2 + 750(n_{12}v_1)^2(n_{12}v_2)^4(v_1v_2) \\ & \quad - 1350(n_{12}v_1)^2(n_{12}v_2)^4v_2^2 + 975(n_{12}v_1)^2(n_{12}v_2)^2v_1^4 - 900(n_{12}v_1)^2(n_{12}v_2)^2v_1^2(v_1v_2) \\ & \quad + 1920(n_{12}v_1)^2(n_{12}v_2)^2v_1^2v_2^2 - 1200(n_{12}v_1)^2(n_{12}v_2)^2(v_1v_2)^2 - 750(n_{12}v_1)^2(n_{12}v_2)^2(v_1v_2)v_2^2 \\ & \quad + 1845(n_{12}v_1)^2(n_{12}v_2)^2v_2^4 + 2250(n_{12}v_1)^2v_1^4(v_1v_2) - 2250(n_{12}v_1)^2v_1^4v_2^2 \\ & \quad + 240(n_{12}v_1)^2v_1^2(v_1v_2)^2 + 1950(n_{12}v_1)^2v_1^2(v_1v_2)v_2^2 - 2190(n_{12}v_1)^2v_1^2v_2^4 \\ & \quad + 60(n_{12}v_1)^2(v_1v_2)^3 + 120(n_{12}v_1)^2(v_1v_2)^2v_2^2 + 1200(n_{12}v_1)^2(v_1v_2)v_2^4 \end{aligned} \right] $$
$$
\begin{aligned}
& -1380(n_{12}v_1)^2 v_2^6 + 245(n_{12}v_1)(n_{12}v_2)^7 - 525(n_{12}v_1)(n_{12}v_2)^5 v_1^2 \\
& + 660(n_{12}v_1)(n_{12}v_2)^5(v_1 v_2) - 1260(n_{12}v_1)(n_{12}v_2)^5 v_2^2 + 615(n_{12}v_1)(n_{12}v_2)^3 v_1^4 \\
& - 400(n_{12}v_1)(n_{12}v_2)^3 v_1^2 (v_1 v_2) + 1600(n_{12}v_1)(n_{12}v_2)^3 v_1^2 v_2^2 - 1000(n_{12}v_1)(n_{12}v_2)^3 (v_1 v_2)^2 \\
& - 1200(n_{12}v_1)(n_{12}v_2)^3 (v_1 v_2) v_2^2 + 2275(n_{12}v_1)(n_{12}v_2)^3 v_2^4 - 615(n_{12}v_1)(n_{12}v_2) v_1^6 \\
& + 1020(n_{12}v_1)(n_{12}v_2) v_1^4 (v_1 v_2) - 1560(n_{12}v_1)(n_{12}v_2) v_1^4 v_2^2 + 840(n_{12}v_1)(n_{12}v_2) v_1^2 (v_1 v_2)^2 \\
& + 960(n_{12}v_1)(n_{12}v_2) v_1^2 (v_1 v_2) v_2^2 - 2025(n_{12}v_1)(n_{12}v_2) v_1^2 v_2^4 + 320(n_{12}v_1)(n_{12}v_2) (v_1 v_2)^3 \\
& + 900(n_{12}v_1)(n_{12}v_2) (v_1 v_2)^2 v_2^2 + 1140(n_{12}v_1)(n_{12}v_2) (v_1 v_2) v_2^4 - 2030(n_{12}v_1)(n_{12}v_2) v_2^6 \\
& + 140(n_{12}v_2)^8 - 270(n_{12}v_2)^6 v_1^2 + 455(n_{12}v_2)^6 (v_1 v_2) \\
& - 860(n_{12}v_2)^6 v_2^2 + 260(n_{12}v_2)^4 v_1^4 - 175(n_{12}v_2)^4 v_1^2 (v_1 v_2) \\
& + 1005(n_{12}v_2)^4 v_1^2 v_2^2 - 450(n_{12}v_2)^4 (v_1 v_2)^2 - 1400(n_{12}v_2)^4 (v_1 v_2) v_2^2 \\
& + 2020(n_{12}v_2)^4 v_1^4 - 290(n_{12}v_2)^2 v_1^6 + 285(n_{12}v_2)^2 v_1^4 (v_1 v_2) \\
& - 810(n_{12}v_2)^2 v_1^4 v_2^2 + 360(n_{12}v_2)^2 v_1^2 (v_1 v_2)^2 + 600(n_{12}v_2)^2 v_1^2 (v_1 v_2) v_2^2 \\
& - 1560(n_{12}v_2)^2 v_1^2 v_2^4 + 200(n_{12}v_2)^2 (v_1 v_2)^3 + 840(n_{12}v_2)^2 (v_1 v_2)^2 v_2^2 \\
& + 1785(n_{12}v_2)^2 (v_1 v_2) v_2^4 - 2460(n_{12}v_2)^2 v_2^6 - 1245 v_1^6 (v_1 v_2) + 1245 v_1^6 v_2^2 \\
& + 90 v_1^4 (v_1 v_2)^2 - 1980 v_1^4 (v_1 v_2) v_2^2 + 1890 v_1^4 v_2^4 + 120 v_1^2 (v_1 v_2)^3 - 2475 v_1^2 (v_1 v_2) v_2^4 \\
& + 2355 v_1^2 v_2^6 + 80 (v_1 v_2)^4 + 100 (v_1 v_2)^3 v_2^3 + 30 (v_1 v_2)^2 v_2^4 - 2870 (v_1 v_2) v_2^6 + 2660 v_2^8) f n_{12}^i \\
& + (-215(n_{12}v_1)^6 (n_{12}v_2) - 360(n_{12}v_1)^5 (n_{12}v_2)^2 - 742(n_{12}v_1)^5 (v_1 v_2) \\
& + 742(n_{12}v_1)^5 v_2^2 - 450(n_{12}v_1)^4 (n_{12}v_2)^3 + 1010(n_{12}v_1)^4 (n_{12}v_2) v_1^2 \\
& - 800(n_{12}v_1)^4 (n_{12}v_2) (v_1 v_2) + 1005(n_{12}v_1)^4 (n_{12}v_2) v_2^2 - 500(n_{12}v_1)^3 (n_{12}v_2)^4 \\
& + 1300(n_{12}v_1)^3 (n_{12}v_2)^2 v_1^2 - 600(n_{12}v_1)^3 (n_{12}v_2)^2 (v_1 v_2) + 1280(n_{12}v_1)^3 (n_{12}v_2)^2 v_2^2 \\
& + 3000(n_{12}v_1)^3 v_1^2 (v_1 v_2) - 3000(n_{12}v_1)^3 v_1^2 v_2^2 + 160(n_{12}v_1)^3 (v_1 v_2)^2 \\
& + 1300(n_{12}v_1)^3 (v_1 v_2) v_2^2 - 1460(n_{12}v_1)^3 v_2^4 - 525(n_{12}v_1)^2 (n_{12}v_2)^5 \\
& + 1230(n_{12}v_1)^2 (n_{12}v_2)^3 v_1^2 - 400(n_{12}v_1)^2 (n_{12}v_2)^3 (v_1 v_2) + 1600(n_{12}v_1)^2 (n_{12}v_2)^3 v_2^2 \\
& - 1845(n_{12}v_1)^2 (n_{12}v_2) v_1^4 + 2040(n_{12}v_1)^2 (n_{12}v_2) v_1^2 (v_1 v_2) - 3120(n_{12}v_1)^2 (n_{12}v_2) v_1^2 v_2^2 \\
& + 840(n_{12}v_1)^2 (n_{12}v_2) (v_1 v_2)^2 + 960(n_{12}v_1)^2 (n_{12}v_2) (v_1 v_2) v_2^2 - 2025(n_{12}v_1)^2 (n_{12}v_2) v_2^4 \\
& - 540(n_{12}v_1)(n_{12}v_2)^6 + 1040(n_{12}v_1)(n_{12}v_2)^4 v_1^2 - 350(n_{12}v_1)(n_{12}v_2)^4 (v_1 v_2) \\
& + 2010(n_{12}v_1)(n_{12}v_2)^4 v_2^2 - 1740(n_{12}v_1)(n_{12}v_2)^2 v_1^4 + 1140(n_{12}v_1)(n_{12}v_2)^2 v_1^2 (v_1 v_2) \\
& - 3240(n_{12}v_1)(n_{12}v_2)^2 v_1^2 v_2^2 + 720(n_{12}v_1)(n_{12}v_2)^2 (v_1 v_2)^2 + 1200(n_{12}v_1)(n_{12}v_2)^2 (v_1 v_2) v_2^2 \\
& - 3120(n_{12}v_1)(n_{12}v_2)^2 v_2^4 - 7470(n_{12}v_1) v_1^4 (v_1 v_2) + 7470(n_{12}v_1) v_1^4 v_2^2 \\
& + 360(n_{12}v_1) v_1^2 (v_1 v_2)^2 - 7920(n_{12}v_1) v_1^2 (v_1 v_2) v_2^2 + 7560(n_{12}v_1) v_1^2 v_2^4 \\
& + 240(n_{12}v_1) (v_1 v_2)^3 - 4950(n_{12}v_1) (v_1 v_2) v_2^4 + 4710(n_{12}v_1) v_2^6 \\
& - 560(n_{12}v_2)^7 + 880(n_{12}v_2)^5 v_1^2 - 552(n_{12}v_2)^5 (v_1 v_2) \\
& + 2597(n_{12}v_2)^5 v_2^2 - 1380(n_{12}v_2)^3 v_1^4 + 1280(n_{12}v_2)^3 v_1^2 (v_1 v_2) \\
& - 3650(n_{12}v_2)^3 v_1^2 v_2^2 + 100(n_{12}v_2)^3 (v_1 v_2)^2 + 2600(n_{12}v_2)^3 (v_1 v_2) v_2^2 \\
& - 5250(n_{12}v_2)^3 v_2^4 + 2660(n_{12}v_2) v_1^6 - 4920(n_{12}v_2) v_1^4 (v_1 v_2) \\
& + 7065(n_{12}v_2) v_1^4 v_2^2 + 60(n_{12}v_2) v_1^2 (v_1 v_2)^2 - 7920(n_{12}v_2) v_1^2 (v_1 v_2) v_2^2 \\
& + 9450(n_{12}v_2) v_1^2 v_2^4 + 160(n_{12}v_2) (v_1 v_2)^3 - 7920(n_{12}v_2) (v_1 v_2) v_2^4 + 8715(n_{12}v_2) v_2^6 v_1^i
\end{aligned}
$$
$$
\begin{aligned}
& + (80(n_{12}v_1)^7 + 130(n_{12}v_1)^6(n_{12}v_2) + 165(n_{12}v_1)^5(n_{12}v_2)^2 - 371(n_{12}v_1)^5v_1^2 \\
& - 76(n_{12}v_1)^5(v_1v_2) - 138(n_{12}v_1)^5v_2^2 + 200(n_{12}v_1)^4(n_{12}v_2)^3 \\
& - 400(n_{12}v_1)^4(n_{12}v_2)v_1^2 - 450(n_{12}v_1)^4(n_{12}v_2)(v_1v_2) - 140(n_{12}v_1)^4(n_{12}v_2)v_2^2 \\
& + 250(n_{12}v_1)^3(n_{12}v_2)^4 - 300(n_{12}v_1)^3(n_{12}v_2)^2v_1^2 - 800(n_{12}v_1)^3(n_{12}v_2)^2(v_1v_2) \\
& - 250(n_{12}v_1)^3(n_{12}v_2)^2v_2^2 + 750(n_{12}v_1)^3v_1^4 + 160(n_{12}v_1)^3v_1^2(v_1v_2) \\
& + 650(n_{12}v_1)^3v_1^2v_2^2 + 60(n_{12}v_1)^3(v_1v_2)^2 + 80(n_{12}v_1)^3(v_1v_2)v_2^2 \\
& + 400(n_{12}v_1)^3v_2^4 + 330(n_{12}v_1)^2(n_{12}v_2)^5 - 200(n_{12}v_1)^2(n_{12}v_2)^3v_1^2 \\
& - 1000(n_{12}v_1)^2(n_{12}v_2)^3(v_1v_2) - 600(n_{12}v_1)^2(n_{12}v_2)^3v_2^2 + 510(n_{12}v_1)^2(n_{12}v_2)v_1^4 \\
& + 840(n_{12}v_1)^2(n_{12}v_2)v_1^2(v_1v_2) + 480(n_{12}v_1)^2(n_{12}v_2)v_1^2v_2^2 + 480(n_{12}v_1)^2(n_{12}v_2)(v_1v_2)^2 \\
& + 900(n_{12}v_1)^2(n_{12}v_2)(v_1v_2)v_2^2 + 570(n_{12}v_1)^2(n_{12}v_2)v_2^4 + 455(n_{12}v_1)(n_{12}v_2)^6 \\
& - 175(n_{12}v_1)(n_{12}v_2)^4v_1^2 - 900(n_{12}v_1)(n_{12}v_2)^4(v_1v_2) - 1400(n_{12}v_1)(n_{12}v_2)^4v_2^2 \\
& + 285(n_{12}v_1)(n_{12}v_2)^2v_1^4 + 720(n_{12}v_1)(n_{12}v_2)^2v_1^2(v_1v_2) + 600(n_{12}v_1)(n_{12}v_2)^2v_1^2v_2^2 \\
& + 600(n_{12}v_1)(n_{12}v_2)^2(v_1v_2)^2 + 1680(n_{12}v_1)(n_{12}v_2)^2(v_1v_2)v_2^2 + 1785(n_{12}v_1)(n_{12}v_2)^2v_2^4 \\
& - 1245(n_{12}v_1)v_1^6 + 180(n_{12}v_1)v_1^4(v_1v_2) - 1980(n_{12}v_1)v_1^4v_2^2 \\
& + 360(n_{12}v_1)v_1^2(v_1v_2)^2 - 2475(n_{12}v_1)v_1^2v_2^4 + 320(n_{12}v_1)(v_1v_2)^3 \\
& + 300(n_{12}v_1)(v_1v_2)^2v_2^2 + 60(n_{12}v_1)(v_1v_2)v_2^4 - 2870(n_{12}v_1)v_2^6 \\
& + 640(n_{12}v_2)^7 - 276(n_{12}v_2)^5v_1^2 - 266(n_{12}v_2)^5(v_1v_2) \\
& - 2968(n_{12}v_2)^5v_2^2 + 320(n_{12}v_2)^3v_1^4 + 100(n_{12}v_2)^3v_1^2(v_1v_2) \\
& + 1300(n_{12}v_2)^3v_1^2v_2^2 + 120(n_{12}v_2)^3(v_1v_2)^2 + 560(n_{12}v_2)^3(v_1v_2)v_2^2 \\
& + 6000(n_{12}v_2)^3v_2^4 - 820(n_{12}v_2)v_1^6 + 30(n_{12}v_2)v_1^4(v_1v_2) - 1980(n_{12}v_2)v_1^4v_2^2 \\
& + 240(n_{12}v_2)v_1^2(v_1v_2)^2 - 3960(n_{12}v_2)v_1^2v_2^4 + 400(n_{12}v_2)(v_1v_2)^3 \\
& + 720(n_{12}v_2)(v_1v_2)^2v_2^2 + 630(n_{12}v_2)(v_1v_2)v_2^4 - 9960(n_{12}v_2)v_2^6v_2^i.
\end{aligned}
\qquad (C6d)
$$

[1] A. Einstein, Sitzber. Preuss. Akad. Wiss. Berlin **1**, 688 (1916).

[2] B. Bertotti, Nuovo Cimento **4**, 898 (1956).

[3] P. Havas, Phys. Rev. **108**, 1351 (1957).

[4] B. Bertotti and J. Plebanski, Ann. Phys. (N.Y.) **11**, 169 (1960).

[5] K. Westpfahl and M. Goller, Lett. Nuovo Cimento **26**, 573 (1979).

[6] P. Havas and J. Goldberg, Phys. Rev. **128**, 398 (1962).

[7] M. Portilla, J. Phys. A **13**, 3677 (1980).

[8] K. Westpfahl, Fortschr. Phys. **33**, 417 (1985).

[9] T. Damour and G. Esposito-Farèse, Phys. Rev. D **53**, 5541 (1996).

[10] K. Thorne and S. Kovàcs, Astrophys. J. **200**, 245 (1975).

[11] R. Crowley and K. Thorne, Astrophys. J. **215**, 624 (1977).

[12] K. Westpfahl and H. Hoyler, Lett. Nuovo Cimento **27**, 581 (1980).

[13] L. Bel, T. Damour, N. Deruelle, J. Ibanez, and J. Martin, Gen. Relativ. Gravit. **13**, 963 (1981).

[14] L. Blanchet, Living Rev. Relativity **17**, 2 (2014).

[15] P. Jaranowski and G. Schäfer, Phys. Rev. D **86**, 061503(R) (2012).

[16] P. Jaranowski and G. Schäfer, Phys. Rev. D **87**, 081503(R) (2013).

[17] T. Damour, P. Jaranowski, and G. Schäfer, Phys. Rev. D **89**, 064058 (2014).

[18] T. Damour, P. Jaranowski, and G. Schäfer, Phys. Rev. D **91**, 084024 (2015).

[19] P. Jaranowski and G. Schäfer, Phys. Rev. D **92**, 124043 (2015).

[20] T. Damour, P. Jaranowski, and G. Schäfer, Phys. Rev. D **93**, 084014 (2016).

[21] L. Bernard, L. Blanchet, A. Bohé, G. Faye, and S. Marsat, Phys. Rev. D **93**, 084037 (2016).
[22] L. Bernard, L. Blanchet, A. Bohé, G. Faye, and S. Marsat, Phys. Rev. D **95**, 044026 (2017).

[23] L. Bernard, L. Blanchet, A. Bohé, G. Faye, and S. Marsat, Phys. Rev. D **96**, 104043 (2017).

[24] T. Marchand, L. Bernard, L. Blanchet, and G. Faye, Phys. Rev. D **97**, 044023 (2018).

[25] L. Bernard, L. Blanchet, G. Faye, and T. Marchand, Phys. Rev. D **97**, 044037 (2018).

[26] S. Foffa and R. Sturani, Phys. Rev. D **87**, 064011 (2013).

[27] S. Foffa and R. Sturani, Phys. Rev. D **87**, 044056 (2013).

[28] C. R. Galley, A. K. Leibovich, R. A. Porto, and A. Ross, Phys. Rev. D **93**, 124010 (2016).

[29] S. Foffa, Phys. Rev. D **89**, 024019 (2014).

[30] T. Ledvinka, G. Schäfer, and J. Bičák, Phys. Rev. Lett. **100**, 251101 (2008).

[31] T. Damour, Phys. Rev. D **94**, 104015 (2016).

[32] T. Damour, Phys. Rev. D **97**, 044038 (2018).

[33] J. Jackson, *Classical Electromagnetism*, 3rd ed. (John Wiley, New York, 1999).

[34] T. Damour and N. Deruelle, C. R. Acad. Sc. Paris **293**, 537 (1981).

[35] T. Damour and G. Schäfer, Gen. Relativ. Gravit. **17**, 879 (1985).

[36] G. Schäfer, Phys. Lett. **A100**, 128 (1984).

[37] V. de Andrade, L. Blanchet, and G. Faye, Classical Quantum Gravity **18**, 753 (2001).

[38] R. Woodard, Lect. Notes Phys. **720**, 403 (2007).

[39] J. Hadamard, *Le Problème de Cauchy et les Équations aux Dérivées Partielles Linéaires Hyperboliques* (Hermann, Paris, 1932).

[40] T. Damour and G. Schäfer, J. Math. Phys. (N.Y.) **32**, 127 (1991).

[41] M. Abramowitz and I. Stegun, *Handbook of Mathematical Functions* (Dover Publications, New York, 1966).