Schäfer (1997), Pati and Will (2002), Königsdörffer et al. (2003), Nissanke and Blanchet (2005), Itoh (2009).

Though the successive post-Newtonian approximations are really a consequence of general relativity, the final equations of motion must be interpreted in a Newtonian-like fashion. That is, once a convenient general-relativistic (Cartesian) coordinate system is chosen, we should express the results in terms of the coordinate positions, velocities, and accelerations of the bodies, and view the trajectories of the particles as taking place in the absolute Euclidean space of Newton. But because the equations of motion are actually relativistic, they must stay manifestly invariant—in harmonic coordinates—when we perform a global post-Newtonian-expanded Lorentz transformation. Furthermore they must possess the correct “perturbative” limit, given by the geodesics of the (post-Newtonian-expanded) Schwarzschild metric, when one of the masses tends to zero, and be conservative, i.e., to admit a Lagrangian or Hamiltonian formulation, when the gravitational radiation reaction is turned off.

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

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

The harmonic-coordinates equations of motion depend on an arbitrary length scale $r'_0$ associated with the logarithms present at the 3PN order. This scale is merely linked with the choice of coordinates—it can be referred to as a “gauge constant”. Initially two gauge scales $r'_1$ and $r'_2$ were introduced by Blanchet and Faye (2001a); however it is useful to choose $r'_1 = r'_2 = r'_0$ for a reason given in Sect. VI of Marchand et al. (2020). One can show that the logarithms appearing in Eq. (339), together with the constant $r'_0$ therein, can be removed by a coordinate transformation, or equivalently at that order by a linear gauge transformation. Therefore, the arbitrariness in the choice of the constant $r'_0$ does not affect the physics, because physical results must be gauge invariant. Indeed we shall verify that $r'_0$ cancels out in the final results.<fnref n="69" />

When retaining the “even” relativistic corrections at the 1PN, 2PN and 3PN orders, and neglecting the “odd” radiation reaction terms at the 2.5PN and 3.5PN orders, we find that the equations of motion admit a conserved energy (and a Lagrangian as well); that energy can be straightforwardly obtained by guess-work starting from Eq. (339). When taking into account the radiation reaction, the energy will enter the left-hand side of the “flux-balance” equation, with the gravitational-wave energy flux in the right-hand side. Several equivalent choices can be made at this stage. The “canonical” choice is to present the flux into manifestly positive definite form, such as Eq. (5) at Newtonian order, or, more generally, Eqs. (179). Then we have to transfer to the left-hand side of the balance equation some odd-parity terms in the form of total time derivatives. Such time derivative represent the analogue of the Schott (1915) terms in electromagnetism. Thus we include a contribution to the energy $E$ at the 2.5PN order, coming from the decision to express the flux into canonical form, i.e., the familiar Einstein quadrupole formula at 2.5PN order. We obtain

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

