The equations of motion for a two point mass binary in harmonic coordinates up to 2.5 PN order, at which the radiation reaction effect first appears, were derived by Damour and Deruelle [52, 46] based on the post-Minkowskian approach . These works used Dirac delta distributions to express the point masses mathematically, therefore they inevitably resorted to a purely mathematical regularization procedure to deal with divergences arising from the nonlinearity of general relativity. Damour  addressed the applicability of the use of a Dirac delta distribution to self-gravitating objects. By investigating the tidal effect exerted by the companion star on the main star, he gave a plausible argument known as the “dominant Schwarzschild condition” which supports that up to 2.5 PN order, the field around the main star is recovered by the energy momentum tensor expressed in terms of a Dirac delta distribution.
Direct validations of the 2.5 PN equations of motion by Damour and Deruelle, where a Dirac delta distribution was not used, have been obtained in several works [87, 95, 113, 130]. Grishchuk and Kopeikin  and Kopeikin  worked on extended but intrinsically spherical bodies with weak internal gravity using the post-Newtonian approximation both inside and outside the stars. They volume-integrated the equations of the conservation of the stress-energy tensor of the matter for an ideal fluid with two compact supports and obtained their equations of motion. The volume integral approach was adopted also by Pati and Will . The present authors and Asada on the other hand derived the 2.5 PN equations of motion  for point particles with arbitrarily strong internal gravity using a regular point particle limit called the strong field point particle limit . These authors also used the local conservation law but adopted a surface integral approach , and they did not specify an explicit form of the stress-energy tensor but assumed that it satisfies some scaling on the initial hypersurface.
Blanchet, Faye, and Ponsot  also derived the 2.5 PN equations of motion using Dirac delta distributions for which Hadamard’s partie finie regularization was employed to handle the divergences due to their use of Dirac delta distributions. In this approach, they have assumed that the two point masses follow regularized geodesic equations. (More precisely, they have assumed that the dynamics of two point masses are described by a regularized action, from which a regularized geodesic equation was shown to be derived.) They also derived the gravitational field up to 2.5 PN order in an explicit form which may help constructing initial data for numerical simulations of compact binaries.
All the works quoted above agree with each other. Namely, our work  shows the applicability of the Damour and Deruelle 2.5 PN equations of motion to a relativistic compact binary consisting of regular stars with strong internal gravity. We mention here that stars consisting of relativistic compact binaries, for which we are searching as gravitational wave sources, have a strong internal gravitational field, and that it is a nontrivial question whether a star follows the same orbit regardless of the strength of its internal gravity.
Currently we have the equations of motion for relativistic compact binaries through the 3.5 PN approximation of general relativity in hand. Actually the 3.5 PN correction to the equations of motion is relatively easily derived [111, 123, 130]. At 3 PN order, an issue on undetermined coefficient associated with the regularization procedures was found which we now briefly discuss.
A 3 PN iteration result was first reported by Jaranowski and Schäfer [98, 99]. There a 3 PN ADM Hamiltonian in the ADM transverse traceless (ADMTT) gauge for two point masses expressed as two Dirac delta distributions was derived based on the ADM canonical approach [125, 135]. However, it was found in [98, 99] that the regularization they had used caused one coefficient to be undetermined in their framework. Moreover, they later found another undetermined coefficient in their Hamiltonian, called . Origins of these two coefficients were attributed to some unsatisfactory features of regularization they had used, such as violation of the Leibniz rule. The former coefficient, which appears as a numerical multiplier of the term that depends on the momenta of the point particles, was then fixed as by a posteriori imposing Poincaré invariance on their 3 PN Hamiltonian . As for the latter coefficient, Damour et al.  succeeded in fixing it as , adopting dimensional regularization1. Moreover, with this method they found the same value of as in , which ensures Lorentz invariance of their Hamiltonian.
On the other hand, Blanchet and Faye have tackled the derivation of the 3 PN equations of motion for two point masses expressed as two Dirac delta distributions in harmonic coordinates [25, 27] based on their previous work . The divergences due to their use of Dirac delta distributions were systematically regularized with the help of Lorentz invariant generalized Hadamard partie finie regularization. They have extended the notion of the Hadamard partie finie regularization to regularize divergent integrals and a singular function which does not permit a power-like expansion near its singularities . Furthermore, the regularization is carefully constructed in  so that it respects Lorentz invariance. Their equations of motion respect the Lorentz invariance in the post-Newtonian pertubative sense and admits a conserved energy of orbital motion modulo the 2.5 PN radiation reaction effect. They found, however, that there exists one and only one undetermined coefficient (which they call ).
Interestingly, the two groups independently constructed a transformation between the two gauges and found that these two results coincide with each other when there exists a relation [55, 64] as  at 3 PN order.
The present authors have derived the 3 PN equations of motion for relativistic compact binaries [91, 92, 93] in harmonic coordinates based on their previous work . Namely, they did not use Dirac delta distributions. As a result, they did not find any undetermined coefficient at all in the equations of motion and found the same value of the parameter as Equation (2). Thus, the issue of the undetermined coefficient problem has been solved.
Physically equivalent 3 PN equations of motion in harmonic coordinates were also completed by Blanchet, Damour, and Esposito-Farèse  based on their previous works [25, 27, 30]. There, they have used the dimensional regularization to overcome the problem with the (generalized) Hadamard partie finie. The physical equivalence between the two results, Blanchet–Damour–Esposito-Farèse’s and our equations of motion, suggests that, at least up to 3 PN order, a particle (with a strong internal gravity) follows a regularized (in some sense) geodesic equation in a dynamical spacetime, a part of whose gravitational field the particle itself generates.
Will and his collaborators [130, 129, 163] have been tackling the 3 PN iteration of the equations of motion where they take into account the density profile of the stars explicitly. Their result may (or may not) give a direct check of the effacement principle  up to 3 PN order which states that the motion of the objects depends only on their masses and not on their internal structures up to the order where the tidal effect comes into play.
This work is licensed under a Creative Commons Attribution-Noncommercial-No Derivative Works 2.0 Germany License.