<footnote n="69">⁶⁹ Notice the dependence upon the irrational number $\pi^2$. Technically, the $\pi^2$ terms arise from non-linear interactions involving some integrals in 3 dimensions such as</footnote>
$$
\begin{aligned}
E = & \frac{m_1 v_1^2}{2} - \frac{G m_1 m_2}{2 r_{12}} \\
& + \frac{1}{c^2} \left\{ \frac{G^2 m_1^2 m_2}{2 r_{12}^2} + \frac{3 m_1 v_1^4}{8} \right. \\
& \left. + \frac{G m_1 m_2}{r_{12}} \left( -\frac{1}{4} (n_{12} v_1) (n_{12} v_2) + \frac{3}{2} v_1^2 - \frac{7}{4} (v_1 v_2) \right) \right\} \\
& + \frac{1}{c^4} \left\{ -\frac{G^3 m_1^3 m_2}{2 r_{12}^3} - \frac{19 G^3 m_1^2 m_2^2}{8 r_{12}^3} + \frac{5 m_1 v_1^6}{16} \right. \\
& \left. + \frac{G m_1 m_2}{r_{12}} \left( \frac{3}{8} (n_{12} v_1)^3 (n_{12} v_2) + \frac{3}{16} (n_{12} v_1)^2 (n_{12} v_2)^2 - \frac{9}{8} (n_{12} v_1) (n_{12} v_2) v_1^2 \right. \right. \\
& \left. - \frac{13}{8} (n_{12} v_2)^2 v_1^2 + \frac{21}{8} v_1^4 + \frac{13}{8} (n_{12} v_1)^2 (v_1 v_2) + \frac{3}{4} (n_{12} v_1) (n_{12} v_2) (v_1 v_2) \right. \\
& \left. - \frac{55}{8} v_1^2 (v_1 v_2) + \frac{17}{8} (v_1 v_2)^2 + \frac{31}{16} v_1^2 v_2^2 \right) \\
& + \frac{G^2 m_1^2 m_2}{r_{12}^2} \left( \frac{29}{4} (n_{12} v_1)^2 - \frac{13}{4} (n_{12} v_1) (n_{12} v_2) + \frac{1}{2} (n_{12} v_2)^2 - \frac{3}{2} v_1^2 + \frac{7}{4} v_2^2 \right) \\
& + \frac{1}{c^5} \left\{ \frac{4 G^2 m_1^2 m_2}{5 r_{12}^2} (n_{12} v_1) \left( v_{12}^2 - \frac{2 G (m_1 - m_2)}{r_{12}} \right) \right\} \\
& + \frac{1}{c^6} \left\{ \frac{35 m_1 v_1^8}{128} \right. \\
& + \frac{G m_1 m_2}{r_{12}} \left( -\frac{5}{16} (n_{12} v_1)^5 (n_{12} v_2) \right. \\
& \left. - \frac{5}{16} (n_{12} v_1)^4 (n_{12} v_2)^2 - \frac{5}{32} (n_{12} v_1)^3 (n_{12} v_2)^3 \right. \\
& + \frac{19}{16} (n_{12} v_1)^3 (n_{12} v_2) v_1^2 + \frac{15}{16} (n_{12} v_1)^2 (n_{12} v_2)^2 v_1^2 + \frac{3}{4} (n_{12} v_1) (n_{12} v_2)^3 v_1^2 \\
& + \frac{19}{16} (n_{12} v_2)^4 v_1^2 - \frac{21}{16} (n_{12} v_1) (n_{12} v_2) v_1^4 - 2 (n_{12} v_2)^2 v_1^4 \\
& \left. + \frac{55}{16} v_1^6 - \frac{19}{16} (n_{12} v_1)^4 (v_1 v_2) - (n_{12} v_1)^3 (n_{12} v_2) (v_1 v_2) \right)
\end{aligned}
$$
$$
\begin{align}
& -\frac{15}{32}(n_{12}v_1)^2(n_{12}v_2)^2(v_1v_2) + \frac{45}{16}(n_{12}v_1)^2v_1^2(v_1v_2) \nonumber \\
& \quad + \frac{5}{4}(n_{12}v_1)(n_{12}v_2)v_1^2(v_1v_2) + \frac{11}{4}(n_{12}v_2)^2v_1^2(v_1v_2) - \frac{139}{16}v_1^4(v_1v_2) \nonumber \\
& \quad -\frac{3}{4}(n_{12}v_1)^2(v_1v_2)^2 + \frac{5}{16}(n_{12}v_1)(n_{12}v_2)(v_1v_2)^2 + \frac{41}{8}v_1^2(v_1v_2)^2 + \frac{1}{16}(v_1v_2)^3 \nonumber \\
& \quad -\frac{45}{16}(n_{12}v_1)^2v_1^2v_2^2 - \frac{23}{32}(n_{12}v_1)(n_{12}v_2)v_1^2v_2^2 + \frac{79}{16}v_1^4v_2^2 - \frac{161}{32}v_1^2(v_1v_2)v_2^2 \nonumber \\
& \quad + \frac{G^2m_1^2m_2}{r_{12}^2} \left( -\frac{49}{8}(n_{12}v_1)^4 + \frac{75}{8}(n_{12}v_1)^3(n_{12}v_2) - \frac{187}{8}(n_{12}v_1)^2(n_{12}v_2)^2 \right. \nonumber \\
& \quad \left. + \frac{247}{24}(n_{12}v_1)(n_{12}v_2)^3 + \frac{49}{8}(n_{12}v_1)^2v_1^2 + \frac{81}{8}(n_{12}v_1)(n_{12}v_2)v_1^2 \right. \nonumber \\
& \quad \left. - \frac{21}{4}(n_{12}v_2)^2v_1^2 + \frac{11}{2}v_1^4 - \frac{15}{2}(n_{12}v_1)^2(v_1v_2) - \frac{3}{2}(n_{12}v_1)(n_{12}v_2)(v_1v_2) \right. \nonumber \\
& \quad \left. + \frac{21}{4}(n_{12}v_2)^2(v_1v_2) - 27v_1^2(v_1v_2) + \frac{55}{2}(v_1v_2)^2 + \frac{49}{4}(n_{12}v_1)^2v_2^2 \right. \nonumber \\
& \quad \left. - \frac{27}{2}(n_{12}v_1)(n_{12}v_2)v_2^2 + \frac{3}{4}(n_{12}v_2)^2v_2^2 + \frac{55}{4}v_1^2v_2^2 - 28(v_1v_2)v_2^2 + \frac{135}{16}v_2^4 \right) \nonumber \\
& \quad + \frac{3G^4m_1^4m_2}{8r_{12}^4} + \frac{G^4m_1^3m_2^2}{r_{12}^4} \left( \frac{9707}{420} - \frac{22}{3}\ln\left(\frac{r_{12}}{r_0'}\right) \right) \nonumber \\
& \quad + \frac{G^3m_1^2m_2^2}{r_{12}^3} \left( \frac{547}{12}(n_{12}v_1)^2 - \frac{3115}{48}(n_{12}v_1)(n_{12}v_2) - \frac{123}{64}(n_{12}v_1)(n_{12}v_12)\pi^2 \right. \nonumber \\
& \quad \left. - \frac{575}{18}v_1^2 + \frac{41}{64}\pi^2(v_1v_12) + \frac{4429}{144}(v_1v_2) \right) \nonumber \\
& \quad + \frac{G^3m_1^3m_2}{r_{12}^3} \left( -\frac{44627}{840}(n_{12}v_1)^2 + \frac{32027}{840}(n_{12}v_1)(n_{12}v_2) + \frac{3}{2}(n_{12}v_2)^2 \right. \nonumber \\
& \quad \left. + \frac{24187}{2520}v_1^2 - \frac{27967}{2520}(v_1v_2) + \frac{5}{4}v_2^2 + 22(n_{12}v_1)(n_{12}v_12)\ln\left(\frac{r_{12}}{r_0'}\right) \right. \nonumber \\
& \quad \left. - \frac{22}{3}(v_1v_12)\ln\left(\frac{r_{12}}{r_0'}\right) \right) \Biggr\} + 1 \leftrightarrow 2 + \mathcal{O}\left(\frac{1}{c^7}\right). \tag{340}
\end{align}
$$

To the terms given above, we must add the same terms but corresponding to the relabelling $1 \leftrightarrow 2$. The time derivative of this energy, as computed by means of the 3PN equations of motion themselves (i.e., by order-reducing all the accelerations), is purely equal to the 2.5PN effect, and agrees with the standard Einstein quadrupole formula (5):

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

where the Newtonian trace-free quadrupole moment reads
$Q_{ij} = m_1 (y_1^i y_1^j - \frac{1}{3} \delta^{ij} y_1^2) + 1 \leftrightarrow 2$. We refer to Iyer and Will (1993, 1995),
Gopakumar et al. (1997) for discussions of the energy balance equation to the next
3.5PN and 4.5PN orders. See also Eq. (177) for the energy balance equation at 4PN order for general fluid systems.

### 3.2.3 The 4PN Lagrangian and Hamiltonian

The conservative part of the equations of motion in harmonic coordinates (339) is derivable from a generalized Lagrangian, depending not only on the positions and velocities of the bodies, but also on their accelerations: $a_1 = dv_1/dt$ and $a_2 = dv_2/dt$, the accelerations occuring already from the 2PN order (Damour and Deruelle 1981a). A general result by Martin and Sanz (1979) states that N-body equations of motion cannot be derived from an ordinary Lagrangian beyond the 1PN level, when the gauge conditions preserve the manifest Lorentz invariance. One can always arrange for the dependence of the Lagrangian on accelerations to be linear, at the price of adding some “multi-zero” terms to the Lagrangian, which do not modify the equations of motion (Damour and Schäfer 1985).<fnref n="70" /> The Lagrangian of two particles in harmonic coordinates at 4PN order is really the Fokker Lagrangian $L \equiv L_F$ of Sect. 3.2.1. It has been derived using the traditional method by Bernard et al. (2016, 2017a, b, 2018), and by Galley et al. (2016), Foffa and Sturani (2012), Foffa et al. (2019) within the EFT approach. The long explicit results are presented in the form of the sequence of PN coefficients,

$$L = L_N + \frac{1}{c^2} L_{1PN} + \frac{1}{c^4} L_{2PN} + \frac{1}{c^6} L_{3PN} + \frac{1}{c^8} L_{4PN} + \mathcal{O}\left(\frac{1}{c^{10}}\right). \quad (342)$$

For convenience we do not write the irrelevant rest-mass contribution. Up to 3PN order the results read

$$L_N = \frac{m_1 v_1^2}{2} + \frac{G m_1 m_2}{2 r_{12}} + 1 \leftrightarrow 2, \quad (343a)$$

$$L_{1PN} = \frac{m_1 v_1^4}{8} - \frac{G^2 m_1^2 m_2}{2 r_{12}^2} + \frac{G m_1 m_2}{r_{12}} \left( -\frac{1}{4}(n_{12}v_1)(n_{12}v_2) + \frac{3}{2}v_1^2 - \frac{7}{4}(v_1 v_2) \right) + 1 \leftrightarrow 2, \quad (343b)$$

<footnote n="70">⁷⁰ The method is a particular application of a general algorithm to eliminate higher-derivative terms in a Lagrangian (Damour and Schäfer 1991).</footnote>
$$
\begin{align}
L_{2\text{PN}} ={}& \frac{m_1 v_1^6}{16} + \frac{G^3 m_1^3 m_2}{2 r_{12}^3} + \frac{19 G^3 m_1^2 m_2^2}{8 r_{12}^3} \nonumber \\
& + \frac{G^2 m_1^2 m_2}{r_{12}^2} \left( \frac{7}{2} (n_{12} v_1)^2 - \frac{7}{2} (n_{12} v_1) (n_{12} v_2) + \frac{1}{2} (n_{12} v_2)^2 \right. \nonumber \\
& \left. + \frac{1}{4} v_1^2 - \frac{7}{4} (v_1 v_2) + \frac{7}{4} v_2^2 \right) \nonumber \\
& + \frac{G m_1 m_2}{r_{12}} \left( \frac{3}{16} (n_{12} v_1)^2 (n_{12} v_2)^2 - \frac{7}{8} (n_{12} v_2)^2 v_1^2 + \frac{7}{8} v_1^4 \right) \tag{343c} \\
& + \frac{3}{4} (n_{12} v_1) (n_{12} v_2) (v_1 v_2) - 2 v_1^2 (v_1 v_2) + \frac{1}{8} (v_1 v_2)^2 + \frac{15}{16} v_1^2 v_2^2 \nonumber \\
& + G m_1 m_2 \left( - \frac{7}{4} (a_1 v_2) (n_{12} v_2) - \frac{1}{8} (n_{12} a_1) (n_{12} v_2)^2 + \frac{7}{8} (n_{12} a_1) v_2^2 \right) \nonumber \\
& + 1 \leftrightarrow 2, \nonumber
\end{align}
$$

$$
\begin{align*}
L_{3\text{PN}} ={}& \frac{5m_1 v_1^8}{128} \\
& + \frac{G^2 m_1^2 m_2}{r_{12}^2} \left( \frac{13}{18} (n_{12}v_1)^4 + \frac{83}{18} (n_{12}v_1)^3 (n_{12}v_2) - \frac{35}{6} (n_{12}v_1)^2 (n_{12}v_2)^2 \right. \\
& \left. - \frac{245}{24} (n_{12}v_1)^2 v_1^2 + \frac{179}{12} (n_{12}v_1) (n_{12}v_2) v_1^2 - \frac{235}{24} (n_{12}v_2)^2 v_1^2 + \frac{373}{48} v_1^4 \right. \\
& \left. + \frac{529}{24} (n_{12}v_1)^2 (v_1v_2) - \frac{97}{6} (n_{12}v_1) (n_{12}v_2) (v_1v_2) - \frac{719}{24} v_1^2 (v_1v_2) \right. \\
& \left. + \frac{463}{24} (v_1v_2)^2 - \frac{7}{24} (n_{12}v_1)^2 v_2^2 - \frac{1}{2} (n_{12}v_1) (n_{12}v_2) v_2^2 + \frac{1}{4} (n_{12}v_2)^2 v_2^2 \right. \\
& \left. + \frac{463}{48} v_1^2 v_2^2 - \frac{19}{2} (v_1v_2) v_2^2 + \frac{45}{16} v_2^4 \right) \\
& + Gm_1 m_2 \left( \frac{3}{8} (a_1v_2) (n_{12}v_1) (n_{12}v_2)^2 + \frac{5}{12} (a_1v_2) (n_{12}v_2)^3 \right. \\
& \left. + \frac{1}{8} (n_{12}a_1) (n_{12}v_1) (n_{12}v_2)^3 + \frac{1}{16} (n_{12}a_1) (n_{12}v_2)^4 + \frac{11}{4} (a_1v_1) (n_{12}v_2) v_1^2 \right. \\
& \left. - (a_1v_2) (n_{12}v_2) v_1^2 - 2(a_1v_1) (n_{12}v_2) (v_1v_2) + \frac{1}{4} (a_1v_2) (n_{12}v_2) (v_1v_2) \right. \\
& \left. + \frac{3}{8} (n_{12}a_1) (n_{12}v_2)^2 (v_1v_2) - \frac{5}{8} (n_{12}a_1) (n_{12}v_1)^2 v_2^2 + \frac{15}{8} (a_1v_1) (n_{12}v_2) v_2^2 \right. \\
& \left. - \frac{15}{8} (a_1v_2) (n_{12}v_2) v_2^2 - \frac{1}{2} (n_{12}a_1) (n_{12}v_1) (n_{12}v_2) v_2^2 \right)
\end{align*}
$$
$$
\begin{aligned}
& -\frac{5}{16}(n_{12}a_1)(n_{12}v_2)^2 v_2^2) \\
& +\frac{G^2m_1^2m_2}{r_{12}}\left(-\frac{235}{24}(a_2v_1)(n_{12}v_1)-\frac{29}{24}(n_{12}a_2)(n_{12}v_1)^2\right. \\
& \qquad \left.-\frac{235}{24}(a_1v_2)(n_{12}v_2)-\frac{17}{6}(n_{12}a_1)(n_{12}v_2)^2+\frac{185}{16}(n_{12}a_1)v_1^2\right. \\
& \qquad \left.-\frac{235}{48}(n_{12}a_2)v_1^2-\frac{185}{8}(n_{12}a_1)(v_1v_2)+\frac{20}{3}(n_{12}a_1)v_2^2\right) \\
& +\frac{Gm_1m_2}{r_{12}}\left(-\frac{5}{32}(n_{12}v_1)^3(n_{12}v_2)^3+\frac{1}{8}(n_{12}v_1)(n_{12}v_2)^3v_1^2+\frac{5}{8}(n_{12}v_2)^4v_1^2\right. \\
& \qquad \left.-\frac{11}{16}(n_{12}v_1)(n_{12}v_2)v_1^4+\frac{1}{4}(n_{12}v_2)^2v_1^4+\frac{11}{16}v_1^6\right. \\
& \qquad \left.-\frac{15}{32}(n_{12}v_1)^2(n_{12}v_2)^2(v_1v_2)+(n_{12}v_1)(n_{12}v_2)v_1^2(v_1v_2)\right. \\
& \qquad \left.+\frac{3}{8}(n_{12}v_2)^2v_1^2(v_1v_2)-\frac{13}{16}v_1^4(v_1v_2)+\frac{5}{16}(n_{12}v_1)(n_{12}v_2)(v_1v_2)^2\right. \\
& \qquad \left.+\frac{1}{16}(v_1v_2)^3-\frac{5}{8}(n_{12}v_1)^2v_1^2v_2^2-\frac{23}{32}(n_{12}v_1)(n_{12}v_2)v_1^2v_2^2+\frac{1}{16}v_1^4v_2^2\right. \\
& \qquad \left.-\frac{1}{32}v_1^2(v_1v_2)v_2^2\right) -\frac{3G^4m_1^4m_2}{8r_{12}^4}+\frac{G^4m_1^3m_2^2}{r_{12}^4}\left(-\frac{9707}{420}+\frac{22}{3}\ln\left(\frac{r_{12}}{r_0'}\right)\right) \\
& \qquad +\frac{G^3m_1^2m_2^2}{r_{12}^3}\left(\frac{383}{24}(n_{12}v_1)^2-\frac{889}{48}(n_{12}v_1)(n_{12}v_2)\right. \\
& \qquad \left.-\frac{123}{64}(n_{12}v_1)(n_{12}v_2)\pi^2-\frac{305}{72}v_1^2+\frac{41}{64}\pi^2(v_1v_12)+\frac{439}{144}(v_1v_2)\right) \\
& \qquad +\frac{G^3m_1^3m_2}{r_{12}^3}\left(-\frac{8243}{210}(n_{12}v_1)^2+\frac{15541}{420}(n_{12}v_1)(n_{12}v_2)+\frac{3}{2}(n_{12}v_2)^2\right. \\
& \qquad \left. +\frac{15611}{1260}v_1^2-\frac{17501}{1260}(v_1v_2)+\frac{5}{4}v_2^2+22(n_{12}v_1)(n_{12}v_2)\ln\left(\frac{r_{12}}{r_0'}\right)\right. \\
& \qquad \left.-\frac{22}{3}(v_1v_2)\ln\left(\frac{r_{12}}{r_0'}\right)\right)+1\leftrightarrow2.
\end{aligned}
\quad (343d)
$$

where the logarithms in harmonic coordinates contain the gauge constant $r'_0$, and the accelerations start appearing at 2PN order [recall the notation, e.g. $(a_1v_2) = a_1 \cdot v_2$]. We refer to de Andrade et al. (2001) for the expressions of the ten conserved quantities corresponding to the integrals of energy [given in Eq. (340)], linear and angular momenta, and center-of-mass position at 3PN order. Note that while it is strictly forbidden to replace the accelerations by the equations of motion in the Lagrangian, this can and *should* be done in the final expressions of the conserved integrals derived from the Lagrangian.

Now we exhibit a transformation of the particles' dynamical variables—or *contact* transformation—which transforms the harmonic-coordinates 3PN Lagrangian (343) into a new Lagrangian, valid in some ADM or ADM-like coordinate system, and such that the associated Hamiltonian coincides with the 3PN Hamiltonian obtained by Jaranowski and Schäfer (1998, 1999). In ADM coordinates the Lagrangian will be ordinary, depending only on the positions and velocities of the bodies. Let this
contact transformation be $\mathbf{Y}_1(t) = \mathbf{y}_1(t) + \delta\mathbf{y}_1(t)$ and $1 \leftrightarrow 2$, where $\mathbf{Y}_1$ and $\mathbf{y}_1$ denote the trajectories in ADM and harmonic coordinates, respectively. For this transformation to be able to remove the accelerations in the Lagrangian (343) at 3PN order, we determine it to be necessarily of the form (de Andrade et al. 2001)

$$ \delta\mathbf{y}_1 = \frac{1}{m_1} \left[ \frac{\partial L}{\partial \boldsymbol{a}_1} + \frac{\partial F}{\partial \boldsymbol{v}_1} + \frac{1}{c^6} \mathbf{X}_1 \right] + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (344) $$

and *idem* for $1 \leftrightarrow 2$, where $F$ is a freely adjustable function of the positions and velocities, made of 2PN and 3PN terms, and where $\mathbf{X}_1$ represents a special correction term, that is purely of order 3PN. The point is that once the function $F$ is specified there is a unique determination of the correction term $\mathbf{X}_1$ for the contact transformation to work (de Andrade et al. 2001). Thus, the freedom we have is entirely encoded into the function $F$, and the work then consists in showing that there exists a unique choice of $F$ for which the Lagrangian (343) is physically equivalent, via the contact transformation (344), to the ADM Hamiltonian of Jaranowski and Schäfer (1998, 1999). An interesting point is that not only the transformation must remove all the accelerations, but it should also cancel out all the logarithms $\ln(r_{12}/r'_0)$, because there are no logarithms in ADM coordinates. As a result, $F$ involves the logarithmic terms$^{71}$

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

together with other non-logarithmic terms (indicated by dots) that are entirely specified by the isometry of the harmonic and ADM descriptions of the motion. The function $F$ and correction term $\mathbf{X}_1$ are given by Eqs. (4.9)–(4.10) of de Andrade et al. (2001). For these choices the ADM Lagrangian becomes

$$ L^{\text{ADM}} = L + \frac{\delta L}{\delta y_1^j} \delta y_1^j + \frac{\delta L}{\delta y_2^j} \delta y_2^j + \frac{dF}{dt} + \mathcal{O}\left(\frac{1}{c^8}\right), \quad (346) $$

and the ADM Hamiltonian is deduced by an ordinary Legendre transformation. The procedure can easily be generalized to 4PN order.

Up to 3PN order all terms were “instantaneous”, as they depend on the system at current time only, see Eqs. (343). However at 4PN order a new feature arises: the non-locality (in time) due to infinite-range temporal correlations of gravitational wave tails (Blanchet and Damour 1988) and computed in the Fokker Lagrangian in Sect. 3.2.1. Accordingly we split the 4PN coefficient as

$$ \eta^\alpha = -\frac{22 G^2 m_1 m_2}{3 c^6} \delta^\alpha \left( \frac{Gm_2}{r_1} + \frac{Gm_1}{r_2} \right) \ln \left( \frac{r_{12}}{r'_0} \right). $$

<footnote>$^{71}$ Hence the logarithms at 3PN order can be removed by the coordinate transformation $x'^\alpha = x^\alpha + \eta^\alpha(x)$ associated with the vector (with $r_a = |\mathbf{x} - \mathbf{y}_a(t)|$ and $\delta^\alpha = \eta^{\alpha\mu}\delta_{\mu\alpha}$):</footnote>
$$L_{4PN} = L_{4PN}^{\text{inst}} + L_{4PN}^{\text{tail}}, \qquad (347)$$

where the tail term reads explicitly

$$L_{4PN}^{\text{tail}} = \frac{G^2 M}{5c^8} M_{ij}^{(3)}(t) \int_{0}^{+\infty} d\tau \ln\left(\frac{c\tau}{2r_{12}}\right) \left[M_{ij}^{(4)}(t-\tau) - M_{ij}^{(4)}(t+\tau)\right]. \quad (348)$$

There is, though, an important difference with Eqs. (336)–(337), is that the Hada-mard cut-off scale $c\tau_0$ has been replaced by (two times) the physical separation between the two bodies, $r_{12} = |\mathbf{y}_1 - \mathbf{y}_2|$. Therefore when varying the tail term (348) [remind that this is a functional variation] one has also to take into account the variation of the “constant” $r_{12}(t)$. See Eq. (4.3) of Bernard et al. (2018) for the result of this variation. Note also that this variation concerns the conservative part of the acceleration; in addition one has to take into account the dissipative part. See Eq. (360) for the complete result in the center-of-mass frame, including both con- servative and dissipative contributions.

As for the instantaneous part of the 4PN term, we split it into non-linear contributions corresponding to increasing powers of G:

$$L_{4PN}^{\text{inst}} = L_{4PN}^{(0)} + G L_{4PN}^{(1)} + G^2 L_{4PN}^{(2)} + G^3 L_{4PN}^{(3)} + G^4 L_{4PN}^{(4)} + G^5 L_{4PN}^{(5)}, \quad (349)$$

where each of the terms reads explicitly

$$L_{4PN}^{(0)} = \frac{7}{256} m_1 v_1^{10} + 1 \leftrightarrow 2, \qquad (350a)$$

$$\begin{aligned} L_{4PN}^{(1)} = & m_1 m_2 \left( \frac{13}{64} (a_2 v_1) (n_{12} v_1)^5 + \frac{5}{128} (a_2 n_{12}) (n_{12} v_1)^6 - \frac{13}{64} (n_{12} v_1)^5 (a_2 v_2) \right. \\ & + \frac{11}{64} (a_2 v_1) (n_{12} v_1)^4 (n_{12} v_2) + \frac{5}{64} (a_2 n_{12}) (n_{12} v_1)^5 (n_{12} v_2) \\ & + \frac{5}{32} (a_2 v_1) (n_{12} v_1)^3 (n_{12} v_2)^2 - \frac{5}{32} (a_1 n_{12}) (n_{12} v_1)^3 (n_{12} v_2)^3 \\ & - \frac{1}{16} (a_2 v_1) (n_{12} v_1)^3 (v_1 v_2) + \frac{11}{64} (a_2 n_{12}) (n_{12} v_1)^4 (v_1 v_2) \\ & + \frac{5}{16} (a_2 n_{12}) (n_{12} v_1)^3 (n_{12} v_2) (v_1 v_2) + \frac{5}{16} (n_{12} v_1)^2 (a_1 v_2) (n_{12} v_2) (v_1 v_2) \\ & - \frac{3}{16} (a_2 v_1) (n_{12} v_1) (v_1 v_2)^2 + \frac{1}{16} (a_1 v_1) (n_{12} v_2) (v_1 v_2)^2 \\ & + \frac{5}{16} (a_1 n_{12}) (n_{12} v_1) (n_{12} v_2) (v_1 v_2)^2 - \frac{77}{96} (a_2 v_1) (n_{12} v_1)^3 v_1^2 \\ & - \frac{27}{128} (a_2 n_{12}) (n_{12} v_1)^4 v_1^2 + \frac{77}{96} (n_{12} v_1)^3 (a_2 v_2) v_1^2 \end{aligned}$$
$$
\begin{aligned}
& -\frac{13}{32}(a_2 v_1)(n_{12} v_1)^2 (n_{12} v_2)^2 v_1^2 - \frac{11}{32}(a_2 n_{12})(n_{12} v_1)^3 (n_{12} v_2)^2 v_1^2 \\
& -\frac{7}{32}(a_2 v_1)(n_{12} v_1)(n_{12} v_2)^2 v_1^2 - \frac{27}{64}(a_2 n_{12})(n_{12} v_1)^2 (n_{12} v_2)^2 v_1^2 \\
& -\frac{19}{32}(a_1 v_1)(n_{12} v_2)^3 v_1^2 - \frac{3}{16}(a_2 v_1)(n_{12} v_1)(v_1 v_2) v_1^2 \\
& -\frac{13}{32}(a_2 n_{12})(n_{12} v_1)^2 (v_1 v_2) v_1^2 + \frac{33}{16}(n_{12} v_1)(a_2 v_2)(v_1 v_2) v_1^2 \\
& -\frac{7}{16}(a_2 n_{12})(n_{12} v_1)(n_{12} v_2)(v_1 v_2) v_1^2 + \frac{1}{16}(a_1 v_2)(n_{12} v_2)(v_1 v_2) v_1^2 \\
& +\frac{123}{64}(a_2 v_1)(n_{12} v_1) v_1^4 + \frac{53}{128}(a_2 n_{12})(n_{12} v_1)^2 v_1^4 - \frac{123}{64}(n_{12} v_1)(a_2 v_2) v_1^4 \\
& +\frac{49}{64}(a_2 v_1)(n_{12} v_2) v_1^4 + \frac{31}{64}(a_2 n_{12})(n_{12} v_1)(n_{12} v_2) v_1^4 \\
& +\frac{49}{64}(a_2 n_{12})(v_1 v_2) v_1^4 - \frac{75}{128}(a_2 n_{12}) v_1^6 + \frac{17}{96}(n_{12} v_1)^3 (a_1 v_2) v_2^2 \\
& -\frac{23}{32}(a_1 v_1)(n_{12} v_1)^2 (n_{12} v_2) v_2^2 + \frac{15}{32}(a_1 n_{12})(n_{12} v_1)^3 (n_{12} v_2) v_2^2 \\
& +\frac{17}{32}(a_1 n_{12})(n_{12} v_1)^2 (v_1 v_2) v_2^2 + \frac{33}{32}(a_2 v_1)(n_{12} v_1) v_1^2 v_2^2 \\
& +\frac{93}{32}(a_1 v_1)(n_{12} v_2) v_1^2 v_2^2 - \frac{23}{32}(a_1 n_{12})(n_{12} v_1)(n_{12} v_2) v_1^2 v_2^2) \\
& +\frac{m_1 m_2}{r_{12}} \left( \frac{35}{128}(n_{12} v_1)^5 (n_{12} v_2)^3 - \frac{35}{256}(n_{12} v_1)^4 (n_{12} v_2)^4 \right. \\
& \qquad \left. -\frac{15}{128}(n_{12} v_1)^4 (n_{12} v_2)^2 (v_1 v_2) + \frac{15}{32}(n_{12} v_1)^3 (n_{12} v_2)^3 (v_1 v_2) \right. \\
& \qquad \left. -\frac{15}{32}(n_{12} v_1)^3 (n_{12} v_2)(v_1 v_2)^2 + \frac{5}{32}(n_{12} v_1)^2 (v_1 v_2)^3 - \frac{3}{16}(n_{12} v_1)(n_{12} v_2)(v_1 v_2)^3 \right. \\
& \qquad \left. +\frac{1}{32}(v_1 v_2)^4 - \frac{5}{32}(n_{12} v_1)^3 (n_{12} v_2)^3 v_1^2 - \frac{5}{16}(n_{12} v_1)(n_{12} v_2)^3 (v_1 v_2) v_1^2 \right. \\
& \qquad \left. +\frac{9}{32}(n_{12} v_1)(n_{12} v_2)(v_1 v_2)^2 v_1^2 - \frac{15}{32}(n_{12} v_2)^2 (v_1 v_2)^2 v_1^2 + \frac{1}{32}(v_1 v_2)^3 v_1^2 \right. \\
& \qquad \left. +\frac{57}{128}(n_{12} v_1)(n_{12} v_2)^3 v_1^4 - \frac{15}{128}(n_{12} v_2)^4 v_1^4 + \frac{39}{128}(n_{12} v_2)^2 (v_1 v_2) v_1^4 + \frac{3}{4}(v_1 v_2)^2 v_1^4 \right. \\
& \qquad \left. -\frac{11}{32}(n_{12} v_2)^2 v_1^6 - \frac{5}{4}(v_1 v_2) v_1^6 + \frac{75}{128} v_1^8 - \frac{75}{128}(n_{12} v_1)^5 (n_{12} v_2) v_2^2 \right. \\
& \qquad \left. -\frac{53}{128}(n_{12} v_1)^4 (v_1 v_2) v_2^2 + \frac{99}{64}(n_{12} v_1)^3 (n_{12} v_2) v_1^2 v_2^2 - \frac{21}{64}(n_{12} v_1)^2 (n_{12} v_2)^2 v_1^2 v_2^2 \right. \\
& \qquad \left. +\frac{11}{64}(n_{12} v_1)^2 (v_1 v_2) v_1^2 v_2^2 + \frac{35}{32}(n_{12} v_1)(n_{12} v_2)(v_1 v_2) v_1^2 v_2^2 - \frac{1}{32}(v_1 v_2)^2 v_1^2 v_2^2 \right. \\
& \qquad \left. -\frac{185}{128}(n_{12} v_1)(n_{12} v_2) v_1^4 v_2^2 + \frac{23}{64}(n_{12} v_2)^2 v_1^4 v_2^2 - \frac{99}{128}(v_1 v_2) v_1^4 v_2^2 + \frac{5}{8} v_1^6 v_2^2 \right. \\
& \qquad \left. +\frac{3}{256} v_1^4 v_2^4 \right) + 1 \leftrightarrow 2,
\end{aligned}
\quad (350b)
$$
$$
\begin{aligned}
L_{\text{4PN}}^{(2)} = & \frac{m_1^2 m_2}{r_{12}} \left( \frac{2099}{288} (a_1 v_1) (n_{12} v_1)^3 + \frac{3341}{480} (a_2 v_1) (n_{12} v_1)^3 + \frac{59}{180} (n_{12} v_1)^3 (a_2 v_2) \right. \\
& \left. + \frac{2197}{240} (a_1 n_{12}) (n_{12} v_1)^3 (n_{12} v_2) + \frac{6661}{720} (a_2 n_{12}) (n_{12} v_1)^3 (n_{12} v_2) \right. \\
& \left. + \frac{10223}{480} (n_{12} v_1)^2 (a_1 v_2) (n_{12} v_2) + \frac{3059}{96} (a_1 v_1) (n_{12} v_1) (n_{12} v_2)^2 \right. \\
& \left. - \frac{3781}{160} (n_{12} v_1) (a_1 v_2) (n_{12} v_2)^2 + \frac{1337}{240} (a_1 n_{12}) (n_{12} v_1) (n_{12} v_2)^3 \right. \\
& \left. + \frac{4621}{480} (a_1 v_2) (n_{12} v_2)^3 \right. \\
& \left. + \frac{1133}{960} (a_1 n_{12}) (n_{12} v_2)^4 + \frac{3613}{48} (a_1 v_1) (n_{12} v_1) (v_1 v_2) \right. \\
& \left. + \frac{4529}{240} (a_2 v_1) (n_{12} v_1) (v_1 v_2) \right. \\
& \left. + \frac{2099}{96} (a_1 n_{12}) (n_{12} v_1)^2 (v_1 v_2) - \frac{23}{60} (a_1 n_{12}) (n_{12} v_1) (n_{12} v_2) (v_1 v_2) \right. \\
& \left. + \frac{6499}{240} (a_1 v_2) (n_{12} v_2) (v_1 v_2) + \frac{247}{30} (a_1 n_{12}) (n_{12} v_2)^2 (v_1 v_2) \right. \\
& \left. + \frac{2503}{96} (a_1 n_{12}) (v_1 v_2)^2 \right. \\
& \left. + \frac{7193}{240} (a_2 v_1) (n_{12} v_1) v_1^2 + \frac{3021}{320} (a_2 n_{12}) (n_{12} v_1)^2 v_1^2 + \frac{4723}{96} (n_{12} v_1) (a_1 v_2) v_1^2 \right. \\
& \left. + \frac{3679}{480} (n_{12} v_1) (a_2 v_2) v_1^2 + \frac{13549}{480} (a_1 v_1) (n_{12} v_2) v_1^2 + \frac{8849}{480} (a_2 v_1) (n_{12} v_2) v_1^2 \right. \\
& \left. + \frac{2063}{96} (a_1 n_{12}) (n_{12} v_1) (n_{12} v_2) v_1^2 + \frac{166}{15} (a_2 n_{12}) (n_{12} v_1) (n_{12} v_2) v_1^2 \right. \\
& \left. + \frac{4621}{320} (a_2 n_{12}) (n_{12} v_2)^2 v_1^2 + \frac{8849}{480} (a_2 n_{12}) (v_1 v_2) v_1^2 + \frac{6943}{384} (a_1 n_{12}) v_1^4 \right. \\
& \left. + \frac{2293}{160} (n_{12} v_1) (a_1 v_2) v_2^2 + \frac{3733}{160} (a_1 v_1) (n_{12} v_2) v_2^2 \right. \\
& \left. + \frac{7}{5} (a_1 n_{12}) (n_{12} v_1) (n_{12} v_2) v_2^2 \right. \\
& \left. - \frac{3733}{160} (a_1 v_2) (n_{12} v_2) v_2^2 - \frac{139}{20} (a_1 n_{12}) (n_{12} v_2)^2 v_2^2 - \frac{5593}{240} (a_1 n_{12}) (v_1 v_2) v_2^2 \right. \\
& \left. + \frac{3613}{192} (a_1 n_{12}) v_1^2 v_2^2 + \frac{2931}{320} (a_1 n_{12}) v_2^4 \right) \\
& + \frac{m_1^2 m_2}{r_{12}^2} \left( - \frac{4027}{800} (n_{12} v_1)^6 + \frac{3227}{800} (n_{12} v_1)^5 (n_{12} v_2) \right. \\
& \left. - \frac{6301}{240} (n_{12} v_1)^4 (n_{12} v_2)^2 \right)
\end{aligned}
$$
$$
\begin{align}
&+\frac{6661}{240}(n_{12}v_1)^3(n_{12}v_2)^3 - \frac{2221}{64}(n_{12}v_1)^4(v_1v_2) + \frac{25267}{720}(n_{12}v_1)^3(n_{12}v_2)(v_1v_2) \nonumber \\
&-\frac{6661}{480}(n_{12}v_1)^2(n_{12}v_2)^2(v_1v_2) - \frac{23401}{480}(n_{12}v_1)^2(v_1v_2)^2 \nonumber \\
&+\frac{4529}{240}(n_{12}v_1)(n_{12}v_2)(v_1v_2)^2 - \frac{8369}{480}(v_1v_2)^3 + \frac{17393}{2880}(n_{12}v_1)^4v_1^2 \nonumber \\
&-\frac{26237}{720}(n_{12}v_1)^3(n_{12}v_2)v_1^2 + \frac{4561}{160}(n_{12}v_1)^2(n_{12}v_2)^2v_1^2 + \frac{691}{240}(n_{12}v_1)(n_{12}v_2)^3v_1^2 \nonumber \\
&+\frac{4621}{240}(n_{12}v_2)^4v_1^2 - \frac{14987}{960}(n_{12}v_1)^2(v_1v_2)v_1^2 + \frac{2601}{160}(n_{12}v_1)(n_{12}v_2)(v_1v_2)v_1^2 \nonumber \\
&+\frac{14649}{320}(n_{12}v_2)^2(v_1v_2)v_1^2 + \frac{97}{5}(v_1v_2)^2v_1^2 - \frac{4879}{320}(n_{12}v_1)^2v_1^4 \nonumber \\
&+\frac{5399}{192}(n_{12}v_1)(n_{12}v_2)v_1^4 + \frac{83}{15}(n_{12}v_2)^2v_1^4 + \frac{749}{128}(v_1v_2)v_1^4 + \frac{20389}{1920}v_1^6 \nonumber \\
&+\frac{107}{180}(n_{12}v_1)^4v_2^2 - \frac{1823}{240}(n_{12}v_1)^3(n_{12}v_2)v_2^2 - \frac{1}{2}(n_{12}v_1)^2(n_{12}v_2)^2v_2^2 \nonumber \\
&+\frac{1021}{120}(n_{12}v_1)^2(v_1v_2)v_2^2 + \frac{1}{2}(n_{12}v_1)(n_{12}v_2)(v_1v_2)v_2^2 + \frac{67}{4}(v_1v_2)^2v_2^2 \nonumber \\
&-\frac{21709}{960}(n_{12}v_1)^2v_1^2v_2^2 - \frac{1873}{480}(n_{12}v_1)(n_{12}v_2)v_1^2v_2^2 - \frac{4621}{320}(n_{12}v_2)^2v_1^2v_2^2 \nonumber \\
&-\frac{42017}{960}(v_1v_2)v_1^2v_2^2 + \frac{11119}{960}v_1^4v_2^2 - \frac{21}{8}(n_{12}v_1)^2v_2^4 - \frac{3}{8}(n_{12}v_1)(n_{12}v_2)v_2^4 \nonumber \\
&+\frac{3}{16}(n_{12}v_2)^2v_2^4 - \frac{105}{8}(v_1v_2)v_2^4 + \frac{105}{16}v_1^2v_2^4 + \frac{115}{32}v_2^6) + 1 \leftrightarrow 2, \tag{350c}
\end{align}
$$

$$
\begin{align*}
L_{4PN}^{(3)} &= \frac{m_1^2 m_2^2}{r_{12}^2} \left( -\frac{1099}{144} (a_2 v_1) (n_{12} v_1) \right. \\
&\quad + \frac{41}{64} \pi^2 (a_2 v_1) (n_{12} v_1) + \frac{2005}{96} (a_2 n_{12}) (n_{12} v_1)^2 \\
&\quad - \frac{123}{128} \pi^2 (a_2 n_{12}) (n_{12} v_1)^2 + \frac{225233}{1800} (a_1 v_1) (n_{12} v_2) - \frac{43}{64} \pi^2 (a_1 v_1) (n_{12} v_2) \\
&\quad - \frac{477941}{3600} (a_1 n_{12}) (v_1 v_2) + \frac{21}{16} \pi^2 (a_1 n_{12}) (v_1 v_2) + \frac{477941}{7200} (a_1 n_{12}) v_1^2 \\
&\quad \left. - \frac{21}{32} \pi^2 (a_1 n_{12}) v_1^2 \right) \\
&\quad + \frac{m_1^2 m_2^2}{r_{12}^3} \left( -\frac{173617}{2880} (n_{12} v_1)^4 - \frac{2155}{1024} \pi^2 (n_{12} v_1)^4 \right. \\
&\quad \left. + \frac{173587}{720} (n_{12} v_1)^3 (n_{12} v_2) + \frac{2155}{256} \pi^2 (n_{12} v_1)^3 (n_{12} v_2) \right)
\end{align*}
$$
$$
\begin{aligned}
& -\frac{85871}{480} (n_{12}v_1)^2 (n_{12}v_2)^2 \\
& \quad -\frac{6465}{1024} \pi^2 (n_{12}v_1)^2 (n_{12}v_2)^2 + \frac{5651}{300} (n_{12}v_1)^2 (v_1v_2) - \frac{939}{256} \pi^2 (n_{12}v_1)^2 (v_1v_2) \\
& \quad -\frac{6851}{300} (n_{12}v_1)(n_{12}v_2)(v_1v_2) + \frac{939}{256} \pi^2 (n_{12}v_1)(n_{12}v_2)(v_1v_2) + \frac{49139}{720} (v_1v_2)^2 \\
& \quad -\frac{195}{512} \pi^2 (v_1v_2)^2 - \frac{3677}{1200} (n_{12}v_1)^2 v_1^2 + \frac{447}{512} \pi^2 (n_{12}v_1)^2 v_1^2 \\
& \quad -\frac{222679}{1200} (n_{12}v_1)(n_{12}v_2)v_1^2 \\
& \quad -\frac{189}{256} \pi^2 (n_{12}v_1)(n_{12}v_2)v_1^2 + \frac{153079}{800} (n_{12}v_2)^2 v_1^2 - \frac{69}{512} \pi^2 (n_{12}v_2)^2 v_1^2 \\
& \quad -\frac{61733}{900} (v_1v_2)v_1^2 - \frac{55}{256} \pi^2 (v_1v_2)v_1^2 + \frac{10337}{320} v_1^4 + \frac{133}{1024} \pi^2 v_1^4 - \frac{116123}{3600} v_1^2 v_2^2 \\
& \quad +\frac{477}{1024} \pi^2 v_1^2 v_2^2 \quad +\frac{m_1^3 m_2}{r_{12}^2} \left(\frac{44023}{720} (a_1v_1)(n_{12}v_1) + \frac{562}{9} (a_2v_1)(n_{12}v_1)\right) \\
& \quad +44 \ln\left(\frac{r_{12}}{r_0'}\right) (n_{12}v_1)(a_1v_2) - 44 \ln\left(\frac{r_{12}}{r_0'}\right) (n_{12}v_1)(a_2v_2) \\
& \quad +\frac{110}{3} \ln\left(\frac{r_{12}}{r_0'}\right) (a_1v_1)(n_{12}v_2) \\
& \quad +\frac{6397}{75} (a_1n_{12})(n_{12}v_1)(n_{12}v_2) + \frac{198097}{4200} (a_2n_{12})(n_{12}v_1)(n_{12}v_2) \\
& \quad +22 \ln\left(\frac{r_{12}}{r_0'}\right) (a_2n_{12})(n_{12}v_1)(n_{12}v_2) + \frac{14377}{280} (a_1v_2)(n_{12}v_2) \\
& \quad -\frac{110}{3} \ln\left(\frac{r_{12}}{r_0'}\right) (a_1v_2)(n_{12}v_2) + \frac{44023}{720} (a_1n_{12})(v_1v_2) \\
& \quad -44 \ln\left(\frac{r_{12}}{r_0'}\right) (a_1n_{12})(v_1v_2) \\
& \quad +44 \ln\left(\frac{r_{12}}{r_0'}\right) (a_1n_{12})v_1^2 + \frac{14377}{560} (a_2n_{12})v_1^2 + \frac{937}{1440} (a_1n_{12})v_2^2) \\
& \quad +\frac{m_1^3 m_2}{r_{12}^3} \left(-\frac{30313}{360} (n_{12}v_1)^4 + \frac{64001}{720} (n_{12}v_1)^3 (n_{12}v_2)\right) \\
& \quad -\frac{185917}{1680} (n_{12}v_1)^2 (n_{12}v_2)^2 \\
& \quad -55 \ln\left(\frac{r_{12}}{r_0'}\right) (n_{12}v_1)^2 (n_{12}v_2)^2 + \frac{179617}{1680} (n_{12}v_1)(n_{12}v_2)^3 \\
& \quad +55 \ln\left(\frac{r_{12}}{r_0'}\right) (n_{12}v_1)(n_{12}v_2)^3 \\
& \quad -\frac{94667}{400} (n_{12}v_1)^2 (v_1v_2) + \frac{338099}{1400} (n_{12}v_1)(n_{12}v_2)(v_1v_2)
\end{aligned}
$$
$$
\begin{align*}
&+ 22 \ln \left( \frac{r_{12}}{r'_{0}} \right) (n_{12}v_{1})(n_{12}v_{2})(v_{1}v_{2}) - \frac{214897}{8400} (n_{12}v_{2})^{2}(v_{1}v_{2}) \\
&- 11 \ln \left( \frac{r_{12}}{r'_{0}} \right) (n_{12}v_{2})^{2}(v_{1}v_{2}) - \frac{737}{18}(v_{1}v_{2})^{2} + \frac{903589}{16800}(n_{12}v_{1})^{2}v_{1}^{2} \\
&- 55 \ln \left( \frac{r_{12}}{r'_{0}} \right) (n_{12}v_{1})^{2}v_{1}^{2} - \frac{96287}{1120}(n_{12}v_{1})(n_{12}v_{2})v_{1}^{2} \\
&+ 55 \ln \left( \frac{r_{12}}{r'_{0}} \right) (n_{12}v_{1})(n_{12}v_{2})v_{1}^{2} \\
&+ \frac{426731}{4200}(n_{12}v_{2})^{2}v_{1}^{2} + 11 \ln \left( \frac{r_{12}}{r'_{0}} \right) (n_{12}v_{2})^{2}v_{1}^{2} + \frac{202687}{3360}(v_{1}v_{2})v_{1}^{2} \tag{350d} \\
&- \frac{55}{3} \ln \left( \frac{r_{12}}{r'_{0}} \right) (v_{1}v_{2})v_{1}^{2} + \frac{22769}{2016}v_{1}^{4} + \frac{55}{3} \ln \left( \frac{r_{12}}{r'_{0}} \right) v_{1}^{4} - \frac{177}{8}(n_{12}v_{1})^{2}v_{2}^{2} \\
&+ 66 \ln \left( \frac{r_{12}}{r'_{0}} \right) (n_{12}v_{1})^{2}v_{2}^{2} - \frac{120397}{4200}(n_{12}v_{1})(n_{12}v_{2})v_{2}^{2} \\
&- 88 \ln \left( \frac{r_{12}}{r'_{0}} \right) (n_{12}v_{1})(n_{12}v_{2})v_{2}^{2} \\
&+ \frac{7}{4}(n_{12}v_{2})^{2}v_{2}^{2} - \frac{43}{2}(v_{1}v_{2})v_{2}^{2} + 22 \ln \left( \frac{r_{12}}{r'_{0}} \right) (v_{1}v_{2})v_{2}^{2} \\
&- \frac{8357}{560}v_{1}^{2}v_{2}^{2} - 22 \ln \left( \frac{r_{12}}{r'_{0}} \right) v_{1}^{2}v_{2}^{2} + \frac{91}{16}v_{2}^{4} \bigg|_{0}^{2} + 1 \leftrightarrow 2,
\end{align*}
$$

$$
\begin{align*}
L_{4\text{PN}}^{(4)} &= \frac{m_1^4 m_2}{r_{12}^4} \left( \frac{282629}{900} (n_{12}v_1)^2 - \frac{880}{3} \ln\left(\frac{r_{12}}{r'_0}\right) (n_{12}v_1)^2 - \frac{283979}{900} (n_{12}v_1)(n_{12}v_2) \right. \\
&\qquad \left. + \frac{880}{3} \ln\left(\frac{r_{12}}{r'_0}\right) (n_{12}v_1)(n_{12}v_2) + \frac{9}{4} (n_{12}v_2)^2 \right. \\
&\qquad \left. + \frac{208529}{3600} (v_1 v_2) - \frac{220}{3} \ln\left(\frac{r_{12}}{r'_0}\right) (v_1 v_2) \right. \\
&\qquad \left. - \frac{211229}{3600} v_1^2 + \frac{220}{3} \ln\left(\frac{r_{12}}{r'_0}\right) v_1^2 + \frac{15}{16} v_2^2 \right) \\
&\qquad + \frac{m_1^3 m_2^2}{r_{12}^4} \left( -\frac{1268557}{50400} (n_{12}v_1)^2 + \frac{659}{96} \pi^2 (n_{12}v_1)^2 - \frac{286}{3} \ln\left(\frac{r_{12}}{r'_0}\right) (n_{12}v_1)^2 \right. \\
&\qquad \left. + \frac{11530469}{25200} (n_{12}v_1)(n_{12}v_2) - \frac{1715}{48} \pi^2 (n_{12}v_1)(n_{12}v_2) \right. \\
&\qquad \left. + 108 \ln\left(\frac{r_{12}}{r'_0}\right) (n_{12}v_1)(n_{12}v_2) \right. \\
&\qquad \left. - \frac{2233689}{5600} (n_{12}v_2)^2 + \frac{2771}{96} \pi^2 (n_{12}v_2)^2 - \frac{82}{3} \ln\left(\frac{r_{12}}{r'_0}\right) (n_{12}v_2)^2 \right. \\
&\qquad \left. - \frac{959797}{8400} (v_1 v_2) + \frac{103}{16} \pi^2 (v_1 v_2) - \frac{202}{3} \ln\left(\frac{r_{12}}{r'_0}\right) (v_1 v_2) \right)
\end{align*}
$$
$$ + \frac{858533}{50400} v_1^2 - \frac{15}{32} \pi^2 v_1^2 + \frac{121}{3} \ln\left(\frac{r_{12}}{r'_0}\right) v_1^2 + \frac{5482669}{50400} v_2^2 \\ - \frac{191}{32} \pi^2 v_2^2 + \frac{70}{3} \ln\left(\frac{r_{12}}{r'_0}\right) v_2^2 \right) + 1 \leftrightarrow 2, \qquad (350e) $$

$$ L_{4PN}^{(5)} = \frac{3 m_1^5 m_2}{8 r_{12}^5} + \frac{m_1^3 m_2^3}{r_{12}^5} \left( \frac{597771}{5600} - \frac{71}{32} \pi^2 - \frac{110}{3} \ln\left(\frac{r_{12}}{r'_0}\right) \right) \\ + \frac{m_1^4 m_2^2}{r_{12}^5} \left( \frac{1734977}{25200} + \frac{105}{32} \pi^2 - \frac{290}{3} \ln\left(\frac{r_{12}}{r'_0}\right) \right) + 1 \leftrightarrow 2. \qquad (350f) $$

By a contact transformation at 4PN order, extending Eq. (344), we can transform the 4PN harmonic-coordinates Lagrangian into a new Lagrangian, valid in some ADM or ADM-like coordinate system, and such that the associated Hamiltonian coincides with the 4PN Hamiltonian obtained with the ADM Hamiltonian formalism of Jaranowski and Schäfer (2013, 2015), Damour et al. (2014, 2015). Denoting by $P_1$ and $P_2$ the canonical momenta in ADM coordinates, and by $R_{12}N_{12} = Y_1 - Y_2$ the relative separation, the 4PN ADM Hamiltonian is [following the same presentation as in (342) and (349)]

$$ H_N = \frac{P_1^2}{2m_1} - \frac{Gm_1m_2}{2R_{12}} + 1 \leftrightarrow 2, \qquad (351a) $$

$$ H_{1PN} = -\frac{P_1^4}{8m_1^3} + \frac{G^2m_1^2m_2}{2R_{12}^2} \\ + \frac{Gm_1m_2}{R_{12}} \left( -\frac{3P_1^2}{2m_1^2} + \frac{(N_{12}P_1)(N_{12}P_2)}{4m_1m_2} + \frac{7(P_1P_2)}{4m_1m_2} \right) + 1 \leftrightarrow 2, \qquad (351b) $$

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

and, for the instantaneous (local-in-time) part at 4PN order,

$$
H_{4PN}^{(0)} = \frac{7P_1^{10}}{256m_1^9} + 1 \leftrightarrow 2, \qquad (352a)
$$
$$
\begin{equation}
\begin{aligned}
H_{4\text{PN}}^{(1)} ={}& \frac{m_1 m_2}{R_{12}} \left\{ \frac{45 P_1^8}{128 m_1^8} \right. \\
& \qquad + \frac{1}{m_1^6 m_2^2} \left( -\frac{9}{64} (N_{12}P_1)^2 (N_{12}P_2)^2 P_1^4 + \frac{15}{64} (N_{12}P_1)^2 P_1^4 P_2^2 \right. \\
& \qquad \left. -\frac{9}{16} (N_{12}P_1)(N_{12}P_2) P_1^4 (P_1 P_2) + \frac{15(N_{12}P_2)^2 P_1^6}{64} \right. \\
& \qquad \left. -\frac{21 P_1^6 P_2^2}{64} - \frac{3 P_1^4 (P_1 P_2)^2}{32} \right) \\
& \qquad + \frac{1}{m_1^5 m_2^3} \left( -\frac{35}{256} (N_{12}P_1)^5 (N_{12}P_2)^3 - \frac{85}{256} (N_{12}P_1)^4 (N_{12}P_2)^2 (P_1 P_2) \right. \\
& \qquad \left. + \frac{25}{128} (N_{12}P_1)^3 (N_{12}P_2)^3 P_1^2 - \frac{7}{128} (N_{12}P_1)^3 (N_{12}P_2) P_1^2 P_2^2 \right. \\
& \qquad \left. + \frac{25}{64} (N_{12}P_1)^3 (N_{12}P_2) (P_1 P_2)^2 - \frac{45}{128} (N_{12}P_1)^2 (N_{12}P_2)^2 P_1^2 (P_1 P_2) \right. \\
& \qquad \left. + \frac{7}{128} (N_{12}P_1)^2 P_1^2 (P_1 P_2) P_2^2 - \frac{3(N_{12}P_1)^2 (P_1 P_2)^3}{64} \right. \\
& \qquad \left. + \frac{33}{256} (N_{12}P_1)(N_{12}P_2)^3 P_1^4 - \frac{25}{256} (N_{12}P_1)(N_{12}P_2) P_1^4 P_2^2 \right. \\
& \qquad \left. + \frac{7}{64} (N_{12}P_1)(N_{12}P_2) P_1^2 (P_1 P_2)^2 - \frac{1}{256} (N_{12}P_2)^2 P_1^4 (P_1 P_2) \right. \\
& \qquad \left. - \frac{7}{256} P_1^4 (P_1 P_2) P_2^2 + \frac{3P_1^2 (P_1 P_2)^3}{64} \right) \\
& \qquad + \frac{1}{m_1^4 m_2^4} \left( -\frac{5}{32} (N_{12}P_1)^2 (N_{12}P_2)^4 P_1^2 + \frac{21}{64} (N_{12}P_1)^2 (N_{12}P_2)^2 P_1^2 P_2^2 \right. \\
& \qquad \left. -\frac{1}{2} (N_{12}P_1)(N_{12}P_2)^3 P_1^2 (P_1 P_2) + \frac{1}{16} (N_{12}P_1)(N_{12}P_2) P_1^2 (P_1 P_2) P_2^2 \right. \\
& \qquad \left. + \frac{7(N_{12}P_2)^4 P_1^4}{32} - \frac{3}{16} (N_{12}P_2)^2 P_1^4 P_2^2 + \frac{1}{8} (N_{12}P_2)^2 P_1^2 (P_1 P_2)^2 \right. \\
& \qquad \left. - \frac{7P_1^4 P_2^4}{128} - \frac{1}{32} P_1^2 (P_1 P_2)^2 P_2^2 \right) \\
& \qquad + \frac{1}{m_1^3 m_2^5} \left( \frac{55}{256} (N_{12}P_1)(N_{12}P_2)^5 P_1^2 - \frac{23}{256} (N_{12}P_2)^4 P_1^2 (P_1 P_2) \right) \right\} \\
& \qquad + 1 \leftrightarrow 2,
\end{aligned}
\tag{352b}
\end{equation}
$$
$$
\begin{aligned}
H_{4\text{PN}}^{(2)} = & \frac{m_1 m_2}{R_{12}^2} \left\{ \frac{13 m_2 P_1^6}{8 m_1^6} \right. \\
& + \frac{1}{m_1^5} \left( \frac{369 (N_{12} P_1)^6}{160} - \frac{889 (N_{12} P_1)^4 P_1^2}{192} + \frac{49 (N_{12} P_1)^2 P_1^4}{16} - \frac{63 P_1^6}{64} \right) \\
& + \frac{1}{m_1^4 m_2} \left( - \frac{549 (N_{12} P_1)^5 (N_{12} P_2)}{128} - \frac{7 (N_{12} P_1)^4 (N_{12} P_2)^2}{4} \right. \\
& \qquad \left. + \frac{1547 (N_{12} P_1)^4 (P_1 P_2)}{256} + \frac{67}{16} (N_{12} P_1)^3 (N_{12} P_2) P_1^2 \right. \\
& \qquad \left. - \frac{7}{2} (N_{12} P_1)^3 (N_{12} P_2) (P_1 P_2) + \frac{7}{6} (N_{12} P_1)^2 (N_{12} P_2)^2 P_1^2 \right. \\
& \qquad \left. - \frac{851}{128} (N_{12} P_1)^2 P_1^2 (P_1 P_2) + \frac{49}{48} (N_{12} P_1)^2 P_1^2 P_2^2 + \frac{21 (N_{12} P_1)^2 (P_1 P_2)^2}{16} \right. \\
& \qquad \left. - \frac{167}{128} (N_{12} P_1) (N_{12} P_2) P_1^4 - \frac{133}{24} (N_{12} P_1) (N_{12} P_2) P_1^2 (P_1 P_2) \right. \\
& \qquad \left. + \frac{197 (N_{12} P_2)^2 P_1^4}{96} + \frac{1099 P_1^4 (P_1 P_2)}{256} - \frac{173 P_1^4 P_2^2}{48} - \frac{77 P_1^2 (P_1 P_2)^2}{96} \right) \\
& + \frac{1}{m_1^3 m_2^2} \left( \frac{3263 (N_{12} P_1)^4 (N_{12} P_2)^2}{1280} - \frac{3571}{320} (N_{12} P_1)^3 (N_{12} P_2) (P_1 P_2) \right. \\
& \qquad \left. + \frac{1067}{480} (N_{12} P_1)^2 (N_{12} P_2)^2 P_1^2 - \frac{1999 (N_{12} P_1)^2 P_1^2 P_2^2}{3840} \right. \\
& \qquad \left. + \frac{4349 (N_{12} P_1)^2 (P_1 P_2)^2}{1280} + \frac{233}{96} (N_{12} P_1) (N_{12} P_2)^3 P_1^2 \right. \\
& \qquad \left. + \frac{3073}{480} (N_{12} P_1) (N_{12} P_2) P_1^2 (P_1 P_2) - \frac{4567 (N_{12} P_2)^2 P_1^4}{3840} \right. \\
& \qquad \left. + \frac{1}{384} (N_{12} P_2)^2 P_1^2 (P_1 P_2) + \frac{2081 P_1^4 P_2^2}{3840} - \frac{3461 P_1^2 (P_1 P_2)^2}{3840} \right) \\
& + \frac{1}{m_1^2 m_2^3} \left( - \frac{13 (N_{12} P_1)^3 (N_{12} P_2)^3}{8} - \frac{19}{384} (N_{12} P_1)^2 (N_{12} P_2)^2 (P_1 P_2) \right)
\end{aligned}
$$
$$ 
\begin{aligned}
& + \frac{191}{192}(N_{12}P_1)(N_{12}P_2)^3 P_1^2 - \frac{47}{32}(N_{12}P_1)(N_{12}P_2)P_1^2 P_2^2 \\
& + \frac{11}{192}(N_{12}P_1)(N_{12}P_2)(P_1P_2)^2 + \frac{1673(N_{12}P_2)^4 P_1^2}{1920} \\
& - \frac{5}{384}(N_{12}P_2)^2 P_1^2(P_1P_2) - \frac{185P_1^2(P_1P_2)P_2^2}{384} + \frac{77(P_1P_2)^3}{96} \\
& + \left\{ \frac{7(N_{12}P_2)^4 P_1^2}{4m_1m_2^4} \right\} + 1 \leftrightarrow 2,
\end{aligned} 
\quad (352c) $$

$$ 
\begin{aligned}
H_{4PN}^{(3)} = & \frac{m_1 m_2}{R_{12}^3} \left\{ \frac{105 m_2^2 P_1^4}{32 m_1^4} \right. \\
& + \frac{m_2}{m_1^3} \left( \frac{375 \pi^2 (N_{12} P_1)^4}{8192} - \frac{23533 (N_{12} P_1)^4}{1280} + \frac{63347 (N_{12} P_1)^2 P_1^2}{1600} \right. \\
& \left. - \frac{1059 \pi^2 (N_{12} P_1)^2 P_1^2}{1024} - \frac{211189 P_1^4}{19200} + \frac{2749 \pi^2 P_1^4}{8192} \right) \\
& + \frac{1}{m_1^2} \left( \frac{5027 (N_{12} P_1)^4}{384} + \frac{2369 (N_{12} P_1)^3 (N_{12} P_2)}{60} \right. \\
& \left. + \frac{35655 \pi^2 (N_{12} P_1)^3 (N_{12} P_2)}{16384} - \frac{22993 (N_{12} P_1)^2 P_1^2}{960} \right. \\
& \left. - \frac{1646983 (N_{12} P_1)^2 (P_1 P_2)}{19200} + \frac{56955 \pi^2 (N_{12} P_1)^2 (P_1 P_2)}{16384} \right. \\
& \left. - \frac{391711 (N_{12} P_1) (N_{12} P_2) P_1^2}{6400} + \frac{43101 \pi^2 (N_{12} P_1) (N_{12} P_2) P_1^2}{16384} \right. \\
& \left. + \frac{26627 (N_{12} P_2)^2 P_1^2}{1600} - \frac{6695 P_1^4}{1152} + \frac{1243717 P_1^2 (P_1 P_2)}{14400} - \frac{40483 \pi^2 P_1^2 (P_1 P_2)}{16384} \right) \\
& + \frac{1}{m_1 m_2} \left( - \frac{3191 (N_{12} P_1)^3 (N_{12} P_2)}{640} - \frac{30383 (N_{12} P_1)^2 (N_{12} P_2)^2}{960} \right. \\
& \left. - \frac{36405 \pi^2 (N_{12} P_1)^2 (N_{12} P_2)^2}{16384} + \frac{8777 (N_{12} P_1)^2 (P_1 P_2)}{384} \right. \\
& \left. + \frac{28561 (N_{12} P_1) (N_{12} P_2) P_1^2}{1920} + \frac{248991 (N_{12} P_1) (N_{12} P_2) (P_1 P_2)}{6400} \right. \\
& \left. - \frac{6153 \pi^2 (N_{12} P_1) (N_{12} P_2) (P_1 P_2)}{2048} + \frac{1411429 (N_{12} P_2)^2 P_1^2}{19200} \right. \\
& \left. - \frac{1059}{512} \pi^2 (N_{12} P_2)^2 P_1^2 + \frac{752969 P_1^2 (P_1 P_2)}{28800} - \frac{2492417 P_1^2 P_2^2}{57600} \right)
\end{aligned}
$$
$$ 
\begin{aligned}
&+ \frac{13723\pi^2 P_1^2 P_2^2}{16384} - \frac{1918349(P_1P_2)^2}{57600} + \frac{10631\pi^2(P_1P_2)^2}{8192} \\
&+ \frac{1}{m_2^2} \left( -\frac{16481}{960}(N_{12}P_1)^2(N_{12}P_2)^2 - \frac{103957(N_{12}P_1)(N_{12}P_2)(P_1P_2)}{2400} \right) \\
&+ \frac{94433(N_{12}P_2)^2 P_1^2}{4800} - \frac{118261P_1^2 P_2^2}{4800} + \frac{791(P_1P_2)^2}{400} \right) \} + 1 \leftrightarrow 2,
\end{aligned}
\quad (352\text{d}) $$

$$ 
\begin{aligned}
H_{4\text{PN}}^{(4)} = & \frac{m_1 m_2}{R_{12}^4} \left\{ \frac{105 m_2^3 P_1^2}{32 m_1^2} \right. \\
& + \frac{m_1^2}{m_2} \left( -\frac{7(N_{12}P_1)(N_{12}P_2)}{2} - \frac{91(P_1P_2)}{8} \right) \\
& + \frac{m_2^2}{m_1} \left( -\frac{28691\pi^2(N_{12}P_1)^2}{24576} + \frac{3200179(N_{12}P_1)^2}{57600} + \frac{282361P_1^2}{19200} - \frac{21837\pi^2 P_1^2}{8192} \right) \\
& + m_1 \left( -\frac{9841(N_{12}P_1)^2}{1600} - \frac{2712013(N_{12}P_1)(N_{12}P_2)}{19200} + \frac{176033\pi^2(P_1P_2)}{24576} \right. \\
& \left. + \frac{63641\pi^2(N_{12}P_1)(N_{12}P_2)}{24576} + \frac{64861P_1^2}{4800} - \frac{2864917(P_1P_2)}{57600} \right) \\
& + m_2 \left( \frac{21745\pi^2(N_{12}P_1)^2}{16384} + \frac{698723(N_{12}P_1)^2}{19200} + \frac{1937033P_1^2}{57600} \right. \\
& \left. - \frac{199177\pi^2 P_1^2}{49152} \right) \left. \right\} + 1 \leftrightarrow 2,
\end{aligned}
\quad (352\text{e}) $$

$$ 
\begin{aligned}
H_{4\text{PN}}^{(5)} = & \frac{m_1 m_2}{R_{12}^5} \left\{ -\frac{m_1^4}{16} + \left( \frac{6237\pi^2}{1024} - \frac{169799}{2400} \right) m_1^3 m_2 \right. \\
& \left. + \left( \frac{44825\pi^2}{6144} - \frac{609427}{7200} \right) m_1^2 m_2^2 \right\} + 1 \leftrightarrow 2.
\end{aligned}
\quad (352\text{f}) $$

The non-local tail term in the Hamiltonian is just the opposite of the tail term in the Lagrangian: $H_{4\text{PN}}^{\text{tail}} = -L_{4\text{PN}}^{\text{tail}}$ (Damour et al. 2014; Bernard et al. 2016).

Arguably, the result given by the ADM-Hamiltonian formalism looks simpler than the harmonic-coordinate counterpart. Indeed, the Lagrangian is ordinary—no accelerations—and there are no logarithms nor associated gauge constant $r'_0$. Of course, one is free to describe the binary motion in whatever coordinates one likes, and the two formalisms, harmonic and ADM describe rigorously the same physics. On the other hand, the higher complexity of the harmonic-coordinates formalism enables one to perform more tests of the computations, notably by inquiring about
the fate of the constant $r'_0$, which must disappear from physical quantities such as the center-of-mass energy and the total gravitational-wave flux. On the other hand, the ADM-Hamiltonian formalism currently provides a limited description of the radiation field, compared to what will be done using harmonic coordinates in Sect. 3.4. Therefore it is important to derive the equations of motion in harmonic coordinates for consistency with the calculation of the radiation field.

We have reviewed in Sect. 1.2.1 the efforts made to extend the equations of motion beyond the current state-of-the-art 4PN order.

### 3.2.4 Equations of motion in the frame of the center of mass

We translate the origin of coordinates to the binary's center-of-mass (CM) by imposing the vanishing of the binary's mass dipole moment: $I_i = 0$ in the notation of Sect. 2.1.2. Actually the dipole moment is computed as the center-of-mass conserved integral $G \equiv I$ associated with the boost symmetry of the equations of motion (Damour et al. 2000d; de Andrade et al. 2001; Blanchet and Iyer 2003; Blanchet and Faye 2019; Bernard et al. 2018). The condition $G = 0$ is solved iteratively with usual order reduction of accelerations using the CM equations of motion. This results in the relationship between the individual positions in the center-of-mass frame $y_1$ and $y_2$, and the relative position $x \equiv y_1 - y_2$ and velocity $v \equiv v_1 - v_2 = dx/dt$ (formerly denoted $y_{12}$ and $v_{12}$). We also use the orbital separation $r \equiv |x|$, together with $n = x/r$ and $\dot{r} = (nv) = n \cdot v$. Mass parameters are: The total mass $m = m_1 + m_2$ (to be distinguished from the conserved ADM mass M); the relative mass difference $\Delta = (m_1 - m_2)/m$; the reduced mass $\mu = m_1 m_2 / m$; the symmetric mass ratio

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

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

$$y_1 = [X_2 + v \Delta \mathcal{P}] x + v \Delta \mathcal{Q} v, \qquad (354a)$$

$$y_2 = [-X_1 + v \Delta \mathcal{P}] x + v \Delta \mathcal{Q} v, \qquad (354b)$$

where all post-Newtonian corrections, beyond Newtonian order, are proportional to the mass ratio $v$ and the mass difference $\Delta$. And, of course, we have $y_1 = -y_2$ for equal masses. Up to the 4PN order the two dimensionless coefficients $\mathcal{P}$ and $\mathcal{Q}$ read [with presentation similar to (342) and (349); the PN coefficients not indicated are zero]
$$P_{1PN} = \frac{v^2}{2} - \frac{Gm}{2r}, \qquad (355a)$$

$$ \begin{aligned} P_{2PN} = & \frac{3 v^4}{8} - \frac{3 v v^4}{2} + \frac{Gm}{r} \left( -\frac{\dot{r}^2}{8} + \frac{3 \dot{r}^2 v}{4} + \frac{19 v^2}{8} + \frac{3 v v^2}{2} \right) \\ & + \frac{G^2 m^2}{r^2} \left( \frac{7}{4} - \frac{v}{2} \right), \end{aligned} \qquad (355b) $$

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

$$ P_{3.5PN} = \frac{G^2 m^2}{r^2} \dot{r} \left( -\frac{8 \dot{r}^2}{5} - \frac{2 v^2}{15} + \frac{68 v^2 v}{35} \right) + \frac{G^3 m^3}{r^3} \dot{r} \left( \frac{52}{15} - \frac{116}{35} v \right), \qquad (355d) $$

$$ P_{4PN}^{(0)} = \left( \frac{35}{128} - \frac{125}{32} v + \frac{145}{8} v^2 - \frac{55}{2} v^3 \right) v^8, \qquad (355e) $$

$$ \begin{aligned} P_{4PN}^{(1)} = & \frac{m}{r} \left( -\frac{5}{128} \dot{r}^6 + \frac{35}{64} v \dot{r}^6 - \frac{125}{64} v^2 \dot{r}^6 + \frac{55}{32} v^3 \dot{r}^6 + \frac{27}{128} \dot{r}^4 v^2 - \frac{115}{64} v \dot{r}^4 v^2 \right. \\ & \left. + \frac{517}{64} v^2 \dot{r}^4 v^2 - \frac{213}{16} v^3 \dot{r}^4 v^2 - \frac{53}{128} \dot{r}^2 v^4 + \frac{3}{2} v \dot{r}^2 v^4 - \frac{95}{8} v^2 \dot{r}^2 v^4 \right. \\ & \left. + 36 v^3 \dot{r}^2 v^4 + \frac{515}{128} v^6 - \frac{749}{32} v v^6 + \frac{91}{4} v^2 v^6 + 42 v^3 v^6 \right), \end{aligned} \qquad (355f) $$

$$ \begin{aligned} P_{4PN}^{(2)} = & \frac{m^2}{r^2} \left( \frac{1133}{960} \dot{r}^4 - \frac{1007}{48} v \dot{r}^4 + \frac{169}{24} v^2 \dot{r}^4 + 9 v^3 \dot{r}^4 - \frac{31}{5} \dot{r}^2 v^2 + 26 v \dot{r}^2 v^2 \right. \\ & \left. - \frac{541}{8} v^2 \dot{r}^2 v^2 - \frac{83}{2} v^3 \dot{r}^2 v^2 + \frac{5631}{320} v^4 - \frac{139}{4} v v^4 + \frac{71}{4} v^2 v^4 - \frac{45}{2} v^3 v^4 \right), \end{aligned} \qquad (355g) $$
$$ P_{4PN}^{(3)} = \frac{m^3}{r^3} \left( -\frac{185497}{8400} r^2 - \frac{64347}{1120} v r^2 - \frac{123}{128} \pi^2 v r^2 + \frac{495}{16} v^2 r^2 + \frac{55}{4} v^3 r^2 \right. \\ \left. - 11 \ln\left(\frac{r}{r_0}\right) r^2 + 44 v \ln\left(\frac{r}{r_0}\right) r^2 + \frac{2737}{1440} v^2 - \frac{87181}{3360} v v^2 + \frac{123}{128} \pi^2 v v^2 \right. \\ \left. - \frac{117}{8} v^2 v^2 + 5 v^3 v^2 - 11 v \ln\left(\frac{r}{r_0}\right) v^2 + 22 \ln\left(\frac{r}{r_0}\right) v^2 \right), \quad (355h) $$

$$ P_{4PN}^{(4)} = \frac{m^4}{r^4} \left( \frac{215279}{3600} + \frac{22043}{720} v - \frac{11}{2} \pi^2 v - \frac{3}{2} v^2 - \frac{1}{2} v^3 \right. \\ \left. - \frac{220}{3} \ln\left(\frac{r}{r_0}\right) + 60 v \ln\left(\frac{r}{r_0}\right) \right). \quad (355i) $$

and

$$ Q_{2PN} = -\frac{7 G m \dot{r}}{4}, \quad (356a) $$

$$ Q_{2.5PN} = \frac{4 G m v^2}{5} - \frac{8 G^2 m^2}{5 r}, \quad (356b) $$

$$ Q_{3PN} = G m \dot{r} \left( \frac{5 \dot{r}^2}{12} - \frac{19 \dot{r}^2 v}{24} - \frac{15 v^2}{8} + \frac{21 v v^2}{4} \right) \\ + \frac{G^2 m^2 \dot{r}}{r} \left( -\frac{235}{24} - \frac{21 v}{4} \right), \quad (356c) $$

$$ Q_{3.5PN} = Gm v^4 \left( \frac{6}{7} - \frac{22 v}{7} \right) + \frac{G^2 m^2}{r} \left( \frac{44}{15} \dot{r}^2 + \frac{6}{35} v^2 + \frac{68}{35} r^2 v + \frac{132}{35} v^2 v \right) \\ + \frac{G^3 m^3}{r^2} \left( -\frac{172}{105} - \frac{64}{35} v \right), \quad (356d) $$

$$ Q_{4PN}^{(1)} = m \left( -\frac{13}{64} \dot{r}^5 + \frac{25}{32} v \dot{r}^5 - \frac{17}{32} v^2 \dot{r}^5 + \frac{77}{96} r^3 v^2 - \frac{187}{48} v \dot{r}^3 v^2 + \frac{19}{4} v^2 r^3 v^2 \right. \\ \left. - \frac{123}{64} r v^4 + \frac{199}{16} v \dot{r} v^4 - 21 v^2 \dot{r} v^4 \right), \quad (356e) $$

$$ Q_{4PN}^{(2)} = \frac{m^2}{r} \left( \frac{4621}{480} \dot{r}^3 + \frac{113}{24} v \dot{r}^3 + \frac{7}{12} v^2 \dot{r}^3 - \frac{3733}{160} \dot{r} v^2 + \frac{95}{4} v \dot{r} v^2 + 28 v^2 \dot{r} v^2 \right), \quad (356f) $$

$$ Q_{4PN}^{(3)} = \frac{m^3}{r^2} \left( \frac{14377}{280} + \frac{71509}{5040} v - \frac{41}{64} \pi^2 v - \frac{49}{4} v^2 + \frac{22}{3} v \ln\left(\frac{r}{r_0}\right) \right. \\ \left. - \frac{110}{3} \ln\left(\frac{r}{r_0}\right) - \frac{44}{3} v \ln\left(\frac{r}{r_0}\right) \right) \dot{r}. \quad (356g) $$
The CM velocities are obtained by differentiating the above expressions—with order reduction of accelerations. There are no tail terms, but the formulas contain logarithmic terms depending on the gauge constant $r'_0$ not affecting physical results.

To proceed in the CM frame we consider first the relative acceleration $\mathbf{a} = \mathbf{a}_1 - \mathbf{a}_2$ and then make the replacements (354). The general form of the order-reduced acceleration (considering first the instantaneous, local part of the acceleration) is

$$ \mathbf{a}_{\text{inst}} = -\frac{Gm}{r^2} [(1 + \mathcal{A}) \mathbf{n} + \mathcal{B} \mathbf{v}], \quad (357) $$

where the various PN coefficients up to 4PN order are given by (Lincoln and Will 1990; Blanchet and Iyer 2003; Bernard et al. 2018)

$$ \mathcal{A}_{1\text{PN}} = -\frac{3 \dot{r}^2 v}{2} + v^2 + 3 v v^2 - \frac{Gm}{r} (4 + 2 v), \quad (358\text{a}) $$

$$ \begin{aligned} \mathcal{A}_{2\text{PN}} = & \frac{15 \dot{r}^4 v}{8} - \frac{45 \dot{r}^4 v^2}{8} - \frac{9 \dot{r}^2 v v^2}{2} + 6 \dot{r}^2 v^2 v^2 + 3 v v^4 - 4 v^2 v^4 \\ & + \frac{Gm}{r} \left( -2 \dot{r}^2 - 25 \dot{r}^2 v - 2 \dot{r}^2 v^2 - \frac{13 v v^2}{2} + 2 v^2 v^2 \right) \\ & + \frac{G^2 m^2}{r^2} \left( 9 + \frac{87 v}{4} \right), \end{aligned} \quad (358\text{b}) $$

$$ \mathcal{A}_{2.5\text{PN}} = \frac{8Gmv}{5r}\dot{r}\left[-\frac{17}{3}\frac{Gm}{r} - 3v^2\right], \quad (358c) $$

$$ \begin{aligned} \mathcal{A}_{3\text{PN}} = & -\frac{35 \dot{r}^6 v}{16} + \frac{175 \dot{r}^6 v^2}{16} - \frac{175 \dot{r}^6 v^3}{16} + \frac{15 \dot{r}^4 v v^2}{2} \\ & - \frac{135 \dot{r}^4 v^2 v^2}{4} + \frac{255 \dot{r}^4 v^3 v^2}{8} - \frac{15 \dot{r}^2 v v^4}{2} + \frac{237 \dot{r}^2 v^2 v^4}{8} \\ & - \frac{45 \dot{r}^2 v^3 v^4}{2} + \frac{11 v v^6}{4} - \frac{49 v^2 v^6}{4} + 13 v^3 v^6 \\ & + \frac{Gm}{r} \left( 79 \dot{r}^4 v - \frac{69 \dot{r}^4 v^2}{2} - 30 \dot{r}^4 v^3 - 121 \dot{r}^2 v v^2 + 16 \dot{r}^2 v^2 v^2 \right. \\ & \left. + 20 \dot{r}^2 v^3 v^2 + \frac{75 v v^4}{4} + 8 v^2 v^4 - 10 v^3 v^4 \right) \\ & + \frac{G^2 m^2}{r^2} \left( \dot{r}^2 + \frac{32573 \dot{r}^2 v}{168} + \frac{11 \dot{r}^2 v^2}{8} - 7 \dot{r}^2 v^3 + \frac{615 \dot{r}^2 v \pi^2}{64} \right. \\ & \left. - \frac{26987 v v^2}{840} + v^3 v^2 - \frac{123 v \pi^2 v^2}{64} \right. \\ & \left. - 110 \dot{r}^2 v \ln\left(\frac{r}{r_0}\right) + 22 v v^2 \ln\left(\frac{r}{r_0}\right) \right) \\ & + \frac{G^3 m^3}{r^3} \left( -16 - \frac{437 v}{4} - \frac{71 v^2}{2} + \frac{41 v \pi^2}{16} \right), \end{aligned} \quad (358\text{d}) $$
$$A_{3.5PN} = \frac{Gmv}{r} \dot{r} \left[ \frac{G^2 m^2}{r^2} \left( \frac{3956}{35} + \frac{184}{5} v \right) + \frac{Gm v^2}{r} \left( \frac{692}{35} - \frac{724}{15} v \right) \right. \\ \left. + v^4 \left( \frac{366}{35} + 12v \right) + \frac{Gm \dot{r}^2}{r} \left( \frac{294}{5} + \frac{376}{5} v \right) \right. \\ \left. - v^2 \dot{r}^2 (114 + 12v) + 112 \dot{r}^4 \right], \qquad (358e)$$

$$A_{4PN}^{(0)} = \left( \frac{315}{128} v - \frac{2205}{128} v^2 + \frac{2205}{64} v^3 - \frac{2205}{128} v^4 \right) \dot{r}^8 \\ + \left( -\frac{175}{16} v + \frac{595}{8} v^2 - \frac{2415}{16} v^3 + \frac{735}{8} v^4 \right) \dot{r}^6 v^2 \\ + \left( \frac{135}{8} v - \frac{1875}{16} v^2 + \frac{4035}{16} v^3 - \frac{1335}{8} v^4 \right) \dot{r}^4 v^4 \\ + \left( -\frac{21}{2} v + \frac{1191}{16} v^2 - \frac{327}{2} v^3 + 99 v^4 \right) \dot{r}^2 v^6 \\ + \left( \frac{21}{8} v - \frac{175}{8} v^2 + 61 v^3 - 54 v^4 \right) v^8, \qquad (358f)$$

$$A_{4PN}^{(1)} = \frac{m}{r} \left( \frac{2973}{40} v \dot{r}^6 + 407 v^2 \dot{r}^6 + \frac{181}{2} v^3 \dot{r}^6 - 86 v^4 \dot{r}^6 + \frac{1497}{32} v^4 \dot{r}^2 \right. \\ \left. - \frac{1627}{2} v^2 \dot{r}^4 v^2 - 81 v^3 \dot{r}^4 v^2 + 228 v^4 \dot{r}^4 v^2 - \frac{2583}{16} v \dot{r}^2 v^4 \right. \\ \left. + \frac{1009}{2} v^2 \dot{r}^2 v^4 + 47 v^3 \dot{r}^2 v^4 - 104 v^4 \dot{r}^2 v^4 + \frac{1067}{32} v v^6 \right. \\ \left. - 58 v^2 v^6 - 44 v^3 v^6 + 58 v^4 v^6 \right), \qquad (358g)$$

$$A_{4PN}^{(2)} = \frac{m^2}{r^2} \left( \frac{2094751}{960} v \dot{r}^4 + \frac{45255}{1024} \pi^2 v \dot{r}^4 + \frac{326101}{96} v^2 \dot{r}^4 - \frac{4305}{128} \pi^2 v^2 \dot{r}^2 \right. \\ \left. - \frac{1959}{32} v^3 \dot{r}^4 - 126 v^4 \dot{r}^4 - 1155 v^2 \ln\left(\frac{r}{r_0}\right) \dot{r}^4 + 385 v \ln\left(\frac{r}{r_0}\right) \dot{r}^4 \right. \\ \left. - \frac{1636681}{1120} v \dot{r}^2 v^2 - \frac{12585}{512} \pi^2 v \dot{r}^2 v^2 - \frac{255461}{112} v^2 \dot{r}^2 v^2 \right. \\ \left. + \frac{3075}{128} \pi^2 v^2 \dot{r}^2 v^2 - \frac{309}{4} v^3 \dot{r}^2 v^2 + 63 v^4 \dot{r}^2 v^2 \right. \\ \left. - 605 v \ln\left(\frac{r}{r_0}\right) \dot{r}^2 v^2 + 825 v^2 \ln\left(\frac{r}{r_0}\right) \dot{r}^2 v^2 + \frac{1096941}{11200} v v^4 \right. \\ \left. + \frac{1155}{1024} \pi^2 v v^4 + \frac{7263}{70} v^2 v^4 - \frac{123}{64} \pi^2 v^2 v^4 + \frac{145}{2} v^3 v^4 - 16 v^4 v^4 \right. \\ \left. + 88 v \ln\left(\frac{r}{r_0}\right) v^4 - 66 v^2 \ln\left(\frac{r}{r_0}\right) v^4 \right), \qquad (358h)$$
$$
\begin{equation}
\begin{aligned}
\mathcal{A}_{4\text{PN}}^{(3)} = & \frac{m^3}{r^3} \left( -2\dot{r}^2 + \frac{1297943}{8400} v\dot{r}^2 - \frac{2969}{16}\pi^2 v\dot{r}^2 + \frac{1255151}{840} v^2 \dot{r}^2 \right. \\
& \left. + \frac{7095}{32}\pi^2 v^2 \dot{r}^2 - 17v^3 \dot{r}^2 - 24v^4 \dot{r}^2 \right. \\
& \left. + 2232v \ln\left(\frac{r}{r_0'}\right) \dot{r}^2 - 3364v^2 \ln\left(\frac{r}{r_0'}\right) \dot{r}^2 \right. \\
& \left. + \frac{1237279}{25200} v v^2 + \frac{3835}{96} \pi^2 v v^2 - \frac{693947}{2520} v^2 v^2 - \frac{229}{8} \pi^2 v^2 v^2 \right. \\
& \left. + \frac{19}{2} v^3 v^2 - \frac{1336}{3} v \ln\left(\frac{r}{r_0'}\right) v^2 + \frac{1616}{3} v^2 \ln\left(\frac{r}{r_0'}\right) v^2 \right),
\end{aligned}
\tag{358i}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
\mathcal{A}_{4\text{PN}}^{(4)} = & \frac{m^4}{r^4} \left( 25 + \frac{6625537}{12600} v - \frac{4543}{96} \pi^2 v + \frac{477763}{720} v^2 + \frac{3}{4} \pi^2 v^2 \right. \\
& \left. + \frac{334}{3} v \ln\left(\frac{r}{r_0'}\right) - \frac{514}{3} v^2 \ln\left(\frac{r}{r_0'}\right) \right).
\end{aligned}
\tag{358j}
\end{equation}
$$

$$
\mathcal{B}_{1\text{PN}} = -4 \dot{r} + 2 \dot{r} v, \tag{359a}
$$

$$
\begin{equation}
\begin{aligned}
\mathcal{B}_{2\text{PN}} = & \frac{9 \dot{r}^3 v}{2} + 3 \dot{r}^3 v^2 - \frac{15 \dot{r} v v^2}{2} - 2 \dot{r} v^2 v^2 \\
& + \frac{Gm}{r} \left( 2 \dot{r} + \frac{41 \dot{r} v}{2} + 4 \dot{r} v^2 \right),
\end{aligned}
\tag{359b}
\end{equation}
$$

$$
\mathcal{B}_{2.5\text{PN}} = \frac{8Gmv}{5r} \left[ 3 \frac{Gm}{r} + v^2 \right], \qquad (359c)
$$

$$
\begin{equation}
\begin{aligned}
\mathcal{B}_{3\text{PN}} = & -\frac{45 \dot{r}^5 v}{8} + 15 \dot{r}^5 v^2 + \frac{15 \dot{r}^5 v^3}{4} + 12 \dot{r}^3 v v^2 \\
& -\frac{111 \dot{r}^3 v^2 v^2}{4} - 12 \dot{r}^3 v^3 v^2 - \frac{65 \dot{r} v v^4}{8} + 19 \dot{r} v^2 v^4 + 6 \dot{r} v^3 v^4 \\
& + \frac{Gm}{r} \left( \frac{329 \dot{r}^3 v}{6} + \frac{59 \dot{r}^3 v^2}{2} + 18 \dot{r}^3 v^3 - 15 \dot{r} v v^2 - 27 \dot{r} v^2 v^2 - 10 \dot{r} v^3 v^2 \right) \\
& + \frac{G^2 m^2}{r^2} \left( -4 \dot{r} - \frac{18169 \dot{r} v}{840} + 25 \dot{r} v^2 + 8 \dot{r} v^3 - \frac{123 \dot{r} v \pi^2}{32} \right. \\
& \left. \qquad + 44 \dot{r} v \ln \left( \frac{r}{r_0'} \right) \right),
\end{aligned}
\tag{359d}
\end{equation}
$$
$$B_{3.5PN} = \frac{Gmv}{r} \left[ \frac{G^2 m^2}{r^2} \left( -\frac{1060}{21} - \frac{104}{5} v \right) + \frac{Gmv^2}{r} \left( \frac{164}{21} + \frac{148}{5} v \right) \right. \\ \left. + v^4 \left( -\frac{626}{35} - \frac{12}{5} v \right) + \frac{Gm\ddot{r}^2}{r} \left( -\frac{82}{3} - \frac{848}{15} v \right) \right. \\ \left. + v^2 \dot{r}^2 \left( \frac{678}{5} + \frac{12}{5} v \right) - 120 \dot{r}^4 \right], \qquad (359e)$$

$$\begin{aligned} B_{4PN}^{(0)} = & \left( \frac{105}{16} v - \frac{245}{8} v^2 + \frac{385}{16} v^3 + \frac{35}{8} v^4 \right) \dot{r}^7 \\ & + \left( -\frac{165}{8} v + \frac{1665}{16} v^2 - \frac{1725}{16} v^3 - \frac{105}{4} v^4 \right) \dot{r}^5 v^2 \\ & + \left( \frac{45}{2} v - \frac{1869}{16} v^2 + 129 v^3 + 54 v^4 \right) \dot{r}^3 v^4 \\ & + \left( -\frac{157}{16} v + 54 v^2 - 69 v^3 - 24 v^4 \right) \dot{r}^6, \end{aligned} \qquad (359f)$$

$$\begin{aligned} B_{4PN}^{(1)} = & \frac{m}{r} \left( -\frac{54319}{160} v \dot{r}^5 - \frac{901}{8} v^2 \dot{r}^5 + 60 v^3 \dot{r}^5 + 30 v^4 \dot{r}^5 + \frac{25943}{48} v \dot{r}^3 v^2 \right. \\ & \left. + \frac{1199}{12} v^2 \dot{r}^3 v^2 - \frac{349}{2} v^3 \dot{r}^3 v^2 - 98 v^4 \dot{r}^3 v^2 - \frac{5725}{32} v \dot{r} v^4 \right. \\ & \left. - \frac{389}{8} v^2 \dot{r} v^4 + 118 v^3 \dot{r} v^4 + 44 v^4 \dot{r} v^4 \right), \end{aligned} \qquad (359g)$$

$$\begin{aligned} B_{4PN}^{(2)} = & \frac{m^2}{r^2} \left( -\frac{9130111}{3360} v \dot{r}^3 - \frac{4695}{256} \pi^2 v \dot{r}^3 - \frac{184613}{112} v^2 \dot{r}^3 + \frac{1845}{64} \pi^2 v^2 \dot{r}^3 \right. \\ & \left. + \frac{209}{2} v^3 \dot{r}^3 + 74 v^4 \dot{r}^3 + 440 v \ln\left(\frac{r}{r_0'}\right) \dot{r}^3 + 550 v^2 \ln\left(\frac{r}{r_0'}\right) \dot{r}^3 \right. \\ & \left. + \frac{8692601}{5600} v \dot{r} v^2 + \frac{1455}{256} \pi^2 v \dot{r} v^2 + \frac{58557}{70} v^2 \dot{r} v^2 - \frac{123}{8} \pi^2 v^2 \dot{r} v^2 \right. \\ & \left. - 70 v^3 \dot{r} v^2 - 34 v^4 \dot{r} v^2 - 154 v \ln\left(\frac{r}{r_0'}\right) \dot{r} v^2 - 264 v^2 \ln\left(\frac{r}{r_0'}\right) \dot{r} v^2 \right), \end{aligned} \qquad (359h)$$

$$\begin{aligned} B_{4PN}^{(3)} = & \frac{m^3}{r^3} \left( 2 - \frac{619267}{525} v + \frac{791}{16} \pi^2 v - \frac{28406}{45} v^2 - \frac{2201}{32} \pi^2 v^2 \right. \\ & \left. + 66 v^3 + 16 v^4 - \frac{1484}{3} v \ln\left(\frac{r}{r_0'}\right) + 1268 v^2 \ln\left(\frac{r}{r_0'}\right) \dot{r} \right). \end{aligned} \qquad (359i)$$
The tail part of the acceleration comes from the variation of the tail part in the Lagrangian (348), taking into account the variation of the scale $r = r_{12}$ entering into the logarithmic kernel, to which one has to add the dissipative contribution. The result is

$$a_{\text{tail}}^{i} = -\frac{8G^2m}{5c^8} x^a \int_{0}^{+\infty} d\tau \ln\left(\frac{c\tau}{2r}\right) Q_{ia}^{(7)}(t-\tau) \\ + \frac{8G^2m}{5c^8} x^a \left[ \left(Q_{ia}^{(3)} \ln r\right)^{(3)} - Q_{ia}^{(6)} \ln r \right] - \frac{2G^2 n^i}{5c^8 v r} \left(Q_{ab}^{(3)}\right)^2, \qquad (360)$$

where to this order the ADM mass M reduces to the rest mass m and the quadrupole moment is the Newtonian one $Q_{ab} = mv^2 x^{(a}x^{b)}$.<fnref n="72" /> The result (360) is made of both conservative and dissipative parts according to the split (see the footnote 68 for the definition of the Hadamard Pf)

$$a_{\text{tail}}^{i} \Big|_{\text{cons}} = -\frac{4G^2m}{5c^8} x^a \text{Pf}_{2r/c} \int_{-\infty}^{+\infty} \frac{dt'}{|t-t'|} Q_{ia}^{(6)}(t') \\ + \frac{8G^2m}{5c^8} x^a \left[ \left(Q_{ia}^{(3)} \ln r\right)^{(3)} - Q_{ia}^{(6)} \ln r \right] - \frac{2G^2 n^i}{5c^8 v r} \left(Q_{ab}^{(3)}\right)^2, \qquad (361a)$$

$$a_{\text{tail}}^{i} \Big|_{\text{diss}} = -\frac{4G^2m}{5c^8} x^a \int_{0}^{+\infty} d\tau \ln\left(\frac{c\tau}{2r}\right) \left[Q_{ia}^{(7)}(t-\tau) + Q_{ia}^{(7)}(t+\tau)\right]. \qquad (361b)$$

The logarithms in (358)–(359), together with the constant $r'_0$ therein, can be removed at 3PN order by a simple coordinate transformation (see the footnote 71), while staying within the class of harmonic coordinates. The modification of the equations of motion affects the 3PN coefficients; we denote the new coefficients as $A_{3PN}^{\text{MH}}$ and $B_{3PN}^{\text{MH}}$ where MH stands for the “modified harmonic” coordinate system, as opposed to the SH (“standard harmonic”) coefficients $A_{3PN}$ and $B_{3PN}$ containing logarithms. Limiting ourselves to 3PN order, we have (Mora and Will 2004; Arun et al. 2008a)

$$Q_{ab}^{(3)} = \frac{2Gm^2v}{r^2} \left[ -4n^{(a}v^{b)} + 3\dot{r}n^{(a}n^{b)} \right], \\ \left(Q_{ab}^{(3)}\right)^2 = \frac{8G^2m^4v^2}{r^4} \left( 4v^2 - \frac{11}{3}\dot{r}^2 \right).$$

<footnote n="72">⁷² We give for reference [in particular the last term in (360) tends to zero when $v \to 0$]</footnote>
$$
\begin{equation}
\begin{aligned}
A_{\text{3PN}}^{\text{MH}} = & -\frac{35r^6 v}{16} + \frac{175r^6 v^2}{16} - \frac{175r^6 v^3}{16} + \frac{15r^4 v v^2}{2} - \frac{135r^4 v^2 v^2}{4} + \frac{255r^4 v^3 v^2}{8} \\
& -\frac{15r^2 v v^4}{2} + \frac{237r^2 v^2 v^4}{8} - \frac{45r^2 v^3 v^4}{2} + \frac{11v v^6}{4} - \frac{49v^2 v^6}{4} + 13v^3 v^6 \\
& + \frac{Gm}{r} \left( 79r^4 v - \frac{69r^4 v^2}{2} - 30r^4 v^3 - 121r^2 v v^2 + 16r^2 v^2 v^2 + 20r^2 v^3 v^2 \right. \\
& \left. + \frac{75v v^4}{4} + 8v^2 v^4 - 10v^3 v^4 \right) \\
& + \frac{G^2 m^2}{r^2} \left( r^2 + \frac{22717r^2 v}{168} + \frac{11r^2 v^2}{8} - 7r^2 v^3 + \frac{615r^2 v \pi^2}{64} - \frac{20827v v^2}{840} \right. \\
& \left. + v^3 v^2 - \frac{123v \pi^2 v^2}{64} \right) \\
& + \frac{G^3 m^3}{r^3} \left( -16 - \frac{1399v}{12} - \frac{71v^2}{2} + \frac{41v \pi^2}{16} \right),
\end{aligned}
\tag{362a}
\end{equation}
$$

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

In the following we shall use the MH coordinate system when constructing the generalized quasi-Keplerian representation of the 3PN motion in Sect. 3.5.2. We now come back to the SH coordinate system.

The conservative part of the CM equations of motion derive from a Lagrangian, which can be constructed as follows. We start from the general-frame Lagrangian (343)–(350), which is a functional of $\mathbf{y}_a$, $\mathbf{v}_a$ and $\mathbf{a}_a$, and admits the CM integral $G[\mathbf{y}_a, \mathbf{v}_a]$ given at 4PN order by Eqs. (B1)–(B3) of Bernard et al. (2018). We perform the change of variables $(\mathbf{y}_1, \mathbf{y}_2) \rightarrow (\mathbf{x}, \mathbf{G})$, where $\mathbf{x} = \mathbf{y}_1 - \mathbf{y}_2$. Since $\mathbf{G} = m_1 \mathbf{y}_1 + m_2 \mathbf{y}_2 + \mathcal{O}(c^{-2})$ to Newtonian order, we find $\mathbf{y}_1 = X_2 \mathbf{x} + \frac{1}{m} \mathbf{G} + \mathcal{O}(c^{-2})$ and $\mathbf{v}_1 =$
$X_2\mathbf{v} + \frac{1}{m}\frac{d\mathbf{G}}{dt} + \mathcal{O}(c^{-2})$ (and $1 \leftrightarrow 2$). Proceeding iteratively, it is easy to see that the old variables $\mathbf{y}_a$ are obtained as functionals of the new variables $(\mathbf{x}, \mathbf{G})$ and their derivatives, up to some differentiation order depending on the looked for PN order. In the process, we do not perform any order reduction of accelerations. Thus we get

$$ \mathbf{y}_a = \mathbf{y}_a[\mathbf{x}, \mathbf{v}, \mathbf{a}, \dots; \mathbf{G}, \frac{d\mathbf{G}}{dt}, \frac{d^2\mathbf{G}}{dt^2}, \dots], \quad (363) $$

which we plug into the original general-frame Lagrangian, still without order reduction of the accelerations. This yields an equivalent Lagrangian, which is now of the “doubly generalized” type

$$ L = L[\mathbf{x}, \mathbf{v}, \mathbf{a}, \dots; \mathbf{G}, \frac{d\mathbf{G}}{dt}, \frac{d^2\mathbf{G}}{dt^2}, \dots]. \quad (364) $$

The ensuing equations of motion read $\frac{\partial L}{\partial \mathbf{x}} = 0$ as usual, together with $\frac{\partial L}{\partial \mathbf{G}} = 0$. Now the latter equation is necessarily equivalent to the conservation law for the CM integral, hence we have

$$ \frac{\partial L}{\partial \mathbf{G}} = 0 \quad \Longleftrightarrow \quad \frac{d^2 \mathbf{G}}{dt^2} = 0. \quad (365) $$

As a result, we can choose $\mathbf{G} = 0$ as a solution of these equations. The CM equations of motion are then given by $\frac{\partial L}{\partial \mathbf{x}} = 0$ in which we pose $\mathbf{G} = 0$ [these equations are nothing but the CM equations of motion (357) with (358)–(359)] and the CM Lagrangian is obtained by setting $\mathbf{G} = 0$ in Eq. (364). Defining as usual the reduced CM Lagrangian as $\mathcal{L} = L/\mu$, we explicitly get (Blanchet and Iyer 2003; Bernard et al. 2018)

$$ \mathcal{L}_N = \frac{v^2}{2} + \frac{Gm}{r}, \quad (366a) $$

$$ \mathcal{L}_{1PN} = \frac{v^4}{8} - \frac{3 v v^4}{8} + \frac{Gm}{r} \left( \frac{\dot{r}^2 v}{2} + \frac{3 v^2}{2} + \frac{v v^2}{2} \right) - \frac{G^2 m^2}{2 r^2}, \quad (366b) $$

$$ \begin{aligned} \mathcal{L}_{2PN} = & \frac{v^6}{16} - \frac{7 v v^6}{16} + \frac{13 v^2 v^6}{16} \\ & + \frac{Gm}{r} \left( \frac{3 \dot{r}^4 v^2}{8} - \frac{\dot{r}^2 (an) v r}{8} + \frac{\dot{r}^2 v v^2}{4} - \frac{5 \dot{r}^2 v^2 v^2}{4} + \frac{7 (an) v r v^2}{8} \right. \\ & \left. \qquad + \frac{7 v^4}{8} - \frac{5 v v^4}{4} - \frac{9 v^2 v^4}{8} - \frac{7 \dot{r} v r (av)}{4} \right) \\ & + \frac{G^2 m^2}{r^2} \left( \frac{\dot{r}^2}{2} + \frac{41 \dot{r}^2 v}{8} + \frac{3 \dot{r}^2 v^2}{2} + \frac{7 v^2}{4} - \frac{27 v v^2}{8} + \frac{v^2 v^2}{2} \right) \\ & + \frac{G^3 m^3}{r^3} \left( \frac{1}{2} + \frac{15 v}{4} \right), \end{aligned} \quad (366c) $$
$$
\begin{equation}
\begin{aligned}
\mathcal{L}_{3\text{PN}} = & \frac{5 v^8}{128} - \frac{59 v v^8}{128} + \frac{119 v^2 v^8}{64} - \frac{323 v^3 v^8}{128} \\
& + \frac{Gm}{r} \left( \frac{5 \dot{r}^6 v^3}{16} + \frac{\dot{r}^4 (an) v r}{16} - \frac{5 \dot{r}^4 (an) v^2 r}{16} - \frac{3 \dot{r}^4 v v^2}{16} \right. \\
& \left. + \frac{7 \dot{r}^4 v^2 v^2}{4} - \frac{33 \dot{r}^4 v^3 v^2}{16} - \frac{3 \dot{r}^2 (an) v r v^2}{16} - \frac{\dot{r}^2 (an) v^2 r v^2}{16} \right. \\
& \left. + \frac{5 \dot{r}^2 v v^4}{8} - 3 \dot{r}^2 v^2 v^4 + \frac{75 \dot{r}^2 v^3 v^4}{16} + \frac{7 (an) v r v^4}{8} \right. \\
& - \frac{7 (an) v^2 r v^4}{2} + \frac{11 v^6}{16} - \frac{55 v v^6}{16} + \frac{5 v^2 v^6}{2} \\
& + \frac{65 v^3 v^6}{16} + \frac{5 \dot{r}^3 v r (av)}{12} - \frac{13 \dot{r}^3 v^2 r (av)}{8} \\
& - \frac{37 \dot{r} v r v^2 (av)}{8} + \frac{35 \dot{r} v^2 r v^2 (av)}{4} \right) \\
& + \frac{G^2 m^2}{r^2} \left( - \frac{109 \dot{r}^4 v}{144} - \frac{259 \dot{r}^4 v^2}{36} + 2 \dot{r}^4 v^3 - \frac{17 \dot{r}^2 (an) v r}{6} \right. \\
& \left. + \frac{97 \dot{r}^2 (an) v^2 r}{12} + \frac{\dot{r}^2 v^2}{4} - \frac{41 \dot{r}^2 v v^2}{6} - \frac{2287 \dot{r}^2 v^2 v^2}{48} \right. \\
& \left. - \frac{27 \dot{r}^2 v^3 v^2}{4} + \frac{203 (an) v r v^2}{12} + \frac{149 (an) v^2 r v^2}{6} \right. \\
& \left. + \frac{45 v^4}{16} + \frac{53 v v^4}{24} + \frac{617 v^2 v^4}{24} - \frac{9 v^3 v^4}{4} \right. \\
& \left. - \frac{235 \dot{r} v r (av)}{24} + \frac{235 \dot{r} v^2 r (av)}{6} \right) \\
& + \frac{G^3 m^3}{r^3} \left( \frac{3 \dot{r}^2}{2} - \frac{12041 \dot{r}^2 v}{420} + \frac{37 \dot{r}^2 v^2}{4} + \frac{7 \dot{r}^2 v^3}{2} - \frac{123 \dot{r}^2 v \pi^2}{64} \right. \\
& \left. + \frac{5 v^2}{4} + \frac{387 v v^2}{70} - \frac{7 v^2 v^2}{4} + \frac{v^3 v^2}{2} + \frac{41 v \pi^2 v^2}{64} \right. \\
& \left. + 22 \dot{r}^2 v \ln\left(\frac{r}{r_0}\right) - \frac{22 v v^2}{3} \ln\left(\frac{r}{r_0}\right) \right) \\
& + \frac{G^4 m^4}{r^4} \left( - \frac{3}{8} - \frac{18469 v}{840} + \frac{22 v}{3} \ln\left(\frac{r}{r_0}\right) \right),
\end{aligned}
\tag{366d}
\end{equation}
$$

$$
\mathcal{L}_{4\text{PN}}^{(0)} = \frac{7}{256} v^{10} - \frac{121}{256} v v^{10} + \frac{785}{256} v^2 v^{10} - \frac{1127}{128} v^3 v^{10} + \frac{2415}{256} v^4 v^{10}, \quad (366e)
$$
$$
\begin{equation}
\begin{aligned}
\mathcal{L}_{\text{4PN}}^{(1)} = & \frac{m}{r} \left( \frac{23}{20} v^2 (av) r \dot{r}^5 - \frac{5}{128} v \dot{r}^8 + \frac{35}{128} v^2 \dot{r}^8 - \frac{35}{64} v^3 \dot{r}^8 + \frac{35}{128} v^4 \dot{r}^8 \right. \\
& \left. + \frac{7}{4} v (av) r \dot{r}^3 v^2 + \frac{361}{24} v^3 (av) r \dot{r}^3 v^2 + \frac{19}{32} v (an) r \dot{r}^4 v^2 \right. \\
& \left. + \frac{85}{32} v^3 (an) r \dot{r}^4 v^2 - \frac{5}{16} v \dot{r}^6 v^2 - \frac{31}{16} v^2 \dot{r}^6 v^2 + \frac{45}{32} v^3 \dot{r}^6 v^2 \right. \\
& \left. - \frac{85}{32} v^4 \dot{r}^6 v^2 + \frac{341}{4} v^2 (av) r \dot{r} v^4 + \frac{245}{32} v^2 (an) r \dot{r}^2 v^4 - \frac{17}{64} v \dot{r}^4 v^4 \right. \\
& \left. - \frac{11}{8} v^2 r^4 v^4 - \frac{193}{16} v^3 r^4 v^4 + \frac{693}{64} v^4 r^4 v^4 + \frac{217}{96} v (an) r v^6 \right. \\
& \left. + \frac{2261}{96} v^3 (an) r v^6 - \frac{11}{48} v \dot{r}^2 v^6 - \frac{29}{4} v^2 \dot{r}^2 v^6 + \frac{1021}{96} v^3 \dot{r}^2 v^6 \right. \\
& \left. - \frac{665}{32} v^4 \dot{r}^2 v^6 + \frac{75}{128} v^8 - \frac{1595}{384} v v^8 + \frac{2917}{128} v^2 v^8 + \frac{493}{192} v^3 v^8 \right. \\
& \left. - \frac{2261}{128} v^4 v^8 \right),
\end{aligned}
\tag{366f}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
\mathcal{L}_{\text{4PN}}^{(2)} = & \frac{m^2}{r^2} \left( \frac{5407}{288} v(av) r \dot{r}^3 - \frac{5531}{3200} v \dot{r}^6 - \frac{487}{40} v^2 \dot{r}^6 + \frac{13}{5} v^3 \dot{r}^6 + \frac{3}{2} v^4 \dot{r}^6 \right. \\
& \left. + \frac{11497}{48} v^2 (av) r \dot{r} v^2 + \frac{469}{4} v^3 (av) r \dot{r} v^2 + \frac{517}{48} v^2 (an) r \dot{r}^2 v^2 \right. \\
& \left. - \frac{61183}{5760} v \dot{r}^4 v^2 + \frac{3079}{96} v^2 \dot{r}^4 v^2 - \frac{161}{8} v^3 \dot{r}^4 v^2 - \frac{27}{2} v^4 \dot{r}^4 v^2 \right. \\
& \left. + \frac{14627}{384} v (an) r v^4 + \frac{3}{16} \dot{r}^2 v^4 - \frac{12091}{640} v \dot{r}^2 v^4 - \frac{22045}{192} v^2 \dot{r}^2 v^4 \right. \\
& \left. + \frac{255}{64} v^3 \dot{r}^2 v^4 + \frac{569}{16} v^4 \dot{r}^2 v^4 + \frac{115}{32} v^6 + \frac{593}{120} v v^6 + \frac{9467}{192} v^2 v^6 \right. \\
& \left. + \frac{599}{64} v^3 v^6 + \frac{195}{16} v^4 v^6 \right),
\end{aligned}
\tag{366g}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
\mathcal{L}_{\text{4PN}}^{(3)} = & \frac{m^3}{r^3} \left( \frac{4937}{1260} v^2 (av) r \dot{r} - \frac{41}{32} \pi^2 v^2 (av) r \dot{r} + \frac{44}{3} v^2 (av) r \ln\left(\frac{r}{r_0'}\right) \dot{r} \right. \\
& \left. + \frac{22}{3} v (av) r \ln\left(\frac{r}{r_0'}\right) \dot{r} - \frac{246373}{2240} v \dot{r}^4 - \frac{2155}{1024} \pi^2 v \dot{r}^4 - \frac{210733}{2016} v^2 \dot{r}^4 \right. \\
& \left. + \frac{205}{128} \pi^2 v^2 \dot{r}^4 + \frac{367}{32} v^3 \dot{r}^4 + \frac{29}{4} v^4 \dot{r}^4 - \frac{55}{3} v \ln\left(\frac{r}{r_0'}\right) \dot{r}^4 \right. \\
& \left. + 55 v^2 \ln\left(\frac{r}{r_0'}\right) \dot{r}^4 + \frac{229319}{6300} v (an) r v^2 - \frac{21}{32} \pi^2 v (an) r v^2 \right. \\
& \left. + \frac{49}{4} v^3 (an) r v^2 + 44 v (an) r \ln\left(\frac{r}{r_0'}\right) v^2 + \frac{44}{3} v^2 (an) r \ln\left(\frac{r}{r_0'}\right) v^2 \right)
\end{aligned}
\tag{366h}
\end{equation}
$$
$$ 
\begin{aligned}
& + \frac{7}{4} \dot{r}^2 v^2 + \frac{516319}{4200} v \dot{r}^2 v^2 + \frac{447}{512} \pi^2 v \dot{r}^2 v^2 + \frac{53099}{560} v^2 \dot{r}^2 v^2 \\
& + \frac{123}{64} \pi^2 v^2 \dot{r}^2 v^2 - \frac{1003}{16} v^3 \dot{r}^2 v^2 - \frac{47}{2} v^4 \dot{r}^2 v^2 - 44 v \ln \left( \frac{r}{r_0'} \right) \dot{r}^2 v^2 \\
& - 110 v^2 \ln \left( \frac{r}{r_0'} \right) \dot{r}^2 v^2 + \frac{91}{16} v^4 - \frac{166703}{20160} v v^4 + \frac{133}{1024} \pi^2 v v^4 \\
& + \frac{10601}{3360} v^2 v^4 - \frac{123}{128} \pi^2 v^2 v^4 + \frac{567}{32} v^3 v^4 - \frac{15}{4} v^4 v^4 \\
& + \frac{55}{3} v \ln \left( \frac{r}{r_0'} \right) v^4 + \frac{77}{3} v^2 \ln \left( \frac{r}{r_0'} \right) v^4,
\end{aligned} 
\quad (366h) $$

$$ 
\begin{aligned}
\mathcal{L}_{\text{4PN}}^{(4)} = & \frac{m^4}{r^4} \left( \frac{9}{4} \dot{r}^2 - \frac{245971}{4200} v \dot{r}^2 + \frac{2771}{96} \pi^2 v \dot{r}^2 - \frac{8089}{140} v^2 \dot{r}^2 - 44 \pi^2 v^2 \dot{r}^2 \right. \\
& \left. + \frac{185}{8} v^3 \dot{r}^2 + \frac{15}{2} v^4 \dot{r}^2 - 339 v \ln \left( \frac{r}{r_0'} \right) \dot{r}^2 + \frac{1616}{3} v^2 \ln \left( \frac{r}{r_0'} \right) \dot{r}^2 \right. \\
& \left. + \frac{15}{16} v^2 + \frac{2039993}{50400} v v^2 - \frac{191}{32} \pi^2 v v^2 - \frac{52907}{1008} v^2 v^2 + 11 \pi^2 v^2 v^2 \right. \\
& \left. - \frac{1}{8} v^3 v^2 + \frac{1}{2} v^4 v^2 + \frac{301}{3} v \ln \left( \frac{r}{r_0'} \right) v^2 - \frac{371}{3} v^2 \ln \left( \frac{r}{r_0'} \right) v^2 \right),
\end{aligned} 
\quad (366i) $$

$$ 
\begin{aligned}
\mathcal{L}_{\text{4PN}}^{(5)} = & \frac{m^5}{r^5} \left( \frac{3}{8} + \frac{1697177}{25200} v + \frac{105}{32} \pi^2 v + \frac{55111}{720} v^2 - 11 \pi^2 v^2 \right. \\
& \left. - \frac{290}{3} v \ln \left( \frac{r}{r_0'} \right) + 120 v^2 \ln \left( \frac{r}{r_0'} \right) \right).
\end{aligned} 
\quad (366j) $$

Remind that one has still to add the conservative 4PN tail term given by

$$ \mathcal{L}_{\text{tail}} = \frac{G^2}{5c^8 v} Q_{ab}^{(3)}(t) \text{Pf}_{2r/c} \int_{-\infty}^{+\infty} \frac{dt'}{|t-t'|} Q_{ab}^{(3)}(t'), \quad (367) $$

whose functional variation yields Eq. (361a).

The above CM Lagrangian in harmonic coordinates depends on accelerations starting at 2PN order, through $(an) \equiv \mathbf{a} \cdot \mathbf{n}$ and $(av) \equiv \mathbf{a} \cdot \mathbf{v}$. Furthermore, it depends on the logarithm of $r/r_0'$. Both the accelerations and the logarithms can be removed by applying a contact transformation, and one can prove that the harmonic-coordinates CM Lagrangian is equivalent to the ADM-coordinates CM Hamiltonian (which has no accelerations nor logarithms). In the CM frame the conjugate variables are the relative separation $X = RN = Y_1 - Y_2$ and the conjugate momentum (per
unit reduced mass) **P** such that $\mu \mathbf{P} = \mathbf{P}_1 = -\mathbf{P}_2$ ($\mathbf{P}_1$ and $\mathbf{P}_2$ are defined in Sect. 3.2.3). Posing $P^2 \equiv \mathbf{P}^2$ and $P_R \equiv N \cdot \mathbf{P}$ and $\mathcal{H} = H/\mu$, we have (Damour et al. 2014, 2015; Schäfer and Jaranowski 2024; Blümlein et al. 2020a)

$$ \mathcal{H}_N = \frac{P^2}{2} - \frac{Gm}{R}, \qquad (368a) $$

$$ \mathcal{H}_{1PN} = \left(-\frac{1}{8} + \frac{3v}{8}\right)P^4 + \frac{Gm}{R}\left(-\frac{P_R^2 v}{2} - \frac{3P^2}{2} - \frac{vP^2}{2}\right) + \frac{G^2m^2}{2R^2}, \qquad (368b) $$

$$ \begin{aligned} \mathcal{H}_{2PN} = & \left(\frac{1}{16} - \frac{5v}{16} + \frac{5v^2}{16}\right)P^6 \\ & + \frac{Gm}{R}\left(-\frac{3P_R^4 v^2}{8} - \frac{P_R^2 P^2 v^2}{4} + \frac{5P^4}{8} - \frac{5vP^4}{2} - \frac{3v^2P^4}{8}\right) \\ & + \frac{G^2m^2}{R^2}\left(\frac{3P_R^2 v}{2} + \frac{5P^2}{2} + 4vP^2\right) \\ & + \frac{G^3m^3}{R^3}\left(-\frac{1}{4} - \frac{3v}{4}\right), \end{aligned} \qquad (368c) $$

$$ \begin{aligned} \mathcal{H}_{3PN} = & \left(-\frac{5}{128} + \frac{35v}{128} - \frac{35v^2}{64} + \frac{35v^3}{128}\right)P^8 \\ & + \frac{Gm}{R}\left(-\frac{5P_R^6 v^3}{16} + \frac{3P_R^4 P^2 v^2}{16} - \frac{3P_R^4 P^2 v^3}{16} + \frac{P_R^2 P^4 v^2}{8} - \frac{3P_R^2 P^4 v^3}{16} - \frac{7P^6}{16} + \frac{21vP^6}{8} - \frac{53v^2P^6}{16} - \frac{5v^3P^6}{16}\right) \\ & + \frac{G^2m^2}{R^2}\left(\frac{5P_R^4 v}{12} + \frac{43P_R^4 v^2}{12} + \frac{17P_R^2 P^2 v}{16} + \frac{15P_R^2 P^2 v^2}{8} - \frac{27P^4}{16} + \frac{17vP^4}{2} + \frac{109v^2 P^4}{16}\right) \\ & + \frac{G^3m^3}{R^3}\left(-\frac{85P_R^2 v}{16} - \frac{7P_R^2 v^2}{4} - \frac{25P^2}{8} - \frac{335vP^2}{48} - \frac{23v^2 P^2}{8} - \frac{3P_R^2 v \pi^2}{64} + \frac{vP^2 \pi^2}{64}\right) \\ & + \frac{G^4m^4}{R^4}\left(\frac{1}{8} + \frac{109v}{12} - \frac{21v \pi^2}{32}\right), \end{aligned} \qquad (368d) $$
$$H_{4PN}^{(0)} = \left( \frac{7}{256} - \frac{63}{256} v + \frac{189}{256} v^2 - \frac{105}{128} v^3 + \frac{63}{256} v^4 \right) P^{10}, \quad (368e)$$

$$H_{4PN}^{(1)} = \frac{m}{R} \left\{ \begin{aligned} & \frac{45}{128} P^8 - \frac{45}{16} P^8 v + \left( \frac{423}{64} P^8 - \frac{3}{32} P_R^2 P^6 - \frac{9}{64} P_R^4 P^4 \right) v^2 \\ & + \left( -\frac{1013}{256} P^8 + \frac{23}{64} P_R^2 P^6 + \frac{69}{128} P_R^4 P^4 - \frac{5}{64} P_R^6 P^2 + \frac{35}{256} P_R^8 \right) v^3 \\ & + \left( -\frac{35}{128} P^8 - \frac{5}{32} P_R^2 P^6 - \frac{9}{64} P_R^4 P^4 - \frac{5}{32} P_R^6 P^2 - \frac{35}{128} P_R^8 \right) v^4 \end{aligned} \right\}, \quad (368f)$$

$$H_{4PN}^{(2)} = \frac{m^2}{R^2} \left\{ \begin{aligned} & \frac{13}{8} P^6 + \left( -\frac{791}{64} P^6 + \frac{49}{16} P_R^2 P^4 - \frac{889}{192} P_R^4 P^2 + \frac{369}{160} P_R^6 \right) v \\ & + \left( \frac{4857}{256} P^6 - \frac{545}{64} P_R^2 P^4 + \frac{9475}{768} P_R^4 P^2 - \frac{1151}{128} P_R^6 \right) v^2 \\ & + \left( \frac{2335}{256} P^6 + \frac{1135}{256} P_R^2 P^4 - \frac{1649}{768} P_R^4 P^2 + \frac{10353}{1280} P_R^6 \right) v^3 \end{aligned} \right\}, \quad (368g)$$

$$H_{4PN}^{(3)} = \frac{m^3}{R^3} \left\{ \begin{aligned} & \frac{105}{32} P^4 \\ & + \left[ \left( \frac{2749}{8192} \pi^2 - \frac{589189}{19200} \right) P^4 + \left( \frac{63347}{1600} - \frac{1059}{1024} \pi^2 \right) P_R^2 P^2 \\ & + \left( \frac{375}{8192} \pi^2 - \frac{23533}{1280} \right) P_R^4 \right] v \\ & + \left[ \left( \frac{18491}{16384} \pi^2 - \frac{1189789}{28800} \right) P^4 \right. \\ & \left. - \left( \frac{127}{3} + \frac{4035}{2048} \pi^2 \right) P_R^2 P^2 + \left( \frac{57563}{1920} - \frac{38655}{16384} \pi^2 \right) P_R^4 \right] v^2 \\ & \left. + \left( -\frac{553}{128} P^4 - \frac{225}{64} P_R^2 P^2 - \frac{381}{128} P_R^4 \right) v^3 \end{aligned} \right\}, \quad (368h)$$

$$H_{4PN}^{(4)} = \frac{m^4}{R^4} \left\{ \frac{105}{32} P^2 \right\} \quad (368i)$$

$$+ \left[ \left( \frac{185761}{19200} - \frac{21837}{8192} \pi^2 \right) P^2 + \left( \frac{3401779}{57600} - \frac{28691}{24576} \pi^2 \right) P_R^2 \right] v \\ + \left[ \left( \frac{672811}{19200} - \frac{158177}{49152} \pi^2 \right) P^2 + \left( -\frac{21827}{3840} + \frac{110099}{49152} \pi^2 \right) P_R^2 \right] v^2 \right\}, \quad (368j)$$

$$H_{4PN}^{(5)} = \frac{m^5}{R^5} \left\{ -\frac{1}{16} + \left( -\frac{169199}{2400} + \frac{6237}{1024} \pi^2 \right) v + \left( -\frac{1256}{45} + \frac{7403}{3072} \pi^2 \right) v^2 \right\},$$
while to this order the tail part is the opposite of that in the Lagrangian: $H_{\text{tail}} = -L_{\text{tail}}$, where $L_{\text{tail}}$ is given by (367) with $r \to R$.

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

Most inspiralling compact binaries will have been circularized by the time they become visible by the detectors LIGO and Virgo; see Sect. 1.1.2 for a discussion. In the case of orbits that are circular—apart from the gradual radiation-reaction inspiral—the complicated equations of motion simplify drastically, since we have $\dot{r} = (nv) = O(1/c^5)$ and, for instance, $\ddot{r} = O(1/c^{10})$. In the case of (quasi-)circular orbits it is convenient to display the successive post-Newtonian corrections employing the useful post-Newtonian parameter

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

The relation between center-of-mass variables and the relative ones still reads (354) but with the simplified coefficients

$$ \begin{aligned} \mathcal{P} = \frac{G^2 m^2}{r^2 c^4} & \left\{ 3 + \left( -\frac{7211}{1260} + \frac{22}{3} \ln\left(\frac{r}{r_0'}\right) - v \right) \gamma \right. \\ & \left. + \left( \frac{31669}{600} - \frac{154}{3} \ln\left(\frac{r}{r_0'}\right) \right) \right. \\ & \left. + \left[ -\frac{13}{180} - \frac{135}{32} \pi^2 + 60 \ln\left(\frac{r}{r_0'}\right) \right] v \right) \gamma^2 \right\} + O\left(\frac{1}{c^9}\right), \end{aligned} \qquad (370a) $$

$$ \mathcal{Q} = \frac{G^2 m^2}{rc^5} \left\{ -\frac{4}{5} + \left( -\frac{316}{105} - \frac{2}{5}v \right) \gamma \right\} + O\left(\frac{1}{c^9}\right). \qquad (370b) $$

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

$$ \mathbf{a} = -\Omega^2 \mathbf{x} - \frac{32}{5} \frac{G^3 m^3 v}{c^5 r^4} \left[ 1 + \gamma \left( -\frac{743}{336} - \frac{11}{4}v \right) + 4\pi\gamma^{3/2} \right] \mathbf{v} + O\left(\frac{1}{c^9}\right), \qquad (371) $$

where $\mathbf{x} \equiv \mathbf{y}_1 - \mathbf{y}_2$ is the relative separation (in harmonic coordinates) and $\Omega$ denotes the angular frequency of the quasi-circular motion. The second term in Eq. (371), opposite to the velocity $\mathbf{v} \equiv \mathbf{v}_1 - \mathbf{v}_2$, represents the radiation reaction force up to 4PN order, which comes from the 2.5PN and 3.5PN coefficients in Eqs. (358)–(359), plus the dissipative part of the 4PN tail contribution, given by Eq. (361b) and which produces the $4\pi$ term in (371). The calculation of the next 4.5PN dissipative
coefficient has been carried out using the EFT by Leibovich et al. (2023). The results at 2.5PN, 3.5PN and 4.5PN orders are in agreement with the flux-balance method (Iyer and Will 1993, 1995; Gopakumar et al. 1997). The radiation-reaction force is responsible for the secular decrease of the separation $r$ and increase of the orbital frequency $\Omega$:

$$ \dot{r} = -\frac{64}{5} \frac{G^3 m^3 v}{r^3 c^5} \left[ 1 + \gamma \left( -\frac{1751}{336} - \frac{7}{4}v \right) + 4\pi\gamma^{3/2} + \mathcal{O}\left(\frac{1}{c^4}\right) \right], \quad (372a) $$

$$ \dot{\Omega} = \frac{96}{5} \frac{Gmv}{r^3} \gamma^{5/2} \left[ 1 + \gamma \left( -\frac{2591}{336} - \frac{11}{12}v \right) + 4\pi\gamma^{3/2} + \mathcal{O}\left(\frac{1}{c^4}\right) \right], \quad (372b) $$

which can be recovered and extended to very high order (namely 2.5PN + 4.5PN) using the phase and frequency evolution described in Sect. 3.4.4.

Concerning conservative effects, the main content of the 4PN equations of motion (371) is the relation between the frequency $\Omega$ and the orbital separation $r$ or equivalently the parameter (369). This is given by the following generalized version of “Kepler’s third law”$^{73}$:

$$ \begin{aligned} \Omega^2 = \frac{Gm}{r^3} & \left\{ 1 + (-3 + v)\gamma + \left(6 + \frac{41}{4}v + v^2\right)\gamma^2 \right. \\ & \left. + \left(-10 + \left[-\frac{75707}{840} + \frac{41}{64}\pi^2 + 22\ln\left(\frac{r}{r_0'}\right)\right]v + \frac{19}{2}v^2 + v^3\right)\gamma^3 \right. \\ & \left. + \left(15 + \left[\frac{19644217}{33600} + \frac{163}{1024}\pi^2 + \frac{256}{5}\gamma_E + \frac{128}{5}\ln(16\gamma) - 290\ln\left(\frac{r}{r_0'}\right)\right]v \right. \right. \\ & \left. \left. + \left[\frac{44329}{336} - \frac{1907}{64}\pi^2 + \frac{1168}{3}\ln\left(\frac{r}{r_0'}\right)\right]v^2 + \frac{51}{4}v^3 + v^4\right)\gamma^4 \right\} \\ & + \mathcal{O}\left(\frac{1}{c^{10}}\right), \end{aligned} \quad (373) $$

where $r_0'$ is the UV type length scale in Eqs. (358)–(359), and $\gamma_E$ is Euler’s constant. As for the circular energy, it is inferred from the CM and then circular-orbit reduction of the general frame expression. We have

$$ \left( \frac{Pf}{2r/c} \int_{-\infty}^{+\infty} \frac{dt'}{|t-t'|} Q_{ij}^{(n)}(t') \right) \Big|_{\text{circ}} = -2Q_{ij}^{(n)} \left[ \ln\left(\frac{2r\Omega}{c}\right) + \gamma_E \right]. $$

<footnote>$^{73}$ In the limit of quasi-circular orbits (in the adiabatic approximation) the tail integral in the equation of motion becomes local and we have</footnote>
$$
\begin{equation}
\begin{aligned}
E = -\frac{\mu c^2 \gamma}{2} & \left\{ 1 + \left(-\frac{7}{4} + \frac{v}{4}\right)\gamma + \left(-\frac{7}{8} + \frac{49}{8}v + \frac{v^2}{8}\right)\gamma^2 \right. \\
& \left. + \left(-\frac{235}{64} + \left[\frac{46031}{2240} - \frac{123}{64}\pi^2 + \frac{22}{3}\ln\left(\frac{r}{r_0}\right)\right]v + \frac{27}{32}v^2 + \frac{5}{64}v^3\right)\gamma^3 \right. \\
& \left. + \left(-\frac{649}{128} + \left[-\frac{3357833}{28800} + \frac{384}{5}\gamma_E + \frac{192}{5}\ln(16\gamma) + \frac{14935}{1024}\pi^2 - 93\ln\left(\frac{r}{r_0}\right)\right]v \right. \right. \\
& \left. \left. + \left[\frac{83959}{8064} - \frac{957}{128}\pi^2 + \frac{371}{3}\ln\left(\frac{r}{r_0}\right)\right]v^2 + \frac{69}{64}v^3 + \frac{7}{128}v^4\right)\gamma^4 \right\} \\
& + \mathcal{O}\left(\frac{1}{c^{10}}\right).
\end{aligned}
\tag{374}
\end{equation}
$$

This expression is that of a physical observable E; however, it depends on the choice
of a coordinate system, as it involves the post-Newtonian parameter $\gamma$ defined from
the harmonic-coordinate separation $r$. But the *numerical* value of E should not
depend on the choice of a coordinate system, so E must admit a frame-invariant
expression, the same in all coordinate systems. To find it we re-express E with the
help of the following frequency-related parameter $x$, instead of the post-Newtonian
parameter $\gamma$ (with $m = m_1 + m_2)^{74}$:

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

We readily obtain from Eq. (373) the expression of $\gamma$ in terms of $x$ as

$$
\begin{equation}
\begin{aligned}
\gamma = x & \left\{ 1 + \left(1 - \frac{v}{3}\right)x + \left(1 - \frac{65}{12}v\right)x^2 \right. \\
& \qquad + \left(1 + \left[-\frac{2203}{2520} - \frac{41}{192}\pi^2 + \frac{22}{3}\ln\left(\frac{xr'_0c^2}{Gm}\right)\right]v + \frac{229}{36}v^2 + \frac{v^3}{81}x^3 \right. \\
& \qquad \left. + \left(1 + \left[-\frac{2067859}{33600} - \frac{256}{15}\gamma_E - \frac{5411}{3072}\pi^2 - 38\ln\left(\frac{xr'_0c^2}{Gm}\right) - \frac{128}{15}\ln(16x)\right]v \right. \right. \\
& \qquad \left. \left. + \left[\frac{153613}{15120} + \frac{6049}{576}\pi^2 + \frac{992}{9}\ln\left(\frac{xr'_0c^2}{Gm}\right)\right]v^2 - \frac{1261}{324}v^3 + \frac{v^4}{243}x^4 \right. \right. \\
& \qquad \left. \left. + \mathcal{O}\left(\frac{1}{c^{10}}\right) \right\},
\end{aligned}
\tag{376}
\end{equation}
$$

<footnote>$^{74}$ This parameter is an invariant in a large class of coordinate systems, namely all those for which the metric becomes asymptotically Minkowskian far away from the matter source, i.e., $g_{z\beta} \to \text{diag}(-1, 1, 1, 1)$.</footnote>
that we substitute back into Eq. (374), making all appropriate post-Newtonian re-
expansions. As a result, we gladly discover that the gauge constant $r'_0$ has cancelled
out. Therefore, the gauge-invariant result is (Damour et al. 2000b; Blanchet and Faye
2000a; Damour et al. 2001a; Blanchet et al. 2004a; Le Tiec et al. 2012a, b; Bini and
Damour 2013; Damour et al. 2014; Bernard et al. 2017a, b; Foffa and Sturani 2012;
Foffa et al. 2019)<fnref n="75" />

$$
\begin{equation}
\begin{aligned}
E = & -\frac{\mu c^2 x}{2} \left\{ 1 + \left(-\frac{3}{4} - \frac{v}{12}\right)x + \left(-\frac{27}{8} + \frac{19}{8}v - \frac{v^2}{24}\right)x^2 \right. \\
& \qquad \left. + \left(-\frac{675}{64} + \left[\frac{34445}{576} - \frac{205}{96}\pi^2\right]v - \frac{155}{96}v^2 - \frac{35}{5184}v^3\right)x^3 \right. \\
& \qquad \left. + \left(-\frac{3969}{128} + \left[-\frac{123671}{5760} + \frac{9037}{1536}\pi^2 + \frac{896}{15}\gamma_E + \frac{448}{15}\ln(16x)\right]v \right. \right. \\
& \qquad \left. \left. + \left[-\frac{498449}{3456} + \frac{3157}{576}\pi^2\right]v^2 + \frac{301}{1728}v^3 + \frac{77}{31104}v^4 \right)x^4 \right. \\
& \qquad \left. + \mathcal{O}\left(\frac{1}{c^{10}}\right) \right\}.
\end{aligned}
\tag{377}
\end{equation}
$$

For quasi-circular orbits one can check that there are no terms of order $\propto x^{9/2}$ in (377), so the result is actually valid up to the 4.5PN order. We shall discuss in Sect. 3.6 how the effects of the spins of the two black holes affect the latter formula. We can write also a similar expression for the angular momentum,

$$
\begin{equation}
\begin{aligned}
J ={}& \frac{G\mu m}{cx^{1/2}} \left\{ 1 + \left(\frac{3}{2} + \frac{v}{6}\right)x + \left(\frac{27}{8} - \frac{19}{8}v + \frac{v^2}{24}\right)x^2 \right. \\
& \qquad \left. + \left(\frac{135}{16} + \left[-\frac{6889}{144} + \frac{41}{24}\pi^2\right]v + \frac{31}{24}v^2 + \frac{7}{1296}v^3\right)x^3 \right. \\
& \qquad \left. + \left(\frac{2835}{128} + \left[\frac{98869}{5760} - \frac{6455}{1536}\pi^2 - \frac{128}{3}\gamma_E - \frac{64}{3}\ln(16x)\right]v \right. \right. \\
& \qquad \left. \left. + \left[\frac{356035}{3456} - \frac{2255}{576}\pi^2\right]v^2 - \frac{215}{1728}v^3 - \frac{55}{31104}v^4 \right)x^4 \right. \\
& \qquad \left. + \mathcal{O}\left(\frac{1}{c^{10}}\right) \right\},
\end{aligned}
\tag{378}
\end{equation}
$$

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

$$
E^{\text{Schw}} = \mu c^2 \left\{ \frac{1 - 2x}{\sqrt{1 - 3x}} - 1 \right\} = - \frac{\mu c^2 x}{2} \sum_{p=0}^{+\infty} \frac{(2p-1)!!(1-2p)}{(p+1)!} \left(\frac{3x}{2}\right)^p .
$$

<footnote n="75">⁷⁵ Note the exact (Schwarzschild) test-mass limit $v \to 0$ of the energy function:</footnote>
$$ \frac{\partial E}{\partial \Omega} = \Omega \frac{\partial J}{\partial \Omega}, \qquad (379) $$

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

We observe in Eqs. (377) and (378) the occurrence of logarithmic terms $\propto \ln x$ starting at the 4PN order. These terms are due to tail effects occurring in the near zone. They contribute to the Fokker action via the generalization of Eq. (348) to include higher mass and current multipole moments (Foffa and Sturani 2020). The 4PN, 5PN and 6PN coefficients of the logarithms in the circular energy function have been explicitly computed using a mix between the traditional approach (Blanchet et al. 2010a; Damour 2010; Le Tiec et al. 2012b) and the EFT method (Foffa and Sturani 2020; Blanchet et al. 2020):

$$ E_{\text{tail}}^{\text{log}} = -\frac{\mu c^2 x}{2} \left\{ \frac{448}{15} v x^4 \ln x + \left( -\frac{4988}{35} - \frac{656}{5} v \right) v x^5 \ln x \right. \\ \left. + \left( -\frac{1967284}{8505} + \frac{914782}{945} v + \frac{32384}{135} v^2 \right) v x^6 \ln x + \mathcal{O}\left(\frac{1}{c^{14}}\right) \right\}. \qquad (380) $$

The 7PN coefficient is also known, however in this case the iterated “tail-of-tail-of-tail” process is also relevant, and the complete 7PN result is not yet known. Furthermore it was shown how to compute, by means of renormalization group (RG) techniques (Goldberger and Ross 2010; Goldberger et al. 2014), the leading logarithmic terms to generic power $n$, which turn out to appear at orders $(3n+1)$PN, with $n=1$ corresponding to the 4PN logarithm in (380), then $n=2$ being the first logarithm square at 7PN order, and so on. We have, up to any power $n$ (Blanchet et al. 2020):

$$ E_{\text{tail}}^{\text{leading-}} (\log)^n = - \frac{\mu c^2 x}{2} \left\{ \frac{64v}{15} \sum_{n=1}^{+\infty} \frac{6n+1}{n!} (4\beta_2^{(m)})^{n-1} x^{3n+1} (\ln x)^n \right\}, \qquad (381a) $$

$$ J_{\text{tail}}^{\text{leading-}} (\log)^n = \frac{G \mu m}{c x^{1/2}} \left\{ - \frac{64v}{15} \sum_{n=1}^{+\infty} \frac{3n+2}{n!} (4\beta_2^{(m)})^{n-1} x^{3n+1} (\ln x)^n \right\}, \qquad (381b) $$

where $\beta_2^{(m)} = -\frac{214}{105}$ is the $\beta$-function coefficient associated with the logarithmic renormalization of the mass-type quadrupole moment ($\ell=2$), and given for an arbitrary multipolarity by Eq. (201). For the linear log-term ($n=1$) one recovers the 4PN coefficient in (380), while further comparison with high-order state-of-the-art self-force calculations by Kavanagh et al. (2015) reveals perfect match with the first terms up to $n \le 7$ in the infinite series (381a). Indeed, Eqs. (381) are first order effects in the test-mass limit $v \to 0$. Note that the comparison necessitates the first law of compact binary dynamics, reviewed in Sect. 3.3.3. The gravitational self-force (GSF) results are expressed in terms of the Detweiler (2008) redshift variable. The first law permits to convey all the information from the redshift variable to the binary’s energy function. The leading logarithms to any power $n$ in the redshift are
given by Eq. (446). All of this shows the great consistency between EFT and RG techniques, the traditional PN method and high-accuracy GSF calculations.

We end this section by giving also the periastron advance in the case of circular orbits to 4PN order (Damour et al. 2015; Bernard et al. 2017b):

$$ K = 1 + 3x + \left(\frac{27}{2} - 7v\right)x^2 + \left(\frac{135}{2} + \left[-\frac{649}{4} + \frac{123}{32}\pi^2\right]v + 7v^2\right)x^3 \\ + \left(\frac{2835}{8} + \left[-\frac{275941}{360} + \frac{48007}{3072}\pi^2 - \frac{1256}{15}\ln x - \frac{592}{15}\ln 2 - \frac{1458}{5}\ln 3\right]v + \frac{2512}{15}\gamma_E\right)x^4 + O\left(\frac{1}{c^{10}}\right), \quad (382) $$

The GSF contribution therein is generally described by the function $\rho(x)$ such that $K^{-2} = 1 - 6x + v\rho(x) + O(v^2)$, i.e.

$$ \rho = 14x^2 + \left(\frac{397}{2} - \frac{123}{16}\pi^2\right)x^3 \\ + \left(-\frac{215729}{180} + \frac{58265}{1536}\pi^2 + \frac{1184}{15}\ln 2 + \frac{2916}{5}\ln 3 + \frac{5024}{15}\gamma_E + \frac{2512}{15}\ln x\right)x^4 \\ + O\left(\frac{1}{c^{10}}\right). \quad (383) $$

The 4PN coefficient is of the type $\rho_{4PN} = a_{4PN} + b_{4PN} \ln x$ and in particular, the numerical value of $a_{4PN}$ (together with $a_{3PN}$ and $a_{2PN}$) is in perfect agreement with the GSF calculations of van de Meent (2017). A comparison between the PN prediction for the periastron advance of circular orbits and numerical calculations based on self-force theory in the small mass ratio limit is done by Le Tiec et al. (2011).

### 3.2.6 The 2.5PN metric in the near zone

The near-zone metric is given by Eqs. (161) for general post-Newtonian matter sources. For point-particles binaries all the potentials $V$, $V_i$, ... parametrizing the metric must be computed and iterated for delta-function sources. This entails UV-type divergencies. Up to the 2.5PN order it is sufficient to cure the divergences by means of the Hadamard self-field regularization (Sect. 3.1.1), equivalent to this order to the dimensional regularization (Sect. 3.1.2).

The post-Newtonian iteration of the metric in closed analytic form is greatly helped—and indeed is made possible at all—by the existence of the Fock (1959) elementary solution of the Poisson equation

$$ \Delta g^{\text{Fock}} = \frac{1}{r_1 r_2}, \quad (384) $$

which takes the interesting closed analytic form
$$g^{\text{Fock}} = \ln S, \qquad (385a)$$

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

with $r_a = |\mathbf{x} - \mathbf{y}_a|$ and $r_{12} = |\mathbf{y}_1 - \mathbf{y}_2|$. Thus $S$ is the Euclidean perimeter of the triangle made of the two source points $\mathbf{y}_a$ and the field point $\mathbf{x}$. An important point is that the Fock function solves the Poisson equation (384) in the sense of distribution theory. The generalization to $d$ dimensions can be found in the Appendix C of Blanchet et al. (2004a).

The Fock function is the starting point of a hierarchy of elementary kernel functions, typically solving Poisson equations of the type $\Delta\varphi = r_1^{2n-1}r_2^{2m-1}$, with $(n, m) \in \mathbb{N}^2$, and their iterations. Such kernels play a crucial role in the post-Newtonian iteration to high order. For definiteness we limit ourselves to the following equations:

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

$$\Delta f = g, \qquad \Delta f_{12} = \frac{r_1}{2r_2}, \qquad \Delta f_{21} = \frac{r_2}{2r_1}, \qquad (386b)$$

$$\Delta h = f, \qquad \Delta h_{12} = \frac{r_1^3}{24r_2}, \qquad \Delta h_{21} = \frac{r_2^3}{24r_1}, \qquad (386c)$$

$$\Delta k = \frac{r_1 r_2}{4}, \qquad \Delta k_{12} = f_{21}, \qquad \Delta k_{21} = f_{12}, \qquad (386d)$$

where the numerical coefficients are chosen for convenience. The general structure of the solutions is made of two parts: an homogeneous regular solution of the Laplace or iterated Laplace operator multiplied by $\ln S$; and a specific polynomial of $r_1, r_2$ and $r_{12}$ (Blanchet et al. 1995a; Jaranowski and Schäfer 1998; Blanchet and Faye 2001a; Jaranowski and Schäfer 2015). If we seek for particular solutions, we are free to add a global homogeneous solution, for example adding a numerical constant to the Fock function. However the final solution should contain the proper homogeneous function which is to be selected by the matching procedure of Sect. 2.2.1, see Marchand et al. (2020) for details. Particular solution of the above equations are

$$g = \ln S + \frac{197}{810}, \qquad (387a)$$

$$f = \frac{1}{12} \left[ (r_1^2 + r_2^2 - r_{12}^2) \left( \ln S - \frac{73}{810} \right) + r_1 r_{12} + r_2 r_{12} - r_1 r_2 \right], \qquad (387b)$$
$$h = \frac{1}{320} \left[ \begin{aligned} & (r_1^4 + r_2^4 - r_{12}^4 - 2r_{12}^2(r_1^2 + r_2^2) + \frac{2}{3}r_1^2r_2^2) \left( \ln S - \frac{37}{81} \right) + \frac{5}{3}r_{12}(r_1^3 + r_2^3) \\ & + r_1r_2(r_{12}^2 - r_1^2 - r_2^2) + r_{12}(r_1 + r_2)(r_1r_2 - r_{12}) + \frac{4}{9}r_1^2r_2^2 \end{aligned} \right], \qquad (387c)$$

$$k = \frac{1}{120} \left[ \begin{aligned} & (r_{12}^4 - 3r_1^4 - 3r_2^4 + 6r_1^2r_2^2 + 2r_{12}^2(r_1^2 + r_2^2)) \ln S + \frac{21}{10}(r_1^4 + r_2^4) \\ & - \frac{r_{12}^4}{30} + 3r_{12}(r_1^3 + r_2^3) + (r_1^2 + r_2^2)\left(3r_1r_2 - \frac{31}{15}r_{12}^2\right) \\ & + r_1r_2r_{12}^2 - \frac{21}{5}r_1^2r_2^2 - r_{12}(r_1 + r_2)(r_{12}^2 - 3r_1r_2) \end{aligned} \right], \qquad (387d)$$

together with $f_{12}$, $h_{12}$ and $k_{12}$ obtained by exchanging the field point $\mathbf{x}$ with the source point $\mathbf{y}_1$, i.e. $f_{12} = f|_{\mathbf{x}\to\mathbf{y}_1}$, $h_{12} = h|_{\mathbf{x}\to\mathbf{y}_1}$ and $k_{12} = k|_{\mathbf{x}\to\mathbf{y}_1}$, and similarly $f_{21}$, $h_{21}$ and $k_{21}$ obtained by exchanging $\mathbf{x}$ and $\mathbf{y}_2$. In addition to satisfying Eqs. (386) the kernels also satisfy

$$\Delta_1 f_{12} = g, \quad \Delta_1 h_{12} = f_{12}, \quad \Delta_1 k = f_{21}, \quad \Delta_1 k_{21} = f, \qquad (388)$$

where $\Delta_a$ is the Laplacian with respect to the source point $\mathbf{y}_a$, together with the equations obtained from 1 $\leftrightarrow$ 2. Note that these equations are true for the specific homogeneous solutions chosen in Eqs. (387). More formulas can be found in the Appendix B of Jaranowski and Schäfer (2015).

To obtain the metric at the 2.5PN order, we need $g$ and also the solutions of more difficult elementary Poisson equations. Namely we meet

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

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

with $a \partial_i$ denoting the partial derivatives with respect to the source point $\mathbf{y}_a$ (and as usual $\partial_i$ being the partial derivative with respect to the field point $\mathbf{x}^i$). It is quite remarkable that the solutions of the latter equations are known in closed analytic form. By combining several earlier results from Carmeli (1965), Ohta et al. (1973), Schäfer (1987), one can write these solutions as (Blanchet et al. 1995a, 1998)

$$K_1 = \left(\frac{1}{2} \Delta - \Delta_1\right) \left[\frac{\ln r_1}{r_2}\right] + \frac{1}{2} \Delta_2 \left[\frac{\ln r_{12}}{r_2}\right] + \frac{r_2}{2r_{12}^2 r_1^2} + \frac{1}{r_{12}^2 r_2^2}, \qquad (390a)$$
$$ H_1 = \frac{1}{2} \Delta_1 \left[ \frac{g}{r_1} + \frac{\ln r_1}{r_{12}} - \Delta_1 \left( \frac{r_1 + r_{12}}{2} g \right) \right] - \frac{r_2}{2r_1^2 r_{12}^2} + \hat{\partial}_i \hat{\partial}_i \left[ \frac{\ln r_{12}}{r_1} + \frac{\ln r_1}{2r_{12}} \right] + \frac{1}{2} \Delta_2 \left[ \frac{\ln r_{12}}{r_1} \right], \quad (390b) $$

or, equivalently, in expanded form,

$$ K_1 = -\frac{1}{r_2^3} + \frac{1}{r_2 r_{12}^2} - \frac{1}{r_1^2 r_2} + \frac{r_2}{2 r_1^2 r_{12}^2} + \frac{r_{12}^2}{2 r_1^2 r_2^3} + \frac{r_1^2}{2 r_2^3 r_{12}^2}, \quad (391a) $$

$$ H_1 = -\frac{1}{2r_1^3} - \frac{1}{4r_{12}^3} - \frac{1}{4r_1^2 r_{12}} - \frac{r_2}{2r_1^2 r_{12}^2} + \frac{r_2}{2r_1^3 r_{12}} + \frac{3r_2^2}{4r_1^2 r_{12}^3} + \frac{r_2^2}{2r_1^3 r_{12}^2} - \frac{r_2^3}{2r_1^3 r_{12}^3}. \quad (391b) $$

Again, those solutions are valid in the sense of distribution theory. Furthermore they tend to zero when $r \to \infty$.

We report the complete expression of the 2.5PN metric in harmonic coordinates, valid at any field point in the near zone. Posing $g_{\alpha\beta} = \eta_{\alpha\beta} + k_{\alpha\beta}$—to distinguish it from the gothic metric deviation $\sqrt{-gg^{\alpha\beta}} = \eta^{\alpha\beta} + h^{\alpha\beta}$—we have (Blanchet et al. 1998)

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

$$
\begin{align}
k_{ij} = & \frac{2Gm_1}{c^2 r_1} \delta^{ij} \nonumber \\
& + \frac{1}{c^4} \left[ \delta_{ij} \left( -\frac{Gm_1}{r_1} (n_1 v_1)^2 + \frac{G^2 m_1^2}{r_1^2} \right. \right. \nonumber \\
& \qquad \left. \left. + G^2 m_1 m_2 \left( \frac{2}{r_1 r_2} - \frac{r_1}{2r_{12}^3} + \frac{r_1^2}{2r_2 r_{12}^3} - \frac{5}{2r_1 r_{12}} + \frac{4}{r_{12} S} \right) \right) \right. \nonumber \\
& \qquad + 4 \frac{Gm_1}{r_1} v_1^i v_1^j + \frac{G^2 m_1^2}{r_1^2} n_1^i n_1^j - 4 G^2 m_1 m_2 n_{12}^i n_{12}^j \left( \frac{1}{S^2} + \frac{1}{r_{12} S} \right) \tag{392c} \\
& + \frac{4G^2 m_1 m_2}{S^2} \left( n_1^{(i} n_2^{j)} + 2 n_1^{(i} n_{12}^{j)} \right) \nonumber \\
& + \frac{G^2 m_1 m_2}{c^5 r_{12}^2} \left( -\frac{2}{3} (n_{12} v_{12}) \delta^{ij} - 6 (n_{12} v_{12}) n_{12}^i n_{12}^j + 8 n_{12}^{(i} v_{12}^{j)} \right) \nonumber \\
& + 1 \leftrightarrow 2 + \mathcal{O}\left(\frac{1}{c^6}\right). \nonumber
\end{align}
$$
To higher order one needs the higher kernels $f, h, \dots$ given by (387) but also the solution of elementary equations still more intricate than (389); the 3PN metric valid in closed form all over the near zone is not currently known.

Finally let us display the 2.5PN metric computed at the location of the particle 1 for instance; to this order the Hadamard self-field regularization, Eq. (255), gives the same result as dimensional regularization. We get

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

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

The metric regularized at the location of the particles has been computed to 3PN order by Blanchet et al. (2010b). At that order it becomes crucial to make full use of dimensional regularization, and poles $\propto 1/\epsilon$ appear in the near-zone metric. The metric regularized at the location of the particles is needed when we compute the redshift observable (427) in Sects. 3.3.3 and 3.3.4.

## 3.3 Conservative dynamics of compact binaries

### 3.3.1 The innermost circular orbit

Having in hand the total relativistic energy for circular orbits (including the rest-mass contribution) at 4PN order, thus $M = m + E$ where $E(x)$ is given by Eq. (377),<fnref n="76" /> we define the innermost circular orbit (ICO) as the minimum, when it exists, of this energy function (Blanchet 2002). Notice that the ICO is not defined as a point of dynamical general-relativistic instability. Hence, we prefer to call this point the ICO rather than, strictly speaking, an innermost stable circular orbit or ISCO. A study of the dynamical stability of circular binary orbits in the post-Newtonian approximation is reported in Sect. 3.3.2.

The previous definition of the ICO is motivated by the comparison with some results of numerical relativity. Indeed we shall confront the prediction of the standard (Taylor-based) post-Newtonian approximation with numerical computations of the energy of binary black holes under the assumptions of conformal flatness for the spatial metric and of exactly circular orbits (Gourgoulhon et al. 2002; Grandclément et al. 2002; Cook and Pfeiffer 2004; Caudill et al. 2006). The latter restriction is implemented by requiring the existence of an “helical” Killing vector (HKV), which is time-like inside the light cylinder associated with the circular motion, and space-like outside. The HKV will be defined in Eq. (424) below. In these numerical approaches there are no gravitational waves; the field is periodic in time, and the gravitational potentials tend to zero at spatial infinity, within a restricted model equivalent to solving five out of the ten Einstein field equations—the Isenberg and Nester (1980), Wilson and Matthews (1989) approximation. Considering an evolutionary sequence of equilibrium configurations the circular-orbit energy $E(\Omega)$

<footnote n="76">⁷⁶ In all of Sect. 3.3 we pose $G = 1 = c$.</footnote>
and the ICO of binary black holes are obtained numerically (see also Bonazzola et al. 1999; Gourgoulhon et al. 2001; Limousin et al. 2005 for related calculations of the ICO in the case of binary neutron stars and strange quark stars; see also Mora and Will 2004).

Since the numerical calculations of Grandclément et al. (2002), Cook and Pfeiffer (2004) have been performed in the case of two corotating black holes, which are essentially spinning with the orbital angular velocity, for the comparison we must include within our PN formalism the effects of spins appropriate to two corotating Kerr black holes. The total relativistic masses of the two Kerr black holes (with $a = 1, 2$ labelling the black holes) are given by

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

We assume the validity of the Christodoulou (1970); Christodoulou and Ruffini (1971) mass formula for each Kerr black holes; i.e., we neglect the influence of the companion. Here $S_a$ is the spin, related to the usual Kerr parameter by $S_a = m_a a_a$, and $\mu_a \equiv m_a^{\text{irr}}$ is the irreducible mass, not to be confused with the reduced mass, and given by $16\pi\mu_a^2 = A_a$ (with $A_a$ the black hole's surface area). The angular velocity of the black hole, defined by the angular velocity of the outgoing photons located at the light-like horizon, is

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

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

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

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

In the limit of slow rotation we get

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

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

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

We now need the relation between the rotation frequency $\omega_a$ of each of the corotating black holes and the orbital frequency $\Omega$ of the binary system. Indeed $\Omega$ is the basic variable describing each equilibrium configuration calculated numerically
by Grandclément et al. (2002), Cook and Pfeiffer (2004), with the irreducible masses held constant along the numerical evolutionary sequences. Here we report the result of an investigation of the condition for corotation based on the first law of mechanics for spinning black holes (Blanchet et al. 2013), which concludes that the corotation condition (at 2PN order) is

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

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

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

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

$$ E_{\text{SO}} = -mv(m\Omega)^{5/3} \left( \frac{m_1^2}{m^2} \left[ \frac{4}{3} + \left( 4 - \frac{31}{18}\nu \right) (m\Omega)^{2/3} \right] \right. \\ \left. + \nu \left[ 1 + \left( \frac{3}{2} - \frac{5}{3}\nu \right) (m\Omega)^{2/3} \right] \frac{S_1}{m_1^2} + 1 \leftrightarrow 2. \right) \quad (399) $$

We shall review in Sect. 3.6 the most up-to-date results for the spin-orbit energy and related quantities; here we are simply employing the formula (592). We immediately infer from Eq. (399) that in the case of corotating black holes the SO effect is equivalent dominantly to a 3PN orbital effect and thus must be retained with the present accuracy. Since our goaled precision is 4PN we have to retain also the next-

<footnote n="77">⁷⁷ Note that this is an iterative process because the masses in Eq. (398) are themselves to be replaced by their expression in terms of the irreducible masses.</footnote>
to-leading correction in (399). With this approximation, the masses in Eq. (399) can be replaced by the irreducible ones. As for the SS interaction (still in the case of aligned spins) it is given by

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

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

$$\Delta E^{\text{corot}} = \mu x_{\mu} \left[ (2 - 6\eta)x_{\mu}^2 + \eta(-10 + 25\eta)x_{\mu}^3 \right. \\ \left. + \eta \left( -21 + \frac{224}{3}\eta - \frac{455}{12}\eta^2 \right) x_{\mu}^4 + \mathcal{O}(x_{\mu}^5) \right]. \qquad (401)$$

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

Figure 6 shows the post-Newtonian results for $E_{\text{ICO}}$ in the case of irrotational and corotational binaries and make comparisons with numerical results. Since $\Delta E^{\text{corot}}$, given by Eq. (401), is at least of order 2PN, the result for $1PN^{\text{corot}}$ is the same as for $1PN$ in the irrotational case; then, obviously, $2PN^{\text{corot}}$ takes into account only the leading 2PN corotation effect, i.e., the spin kinetic energy given by Eq. (397b), while $3PN^{\text{corot}}$ and $4PN^{\text{corot}}$ involve also the corotational SO coupling (401). As we can see the 3PN and 4PN points, and even the 2PN ones, are in good agreement with the numerical values.$^{79}$ The fact that the 2PN, 3PN and 4PN values are so close to each other is a good sign of the convergence of the (Taylor-expanded) series. The agreement found in Fig. 6 is remarkable, all the more if we remember the previous situation, since the first estimations of the ICO in post-Newtonian theory (Kidder et al. 1993a) and numerical relativity (Cook 1994; Pfeiffer et al. 2000; Baumgarte 2000) disagreed with each other, and do not match with the present 3PN/4PN results. The convergence of the standard (Taylor-expanded) post-Newtonian approximation is further discussed by Blanchet (2011a).

<footnote>$^{78}$ This formula extends the works of Blanchet (2002), Blanchet et al. (2013) to the next 4PN order. Note that Blanchet (2002) had initially assumed that the corotation condition was given by the leading-order result $\omega_a = \Omega$; but as was shown by Blanchet et al. (2013) the 1PN correction in Eq. (398) modifies the 3PN terms in (401) with respect to the initial result of Blanchet (2002).</footnote>

<footnote>$^{79}$ As a caveat recall that the numerical works of Gourgoulhon et al. (2002), Grandclément et al. (2002), Cook and Pfeiffer (2004), Caudill et al. (2006) assume that the spatial metric is conformally flat, which is incompatible with the post-Newtonian approximation starting from the 2PN order (see Faye et al. 2004 for a discussion).</footnote>
<figure><img src="image_7.png" /><figcaption>Fig. 6 The binding energy $E_{ICO}$ versus $\Omega_{ICO}$ in the equal-mass case ($v = 1/4$). Top panel: Comparison with the numerical relativity result of Gourgoulhon et al. (2002); Grandclément et al. (2002) valid in the corotating case (marked by a star). Points indicated by $nPN$ are computed from the minimum of Eq. (377), and correspond to irrotational binaries. Points denoted by $nPN^{corot}$ include the contribution (401) to the energy, and describe corotational binaries. Bottom panel: Numerical relativity results of Cook and Pfeiffer (2004), Caudill et al. (2006) for quasi-equilibrium (QE) configurations and various boundary conditions for the lapse function, in the non-spinning (NS), leading-order non spinning (LN) and corotating (CO) cases. The point from Gourgoulhon et al. (2002), Grandclément et al. (2002) (HKV-GGB) is also reported, together with IVP, the initial value approach with effective potential of Cook (1994), Pfeiffer et al. (2000), as well as standard PN predictions from the left panel and non-standard (EOB) ones</figcaption></figure>

### 3.3.2 Dynamical stability of circular orbits

In this section, we investigate the problem of the stability of circular orbits against dynamical perturbations. We propose to use two different methods, one based on a linear perturbation at the level of the center-of-mass equations of motion in harmonic coordinates, the other one consisting of perturbing the center-of-mass Hamiltonian equations in ADM coordinates. We shall find a gauge invariant criterion for the stability of circular orbits (Blanchet and Iyer 2003). We shall check that our two methods agree on the result.

We deal first with the perturbation of the equations of motion, following Kidder et al. (1993a), see their Sect. III.A. We introduce polar coordinates $(r, \varphi)$ in the orbital plane and pose $u \equiv \dot{r}$ and $\Omega \equiv \dot{\varphi}$. In this section we restrict ourselves to the 3PN order. Then the equations of motion are purely instantaneous, as we neglect the non-local tail term (360). Then Eq. (357) yields the system of equations

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

$$ \dot{u} = -\frac{Gm}{r^2} [1 + A + Bu] + r\Omega^2, \qquad (402b) $$
$$ \dot{\Omega} = -\Omega \left[ \frac{Gm}{r^2} \mathcal{B} + \frac{2u}{r} \right], \qquad (402c) $$

where $\mathcal{A}$ and $\mathcal{B}$ are provided by Eqs. (358)–(359), but limited to 3PN order, as functions of $r$, $u$ and $\Omega$ (through $v^2 = u^2 + r^2\Omega^2$).

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

$$ \Omega_0^2 = \frac{Gm}{r_0^3} (1 + \mathcal{A}_0). \qquad (403) $$

Solving iteratively this relation at the 3PN order using the equations of motion, we obtain $\Omega_0$ as a function of the circular-orbit radius $r_0$ in standard harmonic coordinates; the result agrees with the 3PN restriction of Eq. (373).<fnref n="80" />

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

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

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

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

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

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

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

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

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

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

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

<footnote n="80">⁸⁰ Of course one should distinguish $r_0$ from the regularization constant $r'_0$ entering the logarithms at the 3PN and 4PN orders in Eq. (373), and also from $r_0$ in Eq. (72).</footnote>
$$ \gamma_0 = -\Omega_0 \left[ \frac{2}{r_0} + \frac{Gm}{r_0^2} \left( \frac{\partial B}{\partial u} \right)_0 \right]. \qquad (406c) $$

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

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

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

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

where we denote $\gamma_0 = \frac{Gm}{r_0 c^2}$ and we recall that $r_0$ is the radius of the circular orbit in harmonic coordinates.

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

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

This yields a reduced Hamiltonian that is a function of $R, P_R$ and $P_\psi$, and describes the motion in polar coordinates in the orbital plane, hence $\mathcal{H} = \mathcal{H}[R, P_R, P_\psi]$. The Hamiltonian equations then read

$$ \dot{R} = \frac{\partial \mathcal{H}}{\partial P_R}, \qquad \dot{\Psi} = \frac{\partial \mathcal{H}}{\partial P_\psi}, \qquad (410a) $$

$$ \dot{P}_R = -\frac{\partial \mathcal{H}}{\partial R}, \qquad \dot{P}_\psi = 0. \qquad (410b) $$

Evidently the constant $P_\psi$ is nothing but the conserved angular-momentum integral. For circular orbits we have $R = R_0$ (a constant) and $P_R = 0$, so
$$ \frac{\partial \mathcal{H}}{\partial R} [R_0, 0, P_{\psi}^0] = 0, \qquad (411) $$

which gives the angular momentum $P_{\psi}^0$ of the orbit as a function of $R_0$, and

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

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

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

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

We consider now a perturbation of the circular orbit defined by

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

$$ P_{\psi} = P_{\psi}^{0} + \delta P_{\psi}, \qquad (414b) $$

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

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

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

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

$$ \delta \dot{P}_{\psi} = 0, \qquad (415b) $$

$$ \delta \dot{R} = \sigma_0 \delta P_R, \qquad (415c) $$

$$ \delta \Omega = \rho_0 \delta R + \tau_0 \delta P_\psi, \qquad (415d) $$

where the coefficients, which depend on the unperturbed orbit, are given by

$$ \pi_0 = \frac{\partial^2 \mathcal{H}}{\partial R^2} [R_0, 0, P_{\psi}^0], \qquad (416a) $$

$$ \rho_0 = \frac{\partial^2 \mathcal{H}}{\partial R \partial P_{\psi}} [R_0, 0, P_{\psi}^0], \qquad (416b) $$
$$ \sigma_0 = \frac{\partial^2 \mathcal{H}}{\partial P_R^2} [R_0, 0, P_\psi^0], \qquad (416c) $$

$$ \tau_0 = \frac{\partial^2 \mathcal{H}}{\partial P_\psi^2} [R_0, 0, P_\psi^0]. \qquad (416d) $$

By looking to solutions proportional to some $e^{i\sigma t}$ one obtains some real frequencies, and therefore one finds stable circular orbits, if and only if

$$ \hat{C}_0 \equiv \pi_0 \sigma_0 > 0. \qquad (417) $$

Using explicitly the Hamiltonian (368) we readily obtain to 3PN order

$$ \hat{C}_0 = \frac{Gm}{R_0^3} \left\{ \begin{aligned} & 1 + \Gamma_0(-9 + v) + \Gamma_0^2 \left( \frac{117}{4} + \frac{43}{8}v + v^2 \right) \\ & + \Gamma_0^3 \left( -61 + \left[ \frac{4777}{48} - \frac{325}{64}\pi^2 \right]v - \frac{31}{8}v^2 + v^3 \right) + \mathcal{O}\left(\frac{1}{c^8}\right) \end{aligned} \right\}, \qquad (418) $$

where $\Gamma_0 = \frac{Gm}{R_0 c^2}$. This result does not look the same as our previous result (408), but this is simply due to the fact that it depends on the ADM radial separation $R_0$ instead of the harmonic one $r_0$.

Fortunately we have derived in Sect. 3.2.3 the material needed to connect $R_0$ to $r_0$ with the 3PN accuracy. Indeed, with Eqs. (344) [together with (4.9)–(4.10) of de Andrade et al. 2001] we have the relation valid for general orbits in an arbitrary frame between the separation vectors in both coordinate systems. Specializing that relation to circular orbits we find

$$ R_0 = r_0 \left\{ 1 + \gamma_0^2 \left( -\frac{1}{4} - \frac{29}{8}v \right) + \gamma_0^3 \left( \left[ \frac{3163}{1680} + \frac{21}{32}\pi^2 - \frac{22}{3}\ln\left(\frac{r_0}{r'_0}\right) \right]v + \frac{3}{8}v^2 \right) + \mathcal{O}\left(\frac{1}{c^8}\right) \right\}. \qquad (419) $$

Note that the difference between $R_0$ and $r_0$ starts at 2PN order. That relation easily permits to perfectly reconcile both expressions (408) and (418).

Finally let us give to $\hat{C}_0$ an invariant meaning by expressing it with the help of the orbital frequency $\Omega_0$ of the circular orbit, or, more conveniently, of the frequency-related parameter $x_0 \equiv (Gm\Omega_0/c^3)^{2/3}$—cf. Equation (375). This allows us to write the criterion for stability as $C_0 > 0$, where $C_0 = \frac{G^2 m^2}{x_0^3} \hat{C}_0$ admits the gauge-invariant form (Blanchet and Iyer 2003)

$$ C_0 = 1 - 6x_0 + 14v x_0^2 + \left( \left[ \frac{397}{2} - \frac{123}{16}\pi^2 \right]v - 14v^2 \right) x_0^3 + \mathcal{O}(x_0^4). \qquad (420) $$

This form is more interesting than the coordinate-dependent expressions (408) or (418), not only because of its invariant meaning, but also because as we see the 1PN term yields exactly the Schwarzschild result that the innermost stable circular orbit or
ISCO of a test particle (i.e., in the limit $v \to 0$) is located at $x_{\text{ISCO}} = \frac{1}{6}$. Thus we find that, at the 1PN order, but for any mass ratio $v$,

$$x_{\text{ISCO}}^{\text{1PN}} = \frac{1}{6}. \qquad (421)$$

One could have expected that some deviations of the order of $v$ already occur at the 1PN order, but it turns out that only from the 2PN order does one find the occurrence of some non-Schwarzschildian corrections proportional to $v$. At the 2PN order we obtain

$$x_{\text{ISCO}}^{\text{2PN}} = \frac{3}{14v} \left( 1 - \sqrt{1 - \frac{14v}{9}} \right). \qquad (422)$$

For equal masses this gives $x_{\text{ISCO}}^{\text{2PN}} \simeq 0.187$. Notice also that the effect of the finite mass corrections is to increase the frequency of the ISCO with respect to the Schwarzschild result, i.e., to make it more inward$^{81}$:

$$x_{\text{ISCO}}^{\text{2PN}} = \frac{1}{6} \left[ 1 + \frac{7}{18} v + \mathcal{O}(v^2) \right]. \qquad (423)$$

Finally, at the 3PN order, for equal masses $v = \frac{1}{4}$, we find that according to this criterion all the circular orbits are stable. More generally, we find that at the 3PN order all orbits are stable when the mass ratio $v$ is larger than some critical value $v_c \simeq 0.183$.

The stability criterion (420) has been compared by Favata (2011b) to various other stability criteria and shown to perform very well. It has also been generalized to spinning black hole binaries by Favata (2011a). Note that this criterion is based on the physical requirement that a stable perturbation should have a real frequency. It gives an innermost stable circular orbit, when it exists, which differs from the innermost circular orbit or ICO defined in Sect. 3.3.1; see Schäfer (2011) for a discussion on the difference between an ISCO and the ICO in the PN context. Note also that the criterion (420) is based on systematic post-Newtonian expansions, without resorting for instance to Padé approximants. Nevertheless, it performs better than other criteria based on various resummation techniques, as discussed by Favata (2011b).

### 3.3.3 The first law of binary point-particle mechanics

In this section we review an interesting relation for binary systems of point particles modelling black hole binaries and moving on circular orbits, known as the “first law of point-particle mechanics”. This law was obtained using post-Newtonian methods by Le Tiec et al. (2012b), but is actually a particular case of a more general law, valid for systems of black holes and extended fluid balls, derived by Friedman et al. (2002).

<footnote>$^{81}$ This tendency is in agreement with numerical and analytical self-force calculations (Barack and Sago 2009; Le Tiec et al. 2012a).</footnote>
Before tackling the problem it is necessary to make more precise the notion of circular orbits. These are obtained from the conservative part of the dynamics, neglecting the dissipative radiation-reaction force responsible for the gravitational-wave inspiral. In post-Newtonian theory this means neglecting the radiation-reaction force at 2.5PN and 3.5PN orders, i.e., considering only the conservative dynamics at the even-parity 1PN, 2PN and 3PN orders. We have seen in Sects. 2.3.2 and 2.3.4 that this clean separation between conservative even-parity and dissipative odd-parity terms breaks at 4PN order, because of a contribution originating from gravitational-wave tails in the radiation-reaction force. The dissipative part of the 4PN tail term has been given in Eq. (361b). In fact at any high order 4PN, 4.5PN, 5PN, etc. there will be a mixture of conservative and dissipative contributions; here we assume that at any order we can neglect the radiation-reaction dissipation effects.

Consider a system of two compact objects moving on circular orbits. We examine first the case of non-spinning objects. With exactly circular orbits the geometry admits a helical Killing vector (HKV) field $K^\alpha$, satisfying the Killing equation $\nabla^{(\alpha}K^{\beta)} = 0$. Imposing the existence of the HKV is the rigorous way to implement the notion of circular orbits. A Killing vector is only defined up to an overall constant factor. The helical Killing vector $K^\alpha$ extends out to a large distance where the geometry is essentially flat. There,

$$K^\alpha \partial_\alpha = \partial_t + \Omega \partial_\varphi, \qquad (424)$$

in any natural coordinate system which respects the helical symmetry (Sago et al. 2008). We let this equality define the overall constant factor, thereby specifying the Killing vector field uniquely. In Eq. (424) $\Omega$ denotes the angular frequency of the binary's circular motion.

An observer moving with one of the particles (say the particle 1), while orbiting around the other particle, would detect no change in the local geometry. Thus the four-velocity $u_1^\alpha$ of that particle is tangent to the Killing vector $K^\alpha$ evaluated at the location of the particle, which we denote by $K_1^\alpha$. A physical quantity is then defined as the constant of proportionality, say $u_1^T$, between these two vectors, namely

$$u_1^\alpha = u_1^T K_1^\alpha. \qquad (425)$$

The four-velocity of the particle is normalized by $(g_{\mu\nu})_1 u_1^\mu u_1^\nu = -1$, where $(g_{\mu\nu})_1$ denotes the metric at the location of the particle. For a self-gravitating compact binary system, the metric at point 1 is generated by the two particles and has to be regularized according to the self-field regularization discussed in Sect. 3.1. It will in fact be sometimes more convenient to work with the inverse of $u_1^T$, denoted $z_1 \equiv 1/u_1^T$. From Eq. (425) we get

$$z_1 = -(u_1 K_1), \qquad (426)$$

where $(u_1 K_1) = (g_{\mu\nu})_1 u_1^\mu K_1^\nu$ denotes the usual space-time dot product. Thus we can regard $z_1$ as the Killing energy of the particle that is associated with the HKV field $K^\alpha$. The quantity $z_1$ represents also the redshift of light rays emitted from the particle and received on the helical symmetry axis perpendicular to the orbital plane at large
distances from it. In the following we shall refer to $z_1$ as the *redshift observable*; it has been introduced by Detweiler (2008).

If we choose a coordinate system such that Eq. (424) is satisfied everywhere, then in particular $K_1^\alpha = 1$, thus $u_1^T$ simply agrees with $u_1^t$, the $t$-component of the four-velocity of the particle. The Killing vector on the particle is then $K_1^\alpha = u_1^\alpha / u_1^t$, and simply reduces to the particle's ordinary coordinate velocity: $K_1^\alpha = v_1^\alpha$ where $v_1^\alpha = dy_1^\alpha/dt$ and $y_1^\alpha(t) = (t, y_1(t))$ denotes the particle's trajectory in that coordinate system. The redshift observable is then

$$z_1 = \frac{1}{u_1^T} = \sqrt{-(g_{\mu\nu})_1 v_1^\mu v_1^\nu} \qquad (427)$$

It is important to note that for circular orbits this quantity does not depend upon the choice of coordinates; in a perturbative approach in which the perturbative parameter is the particles' mass ratio $v \ll 1$, it does not depend upon the choice of perturbative gauge with respect to the background metric. We shall be interested in the *invariant* scalar function $z_1(\Omega)$, where $\Omega$ is the angular frequency of the circular orbit defined by Eq. (424).

We have obtained in Sect. 3.2.5 the expressions of the post-Newtonian binding energy $E$ and angular momentum $J$ for point-particle binaries on circular orbits. We shall now show that there are some differential and algebraic relations linking $E$ and $J$ to the redshift observables $z_1$ and $z_2$ associated with the two individual particles. Here we prefer to introduce instead of $E$ the total relativistic (ADM) mass of the binary system

$$M = m + E, \qquad (428)$$

where $m$ is the sum of the two post-Newtonian individual masses $m_1$ and $m_2$—those which have been used up to now, for instance in Eq. (339). Note that in the spinning case such post-Newtonian masses acquire some spin contributions given by Eqs. (394)–(397).

For point particles without spins, the ADM mass $M$, angular momentum $J$, and redshifts $z_a$, are functions of three independent variables, namely the orbital frequency $\Omega$ that is imposed by the existence of the HKV, and the individual masses $m_a$. For spinning point particles, we have also the two spins $S_a$ which are necessarily aligned with the orbital angular momentum, and furthermore should satisfy a corotation condition. We first recall that the ADM quantities obey the “thermodynamical” relation already met in Eq. (379),

$$\frac{\partial M}{\partial \Omega} = \Omega \frac{\partial J}{\partial \Omega}. \qquad (429)$$

Such relation is commonly used in post-Newtonian theory, see e.g., Damour et al. (2000b), Blanchet (2002), Blanchet et al. (2020); it states that the gravitational-wave energy and angular momentum *fluxes* are strictly proportional for circular orbits, with $\Omega$ being the coefficient of proportionality. The relation is also used in numerical computations of the binary evolution based on a sequence of quasi-equilibrium
configurations (Gourgoulhon et al. 2002; Grandclément et al. 2002; Cook and Pfeiffer 2004; Caudill et al. 2006), as reviewed in Sect. 3.3.1.

The first law is a generalization of Eq. (429), describing the changes in the ADM quantities not only when the orbital frequency $\Omega$ varies with fixed masses, but also when the individual masses $m_a$ of the particles vary with fixed orbital frequency. That is, one compares together different conservative dynamics with different masses but the same orbital frequency. This situation is answered by the differential equations

$$ \frac{\partial M}{\partial m_a} - \Omega \frac{\partial J}{\partial m_a} = z_a \quad (a = 1, 2). \qquad (430) $$

Finally the relations (429)–(430) can be summarized by the following.

**Theorem 13** (Friedman et al. 2002; LeTiec et al. 2012b) The changes in the ADM mass and angular momentum of a binary system made of point particles on a circular orbit, in response to infinitesimal variations of the individual masses of the point particles, are related together by the first law of binary point-particle mechanics

$$ \delta M - \Omega \delta J = \sum_a z_a \delta m_a. \qquad (431) $$

This law was initially proved in a very general way by Friedman et al. (2002) for systems of black holes and extended bodies under some arbitrary Killing symmetry. The particular form given in Eq. (431) is a specialization to the case of point particle binaries with helical Killing vector. It has been derived directly in this form by Le Tiec et al. (2012b) up to high post-Newtonian order, viz. 3PN and the logarithmic terms at 4PN and 5PN orders. Furthermore, it was proved by Blanchet and Le Tiec (2017) that the first law still holds when taking into account the non-local tail term at 4PN order.

An interesting by-product of the first law (431) is the simple relation

$$ M - 2\Omega J = \sum_a z_a m_a, \qquad (432) $$

which can be seen as a first integral of the differential relation (431). The existence of such a simple algebraic relation between the quantities $z_a$ belonging to each particles, and the globally defined quantities M and J, is remarkable.

The law (431) has been generalized by Le Tiec (2015) for binary systems of point masses moving on generic eccentric orbits, using a natural orbit-averaged definition of the redshift observable (Barack and Sago 2011). An application of the first law is the calculation of the exact gravitational self-force contributions in the binding energy and angular momentum for circular orbits, allowing a coordinate-invariant comparison to numerical relativity (Le Tiec et al. 2012a).

In the case of an asymptotically flat vacuum spacetime with two black holes on circular orbits, we have (Friedman et al. 2002)
$$ \delta M - \Omega \delta J = \sum_a \frac{\kappa_a}{8\pi} \delta A_a, \qquad (433) $$

which can be viewed as a generalization of the celebrated first law of black hole mechanics $\delta M - \omega_H \delta J = \frac{\kappa}{8\pi} \delta A$, which holds for any non-singular, asymptotically flat perturbation of a stationary and axisymmetric black hole of mass M, intrinsic angular momentum (or spin) $J \equiv Ma$, surface area $A$, uniform surface gravity $\kappa$, and angular frequency $\omega_H$ on the horizon (Bardeen et al. 1973; Wald 1973). In the binary black hole case, the role played by the horizon angular velocity $\omega_H$ of a single rotating black hole is taken by the orbital frequency $\Omega$ of the binary.

In order to interpret the first law (431) (and the redshift factor therein) for actual binary black holes, and to reconcile it with (433), we have to remember that the assumption of helical Killing symmetry requires that the two black holes are in co-rotation, and must therefore have non-zero spins. For rotating black holes, the post-Newtonian mass $m_a$ is given by the total mass of the black hole, including the effect of the spin, see Eq. (394). We therefore need the generalization of the first law (431) applicable to systems of point particles with spins (moving on circular orbits). The result will be valid through linear order in the spin of each particle, but will also hold for the quadratic coupling between different spins (interaction spin terms $S_1 \times S_2$ in the language of Sect. 3.6). To be consistent with the HKV symmetry, we must assume that the two spins $S_a$ are aligned or anti-aligned with the orbital angular momentum. We introduce the total (ADM-like) angular momentum $J$ which is related to the orbital angular momentum $L$ by $J = L + \sum_a S_a$ for aligned or anti-aligned spins. The first law now reads

$$ \delta M - \Omega \delta J = \sum_a \left[ z_a \delta m_a + (\Omega_a - \Omega) \delta S_a \right], \qquad (434) $$

where $\Omega_a = |\Omega_a|$ denotes the precession frequency of the spins. This law has been derived by Blanchet et al. (2013) using the canonical Hamiltonian formalism. The spin variables in (434) are the canonical spins $S_a$, that are easily seen to obey, from the algebra satisfied by the canonical variables, the usual Newtonian-looking precession equations $dS_a/dt = \Omega_a \times S_a$. These variables are identical to the “constant-in-magnitude” spins which will be defined and extensively used in Sect. 3.6. Similarly to Eq. (432), we have also a first integral associated with the variational law (434):

$$ M - 2\Omega J = \sum_a \left[ z_a m_a + 2(\Omega_a - \Omega) S_a \right]. \qquad (435) $$

Notice that the relation (434) has been derived for point particles and arbitrary aligned spins. We want now to derive the analogous relation for binary black holes. The key difference is that black holes are extended finite-size objects while point particles have by definition no spatial extension. For point particle binaries the spins can have arbitrary magnitude and still be compatible with the HKV. In this case the law (434) would describe also (super-extremal) naked singularities. For black hole
binaries the HKV constraints the rotational state of each black hole and the binary system must be corotating.

Let us derive, in a heuristic way, the analogue of the first law (434) for general bodies, by introducing some “constitutive relations” $m_a(\mu_a, S_a, \dots)$ specifying the energy content of the bodies, i.e., the relations linking their masses $m_a$ to the spins $S_a$ and to some irreducible masses $\mu_a$. More precisely, we define for each spinning particle the analogue of an irreducible mass $\mu_a \equiv m_a^{\text{irr}}$ via the variational relation $\delta m_a = c_a \delta \mu_a + \omega_a \delta S_a$, in which the “response coefficient” $c_a$ of the body and its proper rotation frequency $\omega_a$ are associated with the internal structure:

$$c_a \equiv \left. \frac{\partial m_a}{\partial \mu_a} \right|_{S_a}, \qquad (436a)$$

$$\omega_a \equiv \left. \frac{\partial m_a}{\partial S_a} \right|_{\mu_a}. \qquad (436b)$$

For instance, using the Christodoulou mass formula (394) for Kerr black holes, we obtain the rotation frequency $\omega_a$ given by Eq. (395). On the other hand, the response coefficient $c_a$ differs from 1 only because of spin effects, and we can check that $c_a = 1 + \mathcal{O}(S_a^2)$.

Within the latter heuristic model a condition for the corotation of black hole binaries has been proposed by Blanchet et al. (2013) as

$$z_a \omega_a + \Omega_a = \Omega. \qquad (437)$$

This condition determines the value of the proper frequency $\omega_a$ of each black hole appropriate to the corotation state. When expanded to 2PN order the condition (437) leads to Eq. (398) that we have already used in Sect. 3.3.1. With Eq. (437) imposed, the first law (434) simplifies considerably:

$$\delta M - \Omega \delta J = \sum_a c_a z_a \delta \mu_a. \qquad (438)$$

This is almost identical to the first law for non-spinning binaries given by Eq. (431); indeed it simply differs from it by the substitutions $c_a \to 1$ and $\mu_a \to m_a$. Since the irreducible mass $\mu_a$ of a rotating black hole is the spin-independent part of its total mass $m_a$, this observation suggests that corotating binaries are very similar to non-spinning binaries, at least from the perspective of the first law. Finally we can easily reconcile the first law (438) for corotating systems with the known first law of binary black hole mechanics given by (433). Indeed, since the surface area is related to the irreducible mass by $\mu_a^2 = A_a/(16\pi)$, it suffices to make the formal identification in Eq. (438):

$$c_a z_a \equiv 4 \mu_a \kappa_a, \qquad (439)$$

where $\kappa_a$ denotes the constant surface gravity. In particular, for two black holes on circular orbits, but far enough apart so that they can be viewed—in first approximation—as isolated, the surface gravity $\kappa_a \to (4\mu_a)^{-1}$, while for the point particles
we have both $z_a \to 1$ and $c_a \to 1$. The nice correspondence (439) shows that the heuristic model based on the constitutive relations (436) is able to capture the physics of corotating black hole binary systems.

### 3.3.4 Post-Newtonian approximation versus gravitational self-force

The high-accuracy predictions from GR we have drawn up to now are well suited to describe the inspiralling phase of compact binaries, when the post-Newtonian parameter $\epsilon$ is small, independently of the mass ratio $q \equiv m_1/m_2$ between the compact bodies. In this section we investigate how well does the post-Newtonian expansion, compared with another very important approximation scheme in general relativity: The gravitational *self-force* approach, based on black-hole perturbation theory, which gives an accurate description of *extreme mass ratio* binaries having $q \ll 1$ or equivalently $v \ll 1$, even in the strong field regime. It is thus expected to provide templates for extreme mass ratio inspirals (EMRI) anticipated to be present in the bandwidth of space-based detectors.

The gravitational self-force analysis was pioneered by De Witt and Brehme (1960), Mino et al. (1997b), Quinn and Wald (1997), Detweiler and Whiting (2003), Gralla and Wald (2008); see Poisson et al. (2011), Detweiler (2011), Barack (2011), Barack and Pound (2018) for reviews. The expansion at leading order in the mass

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Different analytical approximation schemes and numerical techniques to study black hole binaries, depending on the mass ratio $q = m_1/m_2$ and the post-Newtonian parameter $\epsilon^2 \sim v^2/c^2 \sim Gm/(c^2r_{12})$. Post-Newtonian theory and perturbative self-force analysis can be compared in the post-Newtonian regime ($\epsilon \ll 1$ thus $r_{12} \gg Gm/c^2$) of an extreme mass ratio ($m_1 \ll m_2$) binary. Image courtesy A. Le Tiec</figcaption></figure>
ratio reaches impressive PN equivalent levels (Fujita 2012a, b; Kavanagh et al. 2015). The state-of-the-art is currently the second order in the mass ratio, for both the binding energy (Pound et al. 2020) and the energy flux and waveform (Warburton et al. 2021; Wardell et al. 2023).

Consider a system of two (non-spinning) compact objects with mass ratio $q = m_1/m_2 \ll 1$; we shall call the smaller mass $m_1$ the “particle”, and the larger mass $m_2$ the “black hole”. The orbit of the particle around the black hole is supposed to be exactly circular as we neglect the radiation-reaction effects. With circular orbits and no dissipation, we are considering the conservative part of the dynamics, and the geometry admits the HKV field (424). Note that in self-force theory there is a clean split between the dissipative and conservative parts of the dynamics (see e.g., Barack 2009). This split is particularly transparent for an exact circular orbit, since the radial component (along $r$) is the only non-vanishing component of the conservative self-force, while the dissipative part of the self-force are the components along $t$ and $\phi$.

The problem of the comparison between the post-Newtonian and perturbative self-force analyses in their common domain of validity, that of the slow-motion and weak-field regime of an extreme mass ratio binary, is illustrated in Fig. 7. This problem has been tackled by Detweiler (2008), who computed numerically within the self-force (SF) approach the redshift observable $u_1^T$ associated with the particle, and compared it with the 2PN prediction extracted from the post-Newtonian results of Blanchet et al. (1998). This comparison proved to be successful, and therefore was later systematically implemented and extended to higher PN orders by Blanchet et al. (2010a, b) (see also Barack et al. 2010; Barack and Sago 2011 for related comparisons). In this section we review the latter works which have demonstrated an excellent agreement between the analytical post-Newtonian result derived through 3PN order, with inclusion of specific logarithmic terms at 4PN and 5PN orders, and the exact numerical SF result.

For the PN-SF comparison, we require two physical quantities which are precisely defined in the context of each of the approximation schemes. The orbital frequency $\Omega$ of the circular orbit as measured by a distant observer is one such quantity and has been introduced in Eq. (424); the second quantity is the Detweiler (2008) redshift observable $u_1^T$ (or equivalently $z_1 = 1/u_1^T$) associated with the smaller mass $m_1 \ll m_2$ and defined by Eqs. (425) or (426). The truly coordinate and perturbative-gauge independent properties of $\Omega$ and the redshift observable $u_1^T$ play a crucial role in this comparison. In the perturbative self-force approach we use Schwarzschild coordinates for the background, and we refer to “gauge invariance” as a property which holds within the restricted class of gauges for which (424) is a helical Killing vector. In all other respects, the gauge choice is arbitrary. In the post-Newtonian approach we generally work with harmonic coordinates and compute the explicit expression (427) of the redshift observable.

The main difficulty in the post-Newtonian calculation is the control to high PN order of the near-zone metric $(g_{\mu\nu})_1$ entering the definition of the redshift observable (427), and which has to be regularized at the location of the particle by means of dimensional regularization (see Sects. 3.1.2 and 3.1.3). Up to 2.5PN order the Hadamard regularization is sufficient and the regularized metric has been provided in
Eqs. (393). However at 3PN order dimensional regularization plays a crucial role as there are poles $\propto 1/\epsilon$ which develop in the near zone metric. As shown by Blanchet et al. (2010b) the poles are cancelled in fine when the redshift observable is expressed in terms of the gauge invariant post-Newtonian parameter (375). At the next 4PN order we do not know the near zone metric and thus a direct computation of the redshift using Eq. (427) is not possible. Nevertheless we can infer indirectly the 4PN coefficient in the redshift from the known 4PN coefficient in the energy function for circular orbits (377) by using the first law (431). Indeed, as found by Le Tiec et al. (2012a, b) the first law is fully satisfied if and only if there is a precise relationship between the post-Newtonian coefficients in the redshift on the one hand, and those in the conserved energy and angular momentum on the other hand. Finally the complete result for the redshift at 4PN order is

$$ 
\begin{aligned}
z_1 = 1 &+ \left(-\frac{3}{2}X_2 + \frac{\nu}{2}\right)x \\
&+ \left(-\frac{9}{8}X_2 + \left[-\frac{3}{8} - \frac{X_2}{4}\right]\nu + \frac{5}{24}\nu^2\right)x^2 \\
&+ \left(-\frac{27}{16}X_2 + \left[-\frac{27}{16} + \frac{19}{8}X_2\right]\nu + \left[-\frac{19}{16} - \frac{X_2}{16}\right]\nu^2 + \frac{\nu^3}{16}\right)x^3 \\
&+ \left(-\frac{405}{128}X_2 + \left[-\frac{675}{128} + \left(\frac{6889}{192} - \frac{41}{32}\pi^2\right)X_2\right]\nu\right. \\
&\qquad \left. + \left[-\frac{6889}{1152} + \frac{41}{192}\pi^2 - \frac{93}{64}X_2\right]\nu^2\right. \\
&\qquad \left. + \left[\frac{217}{192} - \frac{7}{864}X_2\right]\nu^3 + \frac{91}{10368}\nu^4\right)x^4 \\
&+ \left(-\frac{1701}{256}X_2 + \left[-\frac{3969}{256} + \left(-\frac{24689}{1920} + \frac{128}{5}\gamma_E + \frac{1291}{512}\pi^2 + \frac{64}{5}\ln(16x)\right)X_2\right]\nu\right. \\
&\qquad \left. + \left[\frac{24463}{11520} + \frac{64}{15}\gamma_E + \frac{1291}{3072}\pi^2 + \frac{32}{15}\ln(16x) + \left(-\frac{71207}{768} + \frac{451}{128}\pi^2\right)X_2\right]\nu^2\right. \\
&\qquad \left. + \left[\frac{356035}{6912} - \frac{2255}{1152}\pi^2 + \frac{43}{288}X_2\right]\nu^3 + \left[-\frac{473}{3456} + \frac{55}{20736}X_2\right]\nu^4\right. \\
&\qquad \left. - \left(\frac{187}{62208}\nu^5\right)x^5 + \mathcal{O}(x^6)\right),
\end{aligned} 
\quad (440) $$

where $x$ is defined by (375), $\nu$ is the mass ratio (353), and we recall that $X_2 \equiv m_2/m$ (with 2 being the “black hole”). The redshift of the other particle is deduced by changing $X_2 \to X_1$.

The post-Newtonian result (440) is valid for any mass ratio, and for comparison purpose with the SF calculation we now investigate the small mass ratio regime $q \ll 1$. We introduce a PN parameter appropriate to the small mass limit of the “particle”,
$$y \equiv \left( \frac{G m_2 \Omega}{c^3} \right)^{2/3} = x (1+q)^{-2/3}. \qquad (441)$$

We express the symmetric mass ratio in terms of the asymmetric one: $v = q(1+q)^{-2}$, together with $X_2 = (q+1)^{-1}$. Then Eq. (440) [or, rather, its inverse $u_1^T = 1/z_1$], expanded through first order in $q$, which means including the linear self-force level, reads

$$u_1^T = u_{\text{Schw}}^T + q u_{\text{SF}}^T + \mathcal{O}(q^2). \qquad (442)$$

The Schwarzschildean result is known in closed form as

$$u_{\text{Schw}}^T = \frac{1}{\sqrt{1-3y}}, \qquad (443)$$

and for the SF contribution one obtains$^{82}$

$$u_{\text{SF}}^T(y) = -y - 2y^2 - 5y^3 + \left(-\frac{121}{3} + \frac{41}{32}\pi^2\right)y^4 + \left(-\frac{1157}{15} + \frac{677}{512}\pi^2 - \frac{128}{5}\gamma_E - \frac{64}{5}\ln(16y)\right)y^5 + o(y^5). \qquad (444)$$

The analytic coefficients in (444) were determined up to 2PN order by Detweiler (2008) using the regularized post-Newtonian metric of Blanchet et al. (1998). The 3PN term was computed by Blanchet et al. (2010b) making full use of dimensional regularization. The comparison of the post-Newtonian expansion (444) with the numerical SF data has confirmed with high precision the analytic determination of the 3PN coefficient. Notice that such agreement provides an independent check of the dimensional regularization procedure at the basis of the PN expansion (see Sects. 3.1.2–3.1.3). It is remarkable that such procedure gives results equivalent to those obtained with the subtraction of the singular field in the SF approach (Detweiler and Whiting 2003). The comparison with SF data also revealed the presence of logarithmic contributions at the 4PN and 5PN orders, which have then been computed analytically (Blanchet et al. 2010a, b, ; Damour 2010). Furthermore the PN-SF comparison has permitted to measure the 4PN and 5PN coefficients, with high significant digits (Le Tiec et al. 2012b). The 4PN coefficient in (444) was computed analytically by Bini and Damour (2013).

For reference let us give the higher-order 2SF and 3SF corrections at the 4PN order, i.e. the next two terms in (442), which are straightforward to compute from Eq. (440):

<footnote>$^{82}$ Since there are logarithms in this expansion we use the Landau $o$-symbol for remainders; see the footnote 26.</footnote>
$$u_{2\text{SF}}^{T}(y) = y + 3y^2 + \frac{97}{8}y^3 + \left(\frac{725}{12} - \frac{41}{64}\pi^2\right)y^4 + \left(\frac{674801}{1920} - \frac{9491}{1024}\pi^2 - \frac{64}{5}\gamma_E - \frac{32}{5}\ln(16y)\right)y^5 + o(y^5), \quad (445a)$$

$$u_{3\text{SF}}^{T}(y) = -\frac{28}{27}y - \frac{67}{18}y^2 - \frac{97}{6}y^3 + \left(-\frac{2099}{27} + \frac{41}{72}\pi^2\right)y^4 + \left(-\frac{73417}{135} + \frac{32123}{2304}\pi^2 + \frac{256}{45}\gamma_E + \frac{128}{45}\ln(16y)\right)y^5 + o(y^5). \quad (445b)$$

We report also the leading logarithmic contributions in the redshift variable at any generic power $n$. These have been computed in the energy and angular momentum using renormalization group techniques, see Eqs. (381). We have (at linear order in the mass ratio $q$)

$$u_{\text{SF}}^{T}(y)\Big|_{\text{leading-}} (\log)^n = -\frac{64}{5} \sum_{n=1}^{+\infty} \frac{1}{n!} (4\beta_2^{(m)})^{n-1} y^{3n+2} (\ln y)^n, \quad (446)$$

where we recall that $\beta_2^{(m)} = -\frac{214}{105}$ is the RG coefficient associated with the mass quadrupole moment. We recover the 4PN logarithm in (444), followed by a logarithm square at 7PN order, a logarithm cube at 10PN order, and so on. The agreement with the high-order analytical (linear in mass ratio) self force calculation of Kavanagh et al. (2015) is perfect. This also highlights the predictive power of perturbative self-force calculations, which is obviously due to the fact (illustrated in Fig. 7) that perturbation theory is legitimate in the strong field regime of the coalescence of black hole binary systems, which is inaccessible to the post-Newtonian method.

The accuracy of the numerical computation of the self-force has been drastically improved by Shah et al. (2014). The PN coefficients of the redshift observable were obtained to very high 10.5PN order both numerically and analytically, for a subset of coefficients that are either rational or made of the product of $\pi$ with a rational. The analytical values of the coefficients up to 6PN order have also been obtained from a perturbative self-force calculation by Bini and Damour (2014a, b). Note that the works of Shah et al. (2014) and Bini and Damour (2013, 2014a, b) use alternative techniques (developed by Mano et al. 1996a, b; Mano and Takasugi 1997) with which to represent metric perturbation solutions for black hole space-times.

An interesting feature of the post-Newtonian expansion of the redshift at high order is the appearance of half-integral PN coefficients (i.e., of the type $\frac{p}{2}$PN where $p$ is an odd integer) in the conservative dynamics of binary point particles, i.e. moving on exactly circular orbits. This is interesting because any instantaneous (non-tail or “local”) term at any half-integral PN order in the redshift will be zero for circular orbits, as can be shown by a simple dimensional argument (Blanchet et al. 2014a). Therefore half-integral coefficients can appear only due to truly hereditary (tail or “non-local”) integrals. Using standard post-Newtonian methods it has been proved by Blanchet et al. (2014a, b) that the dominant half-integral PN term in the redshift observable (444) occurs at the 5.5PN order (confirming the earlier finding of
Shah et al. (2014) and originates from the non-linear “tail-of-tail” integrals investigated in Sect. 2.4.2. The results for the 5.5PN coefficient, and also for the next-to-leading 6.5PN and 7.5PN ones due to higher-order tails-of-tails and computed by Blanchet et al. (2014a, b), fully agree with the SF computations. We emphasize that that the latter results (contrary to various analytical and numerical SF calculations by Shah et al. 2014; Bini and Damour 2013, 2014a, b, ; Kavanagh et al. 2015) are achieved with the “traditional” PN approach, which is completely general, i.e. is not tuned to a particular type of source but is applicable to any extended PN source; see the formalism of Sect. 2.

To conclude, the consistency of this cross-cultural comparison between the analytical post-Newtonian and the perturbative self-force approaches confirms the soundness of both approximations in describing the dynamics of compact binaries. Furthermore this interplay between PN and SF efforts is important for the synthesis of template waveforms of EMRs to be analysed by space-based gravitational-wave detectors, and has also an impact on efforts of numerical relativity in the case of comparable masses.

## 3.4 Gravitational waves from compact binaries

We pointed out that the 4PN equations of motion are merely “1.5PN” as regards the radiative aspects of the problem, because the radiation reaction force starts at the 2.5PN order. Because of this difference between conservative and dissipative (radiation reaction) effects, we stop developing high order PN corrections in the radiation reaction force, and instead apply the gravitational wave-generation formalism described in Sect. 2. In this approach, the work done by the radiation reaction force is computed as a flux at future null infinity. The total energy flux $F \equiv (dE/dt)^{GW}$ drives the evolution of the energy in the system through the flux-balance equation

$$ \frac{dE}{dt} = -\mathcal{F}. \qquad (447) $$

Therefore, consistently with the result (377) we obtained for the 4.5PN binary’s center-of-mass energy E, we must now control the total gravitational-wave flux F at the same relative 4.5PN order, which means, beyond the Einstein quadrupole formula (5).

Because the orbit of inspiralling compact binaries is circular, the energy balance equation is sufficient, and there is no need to invoke the angular momentum balance equation.<fnref n="83" /> Furthermore the time average over one orbital period is here irrelevant, and the angular momentum flux $G \equiv (dJ/dt)^{GW}$ is related to the energy flux by $F = \Omega G$, which is a consequence of (429).

We keep in mind that we shall use the flux-balance equation (447) at very high order, where one is missing a complete proof of it, following from first principles in GR. Nevertheless, Eq. (447) has been verified by radiation-reaction calculations for

<footnote n="83">⁸³ By contrast, remind the computation of the evolution of the orbital period $\dot{P}$ and eccentricity $\dot{e}$ of the binary pulsar, see Eqs. (10)–(14).</footnote>
extended post-Newtonian fluids at 1.5PN order (Blanchet 1997b; Blanchet and Faye 2019), and by explicit computations for compact binaries at the 3.5PN (Jaranowski and Schäfer 1997; Pati and Will 2002; Königsdörffer et al. 2003; Nissanke and Blanchet 2005; Itoh 2009) and 4.5PN (Leibovich et al. 2023) levels, fully consistent with the balance equations (Iyer and Will 1993, 1995; Gopakumar et al. 1997).

Obtaining the energy flux $\mathcal{F}$ can be divided into two equally important tasks: computing the source multipole moments $I_L$ and $J_L$ of the compact binary system with due account of self-field and IR regularizations; and controlling the tails and related non-linear effects occurring in the relation between the binary's source moments and the radiative ones $U_L$ and $V_L$ observed at future null infinity (cf. the general formalism of Sect. 2).

### 3.4.1 The binary's multipole moments for circular orbits

The general expressions of the source multipole moments given by Eqs. (7) in Theorem 7 are first worked out explicitly for general fluid systems. For this computation one uses the formula (137), and we insert the metric coefficients (in harmonic coordinates) expressed by means of the retarded-type elementary potentials.<fnref n="84" /> Then we specialize each of the numerous terms to the case of point-particle binaries by inserting the standard expression (324) for the matter stress-energy tensor made out of Dirac delta-functions.

The most important input for the computation of the waveform and flux is the mass quadrupole moment $I_{ij}$, since this moment necessitates the full PN precision. A first computation of the quadrupole moment at the 4PN order was performed by Marchand et al. (2020), using dimensional regularization for the UV divergences but the Hadamard regularization for the IR ones. However, in order to be consistent with the 4PN equations of motion (Bernard et al. 2016, 2017b, a, 2018), it should be computed with dimensional regularization also for the IR divergences. It was thus completed by Larrouturou et al. (2022b), and the expression of the source quadrupole moment, obtained with full dimensional regularization,<fnref n="85" /> has the expected feature of exhibiting poles in $\epsilon = d-3$. Those poles are crucial to cancel the divergences linked with the $d$ dimensional computation of the radiative moment, i.e. $U_{ij}$, as performed by Larrouturou et al. (2022a). Indeed, the source moments are not observables *per se*, but the radiative moments are; therefore, only the latter have to be finite in the $\epsilon \to 0$ limit. However, for the sake of computational simplicity, it was deemed useful to introduce the notion of a “renormalized” source quadrupole moment $I_{ij}$, defined by Eq. (6.2) of Larrouturou et al. (2022a). This renormalized quantity is nothing but the sum of the $d$-dimensional source quadrupole and the corrections, in the form of a “difference” (see Sect. 3.1.3), arising notably from the relation between radiative and canonical moments, involving the non-linear tails-of-tails, tails-of-memory, etc. listed in Sect. 2.4.5. The renormalized quadrupole, now truly called $I_{ij}$, has no poles in $\epsilon$ and when injected into the three dimensional MPM

<footnote n="84">⁸⁴ See Eqs. (159)–(163) and the generalization at 4PN order given in the Appendix A of Marchand et al. (2020).</footnote>

<footnote n="85">⁸⁵ More precisely: the “$B\epsilon$” variant, see Sect. 3.1.3.</footnote>
algorithm, yields the correct radiative quadrupole moment $U_{ij}$, by definition. Last but not least, up to 3PN order, the corrections due to the IR dimensional regularization exactly cancel the ones due to the renormalization, as proven by Larrouturou et al. (2022a). Thus, up to 3PN order, the renormalized source quadrupole coincides with the one computed with Hadamard regularization in the IR, which is why those subtleties did not hit previous computations of the flux. However, this equivalence is no more valid at 4PN order, thus the crucial need for the implementation of dimensional regularization and the previous “renormalization”.

The mass quadrupole moment complete to order 4PN, for non-spinning compact binaries on quasi circular orbits, is

$$I_{ij} = \mu \left\{ A x_{\langle i} x_{j\rangle} + B \frac{r^2}{c^2} v_{\langle i} v_{j\rangle} + \frac{G^2 m^2 v}{c^5 r} C x_{\langle i} v_{j\rangle} \right\} + O\left(\frac{1}{c^9}\right), \quad (448)$$

where $x = y_1 - y_2 = (x_i)$ and $v = v_1 - v_2 = (v_i)$ are the orbital separation and relative velocity. Introducing the post-Newtonian parameter (369), the 4PN coefficients are given by$^{86}$

$$ 
\begin{aligned}
A = 1 &+ \gamma \left(-\frac{1}{42} - \frac{13}{14}v\right) + \gamma^2 \left(-\frac{461}{1512} - \frac{18395}{1512}v - \frac{241}{1512}v^2\right) \\
&+ \gamma^3 \left(\frac{395899}{13200} - \frac{428}{105}\ln\left(\frac{r}{r_0}\right) + \left[\frac{3304319}{166320} - \frac{44}{3}\ln\left(\frac{r}{r_0}\right)\right]v\right. \\
&\left. + \frac{162539}{16632}v^2 + \frac{2351}{33264}v^3\right) \\
&+ \gamma^4 \left(-\frac{1067041075909}{12713500800} + \frac{31886}{2205}\ln\left(\frac{r}{r_0}\right)\right. \\
&\left. + \left[-\frac{85244498897}{470870400} - \frac{2783}{1792}\pi^2 - \frac{128}{7}\gamma_E - \frac{64}{7}\ln(16\gamma)\right]v\right. \\
&\left. - \frac{10886}{735}\ln\left(\frac{r}{r_0}\right) + \frac{8495}{63}\ln\left(\frac{r}{r_0}\right)\right]v \\
&+ \left[\frac{171906563}{4484480} + \frac{44909}{2688}\pi^2 - \frac{4897}{21}\ln\left(\frac{r}{r_0}\right)\right]v^2 \\
&\left. - \frac{22063949}{5189184}v^3 + \frac{71131}{314496}v^4\right),
\end{aligned}
\quad (449a) $$

<footnote>$^{86}$ To 3PN order the mass quadrupole was computed by Blanchet et al. (2002b), Blanchet and Iyer (2005), Blanchet et al. (2005a); to 2PN order it was computed by Blanchet et al. (1995a), Will and Wiseman (1996), Leibovich et al. (2020); to 1PN order by Wagoner and Will (1976), Blanchet and Schäfer (1989).</footnote>
$$
\begin{equation}
\begin{aligned}
B ={}& \frac{11}{21} - \frac{11}{7} v + \gamma \left( \frac{1607}{378} - \frac{1681}{378} v + \frac{229}{378} v^2 \right) \\
& + \gamma^2 \left( -\frac{357761}{19800} + \frac{428}{105} \ln\left(\frac{r}{r_0}\right) - \frac{92339}{5544} v + \frac{35759}{924} v^2 + \frac{457}{5544} v^3 \right) \\
& + \gamma^3 \left( \frac{23006898527}{1589187600} - \frac{4922}{2205} \ln\left(\frac{r}{r_0}\right) \right. \\
& \qquad \left. + \left[ \frac{8431514969}{529729200} + \frac{143}{192} \pi^2 - \frac{64}{7} \gamma_E - \frac{32}{7} \ln(16\gamma) \right] \right. \\
& \qquad \left. - \frac{1266}{49} \ln\left(\frac{r}{r_0}\right) - \frac{968}{63} \ln\left(\frac{r}{r'_0}\right) \right] v \\
& + \left[ \frac{351838141}{5045040} - \frac{41}{24} \pi^2 + \frac{968}{21} \ln\left(\frac{r}{r'_0}\right) \right] v^2 \\
& - \frac{1774615}{81081} v^3 - \frac{3053}{432432} v^4,
\end{aligned}
\tag{449b}
\end{equation}
$$

$$
C = \frac{48}{7} + \gamma \left( -\frac{4096}{315} - \frac{24512}{945} v \right) - \frac{32}{7} \pi \gamma^{3/2}. \quad (449c)
$$

The quadrupole depends on both the scale $r_0$ associated with the Finite Part regularization for the IR, see Eq. (72), and the one $r'_0$ associated with the UV regularization and present in the equations of motion, see Eq. (339). The coefficients $A$ and $B$ represent the conservative part of the quadrupole moment, while $C$ is due to the radiation reaction dissipative effects. As shown by Larrouturou et al. (2022b), the source quadrupole moment contains a non-local tail integral. The 4PN logarithms $\ln(16\gamma)$ in $A$ and $B$ are due to the conservative part of this non-local tail term, while the last term $\propto \pi$ in $C$ corresponds to the dissipative part of the tail term.

Besides the 4PN mass quadrupole, we need also (for the 4PN flux) the mass octupole moment $I_{ijk}$ and current quadrupole moment $J_{ij}$, both of them at the 3PN order. The renormalized versions of these moments coincide at 3PN order with the ones computed with Hadamard regularization for the IR. The octupole mass moment reads (Faye et al. 2015)

$$
\begin{equation}
\begin{split}
I_{ijk} = -\mu\Delta & \left\{ D x_{\langle i} x_{j} x_{k\rangle} + E \frac{r}{c} v_{\langle i} x_{j} x_{k\rangle} + F \frac{r^2}{c^2} v_{\langle i} v_{j} x_{k\rangle} + G \frac{r^3}{c^3} v_{\langle i} v_{j} v_{k\rangle} \right\} \\
& + \mathcal{O}\left(\frac{1}{c^7}\right),
\end{split}
\tag{450}
\end{equation}
$$

where $\Delta = (m_1 - m_2)/m$, and the coefficients are

$$
\begin{equation}
\begin{split}
D ={}& 1 - \gamma v + \gamma^2 \left( -\frac{139}{330} - \frac{11923}{660} v - \frac{29}{110} v^2 \right) \\
& + \gamma^3 \left( \frac{1229440}{63063} + \frac{610499}{20020} v + \frac{319823}{17160} v^2 - \frac{101}{2340} v^3 \right. \\
& \qquad \left. - \frac{26}{7} \ln\left(\frac{r}{r_0}\right) - 22v \ln\left(\frac{r}{r'_0}\right) \right),
\end{split}
\tag{451a}
\end{equation}
$$
$$E = \frac{196}{15} \gamma^2 v, \qquad (451b)$$

$$F = 1 - 2v + \gamma \left( \frac{1066}{165} - \frac{1433}{330} v + \frac{21}{55} v^2 \right) + \gamma^2 \left( -\frac{1130201}{48510} - \frac{989}{33} v + \frac{20359}{330} v^2 - \frac{37}{198} v^3 + \frac{52}{7} \ln\left(\frac{r}{r_0}\right) \right), \qquad (451c)$$

$$G = 0. \qquad (451d)$$

The 3PN current quadrupole moment is (Henry et al. 2021)

$$J_{ij} = -\mu \Delta \left\{ H L^{\langle i} x^{j\rangle} + K \frac{Gm}{c^3} L^{\langle i} v^{j\rangle} \right\} + \mathcal{O}\left(\frac{1}{c^7}\right), \qquad (452)$$

where we denote $L^i \equiv \epsilon_{ijk} x^j v^k$, and where

$$H = 1 + \gamma \left( \frac{67}{28} - \frac{2}{7} v \right) + \gamma^2 \left( \frac{13}{9} - \frac{4651}{252} v - \frac{v^2}{168} \right) \\ + \gamma^3 \left( \frac{2301023}{415800} - \frac{214}{105} \ln\left(\frac{r}{r_0}\right) + \left[ -\frac{243853}{9240} + \frac{123}{128} \pi^2 - 22 \ln\left(\frac{r}{r_0'}\right) \right] v \right. \\ \left. + \frac{44995}{5544} v^2 + \frac{599}{16632} v^3 \right), \qquad (453a)$$

$$K = \frac{188}{35} v \gamma. \qquad (453b)$$

The list of required source moments for the 4PN flux and 3.5PN waveform continues with the 2PN mass 2⁴-pole and current 2³-pole (octupole) moments, and so on. These moments are relatively easy to compute as no regularization subtleties arise. Since the first non-local feature appears at 4PN order, these higher-order source moments are instantaneous. In particular, they cannot contribute to the 4.5PN term in the quasi-circular flux. Here we give the most updated moments (Blanchet et al. 2008; Henry 2023)<fnref n="87" />

<footnote n="87">⁸⁷ The STF projection $\langle \rangle$ applies only on “alive” indices $ijkl...$ but not on the summed up indices $a$ and $b$.</footnote>