1 Introduction

1.1 Gravitational wave detection and post-Newtonian approximation

The motion and associated emission of gravitational waves (GW) of self-gravitating systems have been a main research interest in general relativity. The problem is complicated conceptually as well as mathematically because of the nonlinearity of the Einstein equations. There is no hope in any foreseeable future to have exact solutions describing motions of arbitrarily shaped, massive bodies, so we have to adopt some sort of approximation schemes for solving the Einstein equations to study such problems. In the past years many types of approximation schemes have been developed depending on the nature of the system under consideration. Here we shall focus on a particular scheme called the post-Newtonian (PN) approximation. There are many systems in astrophysics where Newtonian gravity is dominant, but general relativistic gravity plays also an important role in their evolution. For such systems it would be nice to have an approximation scheme which gives a Newtonian description in the lowest order and general relativistic effects as higher order perturbations. The post-Newtonian approximation is perfectly suited for this purpose. Historically Einstein computed first the post-Newtonian effects, e.g. the precession of the perihelion [72]. Studies of the post-Newtonian approximation were made by Lorentz and Droste [117], Einstein, Infeld, and Hoffmann [73], Fock [77], Plebansky and Bazanski [131], and Chandrasekhar and associates [40, 41, 42, 43].

Now it is widely known that the post-Newtonian approximation is important in analyzing a number of relativistic problems, such as the equations of motion of binary pulsars [34, 61, 74, 89], solar-system tests of general relativity [159, 160, 161, 164], and gravitational radiation reaction [39, 42]. Any approximation scheme necessitates one or several small parameters characterizing the nature of the system under consideration. A typical parameter which most of the schemes adopt is the magnitude of the metric deviation from a certain background metric. In particular if the background is Minkowski spacetime and there is no other parameter, the scheme is sometimes called the post-Minkowskian approximation in the sense that the constructed spacetime reduces to Minkowski spacetime in the limit that the parameter tends to zero. This limit is called the weak field limit. In the case of the post-Newtonian approximation the background spacetime is also Minkowski spacetime, but there is another small parameter, that is, the typical velocity of the system divided by the speed of light. We introduce a non-dimensional parameter ϵ to express the “slowness” of the system. These two parameters (the deviation from the flat metric and the velocity) have to have a certain relation in the following sense. As the gravitational field gets weaker, all velocities and forces characteristic of the material systems become smaller, in order to permit the weakening of gravity to remain an important effect in the system’s dynamics. For example in the case of a binary system, the typical velocity would be the orbital velocity υ/cϵ and the deviation from the flat metric would be the Newtonian potential, say Φ. Then these are related by Φ/c2υ2/c2ϵ2 which guarantees that the system is bounded by its own gravity.

In the post-Newtonian approximation, the equations of general relativity take the form of Newton’s equations in an appropriate limit as ϵ → 0. Such a limit is called the Newtonian limit and it will be the basis of constructing the post-Newtonian approximation. However, the limit is not in any sense trivial since it may be thought of as two limits tied together as just described. It is also worth noting that the Newtonian limit cannot be uniform everywhere for all time. For example any compact binary system, no matter how weak the gravity between its components and slow its orbital motion is, will eventually spiral together due to backreaction from the emission of gravitational waves. As the result the effects of its Newtonian gravity will be swamped by those of its gravitational waves. This will mean that higher order effects of the post-Newtonian approximation eventually dominate the lowest order Newtonian dynamics and thus if the post-Newtonian approximation is not carefully constructed, this effect can lead to many formal problems, such as divergent integrals [71]. It has been shown that such divergences may be avoided by carefully defining the Newtonian limit [79]. Moreover, the use of such a limit provides us a strong indication that the post-Newtonian hierarchy is an asymptotic approximation to general relativity [82]. Therefore we shall first discuss in this paper the Newtonian limit and how to construct the post-Newtonian hierarchy before attacking practical problems in later sections.

Before going into the details, we mention the reason for the growth of interest in the post-Newtonian approximation in recent years. Certainly the discovery of the binary neutron star system PSR 1913+16 was a strong reason to have renewed interest in the post-Newtonian approximation, since it is the first system found in which general relativistic gravity plays a fundamental role in its evolution [89]. Particularly the indirect discovery of gravitational waves by the observation of the period shortening led to many fruitful studies of the equations of motion with gravitational radiation emission in 1980s (see [47, 48, 49] for a review). The effect of radiation reaction appears in the form of a potential force at the order of ϵ5 higher than the Newtonian force in the equations of motion. Ehlers and colleagues [71] critically discussed the foundation of the so-called quadrupole formula for the radiation reaction (see also the introduction of [129]). There have been various attempts to show the validity of the quadrupole formula [5, 47, 79, 90, 104, 105, 106, 107, 135, 138, 139, 155, 157, 156]. Namely, Damour [46] proved the formula for compact binaries with the help of the “dominant Schwarzschild condition” [47]. Blanchet and Damour [21] proved it for general fluid systems.

At that time, however, no serious attempts with direct detection of gravitational waves in mind had been made for the study of higher order effects in the equations of motion. The situation changed gradually in the late 1980s because of the increasing expectation of a direct detection of gravitational waves by kilometer-size interferometric gravitational wave detectors, such as LIGO [1, 116], VIRGO [35, 153], GEO [62, 83], and TAMA [114, 149]. Coalescing binary neutron stars are the most promising candidates of sources of gravitational waves for such detectors. The reasons are that (i) we expect, say, the (initial) LIGO to detect the signal of coalescence of binary neutron stars about once per year to once per hundreds of years [38, 103, 102], and (ii) the waveform from coalescing binaries can be predicted with high accuracy compared to other sources [1, 151, 161]. Information carried by gravitational waves tells us not only various physical parameters of neutron stars [45], but also the cosmological parameters [75, 119, 141, 142, 158] if and only if we can make a detailed comparison between the observed signal with theoretical predictions during the epoch of the so-called inspiraling phase where the orbital separation is much larger than the radius of the component stars [44]. This is the place where the post-Newtonian approximation may be applied to make theoretical templates for gravitational waves. The problem is that in order to make any meaningful comparison between theory and observation we need to know the detailed waveforms generated by the motion up to, say, 4 PN order which is of order ϵ8 higher than the Newtonian order [2, 3, 10, 147]. This request from gravitational wave astronomy forces us to construct higher order post-Newtonian equations of motion and waveform templates.

Replying to this request, there have been various works studying the equations of motion for a compact binary system and developing higher order post-Newtonian gravitational waveform templates for such a system. The most systematic among those works that have succeeded in achieving higher order iteration are the ones by Blanchet, Damour, and Iyer who have developed a scheme to calculate the waveform at a higher order, where the post-Minkowskian approximation is used to construct the external field and the post-Newtonian approximation is used to construct the field near the material source. They and their collaborators have obtained the waveform up to 3.5 PN order which is of order ϵ7 higher than the lowest quadrupole wave [23, 24, 29, 31, 32] by using the equations of motion up to that order [22, 91, 93, 111, 123, 130]. The 3.5 PN waveform includes tail terms which manifest nonlinearity of general relativity. Blanchet and Schäfer [33] have investigated a spectral (Fourier) decomposition of the tail and computed the contribution of the tail to the gravitational wave luminosity emitted by a binary system having a general eccentric orbit. Asada and Futamase [12] showed that the dominant part of the tail term originates from the phase shift of the wave due to the Coulomb part of the gravitational field.

In this paper, we mainly discuss the foundation of the Newtonian limit and the post-Newtonian equations of motion for relativistic compact binaries in an inspiralling phase. In the next Section 1.2 we briefly give an historical introduction on the latter topic.

1.2 Post-Newtonian equations of motion

Many authors have derived equations of motion up to the 1 PN order [66, 73, 76, 94, 117, 152], up to 2.5 PN order [30, 47, 50, 51, 52, 46, 87, 95, 113, 130], and up to 3 PN order [22, 91, 93]. The 3.5 PN correction to the equations of motion has been derived by [111, 123, 130]. The post-Newtonian equations of motion are now available in harmonic coordinates up to 3.5 PN order inclusively. See also [17, 85, 86, 96] for the 3.5 PN and the 4.5 PN correction based on the energy balance. Besides the equations of motion, attempts have been made to derive a Hamiltonian for a binary system. The second order post-Newtonian computation of the Hamiltonian was tackled by [125, 126, 127] and completed in [56, 135, 136, 137]. Damour, Jaranowski, and Schä;fer have completed the 3 PN order Hamiltonian in [54]. So far, the post-Newtonian Arnowitt-Deser-Misner (ADM) Hamiltonian is available in ADM transverse-traceless coordinates up to 3.5 PN order inclusively [97].

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 [14]. 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 [47] 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 [87] and Kopeikin [113] 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 [130]. The present authors and Asada on the other hand derived the 2.5 PN equations of motion [95] for point particles with arbitrarily strong internal gravity using a regular point particle limit called the strong field point particle limit [81]. These authors also used the local conservation law but adopted a surface integral approach [73], 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 [30] 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 [95] 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 ωkinetic to be undetermined in their framework. Moreover, they later found another undetermined coefficient in their Hamiltonian, called ωstatic [100]. 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 ωkinetic = 41/24 by a posteriori imposing Poincare invariance on their 3 PN Hamiltonian [53]. As for the latter coefficient, Damour et al. [54] succeeded in fixing it as ωstatic = 0, adopting dimensional regularizationFootnote 1. Moreover, with this method they found the same value of ωkinetic as in [53], 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 [30]. 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 [26]. Furthermore, the regularization is carefully constructed in [28] 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]

$${\omega _{{\rm{static}}}} = - {{11} \over 3}\lambda - {{1987} \over {840}}.$$
(1)

Therefore, it is possible to fix the λ parameter from the result of [54] as

$$\lambda = - {{1987} \over {3080}}.$$
(2)

However, the applicability of mathematical regularization to the current problem is not a trivial issue, but an assumption to be verified, or at least supported by convincing arguments. There is no argument such as the “dominant Schwarzschild condition” [47] 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 [95]. 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 [22] 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 [47] 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.

1.3 Plan of this paper

This article is organized as follows: In Section 2 we introduce the Newtonian limit in general relativity and present how to construct the post-Newtonian hierarchy. There we mention how to avoid divergent integrals which appear at a higher order in the previous treatments. We shall then give a formulation of the dynamics of a Newtonian perfect fluid as the lowest order post-Newtonian approximation in harmonic coordinates.

We focus our attention on the derivation of the 3 PN equations of motion based on the work by the authors of this article from Section 3 to Section 8.

In Section 3, we discuss the basics of our derivation of the post-Newtonian equations of motion for relativistic compact binaries. There we discuss how to incorporate strong internal gravity in the post-Newtonian approximation.

In Section 4, we formulate our method in more detail. As an application, we derive the Newtonian and the 1 PN equations of motion in our formalism.

In Section 5, we briefly explain how to derive the 3 PN gravitational field when possible, and how to derive equations of motion when the gravitational field is unavailable in an explicit closed form.

In Sections 6 and 7, we derive the mass-energy relation and the momentum-velocity relation up to 3 PN order, and finally in Section 8, we derive the 3 PN equations of motion. The resulting equations of motion respect Lorentz invariance, admit a conserved energy (modulo the 2.5 PN radiation reaction effect), and are free from any ambiguity. Section 8 ends with a summary and some remarks on possible future study of the 4 PN order iteration.

In Appendix A, we give a brief sketch of the direct integration of the relaxed Einstein equations (DIRE) method by Will and Wiseman [129, 162, 165] which we have used in this article.

Although in the main body of this article we focus our attention on spherical massive bodies (or point particles), we shall apply our formalism to extended bodies in Appendix B.

In Appendix C, using the strong field point particle limit and the surface integral approach, we discuss that a particle with strong internal gravity moves on a geodesic of some smooth metric part which is produced by the particle itself. This work is done by the authors in the collaboration with Takashi Fukumoto.

Throughout this article, we will use units where G = c = 1 unless otherwise mentioned.

2 Foundation of the Post-Newtonian Approximation

Since the Newtonian limit is the basis of the post-Newtonian approximation, we shall first formulate the Newtonian limit. We follow the formulation by Futamase and Schutz [82]. We will not mention other formulations of the post-Newtonian approximation [63, 70, 134].

2.1 Newtonian limit along a regular asymptotic Newtonian sequence

This formulation is based on the observation that any asymptotic approximation of any theory needs a sequence of solutions of the basic equations of the theory [146, 154]. Namely, if we write the equations in abstract form as

$$E(g) = 0$$
(3)

for an unknown function g, one would like to have a one-parameter (or possibly multi-parameter) family of solutions,

$$E(g(\lambda)) = 0,$$
(4)

where λ is some parameter. Asymptotic approximation then says that a function f(λ) approximates g(λ) to order λp if ∣ f(λ) − g(λ)∣/λp → 0 as λ → 0. We choose the sequence of solutions with appropriate properties in such a way that the properties reflect the character of the system under consideration.

We shall formulate the post-Newtonian approximation according to the general idea just described. As stated in the introduction, we would like to have an approximation which applies to the systems whose motions are described almost by Newtonian theory. Thus we need a sequence of solutions of the Einstein equations parameterized by ϵ (the typical velocity of the system divided by the speed of light) which has Newtonian character as ϵ → 0.

The Newtonian character is most conveniently described by the following scaling law. The Newtonian equations involve six variables, namely the density ρ, the pressure P, the gravitational potential Φ, and the velocity υi, i = 1, 2, 3):

$${\nabla ^2}\Phi - 4\pi \rho = 0,$$
(5)
$${\partial _t}\rho + {\nabla _i}(\rho {v^i}) = 0,$$
(6)
$$\rho {\partial _t}{v^i} + \rho {v^j}{\nabla _j}{v^i} + {\nabla ^i}P + \rho {\nabla ^i}\Phi = 0,$$
(7)

supplemented by an equation of state. For simplicity we have considered a perfect fluid.

It can be seen that the variables {ρ(xi, t), P(xi,t), Φ(xi, t), υi(xj, t)} obeying the above equations satisfy the following scaling law:

$$\begin{array}{*{20}c} {\rho ({x^i},t) \rightarrow {\epsilon ^2}\rho ({x^i},\epsilon t),} \\ {P({x^i},t) \rightarrow {\epsilon ^4}P({x^i},\epsilon t),} \\ {{v^i}({x^k},t) \rightarrow \epsilon {v^i}({x^k},\epsilon t),} \\ {\Phi ({x^i},t) \rightarrow {\epsilon ^2}\Phi ({x^i},\epsilon t).} \\ \end{array}$$
(8)

One can easily understand the meaning of this scaling by noticing that ϵ is the magnitude of the typical velocity (divided by the speed of light). Then the magnitude of the gravitational potential will be of order ϵ2 because of the balance between gravity and the centrifugal force. The scaling of the time variable expresses the fact that the weaker gravity is (ϵ → 0) the longer the time scale is.

Thus we wish to have a sequence of solutions of the Einstein equations which has the above scaling as ϵ → 0. We shall also take the point of view that the sequence of solutions is determined by the appropriate sequence of initial data. This has a practical advantage because there will be no solutions of the Einstein equations which satisfy the above scaling (8) even as ϵ → 0. This is because the Einstein equations are nonlinear in the field variables, so it will not be possible to enforce the scaling everywhere in spacetime. We shall therefore impose it only on the initial data for the solution of the sequence.

Here we first give a general discussion on the formulation of the post-Newtonian approximation independent of any initial value formalism and then present the concrete treatment in harmonic coordinates. The condition is used because of its popularity and some advantages in the generalization to the systems with strong internal gravity.

As the initial data for the matter we take the same data set in the Newtonian case, namely the density ρ, the pressure P, and the coordinate velocity υi. In most of the application, we usually assume a simple equation of state which relates the pressure to the density. The initial data for the gravitational field are gμν, ∂gμν/∂t. Since general relativity is an overdetermined system, there will be constraint equations among the initial data for the field. We shall write the free data for the field as (Qij, Pij) whose explicit forms depend on the coordinate condition one assumes. In any coordinates we shall assume these data for the field vanish since we are interested in the evolution of an isolated system by its own gravitational interaction. It is expected that this choice corresponds to the absence of radiation far away from the source. Thus we choose the following initial data which is indicated by the Newtonian scaling:

$$\begin{array}{*{20}c} {\rho (t = 0,{x^i},\epsilon) = {\epsilon ^2}a({x^i}),} \\ {P(t = 0,{x^i},\epsilon) = {\epsilon ^4}b({x^i}),} \\ {{v^i}(t = 0,{x^k},\epsilon) = \epsilon {c^i}({x^k}),} \\ {{Q_{ij}}(t = 0,{x^i},\epsilon) = 0,\quad \quad \,\,\,\,} \\ {{P_{ij}}(t = 0,{x^i},\epsilon) = 0,\quad \quad \,\,\,} \\ \end{array}$$
(9)

where the functions a, b, and ci are C functions that have compact support contained entirely within a sphere of a finite radius.

Corresponding to the above data, we have a one-parameter set of spacetime parameterized by ϵ. It may be helpful to visualize the set as a fiber bundle, with base space R being the real line (coordinate ϵ) and fiber R4 being the spacetime (coordinates t, xi). The fiber ϵ = 0 is Minkowski spacetime since it is defined by zero data. In the following we shall assume that the solutions are sufficiently smooth functions of ϵ for small ϵ ≠ 0. We wish to take the limit ϵ → 0 along the sequence. The limit is, however, not unique and is defined by giving a smooth nowhere vanishing vector field on the fiber bundle which is nowhere tangent to each fiber [84, 146]. The integral curves of the vector field give a correspondence between points in different fibers, namely events in different spacetimes with different values of ϵ. Remembering the Newtonian scaling of the time variable in the limit, we introduce the Newtonian dynamical time,

$$\tau = \epsilon t,$$
(10)

and define the integral curve as the curve on which τ and xi stay constantFootnote 2. In fact if we take the limit ϵ → 0 along this curve, the orbital period of the binary system with ϵ = 0.01 is 10 times that of the system with ϵ = 0.1 as expected from the Newtonian scaling. This is what we define as the Newtonian limit. Notice that this map never reaches the fiber ϵ = 0 (Minkowski spacetime). There is no pure vacuum Newtonian limit as expected.

In the following we assume the existence of such a sequence of solutions constructed by the initial data satisfying the above scaling with respect to ϵ. We shall call such a sequence a regular asymptotically Newtonian sequence. We have to make further mathematical assumptions about the sequence to make explicit calculations. We will not go into details partly because in order to prove the assumptions we need a deep understanding of the existence and uniqueness properties of the Cauchy problem of the Einstein equations with perfect fluids of compact support which are not available at present.

2.2 Post-Newtonian hierarchy

We shall now define the Newtonian, post-Newtonian, and higher approximations of various quantities as the appropriate higher tangents of the corresponding quantities to the above integral curve at ϵ = 0. For example the hierarchy of approximations for the spacetime metrics can be expressed as follows:

$$\begin{array}{*{20}c} {{g_{\mu \nu}}(\epsilon, \tau, {x^i}) = {g_{\mu \nu}}(0, \tau, {x^i}) + \epsilon ({{\mathcal L}_V}{g_{\mu \nu}})(0, \tau, {x^i})\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {+ {1 \over 2}{\epsilon ^2}({\mathcal L}_V^2{g_{\mu \nu}})(0,\tau,{x^i}) + \cdots + {{{\epsilon ^n}} \over {n!}}({\mathcal L}_V^n{g_{\mu \nu}})(0,\tau,{x^i}) + {R_{n + 1}},} \\ \end{array}$$
(11)

where \({\mathcal L_V}\) is the Lie derivative with respect to the tangent vector of the curves defined above, and the remainder term \(R_{n + 1}^{\mu \nu}\) is

$$R_{n + 1}^{\mu \nu} = {{{\epsilon ^{n + 1}}} \over {(n + 1)!}}\int\nolimits_0^1 d \ell \,{(1 - \ell)^{n + 1}}({\mathcal L}_V^{n + 1}{g_{\mu \nu}})(\ell \epsilon, \tau, {x^i}).$$
(12)

Taylor’s theorem guarantees that the series is an asymptotic expansion about ϵ = 0 under certain assumptions mentioned above. It may be useful to point out that the above definition of the approximation scheme may be formulated purely geometrically in terms of a jet bundle.

The above definition of the post-Newtonian hierarchy gives us an asymptotic series in which each term in the series is manifestly finite. This is based on the ϵ dependence of the domain of dependence of the field point (τ xk). The region is finite with finite values of ϵ, and the diameter of the region increases like ϵ−1 as ϵ → 0. Without this linkage of the region with the expansion parameter ϵ, the post-Newtonian approximation leads to divergences in the higher orders. This is closely related to the retarded expansion. Namely, it is assumed that the slow motion assumption enables one to Taylor expand the retarded integrals in retarded time such as

$$\int d r\,f(\tau - \epsilon r, \ldots) = \int d r\,f(\tau, \ldots) - \epsilon \int d r\,rf{(\tau, \ldots)_{,\tau}} + \ldots,$$
(13)

and assign the second term to a higher order because of its explicit ϵ in front. This is incorrect because rϵ−1 as ϵ → 0 and thus ϵr is not uniformly small in the Newtonian limit. Only if the integrand falls off sufficiently fast, the retardation can be ignored. This happens in the lower order PN terms. But at some higher order there appear many terms which do not fall off sufficiently fast because of the nonlinearity of the Einstein equations. This is the reason that the formal PN approximation produces the divergent integrals. It turns out that such a divergence appears at 4 PN order, indicating a breakdown of the PN approximation in harmonic coordinates [20]Footnote 3. This sort of divergence may be eliminated if we remember that the upper bound of the integral does depend on ϵ as ϵ−1. Thus we would get something like ϵn ln ϵ instead of ϵn ln ∞ in the usual approach. This shows that the asymptotic Newtonian sequence is not differentiable in ϵ at ϵ = 0, but there is no divergence in the expansion and it has still an asymptotic approximation in ϵ that involves logarithms.

Other than the initial value formulation method [79, 82] mentioned above, various methods have been proposed to solve this problem of the divergent integrals. It is known that a higher order post-Newtonian metric does not respect the asymptotically flat condition. This does not mean that the post-Newtonian approximation is useless at such a high order. The problem is related to the fact that a simple post-Newtonian iteration is meaningful only in the near zone — about one wavelength distance away from the material source — and is not useful outside of the near zone, called far zone, where the wave effect (retardation effect) is manifest. So roughly speaking, if a far zone metric satisfying proper boundary conditions at infinity is solved so that we have a boundary condition to the field equations for a post-Newtonian metric in the buffer zone, we can find a post-Newtonian metric which is meaningful in the sense that it respects the correct behaviour at the near zone boundary.

Blanchet and Damour have developed a systematic approach of a matched asymptotic expansion. They solved the far zone metric using a multipolar post-Minkowskinan expansion. The far zone metric satisfies a stationarity condition and is parametrized by radiative multipole moments. On the other hand, they solve a post-Newtonian near zone metric up to a homogeneous solution. They then establish an association between those radiative multipole moments and the source multipole moments that characterize a post-Newtonian near zone metric to fix the homogeneous solution, and find a post-Newtonian metric which satisfies the correct behaviour in the buffer zone.

Will and Wiseman have developed the DIRE method [129, 162, 165] where they split the integral region in the retarded integral into two — one being the near zone and the other being the far zone. The near zone metric is solved by a post-Newtonian expansion. The retarded integral over the far zone is directly evaluated with the assumption of sufficient stationarity of the system in the infinite past.

In fact, both the Blanchet-Damour method and the Will-Wiseman method are proved to give a physically equivalent result [19]. In this paper, for our computation of the 3 PN equations of motion, we will use the Will and Wiseman method to solve the problem of the breakdown of the post-Newtonian approximation in the near zone.

2.3 Explicit calculation in harmonic coordinates

Here we shall use the above formalism to make an explicit calculation in harmonic coordinates. The reduced Einstein equations in the harmonic condition are written as

$${\tilde g^{\alpha \beta}}{\tilde g^{\mu \nu}}{\,_{,\alpha \beta}} = 16\pi {\Theta ^{\mu \nu}} - {\tilde g^{\mu \alpha}}{\,_{,\beta}}{\tilde g^{\nu \beta}}{\,_{,\alpha}},$$
(14)
$${\partial _\mu}[{\tilde g^{\mu \nu}}{\partial _\nu}{x^\alpha}] = 0,$$
(15)

where

$${\tilde g^{\mu \nu}} = {(- g)^{1/2}}{g^{\mu \nu}},$$
(16)
$${\Theta ^{\alpha \beta}} = (- g)({T^{\alpha \beta}} + t_{{\rm{LL}}}^{\alpha \beta}),$$
(17)

where \(t_{LL}^{\mu \nu}\) is the Landau-Lifshitz pseudotensor [115]. In this section we shall choose an isentropic perfect fluid for Tαβ which is enough for most applications,

$${T^{\alpha \beta}} = (\rho + \rho \Pi + P){u^\alpha}{u^\beta} + P{g^{\alpha \beta}},$$
(18)

where ρ is the rest mass density, II the internal energy, P the pressure, and uμ the four-velocity of the fluid with normalization

$${g_{\alpha \beta}}{u^\alpha}{u^\beta} = - 1.$$
(19)

The conservation of energy and momentum is expressed as

$${\nabla _\beta}{T^{\alpha \beta}} = 0.$$
(20)

Defining the gravitational field variable as

$${h^{\mu \nu}} = {\eta ^{\mu \nu}} - {(- g)^{1/2}}{g^{\mu \nu}},$$
(21)

where ημν is the Minkowski metric, the reduced Einstein equations (14) and the gauge condition (15) take the following form:

$$({\eta ^{\alpha \beta}} - {h^{\alpha \beta}}){h^{\mu \nu}}_{,\alpha \beta} = - 16\pi {\Theta ^{\mu \nu}} + {h^{\mu \alpha}}_{,\beta}{h^{\nu \beta}}_{,\alpha},$$
(22)
$${h^{\mu \nu}}_{,\nu} = 0.$$
(23)

Thus the characteristics are determined by the operator (ηαβhαβ)αβ, and thus the light cone deviates from that in the flat spacetime. We may use this form of the reduced Einstein equations in the calculation of the waveform far away from the source because the deviation plays a fundamental role there [12]. However, in the study of the gravitational field near the source it is not necessary to consider the deviation of the light cone from the flat one and thus it is convenient to use the following form of the reduced Einstein equations [5]:

$${\eta ^{\mu \nu}}{h^{\alpha \beta}}_{,\mu \nu} = - 16\pi {\Lambda ^{\alpha \beta}},$$
(24)

where

$${\Lambda ^{\alpha \beta}} = {\Theta ^{\alpha \beta}} + {\chi ^{\alpha \beta \mu \nu}}_{,\mu \nu},$$
(25)
$${\chi ^{\alpha \beta \mu \nu}} = {(16\pi)^{- 1}}({h^{\alpha \nu}}{h^{\beta \mu}} - {h^{\alpha \beta}}{h^{\mu \nu}}).$$
(26)

Equations (23) and (24) together imply the conservation law

$${\Lambda ^{\alpha \beta}}_{,\beta} = 0.$$
(27)

We shall take as our variables the set {ρ, P, υi, hαβ}, with the definition

$${v^i} = {u^i}/{u^0}.$$
(28)

The time component of four-velocity u0 is determined from Equation (19). To make a well-defined system of equations we must add the conservation law for the number density n, which is some function of the density ρ and pressure P:

$${\nabla _\alpha}(n{u^\alpha}) = 0.$$
(29)

Equations (27) and (29) imply that the flow is adiabatic. The role of the equation of state is played by the arbitrary function n(ρ, p).

Initial data for the above set of equations are hαβ, hαβ,0, ρ, P, and υi, but not all these data are independent because of the existence of the constraint equations. Equations (23) and (24) imply the four constraint equations among the initial data for the field,

$$\Delta {h^{\alpha 0}} + 16\pi {\Lambda ^{\alpha 0}} - {\delta ^{ij}}h_{i,j}^\alpha {\,^{,0}} = 0,$$
(30)

where Δ is the Laplacian in the flat space. We shall choose hij and hij,0 as free data and solve Equation (30) for hα0(α = 0, …, 3) and Equation (23) for hα0,0. Of course these constraints cannot be solved explicitly, since Λα0 contains hα0, but they can be solved iteratively as explained below. As discussed above, we shall assume that the free data hij and hij,0 for the field vanish. One can show that such initial data satisfy the O’Murchadha and York criterion for the absence of radiation far away from the source [124].

In the actual calculation, it is convenient to use an expression with explicit dependence of ϵ. The harmonic condition allows us to have such an expression in terms of the retarded integral,

$${h^{\mu \nu}}(\epsilon, \tau, {x^i}) = 4\int\nolimits_{C(\epsilon, \tau, {x^i})} {{d^3}} y\,{\Lambda ^{\mu \nu}}(\tau - \epsilon r,{y^i},\epsilon)/r + h_{\rm{H}}^{\mu \nu}(\epsilon, \tau, {x^i}),$$
(31)

where r = ∣yixi∣ and C(ϵ, τ, xi) is the past flat light cone of the event (τ, xi) in the spacetime given by ϵ, truncated where it intersects with the initial hypersurface τ = 0. \(h_{\rm{H}}^{\mu \nu}\) is the unique solution of the homogeneous wave equation in the flat spacetime,

$$\square h_{\rm{H}}^{\mu \nu} = {\eta ^{\alpha \beta}}h_H^{\mu \nu}{\,_{,\alpha \beta}} = 0.$$
(32)

\(h_{\rm{H}}^{\mu \nu}\) evolves from a given initial data on the τ = 0 initial hypersurface which are subject to the constraint equations (30). The explicit form of \(h_{\rm{H}}^{\mu \nu}\) is available via the Poisson formula (see e.g. [144]),

$$h_{\rm{H}}^{\mu \nu}(\tau, {x^i}) = {\tau \over {4\pi}}\oint\nolimits_{\partial C(\tau, {x^i})} {{h^{\mu \nu}}_{,\tau}} (\tau = 0,{y^i})\,d{\Omega _y} + {1 \over {4\pi}}{\partial \over {\partial \tau}}\left[ {\tau \oint\nolimits_{\partial C(\tau, {x^i})} {{h^{\mu \nu}}} (\tau = 0,{y^i})\,d{\Omega _y}} \right].$$
(33)

We shall henceforth ignore the homogeneous solutions because they play no important role. Because of the ϵ dependence of the integral region, the domain of integral is finite as long as ϵ ≠ 0 and their diameter increases like ϵ−1 as ϵ → 0.

Given the formal expression (31) in terms of initial data (9), we can take the Lie derivative and evaluate these derivatives at ϵ = 0. The Lie derivative is nothing but a partial derivative with respect to ϵ in the coordinate system for the fiber bundle given by (ϵ, τ, xi). Accordingly one should convert all the time indices to τ indices. For example, Tττ = ϵ2Ttt which is of order ϵ4, since Tttρ is of order ϵ2. Similarly the other components of stress-energy tensor Tτi = ϵTti and Tij are of order ϵ4 as well. Thus we expect that the first nonvanishing derivative in Equation (31) will be the forth derivative. In fact we find

$$_4{h^{\tau \tau}}(\tau, {x^i}) = 4\int\nolimits_{{R^3}} {{{2\rho (\tau, {y^i})} \over r}} \,{d^3}y,$$
(34)
$$_4{h^{\tau i}}(\tau, {x^i}) = 4\int\nolimits_{{R^3}} {{{2\rho {{(\tau, {y^k})}_1}{v^i}(\tau, {x^k})} \over r}} \,{d^3}y,$$
(35)
$$_4{h^{ij}}(\tau, {x^k}) = 4\int\nolimits_{{R^3}} {{{2\rho {{(\tau, {y^k})}_1}{v^i}{{(\tau, {y^k})}_1}{v^j}(\tau, {y^k}){+ _4}t_{{\rm{LL}}}^{ij}(\tau, {y^k})} \over r}} \,{d^3}y,$$
(36)

where we have adopted the notation

$$_nf(\tau, {x^i}) = {1 \over {n!}}\underset {\epsilon \rightarrow 0} {\lim} {{{\partial ^n}} \over {\partial {\epsilon ^n}}}f(\epsilon, \tau, {x^i}),$$
(38)

and

$$_4t_{{\rm{LL}}}^{ij} = {1 \over {64\pi}}\left({{}_4{h^{\tau \tau, i}}{}_4{h^{\tau \tau, j}} - {1 \over 2}{\delta ^{ij}}_4{h^{\tau \tau, k}}_4{h^{\tau \tau}}_{,k}} \right).$$
(39)

In the above calculation we have taken the point of view that hμν is a tensor field, defined by giving its components in the assumed harmonic coordinates as the difference between the tensor density \(\sqrt {- g} {g^{\mu \nu}}\) and ημν.

The conservation law (27) also has its first nonvanishing derivatives at this order, which are

$$_2{\rho _{,\tau}} + {{(_2}{\rho _1}{v^i})_{,i}} = 0,$$
(40)
$${{(_2}{\rho _1}{v^i})_{,\tau}} + {{(_2}{\rho _1}{v^i}_1{v^j})_{,j}}{+ _4}{P^{,i}} - {1 \over 4}{\,_2}{\rho _4}{h^{\tau \tau, i}} = 0.$$
(41)

Equations (34), (40), and (41) constitute Newtonian theory of gravity. Thus the lowest nonvanishing derivative with respect to ϵ is indeed Newtonian theory, and the 1 PN and 2 PN equations emerge from the sixth and eighth derivatives, respectively, in the conservation law (27). At the next derivative, the quadrupole radiation reaction term emerges.

3 Post-Newtonian Equations of Motion for Compact Binaries

In this section we review briefly the strong field point particle limit, the surface integral approach for the evaluation of the gravitational force, and scalings of matter and field variables on an initial hypersurface. Also we revisit some elementary but useful equations, the Newtonian velocity momentum relation and the Newtonian equations of motion for extended bodies. These are the basic ideas of our formalism. On the base of our formalism, see also [11, 91, 92, 93, 94, 95, 140]. On the post-Newtonian approximation and related issues from different viewpoints, see [47, 48, 49, 16, 18, 19] and references therein.

3.1 Strong field point particle limit

In a relativistic compact binary system in an inspiralling phase, each star is well approximated by a point particle with a few low order multipole moments. This is because in the inspiralling phase, the stellar size divided by the orbital separation is much smaller than unity.

If we wish to apply the post-Newtonian approximation to the inspiraling phase of binary neutron stars, the strong internal gravity must be taken into account. The usual post-Newtonian approximation explicitly assumes the weakness of the gravitational field everywhere including inside the material source. It is argued by applying the strong equivalence principle that the external gravitational field which governs the orbital motion of the binary system is independent of the internal structure of the components up to tidal interaction. Thus it is expected that the results obtained under the assumption of weak gravity also apply for the case of a neutron star binary. Experimental evidence for the strong equivalence principle is obtained only for systems with weak gravity [159, 160, 161, 164], but at present no experiment is available in a case with strong internal gravity.

In the theoretical aspect, the theory of extended objects in general relativity [69] is still in a preliminary stage for an application to realistic systems. The matched asymptotic expansion technique has been used to treat a system with strong gravity in certain situations [47, 65, 66, 104, 105]. Another way to handle strong internal gravity is by the use of a Dirac delta distribution type source with a fixed mass [73]. However, this makes the Einstein equations mathematically meaningless because of their nonlinearity. Physically, there is no such source in general relativity because of the existence of black holes. Before a body shrinks to a point, it forms a black hole whose size is fixed by its mass. For this reason, it has been claimed that no point particle exists in general relativity.

This conclusion is not correct, however. We can shrink the body keeping the compactness (M/R), i.e. the strength of the internal field fixed. Namely we should scale the mass M just like the radius R. This can be fitted nicely into the concept of the regular asymptotic Newtonian sequence defined in Section 2 because there the mass also scales along the sequence of solutions. In fact, if we take the masses of two stars as M, and the separation between two stars as L, then Φ ∼ M/L. Thus the mass M scales as ϵ2 if we fix the separationFootnote 4. In the above we have assumed that the density scales as ϵ2 to guarantee this scaling for the mass while keeping the size of the body fixed. Now we shrink the size as ϵ2 to keep the compactness of each component. Then the density should scale as ϵ−4. We shall call such a scheme the strong field point particle limit since the limit keeps the strength of internal gravityFootnote 5. The above consideration suggests the following initial data to define a regular asymptotic Newtonian sequence which describes a nearly Newtonian system with strong internal gravity [81]. The initial data are two uniformly rotating fluids with compact spatial support whose stress-energy tensor and size scales as ϵ−4 and ϵ2, respectively. We also assume that each of these fluid configurations would be a stationary equilibrium solution of the Einstein equations if the other were absent. This is necessary for the suppression of irrelevant internal motions of each star. Any remaining motions are tidal effects caused by the other body, which will be of order ϵ6 smaller than the internal self-force. These data allow us to use the Newtonian time τ = ϵt as a natural time coordinate everywhere including the interior region of the stars.

As for multipole moments, the scaling \(R = \mathcal O({\epsilon ^2})\) enables us to incorporate the multipole expansion of the stars into the post-Newtonian approximation. In inspiralling compact binaries the tidally and rotationally induced quadrupole moments are too small to affect the binary orbital motion. The time to coalesce is too short for the binary to be tidally locked and corotate [15, 44, 110]. The phase shift due to the quadrupole-orbit couplings may be negligible in the LIGO bandwidth [44]. However, to detect gravitational waves from inspiralling binaries, we have to have highly accurate prior knowledge about the binary orbital motion, say, 4 PN equations of motion or so. The effect of quadrupole moments on the orbital motion can be of about the same order as that of spin-spin interactions [132], while the spin-spin interactions appear at about 2.5 PN order for slowly rotating stars. Also, at the late inspiralling phase, an effect of extendedness of the stars on the motion will be important. Thus it is important to take the multipole moments of the stars into account in a way that is suitable for compact stars when we derive the equations of motion for an inspiralling compact binary.

3.2 Surface integral approach and body zone

One way to evaluate the gravitational force acting on a star is the volume integral approach. In the Newtonian case, the gravitational force \(F_1^i\) on star 1 becomes

$$F_1^i = \int\nolimits_{{B_1}} {{d^3}} x\,\rho {{\partial \Phi} \over {\partial {x^i}}}.$$
(42)

The integral region B1 covers star 1 but does not cover the star 2 (the companion star). To evaluate the above integral we must know the internal structure, which is generally a difficult task even in the Newtonian dynamics, needless to say in general relativity. By means of the surface integral approach we can put off dealing with the internal structure problem until tidal effects affect the orbital motion. In the Newtonian case, using the Poisson equation and Gauss’s law, we can rewrite the above volume integral into a surface integral,

$$F_1^i = - \oint\nolimits_{\partial {B_1}} d {S_j}\,{t^{ij}},$$
(43)
$${t^{ij}} \equiv {1 \over {4\pi}}\left({{{\partial \Phi} \over {\partial {x^i}}}{{\partial \Phi} \over {\partial {x^j}}} - {{{\delta ^{ij}}} \over 2}{{\partial \Phi} \over {\partial {x^k}}}{{\partial \Phi} \over {\partial {x^k}}}} \right).$$
(44)

Note that the sphere B1 has no intersection, neither with star 1 nor star 2. Thus, we can evaluate the gravitational force acting on star 1 without knowledge about the internal structure of star 1.

The surface integral approach was used by Einstein, Infeld, and Hoffmann in general relativity [73]. They used the vacuum Einstein equations only, and their method can be applied to any object including a black hole. We will take the surface integral approach in this article. But in our formalism, we shall treat only regular objects like a neutron star.

Now let us introduce the body zone to provide the surfaces of the surface integral approach. The scalings of R and m motivate us to define the body zone of star A(A = 1, 2) as \({B_A} \equiv \{{x^i}\vert \vert \vec x - {\vec z_A}(\tau)\vert \; < \epsilon {R_A}\}\) and the body zone coordinates of the star A as \(\alpha _A^{\underset{-}{i}} \equiv {\epsilon ^{- 2}}({x^i} - z_A^i(\tau))\). \(z_A^i(\tau)\) is a representative point of the star A, e.g. the center of the mass of the star A. RA, called the body zone radius, is an arbitrary length scale (much smaller than the orbital separation and not identical to the radius of the star) and constant (i.e. dRA/ = 0). With the body zone coordinates, the star does not shrink when ϵ → 0, while the boundary of the body zone goes to infinity (see Figure 1).

Figure 1
figure 1

Body zone coordinates and near zone coordinates. In the near zone coordinates (τ, xi), both the body zone and the star shrink as ϵ (thin dotted arrow) and ϵ2 (thick dotted arrow) respectively. Both the thin and thick arrows point inside. In the body zone coordinates \((\tau, \alpha _A^{\underset{-}{i}})\), the star does not shrink while the body zone boundary goes to infinity as ϵ−1 (thin dotted line pointing outside).

Then it is appropriate to define the star’s characteristic quantities such as the mass, the spin, and so on with the body zone coordinates. On the other hand the body zone serves us with a surface ∂BA, through which gravitational energy momentum flux flows, and in turn it amounts to the gravitational force acting on star A (see Figure 2). Since the body zone boundary ∂BA is far away from the surface of star A, we can evaluate the gravitational energy momentum flux over ∂BA with the post-Newtonian gravitational field. In fact we shall express our equations of motion in terms of integrals over ∂BA and be able to evaluate them explicitly.

Figure 2
figure 2

Gravitational energy momentum flux through the body zone boundary. The meshed two circles represent stars 1 and 2. Each star is surrounded by the body zone represented here by a striped area. The arrows around star 1 represent the gravitational energy momentum flux flowing through the body zone boundary.

3.3 Scalings on the initial hypersurface

Following [79, 82, 138], we use the initial value formulation to solve the Einstein equations. As initial data for the matter variables and the gravitational field, we take a set of nearly stationary solutions of the exact Einstein equations representing two widely separated fluid balls, each of which rotates uniformly. We assume that these solutions are parametrized by ϵ and that the matter and field variables have the following scalings.

The density scales as ϵ−4 (in the (t, xi) coordinates), implied by the scalings of m and R. The scaling of the density suggests that the natural dynamical time (free fall time) η inside the star may be η = ϵ−2t. Then if we cannot assume almost stationary condition on the stars, it is difficult to use the post-Newtonian approximation [80, 81]. In practice, however, our formalism is still applicable to pulsating stars if the effect of pulsation is not important in the orbital motion.

The velocity of stellar rotation is assumed to be \(\mathcal O(\epsilon)\). In other words, we assume that the star rotates slowly and is pressure supportedFootnote 6. By this assumption, the spin-orbit coupling force appears at 2 PN order rather than the usual 1.5 PN order. The slowly spinning motion assumption is not crucial: In fact, it is straightforward to incorporate a rapidly spinning compact body into our formalism.

From these initial data we have the following scalings of the star A’s stress-energy tensor components \(T_A^{\mu \nu}\) in the body zone coordinates: \(T_A^{\tau \tau} = \mathcal O({\epsilon ^{- 2}}),\;T_A^{\tau {\underset{-}{i}}} = \mathcal O({\epsilon ^{- 4}}),T_A^{i{\underset{-}{j}}} = \mathcal O({\epsilon ^{- 8}})\). Here the underlined indices mean that for any tensor \({A^i},\;{A^{\underset{-}{i}}} = {\epsilon ^{- 2}}{A^i}\). In [94], we have transformed \(T_N^{\mu \nu}\), the components of the stress-energy tensor of the matter in the near zone coordinates, to \(T_A^{\mu \nu}\) using the transformation from the near zone coordinates to the (generalized) Fermi normal coordinates at 1 PN order [13]. It is difficult, however, to construct the (generalized) Fermi normal coordinates at an high post-Newtonian order. Therefore we shall not use it. We simply assume that for \(T_N^{\mu \nu}\) (or rather \(\Lambda _N^{\mu \nu}\), the source term of the relaxed Einstein equations; see Equation (63)),

$$T_N^{\tau \tau} = {\mathcal O}({\epsilon ^{- 2}}),$$
(45)
$$T_N^{\tau \underline i} = {\mathcal O}({\epsilon ^{- 4}}),$$
(46)
$$T_N^{\underline i \underline j} = {\mathcal O}({\epsilon ^{- 8}}),$$
(47)

as their leading scalings.

As for the field variables on the initial hypersurface, we simply assume that the field is of 2.5 PN order except for the field determined by the constraint equations. Note that the radiation reaction effect to the stars first appears at the 2.5 PN order. Futamase showed that even if one takes the field of order 1 PN, initial value of the field does not affect the subsequent motion of the system up to 2.5 PN order [80]. Thus, we expect that the initial value of the field does not affect the orbital motion of the system up to 3 PN order, though a detailed calculation has not been done yet.

It is worth noticing that the initial value formulation has some advantages. First, by using the initial value formulation one can avoid the famous runaway solution problem in a radiation reaction problem. Second, one can construct an initial condition on some spacelike hypersurface rather than at past null infinity. Putting an initial condition for the field in the past null infinity requires a prior knowledge about the spacetime, which is obtained through the time evolution of the field from the initial condition. The initial value formulation can give in a sense a realistic initial condition. In our universe there may be no past null infinity because of the big bang.

An interesting initial condition is the statistical initial condition [138]. Here the binary system is in the background gravitational radiation bath for which we know only its statistical properties. For example, the phase of the radiation is assumed to be random and irrelevant to the motion of the binary. The origins of the radiation are cosmological, or related to the evolution of the system before the initial hypersurface. Then we can evaluate the expected time evolution of the binary system by letting the system evolve from a set of possible initial conditions and taking a statistical ensemble average over the initial conditions.

3.4 Newtonian equations of motion for extended bodies

Before ending this section, we present some equations for Newtonian extended bodies (stars). These equations will give a useful guideline when we develop our formalism.

The basic equations are the equation of continuity, the Euler equation, and the Poisson equation, respectively:

$${{\partial \rho} \over {\partial \tau}} + {{\partial \rho {v^i}} \over {\partial {x^i}}} = 0,$$
(48)
$$\rho {{D{v^i}} \over {D\tau}} = \rho {{\partial \Phi} \over {\partial {x_i}}},$$
(49)
$$\Delta \Phi = 4\pi \rho.$$
(50)

We define the mass, the dipole moment, the quadrupole moment, and the momentum of the star A as

$${m_A} \equiv \int\nolimits_{{B_A}} {{d^3}} x\,\rho,$$
(51)
$$D_A^i \equiv \int\nolimits_{{B_A}} {{d^3}} x\,\rho ({x^i} - z_A^i),$$
(52)
$$I_A^{ij} \equiv \int\nolimits_{{B_A}} {{d^3}} x\,\rho ({x^i} - z_A^i)({x^j} - z_A^j),$$
(53)
$$P_A^i \equiv \int\nolimits_{{B_A}} {{d^3}} x\,\rho {v^i}.$$
(54)

Here \(z_A^i\) is a representative point of the star A. The time derivative of the mass vanishes. Setting the time derivative of the dipole moment to zero gives the velocity momentum relation and a definition of the center of mass,

$${{d{D^i}} \over {d\tau}} = P_A^i - {m_A}v_A^i = 0,$$
(55)

where \(v_A^i = dz_A^i/d\tau\). Using the velocity momentum relation, we calculate the time derivative of the momentum,

$${m_1}{{{d^2}z_1^i} \over {d{\tau ^2}}} = F_1^i,$$
(56)

where \(F_1^i\) is defined by Equation (43). The Newtonian potential can be expressed by the mass and multipole moment as

$$\Phi = \sum\limits_{A = 1,2} {\left[ {{{{m_A}} \over {\vert \vec x - {{\vec z}_A}\vert}} + {1 \over 2}I_A^{ij}{{{\partial ^2}} \over {\partial {x^i}\,\partial {x^j}}}\left({{1 \over {\vert \vec x - {{\vec z}_A}\vert}}} \right)} \right]}.$$
(57)

Substituting Φ into Equation (56), we obtain the equations of motion,

$${m_1}{{{d^2}z_1^i} \over {d{\tau ^2}}} = {\partial \over {\partial z_1^i}}{{{m_1}{m_2}} \over {\vert {{\vec z}_1} - {{\vec z}_2}\vert}} + {1 \over 2}(I_1^{ij} + I_2^{ij}){{{\partial ^3}} \over {\partial z_1^i\partial z_2^j\partial z_2^k}}\left({{1 \over {\vert {{\vec z}_1} - {{\vec z}_2}\vert}}} \right).$$
(58)

Here we ignored the mass multipole moments of the stars that are of higher order than the quadrupole moments.

Actually, it is straightforward to formally include all the Newtonian mass multipole moments in the surface integral approach,

$${m_1}{{{d^2}z_1^i} \over {d{\tau ^2}}} = \sum\limits_{p = 0} {\sum\limits_{q = 0} {{{{{(- 1)}^{p + 1}}(2p + 2q + 1)!!} \over {p!q!}}}} {{I_1^{\langle {M_p}\rangle}I_2^{\langle {N_q}\rangle}} \over {r_{12}^{p + q + 2}}}{N^{\langle i{M_p}{N_q}\rangle}},$$
(59)

where \({\vec r_{12}} = {\vec z_1} - {\vec z_2}\), \({N^i} = r_{12}^i/{r_{12}}\), Mp = m1m2mp is a corrective index, 〈…〉 denotes the symmetric-tracefree operation on the indices between the brackets, and \(I_A^{{M_p}}\) are the Newtonian mass multipole moments of order 2p.

4 Formulation

Following the basic idea explained in the previous Section 3, we develop now our formalism for the derivation of the post-Newtonian equations of motion suitable for relativistic compact binaries. See also [94, 95].

4.1 Field equations

As discussed in the previous Section 3, we express our equations of motion in terms of surface integrals over the body zone boundary where it is assumed that the metric slightly deviates from the flat metric ημν = diag(−ϵ2, 1, 1, 1) (in the near zone coordinates (τ, xi)). Thus we define a deviation field hμν as

$${h^{\mu \nu}} \equiv {\eta ^{\mu \nu}} - \sqrt {- g} \,{g^{\mu \nu}},$$
(60)

where g is the determinant of the metric. Our hμν differs from the corresponding field in [30] in a sign. Indices are raised or lowered by the flat (auxiliary) metric ημν unless otherwise stated.

Now we impose the harmonic coordinate condition on the metric,

$${h^{\mu \nu}}_{,\nu} = 0,$$
(61)

where the comma denotes the partial derivative. In the harmonic gauge, we can recast the Einstein equations into the relaxed form,

$$\square{h^{\mu \nu}} = - 16\pi {\Lambda ^{\mu \nu}},$$
(62)

where □ = ημνμν is the flat d’Alembertian and

$${\Lambda ^{\mu \nu}} \equiv {\Theta ^{\mu \nu}} + {\chi ^{\mu \nu \alpha \beta}}_{,\alpha \beta},$$
(63)
$${\Theta ^{\mu \nu}} \equiv (- g)({T^{\mu \nu}} + t_{{\rm{LL}}}^{\mu \nu}),$$
(64)
$${\chi ^{\mu \nu \alpha \beta}} \equiv {1 \over {16\pi}}({h^{\alpha \nu}}{h^{\beta \mu}} - {h^{\alpha \beta}}{h^{\mu \nu}}).$$
(65)

Here, Tμν and \(t_{{\rm{LL}}}^{\mu \nu}\) denote the stress-energy tensor of the stars and the Landau-Lifshitz pseudotensor [115]. The explicit form of \(t_{{\rm{LL}}}^{\mu \nu}\) in the harmonic gauge is

$$\begin{array}{*{20}c} {(- 16\pi g)t_{{\rm{LL}}}^{\mu \nu} = {g_{\alpha \beta}}{g^{\gamma \delta}}{h^{\mu \alpha}}_{,\gamma}{h^{\nu \beta}}_{,\delta} + {1 \over 2}{g^{\mu \nu}}{g_{\alpha \beta}}{h^{\alpha \gamma}}_{,\delta}{h^{\beta \delta}}_{,\gamma} - 2{g_{\alpha \beta}}{g^{\gamma (\mu}}{h^{\nu)\alpha}}_{,\delta}{h^{\delta \beta}}_{,\gamma}\quad \quad \quad \quad} \\ {+ {1 \over 2}\left({{g^{\mu \alpha}}{g^{\nu \beta}} - {1 \over 2}{g^{\mu \nu}}{g^{\alpha \beta}}} \right)\left({{g_{\gamma \delta}}{g_{\epsilon \zeta}} - {1 \over 2}{g_{\gamma \epsilon}}{g_{\delta \zeta}}} \right){h^{\gamma \epsilon}}_{,\alpha}{h^{\delta \zeta}}_{,\beta}.} \\ \end{array}$$
(66)

χμναβ originates from our use of the flat d’Alembertian instead of the curved space d’Alembertian. In consistency with the harmonic condition, the conservation law is expressed as

$${\Lambda ^{\mu \nu}}_{,\nu} = 0.$$
(67)

Note that the divergence of χμναβ itself vanishes identically due to the symmetry of its indices.

Now we rewrite the relaxed Einstein equations into an integral form,

$${h^{\mu \nu}}(\tau, {x^i}) = 4\int\nolimits_{C(\tau, {x^k})} {{d^3}} y\,{{{\Lambda ^{\mu \nu}}(\tau - \epsilon \vert \vec x - \vec y\vert, {y^k};\epsilon)} \over {\vert \vec x - \vec y\vert}} + h_{\rm{H}}^{\mu \nu}(\tau, {x^i}),$$
(68)

where C(τ, xk) means the past light cone emanating from the event (τ, xk). C(τ, xk) is truncated at the τ = 0 hypersurface. \(h_{\rm{H}}^{\mu \nu}\) is a homogeneous solution of the homogeneous wave equation in the flat spacetime, \(\square h_{\rm{H}}^{\mu \nu} = 0\). \(h_{\rm{H}}^{\mu \nu}\) evolves from a given initial data on the τ = 0 initial hypersurface. The explicit form of \(h_{\rm{H}}^{\mu \nu}\) is available via the Poisson formula (see e.g. [144])

$$h_{\rm{H}}^{\mu \nu}(\tau, {x^i}) = {\tau \over {4\pi}}{\oint\nolimits_{\partial C(\tau, {x^i})} {{h^{\mu \nu}}} _{,\tau}}(\tau = 0,{y^i})\,d{\Omega _y} + {1 \over {4\pi}}{\partial \over {\partial \tau}}\left[ {\tau \oint\nolimits_{\partial C(\tau, {x^i})} {{h^{\mu \nu}}} (\tau = 0,{y^i})\,d{\Omega _y}} \right].$$
(69)

We solve the Einstein equations as follows. First we split the integral region into two zones: the near zone and the far zone.

The near zone is the region containing the gravitational wave source where the wave character of the gravitational radiation is not manifest. In other words, in the near zone the retardation effect on the field is negligible. The near zone covers the whole source system. The size of the near zone is about a little larger than one wave length of the gravitational wave emitted by the source. In this paper we take the near zone as a sphere centered at some fixed point and enclosing the binary system. The radius of this sphere is set to be \(\mathcal R/\epsilon\), where \({\mathcal R}\) is arbitrary but larger than the size of the binary and the wave length of the gravitational radiation. The scaling of the near zone radius is derived from the ϵ dependence of the wavelength of the gravitational radiation emitted due to the orbital motion of the binary. Note that roughly speaking the frequency of such a wave is about twice the Keplerian frequency of the binary. The center of the near zone sphere would be determined, if necessary, for example, to be the center of mass of the near zone. The outside of the near zone is the far zone where the retardation effect of the field is crucial.

For the near zone field point PN, we write the field as

$${h^{\mu \nu}}({{\rm{P}}_{\rm{N}}} \in N) = h_{{P_N}(N)}^{\mu \nu} + h_{{P_N}(F)}^{\mu \nu} + h_{\rm{H}}^{\mu \nu},$$
(70)
$$h_{{P_N}(N)}^{\mu \nu} \equiv 4\int\nolimits_{N = \{y:\vert y\vert \leq R/\epsilon \}} {{d^3}} y\,{{{\Lambda ^{\mu \nu}}(\tau - \epsilon \vert \vec x - \vec y\vert, {y^k};\epsilon)} \over {\vert \vec x - \vec y\vert}},$$
(71)
$$h_{{P_N}(F)}^{\mu \nu} \equiv 4\int\nolimits_{F = \{y:\vert y\vert > R/\epsilon \}} {{d^3}} y\,{{{\Lambda ^{\mu \nu}}(\tau - \epsilon \vert \vec x - \vec y\vert, {y^k};\epsilon)} \over {\vert \vec x - \vec y\vert}}.$$
(72)

\(h_{{P_N}(N)}^{\mu \nu}\) is the near zone integral contribution to the near zone field, and \(h_{{P_N}(F)}^{\mu \nu}\) is the far zone integral contribution to the near zone field.

The far zone contribution can be evaluated with the DIRE method developed in [129]. An explicit calculation shows that apparently there are the far zone contributions to the near zone field at 3 PN order. However, these 3 PN contributions are merely a gauge. Pati and Will showed that the far zone contribution does not affect the equations of motion up to 3 PN order inclusively and that the far zone contribution first appears at 4 PN order. This result is consistent with the earlier result of Blanchet and Damour [20] who used the multipolar-post-Minkowskian formalism. We follow the DIRE method and checke that the far zone contribution does not affect the equations of motion up to 3 PN order in Appendix A. Henceforth we shall focus our attention on the near zone contribution \(h_{{P_N}(N)}^{\mu \nu}\) and do not write down the far zone contribution in the following calculation of the field.

As for the homogeneous solution, we shall ignore it for simplicity. If we take random initial data for the field [138] supposed to be of 1 PN order [79], they are irrelevant to the dynamics of the binary system up to the radiation reaction order [79]. As we have assumed in the previous Section 3.3 that the magnitude of the free data of the gravitational field on the initial hypersurface is 2.5 PN order, we expect that the homogeneous solution does not affect the equations of motion up to 3 PN order. We leave a full implementation of the initial value formulation on the field as future work.

It is worth noticing that when we let τ become sufficiently large, then the condition \(h_{\rm{H}}^{\mu \nu}(\tau, {x^i}) = 0\) corresponds to a no-incoming radiation condition at (Minkowskian) past null infinity (see e.g. [77]). Equation (69) can be written down as Kirchhoff’s formula,

$$h_{\rm{H}}^{\mu \nu}(\tau, {x^i}) = \oint\nolimits_{\partial C(\tau, {x^i})} {{{d{\Omega _y}} \over {4\pi}}} {\left. {\left[ {{\partial \over {\partial \rho}}(\rho {h^{\mu \nu}}(\tau{\prime},{y^i})) + {\partial \over {\partial \tau{\prime}}}(\rho {h^{\mu \nu}}(\tau{\prime},{y^i}))} \right]} \right\vert _{\tau{\prime} = 0,\,\rho = \vert \vec x - \vec y\vert = \tau}}.$$
(73)

Then the no-incoming radiation condition at (Minkowskian) past null-infinity,

$$\underset {\underset {r \rightarrow \infty} {\tau = r,}} {\lim} \left[ {{\partial \over {\partial r}}(r{h^{\mu \nu}}) + {\partial \over {\partial \tau}}(r{h^{\mu \nu}})} \right] = 0,$$
(74)

is a sufficient condition to \(h_{\rm{H}}^{\mu \nu}(\tau, {x^i}) = 0\) when τ goes to infinity.

Now we shall devote ourselves to the evaluation of the near zone contribution to the near zone field,

$${h^{\mu \nu}}(\tau, {x^i}) = 4\int\nolimits_{N(\tau, {x^k})} {{d^3}} y\,{{{\Lambda ^{\mu \nu}}(\tau - \epsilon \vert \vec x - \vec y\vert, {y^k};\epsilon)} \over {\vert \vec x - \vec y\vert}}.$$
(75)

Henceforth we shall omit the subscript PN(N) of the field hμν for notational simplicity.

4.2 Near zone contribution

We shall evaluate the near zone contribution as follows. First, we make a retarded expansion of Equation (75) and change the integral region to a τ = constant spatial hypersurface,

$${h^{\mu \nu}} = 4\sum\limits_{n = 0} {{{{{(- \epsilon)}^n}} \over {n!}}} {\left({{\partial \over {\partial \tau}}} \right)^n}\int\nolimits_N {{d^3}} y\,\vert \vec x - \vec y{\vert ^{n - 1}}\Lambda _N^{\mu \nu}(\tau, {y^k};\epsilon).$$
(76)

Note that the above integral depends on the arbitrary length \({\mathcal R}\) in general. The cancellation between the \({\mathcal R}\) dependent terms in the far zone contribution and those in the near zone contribution was shown by [129] through all the post-Newtonian order. In the following, we shall omit the terms which have negative powers of \({\mathcal R\;({\mathcal R^{- k}};k > 0)}\). In other words, we simply let \({\mathcal R} \rightarrow 0\) whenever it gives a convergent result. On the other hand, we shall retain terms having positive powers of \({\mathcal R}\;({{\mathcal R}^k};k > 0)\) and logarithmic terms \((\ln \;\mathcal R)\) to confirm that the final result, in the end of calculation, is independent of \({\mathcal R}\) (and for logarithmic terms to keep the arguments of logarithm non-dimensional).

Second we split the integral into two parts: a contribution from the body zone BA, and from elsewhere, N/B. Schematically we evaluate the following two types of integrals (we omit indices of the field),

$$\begin{array}{*{20}{c}} {h = {h_B} + {h_{N/B}},\quad \quad \quad \quad \quad \quad \quad \quad \quad } \\ {{h_B} = {^6}\sum\limits_{A = 1,2} {\int_{{B_A}} {{d^3}} } {\alpha _A}{\mkern 1mu} \frac{{f(\tau ,{{\vec z}_A} + {^2}{{\vec \alpha }_A})}}{{|{{\vec r}_A} - {^2}{{\vec \alpha }_A}{|^{1 - n}}}},{\mkern 1mu} {\mkern 1mu} {\mkern 1mu} \quad } \\ {{h_{N/B}} = \int_{N/B} {{d^3}y\frac{{f(\tau ,\vec y)}}{{|\vec x - \vec y{|^{1 - n}}}},} \quad \quad \quad \quad \quad {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} \quad } \end{array}$$
(77)

where \({\vec r_A} \equiv \vec x - {\vec z_A}\). We shall deal with these two contributions successively.

4.2.1 Body zone contribution

As for the body zone contribution, we make a multipole expansion using the scaling of the integrand, i.e. Λμν in the body zone. For example, the n = 0 part in Equation (76), \(h_{B\;n = 0}^{\mu \nu}\) gives

$$h_{B\,n = 0}^{\tau \tau} = 4{\epsilon ^4}\sum\limits_{A = 1,2} {\left({{{P_A^\tau} \over {{r_A}}} + {\epsilon ^2}{{D_A^kr_A^k} \over {r_A^3}} + {\epsilon ^4}{{3I_A^{\langle kl\rangle}r_A^kr_A^l} \over {2r_A^5}} + {\epsilon ^6}{{5I_A^{\langle klm\rangle}r_A^kr_A^lr_A^m} \over {2r_A^7}}} \right)} + {\mathcal O}({\epsilon ^{12}}),$$
(78)
$$h_{B\,n = 0}^{\tau i} = 4{\epsilon ^4}\sum\limits_{A = 1,2} {\left({{{P_A^i} \over {{r_A}}} + {\epsilon ^2}{{J_A^{ki}r_A^k} \over {r_A^3}} + {\epsilon ^4}{{3J_A^{kli}r_A^kr_A^l} \over {2r_A^5}}} \right)} + {\mathcal O}({\epsilon ^{10}}),$$
(79)
$$h_{B\,n = 0}^{ij} = 4{\epsilon ^2}\sum\limits_{A = 1,2} {\left({{{Z_A^{ij}} \over {{r_A}}} + {\epsilon ^2}{{Z_A^{kij}r_A^k} \over {r_A^3}} + {\epsilon ^4}{{3Z_A^{\langle kl\rangle ij}r_A^kr_A^l} \over {2r_A^5}} + {\epsilon ^6}{{5Z_A^{\langle klm\rangle ij}r_A^kr_A^lr_A^m} \over {2r_A^7}}} \right)} + {\mathcal O}({\epsilon ^{10}}).$$
(80)

Here the operator 〈…〉 denotes a symmetric and tracefree (STF) operation on the indices in the brackets. See [129, 150, 165] for some useful formulas of STF tensors. Also we define \({r_A} \equiv \vert {\vec r_A}\vert\). To derive the 3 PN equations of motion, hττ up to \({\mathcal O}({\epsilon ^{10}})\) and hτi as well as hij up to \({\mathcal O}({\epsilon ^8})\) are required.

In the above equations we define the multipole moments as

$$I_A^{{K_l}} \equiv {\epsilon ^2}\int\nolimits_{{B_A}} {{d^3}} {\alpha _A}\,\Lambda _N^{\tau \tau}\,\alpha _A^{{{\underline K}_l}},$$
(81)
$$J_A^{{K_l}i} \equiv {\epsilon ^4}\int\nolimits_{{B_A}} {{d^3}} {\alpha _A}\,\Lambda _N^{\tau \underline i}\,\alpha _A^{{{\underline K}_l}},$$
(82)
$$Z_A^{{K_l}ij} \equiv {\epsilon ^8}\int\nolimits_{{B_A}} {{d^3}} {\alpha _A}\,\Lambda _N^{\underline i \underline j}\,\alpha _A^{{{\underline K}_l}},$$
(83)

where we introduced a collective multi-index \({I_l} \equiv {i_1}{i_2} \ldots {i_l}\) and \(\alpha _A^{{{\underset{-}{I}}_l}} \equiv \alpha _A^{{{\underset{-}{i}}_1}}\alpha _A^{{{\underset{-}{i}}_2}} \ldots \alpha _A^{{{\underset{-}{i}}_l}}\). Then \(P_A^\tau \equiv I_A^{{K_0}}\), \(D_A^{{k_1}} \equiv I_A^{{K_1}}\), and \(P_A^{{k_1}} \equiv J_A^{{K_1}}\). We simply call \(P_A^\mu\) the four-momentum of the star A, \(P_A^i\) the momentum, and \(P_A^\tau\) the energyFootnote 7. Also we call \(D_A^k\) the dipole moment of the star A, and \(I_A^{kl}\) the quadrupole moment of the star A.

Then we transform these moments into more convenient forms. By the conservation law (67), we have

$$\Lambda _N^{\tau i} = {(\Lambda _N^{\tau \tau}y_A^i)_{,\tau}} + {(\Lambda _N^{\tau j}y_A^i)_{,j}} + v_A^i\Lambda _N^{\tau \tau},$$
(84)
$$\Lambda _N^{ij} = {(\Lambda _N^{\tau (i}y_A^{j)})_{,\tau}} + {(\Lambda _N^{k(i}y_A^{j)})_{,k}} + v_A^{(i}\Lambda _N^{j)\tau},$$
(85)

where \(v_A^i \equiv \dot z_A^i\) an overdot denotes a time derivative with respect to τ, and \({\vec y_A} \equiv \vec y - {\vec z_A}\). Using these equations and noticing that the body zone remains unchanged (in the near zone coordinates), i.e. A = 0, we have

$$P_A^i = P_A^\tau v_A^i + Q_A^i + {\epsilon ^2}{{dD_A^i} \over {d\tau}},$$
(86)
$$J_A^{ij} = {1 \over 2}\left({M_A^{ij} + {\epsilon ^2}{{dI_A^{ij}} \over {d\tau}}} \right) + v_A^{(i}D_A^{j)} + {1 \over 2}{\epsilon ^{- 2}}Q_A^{ij},$$
(87)
$$Z_A^{ij} = {\epsilon ^2}P_A^\tau v_A^iv_A^j + {1 \over 2}{\epsilon ^6}{{{d^2}I_A^{ij}} \over {d{\tau ^2}}} + 2{\epsilon ^4}v_A^{(i}{{dD_A^{j)}} \over {d\tau}} + {\epsilon ^4}{{dv_A^{(i}} \over {d\tau}}D_A^{j)} + {\epsilon ^2}Q_A^{(i}v_A^{j)} + {\epsilon ^2}R_A^{(ij)} + {1 \over 2}{\epsilon ^2}{{dQ_A^{ij}} \over {d\tau}},$$
(88)
$$Z_A^{kij} = {3 \over 2}A_A^{kij} - A_A^{(ij)k},$$
(89)

where

$$M_A^{ij} \equiv 2{\epsilon ^4}\int\nolimits_{{B_A}} {{d^3}} {\alpha _A}\,\alpha _A^{\left[ {\underline i} \right.}\Lambda _N^{{}^{\left. {\underline j} \right]}\tau},$$
(90)
$$Q_A^{{K_l}i} \equiv {\epsilon ^{- 4}}\oint\nolimits_{\partial {B_A}} d {S_k}\,\left({\Lambda _N^{\tau k} - v_A^k\Lambda _N^{\tau \tau}} \right)y_A^{{K_l}}y_A^i,$$
(91)
$$R_A^{{K_l}\,i\,j} \equiv {\epsilon ^{- 4}}\oint\nolimits_{\partial {B_A}} d {S_k}\,\left({\Lambda _N^{k\,j} - \upsilon _A^k\Lambda _N^{\tau \,j}} \right)y_A^{{K_l}}y_A^i,$$
(92)

and

$$A_A^{kij} \equiv {\epsilon ^2}J_A^{k(i}v_A^{j)} + {\epsilon^2}v_A^kJ_A^{(ij)} + R_A^{k(ij)} + {\epsilon ^4}{{dJ_A^{k(ij)}} \over {d\tau}}.$$
(93)

The operators [ ] and () attached to the indices denote anti-symmetrization and symmetrization. \(M_A^{ij}\) is the spin of the star A and Equation (86) is the momentum-velocity relation. Thus our momentum-velocity relation is a direct analog of the Newtonian momentum-velocity relation (see Section 3.4). In general, we have

$$J_A^{{K_l}i} = J_A^{({K_l}i)} + {{2l} \over {l + 1}}J_A^{({K_{l - 1}}[{k_l})i]},$$
(94)
$$Z_A^{{K_l}ij} = {1 \over 2}\left[ {Z_A^{({K_l}i)j} + {{2l} \over {l + 1}}Z_A^{({K_{l - 1}}[{k_l})i]j} + Z_A^{({K_l}j)i} + {{2l} \over {l + 1}}Z_A^{({K_{l - 1}}[{k_l})j]i}} \right],$$
(95)

and

$$J_A^{({K_l}i)} = {1 \over {l + 1}}{\epsilon ^2}{{dI_A^{{K_l}i}} \over {d\tau}} + v_A^{(i}I_A^{{K_l})} + {1 \over {l + 1}}{\epsilon ^{- 2l}}Q_A^{{K_l}i},$$
(96)
$$Z_A^{({K_l}i)j} + Z_A^{({K_l}j)i} = {\epsilon ^2}v_A^{(i}J_A^{{K_l})j} + {\epsilon ^2}v_A^{(j}J_A^{{K_l})i} + {2 \over {l + 1}}{\epsilon ^4}{{dJ_A^{{K_l}(ij)}} \over {d\tau}} + {\epsilon ^{- 2l + 2}}R_A^{{K_l}(ij)},$$
(97)

where l is the number of indices in the multi-index Kl.

Now, from the above equations, especially Equation (88), we find that the body zone contributions, \(h_{B\;n = 0}^{\mu \nu}\), are of order \({\mathcal O}({\epsilon ^4})\). Note that if we can not or do not assume the (nearly) stationarity of the initial data for the stars, then, instead of Equation (88) we have

$$Z_A^{ij} = {\epsilon ^2}P_A^\tau v_A^iv_A^j + {1 \over 2}{{{d^2}I_A^{ij}} \over {d{\eta ^2}}} + \ldots,$$
(98)

where we used the dynamical time η (see Section 3.3). In this case the lowest order metric differs from the Newtonian form. From our (almost) stationarity assumption the remaining motion inside a star, apart from the spinning motion, is caused only by the tidal effect by the companion star and from Equation (88); it appears at 3 PN order [81].

To obtain the lowest order \(h_{B\;n = 0}^{\mu \nu}\), we have to evaluate the surface integrals \(Q_A^{{K_l}i}\), \(R_A^{{K_l}ij}\). Generally, in \(h_{B\;n = 0}^{\mu \nu}\) the moments \(J_A^{{K_l}i}\) and \(Z_A^{{K_l}ij}\) appear formally at the order \({\epsilon ^{2l + 4}}\) and \({\epsilon ^{2l + 2}}\). Thus \(Q_A^{{K_l}i}\) and \(R_A^{{K_l}ij}\) appear as

$$h_{B\,n = 0}^{\tau i} \sim \cdots + {\epsilon ^4}{{r_A^{{K_l}}Q_A^{\langle {K_l}\rangle i}} \over {r_A^{2l + 1}}} + \ldots,$$
(99)
$$h_{B\,n = 0}^{ij} \sim \cdots + {\epsilon ^4}{{r_A^{{K_l}}R_A^{\langle {K_l}\rangle (ij)}} \over {r_A^{2l + 1}}} + \ldots,$$
(100)

where we omitted irrelevant terms and numerical coefficients. Thus one may expect that \(Q_A^{{K_l}i}\) and \(R_A^{{K_l}ij}\) appear at the order ϵ4 for any l and we have to calculate an infinite number of moments. In fact, this is not the case and it was shown in [91] that only l = 0, 1 terms of \(R_A^{{K_l}ij}\) contribute to the 3 PN equations of motion. The important thing here is that \({\epsilon ^4}Q_A^{{K_l}i}\) and \({\epsilon ^4}R_A^{{K_l}ij}\) are at most \({\mathcal O}({\epsilon ^4})\) in \(h_{B\;n = 0}^{\mu \nu}\).

Finally, since the order of \(h_{Bn}^{\mu \nu}(n \geq 1)\) is higher than that of \(h_{B\;n = 0}^{\mu \nu}\), we conclude that \(h_B^{\mu \nu} = {\mathcal O}({\epsilon ^4})\).

4.2.2 N/B contribution

As far as the N/B contribution is concerned, since the integrand \(\Lambda _N^{\mu \nu} = - gt_{{\rm{LL}}}^{\mu \nu} + {\chi ^{\mu \nu \alpha \beta}}_{,\alpha \beta}\) is at least quadratic in the small deviation field hμν, we make the post-Newtonian expansion in the integrand. Then, basically, with the help of (super-)potentials \(g(\vec x)\) which satisfy \(\Delta g(\vec x) = f(\vec x)\), Δ denoting the Laplacian, we have for each integral (see e.g. the n = 0 term in Equation (77))

$$\int\nolimits_{N/B} {{d^3}} y{{f(\vec y)} \over {\vert \vec x - \vec y\vert}} = - 4\pi g(\vec x) + \oint\nolimits_{\partial (N/B)} d {S_k}\,\left[ {{1 \over {\vert \vec x - \vec y\vert}}{{\partial g(\vec y)} \over {\partial {y^k}}} - g(\vec y){\partial \over {\partial {y^k}}}\left({{1 \over {\vert \vec x - \vec y\vert}}} \right)} \right].$$
(101)

Equation (101) can be derived without using Dirac delta distributions (see Appendix B of [91]). For the n ≥ 1 terms in Equation (77), we use potentials many times to convert all the volume integrals into surface integrals and “\(- 4\pi g(\vec x)\)” termsFootnote 8.

In fact finding the super-potentials is one of the most formidable tasks especially when we proceed to high post-Newtonian orders. Fortunately, up to 2.5 PN order, all the required super-potentials are available. At 3 PN order, there appear many integrands for which we could not find the required super-potentials. To obtain the 3 PN equations of motion, we devise an alternative method similar to the method employed by Blanchet and Faye [28]. The details of the method will be explained later.

Now the lowest order integrands can be evaluated with the body zone contribution \(h_B^{\mu \nu}\), and since \(h_B^{\mu \nu}\) is \({\mathcal O}({\epsilon ^4})\), we find

$$\Lambda _N^{\tau \tau} = {\mathcal O}({\epsilon ^6}),$$
(102)
$$\Lambda _N^{\tau i} = {\mathcal O}({\epsilon ^6}),$$
(103)
$$\Lambda _N^{ij} = (- g)t_{{\rm{LL}}}^{ij} + {\mathcal O}({\epsilon ^8}) = {\epsilon ^4}{1 \over {64\pi}}\left({{\delta ^i}_k{\delta ^j}_l - {1 \over 2}{\delta ^{ij}}{\delta _{kl}}} \right){\,_4}h_B^{\tau \tau, k}\,{}_4h_B^{\tau \tau, l} + {\mathcal O}({\epsilon ^5}),$$
(104)

where we expanded \(h_B^{\mu \nu}\) in an ϵ series,

$$h_B^{\mu \nu} = \sum\limits_{n = 0} {{\epsilon ^{4 + n}}_n} h_B^{\mu \nu}.$$
(105)

Similarly, in the following we expand hμν in an ϵ series. From these equations we find that the deviation field in N/B, hμν, is \({\mathcal O}({\epsilon ^4})\). (It should be noted that in the body zone hμν is assumed to be of order unity and within our method we can not calculate hμν there explicitly. To obtain in the body zone, we have to know the internal structure of the star.)

4.3 Lorentz contraction and multipole moments

In Equations (81, 82, 83), we have defined multipole moments of a star. The definition of those multipole moments are operational ones and are not necessarily equal to “intrinsic” multipole moments of the star. This is clear, for example, if we remember that a moving ball has spurious multipole moments due to Lorentz contraction.

We have adopted the generalized Fermi normal coordinates (GFC) [13] as the star’s reference coordinates to address this problem. A question specific to our formalism is whether the differences between the multipole moments defined in Equations (81, 82, 83) and the multipole moments in GFC give purely monopole terms. If so, we of course have to take into account such terms in the field to compute the equations of motion for two intrinsically spherical star binaries. This problem is addressed in the Appendix C of [91] and the differences are mainly attributed to the shape of the body zone. The body zone BA which is spherical in the near zone coordinates (NZC) is not spherical in the GFC mainly because of a kinematic effect (Lorentz contraction). To derive the 3 PN equations of motion, it turns out that it is sufficient to compute the difference in the STF quadrupole moment up to 1 PN order. The result is

$$\begin{array}{*{20}c} {\delta I_A^{\langle ij\rangle} \equiv I_{A,{\rm{NZC}}}^{\langle ij\rangle} - I_{A,{\rm{GFC}}}^{\langle ij\rangle}\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {= {\epsilon ^{- 6}}{1 \over 2}v_A^kv_A^l\oint\nolimits_{\partial {B_A}} d {S_k}y_A^l\,y_A^{\langle i}y_A^{j\rangle}{\Lambda ^{\tau \tau}} + {\mathcal O}({\epsilon ^3})} \\ {= - {\epsilon ^2}{{4m_A^3} \over 5}v_A^{\langle i}v_A^{j\rangle} + {\mathcal O}({\epsilon ^3}),\quad \quad \quad \quad \quad \quad} \\ \end{array}$$
(106)

where \(I_{A,{\rm{NZC}}}^{ij} \equiv I_A^{ij} \cdot I_{A,{\rm{GFC}}}^{ij}\) are the quadrupole moments defined in the generalized Fermi normal coordinates. Note that the difference \(\delta I_A^{\langle ij\rangle}\) is expressed in a surface integral form.

As is obvious from Equation (106), this difference appears even if the companion star does not exist. We note that we could derive the 3 PN metric for an isolated star A moving at a constant velocity using our method explained in this section by simply letting the mass of the companion star be zero. Actually, the \(\delta I_A^{\langle ij\rangle}\) above is a necessary term which makes the so-obtained 3 PN metric equal to the Schwarzschild metric boosted at the constant velocity \({\vec v_A}\) in harmonic coordinates.

The reason why the influence of a body’s Lorentz contraction appears starting from 3 PN order is as follows. The body’s Lorentz contraction appears as an apparent deformation of the body, namely, apparent quadrupole moment as the leading order in the frame where the body is moving. The body’s (real) quadrupole affects the field from 2 PN order through (see Equation (78))

$$h_{{\rm{quad}}}^{\tau \tau} = {\epsilon ^8}\sum\limits_A {{{3I_A^{\langle ij\rangle}r_A^ir_A^j} \over {2r_A^5}}}.$$
(107)

The radius of a compact body A is of order of its mass mA, and the Lorentz contraction deforms the body so that the radius of the body changes by the amount \({\epsilon ^2}{m_A}v_A^2 + \mathcal O({\epsilon ^4})\). So the apparent quadrupole moments due to Lorentz contraction is of order of

$$I_{A\,{\rm{apparent}}}^{\langle ij\rangle}\sim{\epsilon ^2}m_A^3v_A^2.$$
(108)

This field is the 3 PN field and thus the Lorentz contraction of a body affects the equations of motion starting from 3 PN order.

4.4 General form of the equations of motion

From the definition of the four-momentum,

$$P_A^\mu (\tau) = {\epsilon ^2}\int\nolimits_{{B_A}} {{d^3}} {\alpha _A}\Lambda _N^{\tau \mu},$$
(109)

and the conservation law (67), we have the evolution equations for the four-momentum:

$${{dP_A^\mu} \over {d\tau}} = - {\epsilon ^{- 4}}\oint\nolimits_{\partial {B_A}} d {S_k}\,\Lambda _N^{k\mu} + {\epsilon ^{- 4}}v_A^k\oint\nolimits_{\partial {B_A}} d {S_k}\,\Lambda _N^{\tau \mu}.$$
(110)

Here we used the fact that the size and the shape of the body zone are defined to be fixed (in the near zone coordinates), while the center of the body zone moves at the velocity of the star’s representative point.

Substituting the momentum-velocity relation (86) into the spatial components of Equation (110), we obtain the general form of equations of motion for star A,

$$\begin{array}{*{20}c} {P_A^\tau {{dv_A^i} \over {d\tau}} = - {\epsilon ^{- 4}}\oint\nolimits_{\partial {B_A}} d {S_k}\,\Lambda _N^{ki} + {\epsilon ^{- 4}}v_A^k\oint\nolimits_{\partial {B_A}} d {S_k}\,\Lambda _N^{\tau i}\quad \quad \quad} \\ {+ {\epsilon ^{- 4}}v_A^i\left({\oint\nolimits_{\partial {B_A}} d {S_k}\,\Lambda _N^{k\tau} - v_A^k\oint\nolimits_{\partial {B_A}} d {S_k}\,\Lambda _N^{\tau \tau}} \right)} \\ {- {{dQ_A^i} \over {d\tau}} - {\epsilon ^2}{{{d^2}D_A^i} \over {d{\tau ^2}}}.\quad \quad \quad \quad \quad \quad \quad \quad} \\ \end{array}$$
(111)

All the right hand side terms in Equation (111) except the dipole moment are expressed as surface integrals. We can specify the value of \(D_A^i\) freely to determine the representative point \(z_A^i(\tau)\) of star A. Up to 2.5 PN order we take \(D_A^i = 0\) and simply call \(z_A^i\) the center of mass of star A. Note that in order to obtain the spin-orbit coupling force in the same form as in previous works [46, 108, 152], we have to make another choice for \(z_A^i\) (see [94] and Appendix B.1). At 3 PN order, yet another choice of the value of the dipole moment \(D_A^i\) shall be examined (see Sections 7 and 8.2).

In Equation (111), \(P_A^\tau\) rather than the mass of star A appears. Hence we have to derive a relation between the mass and \(P_A^\tau\). We shall derive that relation by solving the temporal component of the evolution equations (110) functionally.

Then, since all the equations are expressed with surface integrals except \(D_A^i\) to be specified, we can derive the equations of motion for a strongly self-gravitating star using the post-Newtonian approximation.

4.5 On the arbitrary constant RA

Since we have introduced the body zone by hand, the arbitrary body zone size RA seems to appear in the metric, the multipole moments of the stars, and the equations of motion. More specifically RA appears in them because of (i) the splitting of the deviation field into two parts (i.e. B and N/B contributions), the definition of the moments, and (ii) the surface integrals that we evaluate to derive the equations of motion.

4.5.1 RA dependence of the field

B and N/B contributions to the field depend on the body zone boundary ϵRA. But hμν itself does not depend on ϵRA. Thus it is natural to expect that there are renormalized multipole moments which are independent of RA since we use nonsingular matter sources. This renormalization would absorb the ϵRA dependence occuring in the computation of the N/B field (see Section 4.8 for an example of such a renormalization). One possible practical obstacle for this expectation might be the ln(ϵRA) dependence of multipole moments. Although at 3 PN order there appear such logarithmic terms, it is found that we could remove them by rechoosing the value of the dipole moment \(D_A^i\) of the star.

Though we use the same symbol for the moments henceforth as before for notational simplicity, it should be understood that they are the renormalized ones. For instance, we use the symbol “\(P_A^\mu\)” for the renormalized \(P_A^\mu\).

4.5.2 RA dependence of the equations of motion

Since we compute integrals over the body zone boundary, in general the resulting equations of motion seem to depend on the size of the body zone boundary, ϵRA. Actually this is not the case.

In the derivation of Equation (111), if we did not use the conservation law (67) until the final step, we have

$$\begin{array}{*{20}c} {P_A^\tau {{dv_A^i} \over {d\tau}} + {\epsilon ^{- 4}}\oint\nolimits_{\partial {B_A}} d {S_k}\,\Lambda _N^{ki} - {\epsilon ^{- 4}}v_A^k\oint\nolimits_{\partial {B_A}} d {S_k}\,\Lambda _N^{\tau i} - {\epsilon ^{- 4}}v_A^i\left({\oint\nolimits_{\partial {B_A}} d {S_k}\,\Lambda _N^{k\tau} - v_A^k\oint\nolimits_{\partial {B_A}} d {S_k}\,\Lambda _N^{\tau \tau}} \right) +} \\ {{{dQ_A^i} \over {d\tau}} + {\epsilon ^2}{{{d^2}D_A^i} \over {d{\tau ^2}}}} \\ {= {\epsilon ^{- 4}}\int\nolimits_{{B_A}} {{d^3}} y\,\Lambda _{N,\,\nu}^{i\nu} - {\epsilon ^{- 4}}v_A^i\int\nolimits_{{B_A}} {{d^3}} y\,\Lambda _{N,\,\nu}^{\tau \nu} + {\epsilon ^{- 4}}{d \over {d\tau}}\left({\int\nolimits_{{B_A}} {{d^3}} y\,\Lambda _{N,\,\nu}^{\tau \nu}y_A^i} \right).} \\ \end{array}$$
(112)

Now the conservation law is satisfied for whatever value we take for RA, then the right hand side of the above equation is zero for any RA. Hence the equations of motion (111) do not depend on RA (a similar argument can be found in [73]).

Along the same line, the momentum-velocity relation (86) does not depend on RA.

In Section 4.8 we shall explicitly show the irrelevance of the field and the equations of motion to RA by checking the cancellation among the RA dependent terms up to 0.5 PN order.

4.6 Newtonian equations of motion

We first derive the Newtonian equations of motion and the 1 PN correction to the acceleration. The derivation of the 1 PN equations of motion includes some essences of our formalism and shows how it works properly. Thus we shall give a detailed explanation about the derivation, though the calculation is elementary and the resulting equations of motion are well-known.

Now let us derive the Newtonian mass-energy relation first. From Equations (102, 103) and the time component of Equation (110), we have

$${{dP_A^\tau} \over {d\tau}} = {\mathcal O}({\epsilon ^2}).$$
(113)

then we define the mass of star A as the integrating constant of this equation,

$${m_A} \equiv \underset{\epsilon \rightarrow 0} {\lim} P_A^\tau.$$
(114)

here mA is the ADM mass that star A had if it were isolated. (We took ϵ zero limit in Equation (114) to ensure that the mass defined above does not include the effect of the companion star and the orbital motion of the star itself. By this limit we ensure that this mass is the integrating constant of Equation (113).) By definition mA is constant. Then we obtain the lowest order hττ:

$$_4{h^{\tau \tau}}{= _4}h_B^{\tau \tau} = 4 \sum\limits_{A = 1,2} {{{{m_A}} \over {{r_A}}}}.$$
(115)

Second, from Equation (91) with Equations (102) and (103) we obtain \(Q_A^i = 0\) at the lowest order. Thus we have the Newtonian momentum-velocity relation \(P_A^i = {m_A}v_A^i + \mathcal O({\epsilon ^2})\) from Equation (86) (we set \(D_A^i = 0\)).

Substituting Equation (104) with the lowest order hττ into the general form of equations of motion (111) and using the Newtonian momentum-velocity relation, we have for star 1:

$$\begin{array}{*{20}c} {{m_1}{{dv_1^i} \over {d\tau}} = - \oint\nolimits_{\partial {B_1}} d {S_{k(4)}}\,[(- g)t_{{\rm{LL}}}^{ik}]\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {= - {1 \over {4\pi}}\left({\delta _l^i\delta _m^k - {1 \over 2}{\delta ^{ik}}{\delta _{lm}}} \right)\sum\limits_{A = 1,2} {\sum\limits_{B = 1,2} {\oint\nolimits_{\partial {B_1}} d}} {S_k}\,{{{m_A}{m_B}y_A^ly_B^m} \over {y_A^3y_B^3}},\quad \quad \quad \quad \quad \quad} \\ {= - {1 \over {4\pi}}\left({\delta _l^i\delta _m^k - {1 \over 2}{\delta ^{ik}}{\delta _{lm}}} \right)\oint d {\Omega _{{{\bf{n}}_{\bf{1}}}}}\,n_1^k({{m_1^2n_1^ln_1^m} \over {{{(\epsilon {R_1})}^2}}} + {{2{m_1}{m_2}n_A^{(l}(r_{12}^{m)} + \epsilon {R_1}n_1^{m)})} \over {\vert {{\vec r}_{12}} + \epsilon {R_1}{{\vec n}_1}{\vert ^3}}}} \\ {\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad + {{{{(\epsilon {R_1})}^2}m_2^2(r_{12}^l + \epsilon {R_1}n_1^l)(r_{12}^m + \epsilon {R_1}n_1^m)} \over {\vert {{\vec r}_{12}} + \epsilon {R_1}{{\vec n}_1}{\vert ^3}}}),} \\ \end{array}$$
(116)

where \(y_A^i\) is the integral variable and \({\vec y_A} = \vec y - {\vec z_A}\). We defined \({{\vec n}_1}\) as the spatial unit vectorFootnote 9 emanating from \({{\vec z}_1}\), \({{\vec r}_{12}} \equiv {{\vec z}_1} - {{\vec z}_2}\), and \(n_{12}^i \equiv r_{12}^i/{r_{12}}\). We used \({\vec r_1} = \epsilon {R_1}{\vec n_1}\) and \({\vec r_2} = {\vec r_{12}} + \epsilon {R_1}{\vec n_1}\). For details see Figure 3).

Figure 3
figure 3

The vectors used in the surface integral over the boundary of the body zone 1.

In the above equation, by virtue of the angular integral the first term (which is singular when the ϵ zero limit is taken) vanishes. The third term vanishes by letting ϵ go to zero. Only the second term survives and gives the Newtonian equations of motion as expected. This completes the Newtonian order calculations.

4.7 First post-Newtonian equations of motion

Next, at 1 PN order, we need 6hττ and 4hμν. The n = 1 term in the retardation expansion series of hτ τ, Equation (76), gives no contribution at the 1 PN order by the constancy of the mass mA, i.e. 5hττ = 0.

Now we obtain 4hτ i from the Newtonian momentum-velocity relation:

$$_4{h^{\tau i}} = 4\sum\limits_{A = 1,2} {{{{m_A}v_A^i} \over {{r_A}}}}.$$
(117)

We evaluate the surface integrals in the evolution equation for \(P_A^\tau\) at 1 PN order. The result for star 1 is

$${{dP_1^\tau} \over {d\tau}} = {\epsilon ^2}{m_1}{d \over {d\tau}}\left({{1 \over 2}v_1^2 + {{3{m_2}} \over {{r_{12}}}}} \right) + {\mathcal O}({\epsilon ^3}),$$
(118)

where we used the Newtonian equations of motion. From this equation we have the mass-energy relation at 1 PN order,

$$P_1^\tau = {m_1}\left[ {1 + {\epsilon ^2}\left({{1 \over 2}v_1^2 + {{3{m_2}} \over {{r_{12}}}}} \right)} \right] + {\mathcal O}({\epsilon ^3}).$$
(119)

Then we have to calculate the 1 PN order \(Q_A^i\). The result is \(Q_A^i = {\epsilon ^2}m_A^2v_A^i/(6\epsilon {R_A})\). As \(Q_A^i\) depends on RA, we ignore it (see Section 4.5). As a result we obtain the momentum-velocity relation at 1 PN order, \(P_A^i = P_A^\tau v_A^i + \mathcal O({\epsilon ^3})\) from Equation (86).

Now as for 4hij, we first calculate the surface integrals \(Q_A^{ij}\), \(R_A^{ij}\), and \(R_A^{kij}\) from Equations (91) and (92). We then find that they depend on RA, hence we ignore them and obtain

$$h_B^{ij} = 4{\epsilon ^4}\sum\limits_{A = 1,2} {{{{m_A}v_A^iv_A^j} \over {{r_A}}}} + {\mathcal O}({\epsilon ^5}).$$
(120)

To derive 6hττ and 4hij, we have to evaluate non-compact support integrals for \(_6h_{N/B}^{\tau \tau}\) and \(_4h_{N/B}^{ij}\), and the n = 2 term in Equation (76) for hττ:

$${h^{\tau \tau}} = 4{\epsilon ^4}\sum\limits_{A = 1,2} {{{P_A^\tau} \over {{r_A}}}} + 4{\epsilon ^6}\int\nolimits_{N/B} {{{{d^3}y} \over {\vert \vec x - \vec y{\vert ^6}}}} [(- g)t_{{\rm{LL}}}^{\tau \tau}] + 2{\epsilon ^6}{\partial \over {\partial {\tau ^2}}}\sum\limits_{A = 1,2} {P_A^\tau} {r_A},$$
(121)
$${h^{ij}} = 4{\epsilon^4}\sum\limits_{A = 1,2} {{{{m_A}\upsilon _A^i\upsilon _A^j} \over {{r_A}}}} + 4{\epsilon^4}\int\nolimits_{N/B} {{{{d^3}y} \over {\vert \vec x - \vec y\vert}}\,{}_4} [(- g)t_{{\rm{LL}}}^{ij}],$$
(122)

with

$$_6[(- 16\pi g)t_{{\rm{LL}}}^{\tau \tau}] = - {7 \over 8}{\,_4}{h^{\tau \tau, k}}_4{h^{\tau \tau}}_{,k},$$
(123)
$$_4[(- 16\pi g)t_{{\rm{LL}}}^{ij}] = {1 \over 4}\left({{\delta ^i}_k{\delta ^i}_k - {1 \over 2}{\delta ^{ij}}{\delta _{kl}}} \right){\,_4}{h^{\tau \tau, k}}_4{h^{\tau \tau, l}}.$$
(124)

The evaluation of the twice retardation expansion term (the last term in Equation (121)) is straightforward. For the non-compact support integrals, it is sufficient to consider the following integral:

$$\int\nolimits_{N/B} {{{{d^3}y} \over {16\pi |\vec x - \vec y|}}\,{}_4} {h^{\tau \tau ,i}}_4{h^{\tau \tau ,j}} = {1 \over \pi}\sum\limits_{A = 1,2} {P_A^{\tau 2}} \int\nolimits_{N/B} {{{{d^3}y} \over {|\vec x - \vec y|}}} {{r_A^ir_A^j} \over {r_A^6}} + {2 \over \pi}P_1^\tau P_2^\tau \int\nolimits_{N/B} {{{{d^3}y} \over {|\vec x - \vec y|}}} {{r_1^{(i}r_2^{j)}} \over {r_{\,\,1}^3r_{\,\,2}^3}}.$$
(125)

To evaluate the above Poisson integrals, we use Equation (101), thus we need to find superpotentials for the integrands. For this purpose, it is convenient to transform the tensorial integrands into scalar integrands with differentiation operators,

$${{r_A^ir_A^j} \over {r_A^6}} = {1 \over 8}\left({{\delta ^i}_k{\delta ^j}_l + {\delta ^{ij}}{\delta _{kl}}} \right){{{\partial ^2}} \over {\partial {y^k}{y^l}}}\left({{1 \over {r_A^2}}} \right),$$
(126)
$${{r_1^ir_2^j} \over {r_1^3r_2^3}} = {\partial \over {\partial z_1^i}}{\partial \over {\partial z_2^i}}\left({{1 \over {{r_1}{r_2}}}} \right).$$
(127)

then it is relatively easy to find the super-potentials for these scalars. The results are \(\Delta \ln {r_1} = 1/r_1^2\) and Δ ln S = 1/(r1r2) where S = r1 + r2 + r12 [77]. Equation (101) for each integrand then becomes

$$\int\nolimits_{N/B} {{{{d^3}y} \over {\vert \vec x - \vec y\vert}}} {{r_A^ir_A^j} \over {r_A^6}} = - 4\pi \left({{{{\delta ^{ij}}} \over {4r_A^2}} - {{r_A^ir_A^j} \over {4r_A^4}}} \right) + {{4\pi {\delta ^{ij}}} \over {3\epsilon {R_1}{r_1}}} + {\mathcal O}(\epsilon {R_A}),$$
(128)
$$\begin{array}{*{20}c} {\int\nolimits_{N/B} {{{{d^3}y} \over {\vert \vec x - \vec y\vert}}} {{r_1^ir_2^j} \over {r_1^3r_2^3}} = - 4\pi \left({- {{{\delta ^{ij}}} \over {{r_{12}}S}} - {{r_1^ir_{12}^j} \over {{r_1}{r_{12}}{S^2}}} + {{r_{12}^ir_{12}^j} \over {r_{12}^2{S^2}}} + {{r_{12}^ir_{12}^j} \over {r_{12}^3S}} - {{r_1^ir_{12}^j} \over {{r_1}{r_2}{S^2}}} + {{r_{12}^ir_2^j} \over {{r_{12}}{r_2}{S^2}}}} \right)} \\ {+ {\mathcal O}(\epsilon {R_A}){.}\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ \end{array}$$
(129)

The second last term of Equation (128) and the terms abbreviated as \({\mathcal O}(\epsilon {R_A})\) in the above two equations arise from the surface integrals in Equation (101). Since they depend on ϵRA, we ignore it (see Section 4.8). Substituting Equations (128) and (129) back into Equation (125), we can compute the non-compact support integrals.

Using the above results and the Newtonian equations of motion for the twice retardation expansion term, we finally obtain

$$\begin{array}{*{20}c} {{h^{\tau \tau}} = 4{\epsilon ^4}\sum\limits_{A = 1,2} {{{P_A^\tau} \over {{r_A}}}\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad}} \\ {+ \,{\epsilon ^6} \left[ {- 2\sum\limits_{A = 1,2} {{{{m_A}} \over {{r_A}}}\{{{({{\vec n}_A}\cdot{{\vec v}_A})}^2} - v_A^2\} + 2{{{m_1}{m_2}} \over {r_{12}^2}}{{\vec n}_{12}}\cdot({{\vec n}_1} - {{\vec n}_2})}} \right.} \\ {\quad \quad \quad + \left. {7\sum\limits_{A = 1,2} {{{m_A^2} \over {r_A^2}}} + 14{{{m_1}{m_2}} \over {{r_1}{r_2}}} - 14{{{m_1}{m_2}} \over {{r_{12}}}}\sum\limits_{A = 1,2} {{1 \over {{r_A}}}}} \right] + {\mathcal O}({\epsilon ^7}),} \\ \end{array}$$
(130)
$$\begin{array}{*{20}c} {{h^{ij}} = 4{\epsilon ^4}\sum\limits_{A = 1,2} {{{{m_A}v_A^iv_A^j} \over {{r_A}}}} \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {+ \,{\epsilon ^4}\left[ {\sum\limits_{A = 1,2} {{{m_A^2} \over {r_A^2}}} n_A^in_A^j - {{8{m_1}{m_2}} \over {{r_{12}}S}}n_{12}^in_{12}^j\quad \quad \quad \quad \quad \quad} \right.\quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {- \left. {8\left({{\delta ^i}_k{\delta ^j}_l - {1 \over 2}{\delta ^{ij}}{\delta _{kl}}} \right){{{m_1}{m_2}} \over {{S^2}}}{{({{\vec n}_{12}} - {{\vec n}_1})}^{(k}}{{({{\vec n}_{12}} + {{\vec n}_2})}^{l)}}} \right] + {\mathcal O}({\epsilon ^5}),} \\ \end{array}$$
(131)

where \({\vec n_A} \equiv {\vec r_A}/{r_A}\).

Evaluating the surface integrals in Equation (111) as in the Newtonian case, we obtain the 1 PN equations of motion,

$$\begin{array}{*{20}c} {{m_1}{{dv_1^i} \over {d\tau}} = - {{{m_1}{m_2}} \over {r_{12}^2}}n_{12}^i\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {+ \,{\epsilon ^2}{{{m_1}{m_2}} \over {r_{12}^2}}\left[ {n_{12}^i\left({- v_1^2 - 2v_2^2 + {3 \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} + 4({{\vec v}_1}\cdot{{\vec v}_2}) + {{5{m_1}} \over {{r_{12}}}} + {{4{m_2}} \over {{r_{12}}}}} \right)} \right.} \\ {\left. {+ {V^i}\left({4({{\vec n}_{12}}\cdot{{\vec v}_1}) - 3({{\vec n}_{12}}\cdot{{\vec v}_2})} \right)} \right],\quad \quad \quad \quad \quad} \\ \end{array}$$
(132)

where we defined the relative velocity as \(\vec V \equiv {\vec v_1} - {\vec v_2}\) and we used the Newtonian equations of motion as well as Equation (119).

Finally let us give a summary of our procedure (see Figure 4). With the n PN order equations of motion and hμν in hand, we first derive the n + 1 PN evolution equation for \(P_A^\tau\). Then we solve it functionally and obtain the mass-energy relation at n + 1 PN order. Next we calculate \(Q_A^i\) at n + 1 PN order and derive the momentum-velocity relation at n + 1 PN order. Then we calculate \(Q_A^{{K_l}i}\) and \(R_A^{{K_l}ij}\). With the n + 1 PN mass-energy relation, the n + 1 PN momentum-velocity relation, \(Q_A^{{K_l}i}\) and \(R_A^{{K_l}ij}\) we next derive the n + 1 PN deviation field hμν. Finally we evaluate the surface integrals which appear in the right hand side of Equation (111) and obtain the n + 1 PN equations of motion. In the above calculations we use the n PN order equations of motion to reduce the order of the equations of motion whenever an acceleration appears in the right hand of the resulting equations of motion. For instance, when we meet \({\epsilon ^2}dv_1^i/d\tau\) in the right hand side of the equations motion and we have to evaluate this up to ϵ2, then using the Newtonian equations of motion, we replace it by \(- {\epsilon ^2}{m_2}r_{12}^i/r_{12}^3\). Basically we shall derive the 3 PN equations of motion with the procedure as described above.

Figure 4
figure 4

Flowchart of the post-Newtonian iteration.

4.8 Body zone boundary dependent terms

As explained in Section 4.5 we discard the body zone boundary dependent terms in the field, since we expect that they cancel out between the body zone contribution and the N/B contribution. Before moving on to the higher order calculations, however, it is instructive to see that such a cancellation really occurs in the field and consequently the equations of motion up to 0.5 PN order.

First, we show the independence of the body zone radius RA in the 0.5 PN field. Returning back to the derivation of the 1 PN hττ with care for the RA dependence (see Equation (128)), we get

$$\begin{array}{*{20}c} {{h^{\tau \tau}} = 4{\epsilon ^4}\sum\limits_{A = 1,2} {{{P_A^\tau} \over {{r_A}}}} + 4{\epsilon ^6}\int\nolimits_{C/B} {{d^3}} y\,{{{}_6\Lambda _N^{\tau \tau}(\tau, {y^k})} \over {\vert \vec x - \vec y\vert}} + 2{\epsilon ^6}{{{\partial ^2}} \over {\partial {\tau ^2}}}\sum\limits_{A = 1,2} {P_A^\tau} {r_A} + {\mathcal O}({\epsilon ^7})} \\ {= 4{\epsilon ^4}\sum\limits_{A = 1,2} {{1 \over {{r_A}}}} \left({P_A^\tau - {\epsilon ^2}{{7m_A^2} \over {2\epsilon {R_A}}}} \right) + {\mathcal O}({\epsilon ^6}).\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ \end{array}$$
(133)

Now we split \(P_A^\tau\) as \(P_A^\tau = \bar P_A^\tau + \tilde P_A^\tau\) where \(\bar P_A^\tau\) does not depend on RA while \(\tilde P_A^\tau\) does. In order to evaluate the RA dependent part of \(P_A^\tau\), \(\tilde P_A^\tau\) we use the fact that the integral of \(\Lambda _N^{\tau \tau}\) over the near zone does not depend on the size of the body zone. (Notice that \(P_A^\tau\) is defined as the volume integral of \(\Lambda _N^{\tau \tau}\) over BA.) In fact, from the relevant expression of the pseudotensor, we obtain at the lowest order

$$\begin{array}{*{20}c} {{\epsilon ^2}\int\nolimits_{N/B} {{d^3}} \alpha {\epsilon ^6}_6\Lambda _N^{\tau \tau} = {\epsilon ^2}\int\nolimits_{N/B} {{d^3}} {y_6}\Lambda _N^{\tau \tau}\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {= - {\epsilon ^2}{7 \over 2}\sum\limits_{A = 1,2} {{{m_A^2} \over {\epsilon {R_A}}}\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \quad \quad} \\ {\quad + ({\rm{terms}}\,{\rm{independent}}\,{\rm{of}}\,{R_A},\,{\rm{or}}\,{\rm{terms}}\,{\rm{having}}\,{\rm{positive}}\,{\rm{power}}({\rm{s}})\,{\rm{of}}\,{R_A}){.}} \\ \end{array}$$
(134)

Hence we find \(\tilde P_A^\tau = {\epsilon ^2}7m_A^2/(2\epsilon {R_A}) + {\mathcal O}({\epsilon ^2})\). The above equation shows us that the 0.5 PN field is independent of RA and fully expressed by the RA independent energy \(\bar P_A^\tau\), or mass up to 0.5 PN order.

In the similar manner we split \(P_A^i\) as \(P_A^i = \bar P_A^i + \tilde P_A^i\) and obtain \(\tilde P_A^i = {\epsilon ^2}11m_A^2v_A^i/(3\epsilon {R_A}) + {\mathcal O}({\epsilon ^2})\) from the definition of \(P_A^i\). Thus from the fact that \(Q_A^i = {\epsilon ^2}m_A^2v_A^i/(6\epsilon {R_A}) + {\mathcal O}({\epsilon ^2})\), we find that the 0.5 PN momentum-velocity relation does not depend on \({R_A}:\;\bar P_A^i = \bar P_A^\tau v_A^i + {\mathcal O}({\epsilon ^2})\). Finally evaluating the surface integrals in the general form of equations of motion using the “renormalized” (barred) moments, we find that the equations of motion are independent of RA as was expected.

As remarked in Section 4.5 from now on we use the same symbol for the renormalized moments as for the “bare” moments henceforth as before for notational simplicity.

5 Third Post-Newtonian Gravitational Field

In the post-Newtonian approximation, we need to solve Poisson equations to find the metric. Up to 2.5 PN order, explicit forms of the metric have been obtained in [30]. However, it seems impossible to derive the 3 PN accurate gravitational field in harmonic coordinates in a closed form completely. The problem is that it seems difficult (if at all possible) to find a particular solution of the Poisson equations for non-compact sources. The works so far overcome this problem by not solving the Poisson equations but keeping the Poisson integral unevaluated. Then to derive the equations of motion, we basically interchange the order of operations; first we evaluate surface integrals with the Poisson integrals as integrands (in the surface integral approach [91]) or compute derivatives of the Poisson integrals (when one adopts the geodesic equation [27]), and then we evaluate the remaining volume integrals. We first explain the usual method and a method to derive a field just around the star, and then explain the method mentioned above.

5.1 Super-potential method

Up to 2.5 PN order, we have solved all the Poisson equations necessary to derive the 2.5 PN gravitational field. At 3 PN order, we have found a part of the solutions of the Poisson equations, which we call (super-)potentialsFootnote 10. For example,

$${{r_1^i\,r_1^j\,r_1^k\,r_2^l} \over {r_1^5\,r_2^3}} = \Delta \left[ {- {1 \over 3}{{{\partial ^4}} \over {\partial z_1^i\,\partial z_1^j\,\partial z_1^k\,\partial z_2^l}}{f^{(1, - 1)}} + {1 \over 3}({\delta ^{ij}}{\partial _{z_1^k}} + {\delta ^{ik}}{\partial _{z_1^j}} + {\delta ^{jk}}{\partial _{z_1^i}}){\partial _{z_2^l}}\ln S} \right],$$
(135)

where S = r1 + r2 + r12, and f(1,−1) which satisfies Δf(1,−1) = r1/r2 is given in [99] as

$${f^{(1, - 1)}} = {1 \over {18}}(- r_1^2 - 3{r_1}{r_{12}} - r_{12}^2 + 3{r_1}{r_2} + 3{r_{12}}{r_2} + r_2^2) + {1 \over 6}(- r_1^2 + r_{12}^2 + r_2^2)\ln S.$$
(136)

It is possible to add any homogeneous solution to super-potentials. In our formalism, the only place where we use super-potentials is Equation (101). In the case above, we could add, say, 1/r1 to f(1,−1). (Note that to evaluate the surface integrals in the general form of equations of motion (111) we need super-potentials in the spatial region N/B which do not include any singularity due to the point particle limit.) It is easy to see that contribution from a possible additional homogeneous solution cancels out between the “\(- 4\pi g(\vec x)\)” term and the surface integral in Equation (101).

Useful super-potentials for derivation of the 3 PN field are given in [30, 27, 91, 99].

5.2 Super-potential-in-series method

As all what we need to do is to evaluate the surface integrals in the general form of equations of motion (111), we need an expression for the gravitational field only around the star. In fact, we have developed such a method in [91] for the source term of the following form

$${\partial _{z_A^i}}{\partial _{z_{A^{\prime}}^j}}g(\vec x) \equiv {\partial _{z_A^i}}{\partial _{z_{A^{\prime}}^j}}\left({{{{{(\ln {r_1})}^p}{{(\ln {r_2})}^q}} \over {r_1^ar_2^b}}} \right),$$
(137)

where a and b are integers and p = 0, 1, q = 0, 1. Note that A, A′ = 1, 2. Then, we take spatial derivatives out of the Poisson integral,

$$\int\nolimits_{N/B} {{{{d^3}y} \over {\vert \vec x - \vec y\vert}}} {\partial _{z_A^i}}{\partial _{z_{A^{\prime}}^j}}g(\vec y) = {\partial _{z_A^i}}{\partial _{z_{A^{\prime}}^j}}\int\nolimits_{N/B} {d^3}y\,{{g(\vec y)} \over {\vert \vec x - \vec y\vert}} + {\partial _{z_A^i}}\oint\nolimits_{\partial {B_{A^{\prime}}}} d{S_j}\,{{g(\vec y)} \over {\vert \vec x - \vec y\vert}} + \oint\nolimits_{\partial {B_A}} d{S_i}\,{{{\partial _{z_{A^{\prime}}^j}}g(\vec y)} \over {\vert \vec x - \vec y\vert}}.$$
(138)

Note that the integration region is N/B and therefore \(g(\vec x)\) is nonsingular in N/B. For this kind of source term, we have given a method in [91] to find a field \(F_{[A,c]}^{(m,n)}\) in the neighborhood of star A in the following sense:

$$\Delta F_{[A,c]}^{(p,m,q,n)} - {(\ln {r_1})^p}r_1^m{(\ln {r_2})^q}r_2^n = {\mathcal O}\left({{{r_A^{c + 1}} \over {r_{12}^{c + 1}}}} \right)\qquad {\rm{as}}\;{r_A} \rightarrow 0.$$
(139)

We have checked at 3 PN order that the resulting field from this method is equal to the field obtained from the usual (super-potential) method whenever the super-potentials are available. Unfortunately, however, this method is not perfect and we need another method to derive the equations of motion which we explain now.

5.3 Direct-integration method

In the surface integral approach, we need to evaluate a surface integral of the Landau-Lifshitz pseudotensor which has a form hμν,αhλσ,β. From an order counting, it would be clear that we need to evaluate a surface integral of the form “Newtonian potential” × “3 PN potential” to find the 3 PN equations of motion, where it seems impossible to find the “3 PN potential” in a closed form, as mentioned above. The surface integral mentioned here thus has the form

$$\oint\nolimits_{\partial {B_1}} d {S_k}{{r_A^l} \over {r_A^3}}\underset {\epsilon {R_A}} {{\rm{disc}}} \int\nolimits_{N/B} {{d^3}} y\,{{f(\vec y)} \over {\vert \vec x - \vec y\vert}}$$
(140)

for star 1. Here the operator \({\rm{dis}}{{\rm{c}}_{\epsilon {R_A}}}\) means to discard all the ϵRA dependent terms other than logarithms of ϵRA [91]. The method to evaluate this type of integral is to exchange the order of integrals, first calculate the surface integral, and then calculate the volume integral. One caveat is that we can not simply exchange the order of integrals, and we put an operation discϵRA in front of the Poisson integral as in Equation (140) above.

We first exchange the order of integrals,

$$\oint\nolimits_{\partial {B_1}} d {S_k}\,{{r_A^l} \over {r_A^3}}\underset {\epsilon {R_A}} {{\rm{disc}}} \int\nolimits_{N/B} {{d^3}} y\,{{f({{\vec y}_1})} \over {\vert {{\vec r}_1} - {{\vec y}_1}\vert}} = \underset {{r_1}\prime \to \epsilon {R_1}} {\lim} \underset {\epsilon {R_A}} {{\rm{disc}}} \int\nolimits_{N/B} {{d^3}} {y_1}\,f({\vec y_1})\,\oint\nolimits_{\partial {{B\prime}_1}} d {S_k}\,{1 \over {\vert {{\vec y}_1} - r_1\prime {{\vec n}_1}\vert}}{\partial _{z_A^l}}{1 \over {{r_A}}},$$
(141)

where we defined a sphere \({B_1^\prime}\) whose center is \({\vec z_1}\) and that has a radius \({r_1^\prime}\) which is a constant slightly larger than ϵRA for any (small) \(\epsilon \; (\epsilon {R_1} < r_1^{\prime} \ll {r_{12}})\).

The reason we introduced \({r_1^\prime}\) is as follows. Suppose that we treat an integrand for which the super-potential is available. By calculating the Poisson integral, we have a piece of field corresponding to the integrand. The piece generally depends on ϵRA, however we reasonably discard such RA-dependent terms (other than logarithmic dependence) as explained in Section 4.5. Using the so-obtained RA-independent field, we evaluate the surface integrals in the general form of the 3 PN equations of motion by discarding the ϵRA dependence emerging from the surface integrals, and obtain the equations of motion. Thus the “discarding-ϵRA” procedure must be employed each time when the field is derived and also when equations of motion are derived, not just once. Thus \({r_1^\prime}\) was introduced to distinguish the two kinds of ϵRA dependence and to discard the ϵRA dependence in the right order. We show here a simple example. Let us consider the following integral:

$$\oint\nolimits_{\partial {B_1}} d {S_k}\,{{r_1^k} \over {r_1^3}}\int\nolimits_{N/B} {{{{d^3}y} \over {\vert \vec x - \vec y\vert}}} {1 \over {y_1^2}}.$$
(142)

Using \(\Delta \ln {r_1} = 1/r_1^2\), we can integrate the Poisson integral and obtain the “field”,

$$\int\nolimits_{N/B} {{{{d^3}y} \over {\vert \vec x - \vec y\vert}}} {1 \over {y_1^2}} = - 4\pi \ln \left({{{{r_1}} \over {R/\epsilon}}} \right) + 4\pi - 4\pi {{\epsilon {R_1}} \over {{r_1}}} + {\mathcal O}\left({{{\left({\epsilon {R_A}} \right)}^2}} \right).$$
(143)

Since the “body zone contribution” must have an ϵR1 dependence hidden in the “moments” as \(4\pi \epsilon {R_1}/{r_1}[ + {\mathcal O}({(\epsilon {R_A})^2})]\) (see Section 4.5), the last term should be discarded before we evaluate the the surface integral in Equation (142) using this “field”. The surface integral gives the “equations of motion”,

$$16{\pi ^2}\left({\ln \left({{{R/\epsilon} \over {\epsilon {R_1}}}} \right) + 1} \right).$$
(144)

On the other hand, we can derive the “equations of motion” by first evaluating the surface integral over \(\partial B_1^{\prime}\),

$$\begin{array}{*{20}c} {\oint\nolimits_{\partial {B_1}} d {S_k}\,{{r_1^k} \over {r_1^3}}\int\nolimits_{N/B} {{{{d^3}y} \over {\vert \vec x - \vec y\vert}}} {1 \over {y_1^2}} = \int\nolimits_{N/B} {{{{d^3}y} \over {y_1^2}}} \oint\nolimits_{\partial B_1^\prime} d {S_k}{{r_1^k} \over {r_1^3}}{1 \over {\vert {{\vec r}_1} - {{\vec y}_1}\vert}}} \\ {= 16{\pi ^2}\left[ {\int\nolimits_{\epsilon {R_1}}^{r_1^\prime} {{{dy} \over {r_1^\prime}}} + \int\nolimits_{r_1^\prime}^{R/\epsilon} {{{dy} \over {{y_1}}}}} \right]} \\ {= 16{\pi ^2}\left({\ln \left({{{{\mathcal R}/\epsilon} \over {r_1^\prime}}} \right) + 1 - {{\epsilon {R_1}} \over {r_1^\prime}}} \right).} \\ \end{array}$$
(145)

Thus, if we take ∂B1 as the integral region instead of \(\partial B_1^{\prime}\) in the first equality in Equation (145), or if we take \(\lim \nolimits _{r_1^{\prime}\rightarrow \epsilon R_{1}}\) without employing \({\rm{dis}}{{\rm{c}}_{\epsilon {R_A}}}\) beforehand, we will obtain an incorrect result,

$$16{\pi ^2}\ln \left({{{{\mathcal R}/\epsilon} \over {\epsilon {R_1}}}} \right),$$
(146)

which disagrees with Equation (144).

With this caution in mind, it is straightforward (though tedious) to evaluate the surface integrals and then the volume integrals, and thus evaluate all the necessary integrals for our derivation of the 3 PN equations of motion.

When we have derived the 3 PN equations of motion, we have used the super-potential method whenever possible, and used a combination of the above three methods when necessary. In fact, for a computational check, we have used the direct-integration method to evaluate the contributions to the equations of motion from all of the 3 PN N/B nonretarded field, \(_8h_{N/B\;n = 0}^{\tau i}\) and \(_{10}h_{N/B\;n = 0}^{\tau \tau}{+ _8}h_{N/B\;n = {0^l}}^l\). As expected, we obtain the same result from two computations; the result from the direct-integration method agrees with that from the combination of the three methods: the direct-integration method, the super-potential method, and the super-potential-in-series method.

6 Third Post-Newtonian Mass-Energy Relation

It was found that the direct-integration part does not play any role in the evaluation of the evolution equation of \(P_A^\tau\) at 3 PN order. Thus we use the same method as in the evaluation of the 2.5 PN equations of motion. Evaluating the surface integrals in Equation (110), we obtain the evolution equation of \(P_A^\tau\) as

$$\begin{array}{*{20}c} {{{dP_{1\Theta}^\tau} \over {d\tau}} = - {\epsilon ^2}{{{m_1}{m_2}} \over {r_{12}^2}}\left[ {4({{\vec n}_{12}}\cdot{{\vec v}_1}) - 3({{\vec n}_{12}}\cdot{{\vec v}_2})} \right]\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {+ {\epsilon ^4}{{{m_1}{m_2}} \over {r_{12}^2}}\left[ {- {9 \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3} + {1 \over 2}v_1^2({{\vec n}_{12}}\cdot{{\vec v}_2}) + 6({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}\quad} \right.} \\ {\quad \quad \quad \quad \quad \quad \quad \quad - 2v_1^2({{\vec n}_{12}}\cdot{{\vec v}_1}) + 4({{\vec v}_1}\cdot{{\vec v}_2})({{\vec n}_{12}}\cdot\vec V) + 5v_2^2({{\vec n}_{12}}\cdot{{\vec v}_2}) - 4v_2^2({{\vec n}_{12}}\cdot{{\vec v}_1})\quad} \\ {\quad \quad \quad \quad \quad \quad \quad + \left. {{{{m_1}} \over {{r_{12}}}}(- 4({{\vec n}_{12}}\cdot{{\vec v}_2}) + 6({{\vec n}_{12}}\cdot{{\vec v}_1})) + {{{m_2}} \over {{r_{12}}}}(- 10({{\vec n}_{12}}\cdot{{\vec v}_1}) + 11({{\vec n}_{12}}\cdot{{\vec v}_2}))} \right]} \\ {\quad \quad \, + {\epsilon ^6}{{{m_1}{m_2}} \over {{r_{12}}}}\left[ {- \left({{3 \over 2}v_1^4 + 2v_1^2v_2^2 + 4v_2^4} \right)({{\vec n}_{12}}\cdot{{\vec v}_1}) + \left({{5 \over 8}v_1^4 + {3 \over 2}v_1^2v_2^2 + 7v_2^4} \right)({{\vec n}_{12}}\cdot{{\vec v}_2})} \right.} \\ {\quad \quad \quad \quad + \left({2v_1^2 + 4v_2^2} \right)({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec v}_1}\cdot{{\vec v}_2}) - \left({2v_1^2 + 8v_2^2} \right)({{\vec n}_{12}}\cdot{{\vec v}_2})({{\vec v}_1}\cdot{{\vec v}_2})} \\ {\quad \quad \quad \quad + \left({3v_1^2 + 12v_2^2} \right)({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} - \left({{3 \over 4}v_1^2 + 12v_2^2} \right){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3}} \\ {\quad \quad \quad \quad \quad \quad \quad + 2({{\vec n}_{12}}\cdot{{\vec v}_2}){{({{\vec v}_1}\cdot{{\vec v}_2})}^2} - 6({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}({{\vec v}_1}\cdot{{\vec v}_2}) + 6{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3}({{\vec v}_1}\cdot{{\vec v}_2})} \\ {- {{15} \over 2}({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^4} - {{45} \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^5}\quad \quad \quad \quad} \\ {+ {{{m_1}} \over {{r_{12}}}}\left\{{\left({- 42v_1^2 - {{117} \over 4}v_2^2} \right)({{\vec n}_{12}}\cdot{{\vec v}_1}) + 60{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^3}} \right.} \\ {\quad \quad \quad \quad \quad \quad + \left({{{137} \over 4}v_1^2 + {{37} \over 2}v_2^2} \right)({{\vec n}_{12}}\cdot{{\vec v}_2}) + {{297} \over 4}({{\vec n}_{12}} \cdot {{\vec v}_1})(\vec v \cdot {{\vec v}_2})} \\ {\quad \quad \quad \quad - {{219} \over 4}({{\vec n}_{12}} \cdot {{\vec v}_2})(\vec v \cdot {{\vec v}_2}) - 151{{({{\vec n}_{12}} \cdot {{\vec v}_1})}^2}({{\vec n}_{12}}\cdot{{\vec v}_2})} \\ {\left. {\quad + 109({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} - 23{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3}} \right\}} \\ {\quad \quad \quad \quad + {{{m_2}} \over {{r_{12}}}}\left\{{- \left({13v_1^2 + 18v_2^2} \right)({{\vec n}_{12}}\cdot{{\vec v}_1}) + \left({17v_1^2 + 25v_2^2} \right)({{\vec n}_{12}}\cdot{{\vec v}_2})} \right.} \\ {\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad + 26({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec v}_1}\cdot{{\vec v}_2}) - 28({{\vec n}_{12}}\cdot{{\vec v}_2})({{\vec v}_1}\cdot{{\vec v}_2}) + 2{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}({{\vec n}_{12}}\cdot{{\vec v}_2})} \\ {\quad + \left. {16({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} - 20{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3}} \right\}} \\ {\qquad \qquad \quad \,\quad \quad \quad + {{m_1^2} \over {r_{12}^2}}\left({{{33} \over 4}({{\vec n}_{12}}\cdot{{\vec v}_1}) - {{13} \over 2}({{\vec n}_{12}}\cdot{{\vec v}_2})} \right) - {{{m_1}{m_2}} \over {r_{12}^2}}\left({{{35} \over 4}({{\vec n}_{12}}\cdot{{\vec v}_1}) + {{17} \over 4}({{\vec n}_{12}}\cdot{{\vec v}_2})} \right)} \\ {\, + \left. {{{m_2^2} \over {r_{12}^2}}\left({- 12({{\vec n}_{12}}\cdot{{\vec v}_1}) + {{23} \over 2}({{\vec n}_{12}}\cdot{{\vec v}_2})} \right)} \right]\quad \quad \quad \quad} \\ {+ {\mathcal O}({\epsilon ^7}).\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ \end{array}$$
(147)

Remarkably, we can integrate Equation (147) functionally:

$$P_{1\Theta}^\tau = {m_1}\left({1 + {\epsilon ^2}_2{\Gamma _1} + {\epsilon ^4}_4{\Gamma _1} + {\epsilon ^6}_6{\Gamma _1}} \right) + {\mathcal O}({\epsilon ^7}),$$
(148)

with

$$_2{\Gamma _1} = {1 \over 2}v_1^2 + {{3{m_2}} \over {{r_{12}}}},$$
(149)
$$_4{\Gamma _1} = - {{3{m_2}} \over {2{r_{12}}}}{({\vec n_{12}}\cdot{\vec v_2})^2} + {{2{m_2}} \over {{r_{12}}}}v_2^2 + {{7{m_2}} \over {2{r_{12}}}}v_1^2 - {{4{m_2}} \over {{r_{12}}}}({\vec v_1}\cdot{\vec v_2}) + {3 \over 8}v_1^4 + {{7m_2^2} \over {2r_{12}^2}} - {{5{m_1}{m_2}} \over {2r_{12}^2}},$$
(150)
$$\begin{array}{*{20}c} {_6{\Gamma _1} = {{m_1^2{m_2}} \over {2r_{12}^3}} + {{21{m_1}m_2^2} \over {4r_{12}^3}} + {{5m_2^3} \over {2r_{12}^3}} + {5 \over {16}}v_1^6\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {\,\, + {{m_2^2} \over {r_{12}^2}}\left({{{45} \over 4}v_1^2 + {{19} \over 2}v_2^2 + {1 \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2} - 19({{\vec v}_1}\cdot{{\vec v}_2}) - ({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2}) - 3{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}} \right)\quad \quad} \\ {\quad \quad + {{{m_1}{m_2}} \over {r_{12}^2}}\left({{{43} \over 8}v_1^2 + {{53} \over 8}v_2^2 - {{69} \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2} - {{53} \over 4}({{\vec v}_1}\cdot{{\vec v}_2}) + {{85} \over 4}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2}) - {{69} \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}} \right)} \\ {+ {{{m_2}} \over {{r_{12}}}}\left({{{33} \over 8}v_1^4 + {3 \over 2}v_1^2v_2^2 + v_2^4 - 6v_1^2({{\vec v}_1}\cdot{{\vec v}_2}) - 4v_2^2({{\vec v}_1}\cdot{{\vec v}_2}) - {7 \over 4}v_1^2{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}\quad \quad \quad \quad \quad} \right.} \\ {\quad - \left. {{5 \over 2}v_2^2{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} + 2{{({{\vec v}_1}\cdot{{\vec v}_2})}^2} + 2{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}({{\vec v}_1}\cdot{{\vec v}_2}) + {9 \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^4}} \right).\quad \quad \quad} \\ \end{array}$$
(151)

Equation (148) together with Equations (149, 150, 151) gives the 3 PN order mass-energy relation.

6.1 Meaning of \(P_{A\Theta}^\tau\)

In this section we explain the meaning of \(P_{A\Theta}^\tau\). First of all, we expand in an ϵ series the four-velocity of star A normalized as \({g_{\mu \nu}}u_A^\mu u_A^\nu = - {\epsilon ^{- 2}}\), where \(u_A^i = u_A^\tau v_A^i\). The result is

$$\begin{array}{*{20}c} {u_A^\tau = 1 + {\varepsilon ^2}\left[ {{1 \over 2}v_A^2 + {1 \over 4}\,{}_4{h^{\tau \tau}}} \right]\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {+ {\varepsilon ^4}\left[ {{1 \over 4}\,{}_6{h^{\tau \tau}} + {1 \over 4}\,{}_4{h^k}_k - {3 \over {32}}{{{(_4}{h^{\tau \tau}})}^2} + {5 \over 8}\,{}_4{h^{\tau \tau}}v_A^2 - \,{}_4{h^\tau}_kv_A^k + {3 \over 8}v_A^4} \right]\quad \quad \quad} \\ {+ {\varepsilon ^5}{1 \over 4}\left[ {{}_7{h^{\tau \tau}} + \,{}_5{h^k}_k} \right]\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {\quad \quad + {\varepsilon ^6}\left[ {{1 \over 4}\,{}_8{h^{\tau \tau}} + {1 \over 4}\,{}_6{h^k}_k + {1 \over {16}}\,{}_4{h^{\tau \tau}}_4{h^k}_k - {3 \over {16}}\,{}_4{h^{\tau \tau}}_6{h^{\tau \tau}} + {7 \over {128}}{{{(_4}{h^{\tau \tau}})}^3} + {1 \over 4}\,{}_4{h^{\tau k}}_4{h^\tau}_k} \right.} \\ {\quad \quad {- _6}{h^\tau}_kv_A^k - {1 \over 4}\,{}_4{h^{\tau \tau}}_4{h^\tau}_kv_A^k + {1 \over 2}\,{}_4{h_{kl}}v_A^kv_A^l + {1 \over 8}\,{}_4{h^k}_kv_A^2 + {5 \over {64}}{{{(_4}{h^{\tau \tau}})}^2}v_A^2} \\ {\left. {+ {5 \over 8}\,{}_6{h^{\tau \tau}}v_A^2 - {3 \over 2}\,{}_4{h^\tau}_kv_A^kv_A^2 + {{27} \over {32}}\,{}_4{h^{\tau \tau}}v_A^4 + {5 \over {16}}v_A^6} \right]\quad \quad \quad \quad \quad} \\ {\quad + {\mathcal O}({\varepsilon ^7}).\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ \end{array}$$
(152)

The field should be evaluated somehow at \({\vec z_A}\). This is a formal series since the metric derived via the point particle description diverges at \({\vec z_A}\).

Now let us regularize this equation with the Hadamard partie finie regularization (see [88, 143] and for example, [26, 30] in the literature of the post-Newtonian approximation). Consider a function \(f(\vec x)\) which can be expanded around \({\vec z_1}\) in the form

$$f({\vec r_1}) = \sum\limits_{p = {p_0}} {{1 \over {r_1^p}}} \underset {p}{\hat f} ({\vec N_1}).$$
(153)

then the Hadamard partie finie at \({\vec z_1}\) of the function \(f(\vec x)\) is defined by

$$\oint {{{d{\Omega _{{{\bf{n}}_{\bf{1}}}}}} \over {4\pi}}} \underset {0}{\hat f}({\vec n_1}).$$
(154)

For example, by this procedure hττ becomes (see Equation (130))

$$[{h^{\tau \tau}}]_1^H = 4{\epsilon ^4}{{P_2^\tau} \over {{r_{12}}}} + {\epsilon ^6}\left[ {- 2{{{m_2}} \over {{r_{12}}}}\{{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} - v_2^2\} - 16{{{m_1}{m_2}} \over {r_{12}^2}} + 7{{m_2^2} \over {r_{12}^2}}} \right] + {\mathcal O}({\epsilon ^7})$$
(155)

for star 1. In the above equation, \([f]_A^{\rm{H}}\) means that we regularize the quantity f at star A by the Hadamard partie finie. Evaluating Equation (152) and \([\sqrt {- g} ]_A^{\rm{H}}\) by this procedure, then comparing the result to Equation (148) combined with Equations (149, 150, 151), we find at least up to 3 PN order:

$$P_{A\Theta}^\tau = {m_A}[\sqrt {- g} u_A^\tau ]_A^{\rm{H}}.$$
(156)

It is important that even after regularizing all the divergent terms, there remains a nonlinear effect. Equation (156) is natural. And note that we have never assumed this relation in advance. This relation has been derived by solving the evolution equation for \(P_{A\Theta}^\tau\) functionally. In this regard, it is worth mentioning that the naturality of Equation (156) supports the use of the Hadamard partie finie regularization (or any regularization procedures if we can reproduce Equation (156) with them) to derive the 3 PN mass-energy relation to deal with divergences when one uses Dirac delta distributionsFootnote 11.

Finally, we note that up to 3 PN order

$$[\sqrt {- g} u_A^\tau ]_A^{\rm{H}} = [\sqrt {- g} ]_A^{\rm{H}}[u_A^\tau ]_A^{\rm{H}} + {\mathcal O}({\epsilon ^7})$$
(157)

is satisfied if we use the Hadamard partie finie regularization explained above.

7 Third Post-Newtonian Momentum-Velocity Relation

We now derive the 3 PN momentum-velocity relation by calculating the \(Q_A^i\) integral at 3 PN order. From the definition of the \(Q_A^i\) integral, Equation (91),

$$Q_A^i = {\epsilon ^6}\oint\nolimits_{\partial {B_A}} {d{S_k}({}_{10}\Lambda _N^{\tau k} - v_A^k{}_{10}\Lambda _N^{\tau \tau})y_A^i,}$$
(158)

we find that the calculation required is almost the same as that in the equation for \(dP_A^\tau/d\tau\). Namely, it turns out that we do not need to use the direct-integration method to compute the \(Q_A^i\) integral. Therefore it is straightforward to evaluate the surface integrals in the definition of \(Q_A^i\). Here we split \(Q_A^i\) into \(Q_{A\Theta}^i\) and \(Q_{A\chi}^i\) as

$$Q_A^i = Q_{A\Theta}^i + Q_{A\chi}^i,$$
(159)

with

$$Q_{A\Theta}^i = {\epsilon ^{- 4}}\oint\nolimits_{\partial {B_A}} d {S_k}\left({\Theta _N^{\tau k} - v_A^k\Theta _N^{\tau \tau}} \right)y_A^i,$$
(160)
$$Q_{A\chi}^i = {\epsilon ^{- 4}}\oint\nolimits_{\partial {B_A}} d {S_k}\left({\chi _N^{\tau k\alpha \beta}{}_{,\alpha \beta} - v_A^k\chi _N^{\tau \tau \alpha \beta}{}_{,\alpha \beta}} \right)y_A^i,$$
(161)

and we show only \(Q_{A\Theta}^i\):

$${}_{\leq 6}Q_{1\Theta}^i = - {\epsilon ^6}{{m_1^3\,{m_2}\,n_{12}^{\langle ij\rangle}v_{12}^j} \over {2r_{12}^3}} = - {\epsilon ^6}{d \over {d\tau}}\left({{{m_1^3{m_2}} \over {6r_{12}^3}}r_{12}^i} \right) = {\epsilon ^6}{d \over {d\tau}}\left({{1 \over 6}m_1^3a_1^i} \right),$$
(162)

where nf is the quantity f up to order n inclusively. Here it should be understood that \(a_A^i\) in the last expression is evaluated with the Newtonian acceleration.

\(Q_A^i\) of \({\mathcal O}({\epsilon ^6})\) appears at the 4 PN or higher order field. Thus up to 3 PN order, \(_6Q_A^i\) affects the equations of motion only through the 3 PN momentum-velocity relation. For this reason, not \(Q_{A\chi}^i\) but \(Q_{A\Theta}^i\) is necessary to derive the 3 PN equations of motion. The explicit expression for \(Q_{A\chi}^i\) is given in [91].

Now with \(Q_{A\Theta}^i\) in hand, we obtain the momentum-velocity relation. It turns out that the χ part of the momentum velocity relation is a trivial identityFootnote 12. Thus, defining the Θ parts of \(P_A^\mu\) and \(D_A^i\) in the same way as for \(Q_A^i\), we obtain

$$P_{1\Theta}^i = P_{1\Theta}^\tau v_1^i + {\epsilon ^6}{d \over {d\tau}}\left({{1 \over 6}m_1^3a_1^i} \right) + {\epsilon ^2}{{dD_{1\Theta}^i} \over {d\tau}}.$$
(163)

As explained in the previous Sections 3.4 and 4.4, we define the representative point \(z_A^i\) of star A by choosing the value of \(D_A^i\). In other words, one can freely choose \(D_A^i\) in principleFootnote 13. One may set \(D_A^i\) equal to zero up to 2.5 PN order. Alternatively, one may find it “natural” to see a three-momentum proportial to a three-velocity and take another choice,

$$D_{A\Theta}^i = - {\epsilon ^4}{1 \over 6}m_A^3a_A^i = {\epsilon ^4}\delta _{A\Theta}^i.$$
(164)

Henceforth, we shall define \(z_A^i\) by this equation.

Finally, it is important to realize that the nonzero dipole moment \(D_A^i\) of order ϵ4 affects the 3 PN field and the 3 PN equations of motion in essentially the same manner as the Newtonian dipole moment affects the Newtonian field and equations of motion. From Equations (78, 79, 80) we see that \(\delta _{A\Theta}^i\) appears only at 10hττ as

$${h^{\tau \tau}}{\vert _{{\delta _{A\Theta}}}} = 4{\epsilon ^{10}}\sum\limits_{A = 1,2} {{{\delta _{A\Theta}^kr_A^k} \over {r_A^3}}} + {\mathcal O}({\epsilon ^{11}}){.}$$
(165)

then the corresponding acceleration becomes

$${m_1}a_1^i{\vert _{{\delta _{A\Theta}}}} = - {\epsilon ^6}{{3{m_1}\delta _{2\Theta}^i} \over {r_{12}^3}}n_{12}^{\langle ik\rangle} + {\epsilon ^6}{{3{m_2}\delta _{1\Theta}^i} \over {r_{12}^3}}n_{12}^{\langle ik\rangle} - {\epsilon ^6}{{{d^2}\delta _{1\Theta}^i} \over {d{\tau ^2}}}.$$
(166)

The last term compensates the \(Q_A^i\) integral contribution appearing through the momentum-velocity relation (163).

Note that this change of the acceleration does not affect the existence of the conservation of the (Newtonian-sense) energy,

$${m_1}a_1^i{\vert _{{\delta _{A\Theta}}}}v_1^i + {m_2}a_2^i{\vert _{{\delta _{A\Theta}}}}v_2^i = {\epsilon ^6}{d \over {d\tau}}\sum\limits_{A = 1,2} {\left[ {\delta _{A\Theta}^k{{dv_A^k} \over {d\tau}} - v_A^k{d \over {d\tau}}\delta _{A\Theta}^k} \right]}.$$
(167)

8 Third Post-Newtonian Equations of Motion

8.1 Third post-Newtonian equations of motion with logarithmic terms

To derive the 3 PN equations of motion, we evaluate the surface integrals in the general form of the equations of motion (111) using the field 8hττ, the field ≤6hμν, the 3 PN body zone contributions, and the 3 PN N/B contributions corresponding to the results from the super-potential method and the super-potential-in-series method. We then combine the result with the terms from the direct-integration method.

From 3 PN order, the effects of the \(Q_A^{{K_l}i}\) and \(R_A^{{K_l}ij}\) integrals appearing in the 3 PN field in \(h_B^{\mu \nu}\) contribute to the 3 PN equations of motion. \(_6Q_{A\Theta}^i\) given in Equation (162) affects the 3 PN equations of motion through the 3 PN momentum-velocity relation. Since we define the representative points of the stars via Equation (164), we add the corresponding acceleration given by Equation (166). Furthermore, our choice of the representative points of the stars makes \(D_{A\chi}^i\) appear independently of \(D_{A\Theta}^i\) in the field, and hence \(_4D_{A\chi}^i\) affects the 3 PN equations of motion. In summary, the \(_{\leq 4}Q_A^{{K_l}i}{,_{\; \leq 4}}R_A^{{K_l}ij},{\;_6}Q_{A\Theta}^i,\;\delta _{A\Theta}^i\) and \(_4D_{A\chi}^i\) contributions to the 3 PN field can be written as

$${}_{10}{h^{\tau \tau}} + {}_8{h^k}_k = 4\sum\limits_{A = 1,2} {{{r_A^k} \over {r_A^3}}} \left({\delta _{A\Theta}^k + {}_4D_{A\chi}^k + {}_4R_A^{kll} - {1 \over 2}{}_4R_A^{llk}} \right) + \ldots,$$
(168)

where “…” are other contributions. On the other hand, \(_6Q_{A\Theta}^i\) and \(\delta _{A\Theta}^i\) affect the equations of motion through the momentum-velocity relation in Equation (111),

$${m_1}{\left({{{dv_1^i} \over {d\tau}}} \right)_{\leq 3{\rm{PN}}}} = - {\epsilon ^6}{{{d_6}Q_{1\Theta}^i} \over {d\tau}} - {\epsilon ^6}{{{d^2}\delta _{1\Theta}^i} \over {d{\tau ^2}}} + \ldots,$$
(169)

but they cancel each other out, since we choose Equation (164). Then these contributions to a 3 PN acceleration can be summarized into

$$\left({\begin{array}{*{20}c} {{\rm{the}}\;{\rm{contribution}}\;{\rm{to}}\;{m_1}a_1^i\;{\rm{from}}\;}\\ {{}_{\leq 4}Q_A^{{K_l}i},{}_{\leq 4}R_A^{{K_l}ij},{}_6Q_{A\Theta}^i,\delta _{A\Theta}^i,{}_4D_{A\chi}^i}\\ \end{array}} \right) = {\epsilon ^6}{{118} \over 9}{{m_1^3m_2^2} \over {r_{12}^6}}r_{12}^i + {\epsilon ^6}{{118} \over 9}{{m_1^2m_2^3} \over {r_{12}^6}}r_{12}^i.$$
(170)

Collecting these contributions mentioned above, we obtain the 3 PN equations of motion. However, we found that logarithmic terms having the arbitrary constants ϵRA in their arguments survive,

$$\begin{array}{*{20}c} {{m_1}{{\left({{{dv_1^i} \over {d\tau}}} \right)}^{{\rm{with}}\,{\rm{log}}}} = {m_1}{{\left({{{dv_1^i} \over {d\tau}}} \right)}_{\le 2.5{\rm{PN}}}}\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {+ {\epsilon ^6}{{m_1^2{m_2}} \over {r_{12}^3}}\left[ {{{44m_1^2} \over {3r_{12}^2}}n_{12}^i\ln \left({{{{r_{12}}} \over {\epsilon {R_1}}}} \right) - {{44m_2^2} \over {3r_{12}^2}}n_{12}^i\ln \left({{{{r_{12}}} \over {\epsilon {R_2}}}} \right)\quad \quad \quad \quad \quad \quad} \right.\quad \quad \quad \quad \quad} \\ {+ \left. {{{22{m_1}} \over {{r_{12}}}}\left({5{{({{\vec n}_{12}}\cdot\vec V)}^2}n_{12}^i - {V^2}n_{12}^i - 2({{\vec n}_{12}}\cdot\vec V){V^i}} \right)\ln \left({{{{r_{12}}} \over {\epsilon {R_1}}}} \right)} \right]} \\ {+ \ldots + {\mathcal O}({\epsilon ^7}),\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ \end{array}$$
(171)

where the acceleration through 2.5 PN order, \({(dv_1^i/d\tau)_{\leq 2.5{\rm{PN}}}}\), is the Damour and Deruelle 2.5 PN acceleration. In our formalism, we have computed it in [95]. The “…” stands for the terms that do not include any logarithms.

Since this equation contains two arbitrary constants, the body zone radii RA, at first sight its predictive power on the orbital motion of the binary seems to be limited. In the next Section 8.2, we shall show that by a reasonable redefinition of the representative points of the stars, we can remove RA from our equations of motion. There, we show the explicit form of the 3 PN equations of motion we obtained.

8.2 Arbitrary constant ϵRA

The reason logarithms appear in the 3 PN equations of motion (171) is in a sense easy to understand. Since the post-Newtonian approximation is a weak field expansion, at some level of iteration in the evaluation of field, we have inevitably logarithms in the field through volume integrals such as

$$\int\nolimits_{N/B} {{{{d^3}y} \over {\vert \vec x - \vec y\vert}}} {\left({{m \over y}} \right)^n},$$
(172)

where m and y are typical the mass and length scale in the orbital motion, respectively, and n is a positive integer. For instance consider \(m_1^3({\vec r_1} \cdot {\vec a_1})/r_1^5\) as an integrand. Then we find

$$\int\nolimits_{N/B} {{{{d^3}y} \over {\vert \vec x - \vec y\vert}}} {{m_1^3({{\vec y}_1}\cdot{{\vec a}_1})} \over {y_1^5}} = {{m_1^3({{\vec r}_1}\cdot{{\vec a}_1})} \over {3r_1^3}}\ln \left({{{{r_1}} \over {\epsilon {R_1}}}} \right) + \ldots.$$
(173)

Actually, we could remove the ϵRA dependence in our 3 PN equations of motion via an alternative choice of the center of mass. The following alternative choice of the representative point of star A removes the ϵRA dependence in Equation (171):

$$D_{A\Theta, {\rm{new}}}^i(\tau) = {\epsilon ^4}\delta _{A\Theta}^i(\tau) - {\epsilon ^4}{{22} \over 3}m_A^3a_A^i\ln \left({{{{r_{12}}} \over {\epsilon {R_A}}}} \right) \equiv {\epsilon ^4}\delta _{A\Theta}^i(\tau) + {\epsilon ^4}\delta _{A\ln}^i(\tau) \equiv {\epsilon ^4}\delta _A^i(\tau){.}$$
(174)

Note that this redefinition of the center of mass does not affect the existence of the energy conservation as was shown by Equation (167). We can examine the effect of this redefinition on the equations of motion using Equation (166) (using \(\delta _{A\;\ln}^i\) instead of δAΘ). Then we have

$$\begin{array}{*{20}c} {{m_1}a_1^i{\vert _{{\delta _{A\ln}}}} = - {\epsilon ^6}{{3{m_1}\delta _{2\ln}^i} \over {r_{12}^3}}n_{12}^{\langle ik\rangle} + {\epsilon ^6}{{3{m_2}\delta _{1\ln}^i} \over {r_{12}^3}}n_{12}^{\langle ik\rangle} - {\epsilon ^6}{{{d^2}\delta _{1\ln}^i} \over {d{\tau ^2}}}\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad}\\ {= - {{44m_1^3m_2^2} \over {3r_{12}^5}}n_{12}^i\ln \left({{{{r_{12}}} \over {\epsilon {R_1}}}} \right) + {{44m_1^2m_2^3} \over {3r_{12}^5}}n_{12}^i\ln \left({{{{r_{12}}} \over {\epsilon {R_2}}}} \right)\quad \quad \quad \quad \quad \quad \quad \quad}\\ {- {{22m_1^3{m_2}} \over {r_{12}^4}}\left({5{{({{\vec n}_{12}}\cdot\vec V)}^2}n_{12}^i - {V^2}n_{12}^i - 2({{\vec n}_{12}}\cdot\vec V){V^i}} \right)\ln \left({{{{r_{12}}} \over {\epsilon {R_1}}}} \right)\quad \quad \quad \quad}\\ {+ {{22m_1^3{m_2}} \over {3r_{12}^4}}\left({{{{m_1}} \over {{r_{12}}}}n_{12}^i + {{{m_2}} \over {{r_{12}}}}n_{12}^i - {V^2}n_{12}^i + 8{{({{\vec n}_{12}}\cdot\vec V)}^2}n_{12}^i - 2{{({{\vec n}_{12}}\cdot\vec V)}^2}{V^i}} \right).}\\ \end{array}$$
(175)

Comparing the above equations with Equation (171), we easily conclude that the representative point \(z_A^i\) of star A defined by

$$D_{A\Theta, {\rm{new}}}^i(\tau) = {\epsilon ^{- 6}}\int\nolimits_{{B_A}} {{d^3}} y\,({y^i} - z_A^i(\tau))\,\Theta _N^{\tau \tau}(\tau, {y^k}) = {\epsilon ^4}\delta _A^i(\tau)$$
(176)

obeys the equations of motion free from any logarithmic term and hence free from any ambiguity up to 3 PN order inclusively. We note that in our formalism \(z_A^i\) is defined by the value of \(D_A^i\) and in turn we have a freedom to assign to \(D_A^i\) any value as we like (though it may be natural to set the value of \(D_A^i\) such that \(z_A^i\) resides inside star A). We also note that we define \(z_A^i\) order by order.

We mention here that Blanchet and Faye [27] have already noticed that in their 3 PN equations of motion a suitable coordinate transformation removes (parts of) the logarithmic dependence of arbitrary parameters corresponding (roughly) to our body zone radiiFootnote 14. It is well-known that choosing different values of dipole moments corresponds to a coordinate transformation.

8.3 Consistency relation

Our new choice of the dipole moments of the stars is in a sense natural. To see this, let us consider the harmonic condition:

$$0 = {h^{\tau \mu}}_{,\mu} = 4{\epsilon ^4}\sum\limits_{A = 1,2} {\left[ {{{\dot P_A^\tau} \over {{r_A}}} + {{r_A^i} \over {r_A^3}}\left({P_A^\tau v_A^i + {\epsilon ^2}\dot D_A^i - P_A^i} \right)} \right] + \sum\limits_{A = 1,2} {\oint\nolimits_{\partial {B_A}} {{{d{S_i}} \over {\vert \vec x - \vec y\vert}}(\Lambda _N^{\tau i} - v_A^i\Lambda _N^{\tau \tau}) \ldots,}}}$$
(177)
$$0 = {h^{i\mu}}_{,\mu} = 4{\epsilon ^4}\sum\limits_{A = 1,2} {{{\dot P_A^i} \over {{r_A}}} + \sum\limits_{A = 1,2} {\oint\nolimits_{\partial {B_A}} {{{d{S_j}} \over {\vert \vec x - \vec y\vert}}\left(\Lambda _N^{ij} - v_A^j\Lambda _N^{\tau i}\right) + \ldots,}}}$$
(178)

where “…” are irrelevant terms. These equations are a manifestation of the fact that the harmonic condition is consistent with the evolution equation for \(P_A^\tau\), the momentum-velocity relation, and the equations of motion (and relations among higher multipole moments, hidden in “…”). Thus if the logarithmic dependence of ϵRA arises from the second term of Equation (178), \(P_A^i\) must have the same logarithmic dependence (times minus sign) to ensure harmonicity. This and the momentum velocity relation in turn mean \(P_A^\tau\), \(v_A^i = \dot z_A^i\) or \(D_A^i\) have corresponding logarithmic dependence. Since we already know that the \(P_A^\tau\) have no logarithm up to 3 PN order, \(z_A^i\) or \(D_A^i\) should have logarithms. This is consistent with the fact that a choice of \(D_A^i\) determines \(z_A^i\). \(z_A^i\) depends on logarithms if the old choice is taken, while it does not if our new choice is taken.

There is yet another fact which supports our interpretation. Let us retain \(D_A^i \neq 0\) for a while. We then find that the near zone dipole moment \(D_N^i\) defined by a volume integral of \(\Lambda _N^{\tau \tau}{y^i}\) becomes

$${\epsilon ^2}D_N^i \equiv {\epsilon ^{- 4}}\int\nolimits_N {{d^3}} y\,\Lambda _N^{\tau \tau}{y^i} = \sum\limits_{A = 1,2} {P_A^\tau} z_A^i + {\epsilon ^2}\sum\limits_{A = 1,2} {D_A^i} + {\epsilon ^{- 4}}\int\nolimits_{N/B} {{d^3}} y\,\Lambda _N^{\tau \tau}{y^i}.$$
(179)

then if we take the old choice of \(D_A^i\), the volume integral becomes

$$\int\nolimits_{N/B} {{d^3}} y\,\Lambda _N^{\tau \tau}{y^i} = {\epsilon ^4}{{22} \over 3}\sum\limits_{A = 1,2} m_A^3a_A^i\ln \left({{{{r_{12}}} \over {\epsilon {R_A}}}} \right) + \ldots,$$
(180)

where terms denoted by “…” have no logarithmic dependence. Notice that the near zone dipole moment can be freely determined, say, \(D_N^i = 0\) because we can define the origin of the near zone freely. By taking temporal derivatives twice of \(D_N^i\) we see that \(D_{A,{\rm{new}}}^i\) gives a natural definition of the center of the mass in terms of which the 3 PN equations of motion are independent of ϵRA.

8.4 Third post-Newtonian equations of motion

By adding \({m_1}a_1^i{\vert _{{\delta _{A\;\ln}}}}\) to Equation (171), we obtain our 3 PN equations of motion for two spherical compact stars whose representative points are defined by Equation (176),

$$\begin{array}{*{20}c} {{m_1}{{dv_1^i} \over {d\tau}} = - {{{m_1}{m_2}} \over {r_{12}^2}}n_{12}^i\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {+ {\epsilon ^2}{{{m_1}{m_2}} \over {r_{12}^2}}n_{12}^i\left[ {- v_1^2 - 2v_2^2 + 4({{\vec v}_1}\cdot{{\vec v}_2}) + {3 \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} + {{5{m_1}} \over {{r_{12}}}} + {{4{m_2}} \over {{r_{12}}}}} \right]\quad} \\ {+ {\epsilon ^2}{{{m_1}{m_2}} \over {r_{12}^2}}{V^i}\left[ {4({{\vec n}_{12}}\cdot{{\vec v}_1}) - 3({{\vec n}_{12}}\cdot{{\vec v}_2})} \right]\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {\quad \quad \quad \quad + {\epsilon ^4}{{{m_1}{m_2}} \over {r_{12}^2}}n_{12}^i\left[ {- 2v_2^4 + 4v_2^2({{\vec v}_1}\cdot{{\vec v}_2}) - 2{{({{\vec v}_1}\cdot{{\vec v}_2})}^2} + {3 \over 2}v_1^2{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} + {9 \over 2}v_2^2{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}} \right.} \\ {\quad \quad \quad \quad \quad \quad \quad \quad \quad - 6({{\vec v}_1}\cdot{{\vec v}_2}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} - {{15} \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^4} - {{57} \over 4}{{m_1^2} \over {r_{12}^2}} - 9{{m_2^2} \over {r_{12}^2}} - {{69} \over 2}{{{m_1}{m_2}} \over {r_{12}^2}}} \\ {\quad \quad \quad + {{{m_1}} \over {{r_{12}}}}\left({- {{15} \over 4}v_1^2 + {5 \over 4}v_2^2 - {5 \over 2}({{\vec v}_1}\cdot{{\vec v}_2}) + {{39} \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}} \right.} \\ {\quad \quad \quad \quad - \left. {- 39({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2}) + {{17} \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}} \right)} \\ {\left. {\quad \quad \quad \quad \quad \quad \quad + {{{m_2}} \over {{r_{12}}}}\left({4v_2^2 - 8({{\vec v}_1}\cdot{{\vec v}_2}) + 2{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2} - 4({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2}) - 6{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}} \right)} \right]} \\ {+ {\epsilon ^4}{{{m_1}{m_2}} \over {r_{12}^2}}{V^i}\left[ {{{{m_1}} \over {{r_{12}}}}\left({- {{63} \over 4}({{\vec n}_{12}}\cdot{{\vec v}_1}) + {{55} \over 4}({{\vec n}_{12}}\cdot{{\vec v}_2})} \right) + {{{m_2}} \over {{r_{12}}}}\left({- 2({{\vec n}_{12}}\cdot{{\vec v}_1}) - 2({{\vec n}_{12}}\cdot{{\vec v}_2})} \right)} \right.} \\ {\quad \quad \quad \quad + v_1^2({{\vec n}_{12}}\cdot{{\vec v}_2}) + 4v_2^2({{\vec n}_{12}}\cdot{{\vec v}_1}) - 5v_2^2({{\vec n}_{12}}\cdot{{\vec v}_2}) - 4({{\vec v}_1}\cdot{{\vec v}_2})({{\vec n}_{12}}\cdot\vec V)} \\ {- \left. {6({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} + {9 \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3}} \right]\quad \quad \quad \quad \quad} \\ {+ {\epsilon ^5}{{4m_1^2{m_2}} \over {5r_{12}^3}}\left[ {n_{12}^i({{\vec n}_{12}}\cdot\vec V)\left({- 6{{{m_1}} \over {{r_{12}}}} + {{52} \over 3}{{{m_2}} \over {{r_{12}}}} + 3{V^2}} \right) + {V^i}\left({2{{{m_1}} \over {{r_{12}}}} - 8{{{m_2}} \over {{r_{12}}}} - {V^2}} \right)} \right]} \\ {+ {\epsilon ^6}{{{m_1}{m_2}} \over {r_{12}^2}}n_{12}^i\left[ {{{35} \over {16}}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^6} - {{15} \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^4}v_1^2 + {{15} \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^4}({{\vec v}_1}\cdot{{\vec v}_2})} \right.\quad \quad \quad} \\ {+ 3{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}{{({{\vec v}_1}\cdot{{\vec v}_2})}^2} - {{15} \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^4}v_2^2 + {3 \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}v_1^2v_2^2} \\ {- 12{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}({{\vec v}_1}\cdot{{\vec v}_2})v_2^2 - 2{{({{\vec v}_1}\cdot{{\vec v}_2})}^2}v_2^2 + {{15} \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}v_2^4} \\ {+ 4({{\vec v}_1}\cdot{{\vec v}_2})v_2^4 - 2v_2^6\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {+ {{{m_1}} \over {{r_{12}}}}\left({- {{171} \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^4} + {{171} \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^3}({{\vec n}_{12}}\cdot{{\vec v}_2})\quad \quad \quad} \right.} \\ {\quad \quad \quad - {{723} \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} + {{383} \over 2}({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3}} \\ {\quad \quad \quad \quad \quad \quad - {{455} \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^4} + {{229} \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}v_1^2 - {{205} \over 2}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2})v_1^2} \\ {\,\quad \quad + {{191} \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}v_1^2 - {{91} \over 8}v_1^4 - {{229} \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}({{\vec v}_1}\cdot{{\vec v}_2})} \\ {\quad \quad \quad + 244({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2})({{\vec v}_1}\cdot{{\vec v}_2}) - {{225} \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}({{\vec v}_1}\cdot{{\vec v}_2})} \\ {\,\quad \quad + {{91} \over 2}v_1^2({{\vec v}_1}\cdot{{\vec v}_2}) - {{177} \over 4}{{({{\vec v}_1}\cdot{{\vec v}_2})}^2} + {{229} \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}v_2^2} \\ {\quad \quad \quad \quad - {{283} \over 2}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2})v_2^2 + {{259} \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}v_2^2 - {{91} \over 4}v_1^2v_2^2} \\ {\, + \left. {43({{\vec v}_1}\cdot{{\vec v}_2})v_2^2 - {{81} \over 8}v_2^4} \right)\quad \quad \quad \quad \quad \quad \quad \quad} \\ {+ {{{m_2}} \over {{r_{12}}}}\left({- 6{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} + 12({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3}\quad \quad \quad \quad} \right.} \\ {+ 6{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^4} + 4({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2})({{\vec v}_1}\cdot{{\vec v}_2})\quad \quad \quad} \\ {+ 12{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}({{\vec v}_1}\cdot{{\vec v}_2}) + 4{{({{\vec v}_1}\cdot{{\vec v}_2})}^2}\quad \quad \quad \quad \quad} \\ {\quad \quad \quad \quad \quad - \left. {4({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2})v_2^2 - 12{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}v_2^2 - 8({{\vec v}_1}\cdot{{\vec v}_2})v_2^2 + 4v_2^4} \right)} \\ {+ {{m_2^2} \over {r_{12}^2}}\left({- {{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2} + 2({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2}) + {{43} \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}} \right.\quad \quad \quad} \\ {\left. {+ 18({{\vec v}_1}\cdot{{\vec v}_2}) - 9v_2^2} \right)} \\ {+ {{{m_1}{m_2}} \over {r_{12}^2}}\left({{{415} \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2} - {{375} \over 4}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2}) + {{1113} \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}} \right.} \\ {\quad \quad \quad \quad - \left. {{{615{\pi ^2}} \over {64}}{{({{\vec n}_{12}}\cdot\vec V)}^2} + 18v_1^2 + {{123{\pi ^2}} \over {64}}{V^2} + 33({{\vec v}_1}\cdot{{\vec v}_2}) - {{33} \over 2}v_2^2} \right)} \\ {\quad + {{m_1^2} \over {r_{12}^2}}\left({- {{2069} \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2} + 543({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2}) - {{939} \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}} \right.} \\ {\quad + \left. {{{471} \over 8}v_1^2 - {{357} \over 4}({{\vec v}_1}\cdot{{\vec v}_2}) + {{357} \over 8}v_2^2} \right)\quad \quad \quad \quad \quad} \\ {\quad + \left. {{{16m_2^3} \over {r_{12}^3}} + {{m_1^2{m_2}} \over {r_{12}^3}}\left({{{547} \over 3} - {{41{\pi ^2}} \over {16}}} \right) - {{13m_1^3} \over {12r_{12}^3}} + {{{m_1}m_2^2} \over {r_{12}^3}}\left({{{545} \over 3} - {{41{\pi ^2}} \over {16}}} \right)} \right]} \\ {+ {\epsilon ^6}{{{m_1}{m_2}} \over {r_{12}^2}}{V^i}\left[ {{{15} \over 2}({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^4} - {{45} \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^5} - {3 \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3}v_1^2\quad \quad \quad \quad \quad \quad \quad \quad \quad} \right.} \\ {+ 6({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}({{\vec v}_1}\cdot{{\vec v}_2}) - 6{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3}({{\vec v}_1}\cdot{{\vec v}_2})\quad \quad \quad \quad \quad} \\ {- 2({{\vec n}_{12}}\cdot{{\vec v}_2}){{({{\vec v}_1}\cdot{{\vec v}_2})}^2} - 12({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}v_2^2 + 12{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3}v_2^2} \\ {+ ({{\vec n}_{12}}\cdot{{\vec v}_2})v_1^2v_2^2 - 4({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec v}_1}\cdot{{\vec v}_2})v_2^2 + 8({{\vec n}_{12}}\cdot{{\vec v}_2})({{\vec v}_1}\cdot{{\vec v}_2})v_2^2\quad} \\ {+ 4({{\vec n}_{12}}\cdot{{\vec v}_1})v_2^4 - 7({{\vec n}_{12}}\cdot{{\vec v}_2})v_2^4\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {+ {{{m_2}} \over {{r_{12}}}}\left({- 2{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}({{\vec n}_{12}}\cdot{{\vec v}_2}) + 8({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} + 2{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3}} \right.} \\ {+ 2({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec v}_1}\cdot{{\vec v}_2}) + 4({{\vec n}_{12}}\cdot{{\vec v}_2})({{\vec v}_1}\cdot{{\vec v}_2})\quad \quad \quad \quad} \\ {- \left. {2({{\vec n}_{12}}\cdot{{\vec v}_1})v_2^2 - 4({{\vec n}_{12}}\cdot{{\vec v}_2})v_2^2} \right)\quad \quad \quad \quad \quad \quad \quad} \\ {+ {{{m_1}} \over {{r_{12}}}}\left({- {{243} \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^3} + {{565} \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}({{\vec n}_{12}}\cdot{{\vec v}_2})} \right.\quad \quad \quad \quad \quad \quad} \\ {\quad \quad \quad - {{269} \over 4}({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} - {{95} \over {12}}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3} + {{207} \over 8}({{\vec n}_{12}}\cdot{{\vec v}_1})v_1^2} \\ {\quad \quad \quad \quad - {{137} \over 8}({{\vec n}_{12}}\cdot{{\vec v}_2})v_1^2 - 36({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec v}_1}\cdot{{\vec v}_2}) + {{27} \over 4}({{\vec n}_{12}}\cdot{{\vec v}_2})({{\vec v}_1}\cdot{{\vec v}_2})} \\ {+ \left. {{{81} \over 8}({{\vec n}_{12}}\cdot{{\vec v}_1})v_2^2 + {{83} \over 8}({{\vec n}_{12}}\cdot{{\vec v}_2})v_2^2} \right)\quad \quad \quad \quad \quad} \\ {+ {{m_2^2} \over {r_{12}^2}}\left({4({{\vec n}_{12}}\cdot{{\vec v}_1}) + 5({{\vec n}_{12}}\cdot{{\vec v}_2})} \right)\quad \quad \quad \quad \quad \quad} \\ {\quad \quad \quad \quad + {{{m_1}{m_2}} \over {r_{12}^2}}\left({- {{307} \over 8}({{\vec n}_{12}}\cdot{{\vec v}_1}) + {{479} \over 8}({{\vec n}_{12}}\cdot{{\vec v}_2}) + {{123{\pi ^2}} \over {32}}({{\vec n}_{12}}\cdot\vec V)} \right)} \\ {+ \left. {{{m_1^2} \over {r_{12}^2}}\left({{{311} \over 4}({{\vec n}_{12}}\cdot{{\vec v}_1}) - {{357} \over 4}({{\vec n}_{12}}\cdot{{\vec v}_2})} \right)} \right]\quad \quad \quad} \\ {+ {\mathcal O}({\epsilon ^7}),\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ \end{array}$$
(181)

in the harmonic gauge.

Now we list some features of our 3 PN equations of motion. In the test-particle limit, our 3 PN equations of motion coincide with a geodesic equation for a test-particle in the Schwarzschild metric in harmonic coordinates (up to 3 PN order). Suppose that star 1 is a test particle, star 2 is represented by the Schwarzschild metric, and \({\vec v_2} = 0\). Then the geodesic equation for star 1 becomes

$$\begin{array}{*{20}c} {a_1^i = - {{{m_2}} \over {r_{12}^2}}n_{12}^i + {\epsilon^2}{{{m_2}} \over {r_{12}^2}}\left({- v_1^2n_{12}^i + {{4{m_2}} \over {{r_{12}}}}n_{12}^i + 4({{\vec n}_{12}}\cdot{{\vec v}_1})v_1^i} \right)\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad}\\ {+ {\epsilon^4}{{m_2^2} \over {r_{12}^3}}\left({2{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}n_{12}^i - 9{{{m_2}} \over {{r_{12}}}}n_{12}^i - 2({{\vec n}_{12}}\cdot{{\vec v}_1})v_1^i} \right)\quad \quad \quad \quad}\\ {+ {\epsilon^6}{{m_2^3} \over {r_{12}^4}}\left({- {{({{\vec n}_{12}}\cdot{{\vec v}_1})}^3}n_{12}^i + {{16{m_2}} \over {{r_{12}}}}n_{12}^i + 4({{\vec n}_{12}}\cdot{{\vec v}_1})v_1^i} \right) + {\mathcal O}({\epsilon^7})}\\ \end{array}$$
(182)

in the harmonic gauge. Thus, in the test particle limit Equation (181) coincides with the geodesic equation for a test particle in the Schwarzschild metric up to 3 PN order.

With the help of the formulas developed in [28], we have checked the Lorentz invariance of Equation (181) (in the post-Newtonian perturbative sense). Also, we have checked that our 3 PN acceleration admits a conserved energy of the binary orbital motion (modulo the 2.5 PN radiation reaction effect). In fact, the energy E of the binary associated with Equation (181) is

$$\begin{array}{*{20}c} {E = {1 \over 2}{m_1}v_1^2 - {{{m_1}{m_2}} \over {2{r_{12}}}}\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {\quad \quad \quad \quad + {\epsilon ^2}\left[ {{3 \over 8}{m_1}v_1^4 + {{m_1^2{m_2}} \over {2r_{12}^2}} + {{{m_1}{m_2}} \over {2{r_{12}}}}\left({3v_1^2 - {7 \over 2}({{\vec v}_1}\cdot{{\vec v}_2}) - {1 \over 2}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2})} \right)} \right]\quad \quad \quad \quad} \\ {+ {\epsilon ^4}\left[ {{5 \over {16}}{m_1}v_1^6 - {{m_1^3{m_2}} \over {2r_{12}^3}} - {{19m_1^2m_2^2} \over {8r_{12}^3}}\quad} \right.\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {\quad \quad \quad \quad \quad + {{m_1^2{m_2}} \over {2r_{12}^2}}\left({- 3v_1^2 + {7 \over 2}v_2^2 + {{29} \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2} - {{13} \over 2}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2}) + {{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}} \right)} \\ {\quad \quad \quad \quad \quad \quad + {{{m_1}{m_2}} \over {4{r_{12}}}}\left({{3 \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^3}({{\vec n}_{12}}\cdot{{\vec v}_2}) + {3 \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} - {9 \over 2}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2})v_1^2} \right.} \\ {\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad - {{13} \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}v_1^2 + {{21} \over 2}v_1^4 + {{13} \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}({{\vec v}_1}\cdot{{\vec v}_2})\quad \quad \quad \quad \quad \quad \quad \quad} \\ {\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \,\,\,+ \left. {\left. {3({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2})({{\vec v}_1}\cdot{{\vec v}_2}) - {{55} \over 2}v_1^2({{\vec v}_1}\cdot{{\vec v}_2}) + {{17} \over 2}{{({{\vec v}_1}\cdot{{\vec v}_2})}^2} + {{31} \over 4}v_1^2v_2^2} \right)} \right]} \\ {+ {\epsilon ^6}\left[ {{{35} \over {128}}{m_1}v_1^8 + {{3m_1^4{m_2}} \over {8r_{12}^4}} + {{469m_1^3m_2^2} \over {18r_{12}^4}}\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \right.} \\ {\quad \quad \quad \quad + {{m_1^2m_2^2} \over {2r_{12}^3}}\left({{{547} \over 6}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2} - {{3115} \over {24}}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2}) - {{123{\pi ^2}} \over {32}}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot\vec V)} \right.} \\ {- \left. {{{575} \over 9}v_1^2 + {{41{\pi ^2}} \over {32}}(\vec V\cdot{{\vec v}_2}) + {{4429} \over {72}}({{\vec v}_1}\cdot{{\vec v}_2})} \right)\quad} \\ {\quad \quad + {{m_1^3{m_2}} \over {2r_{12}^3}}\left({- {{437} \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2} + {{317} \over 4}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2}) + 3{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} + {{301} \over {12}}v_1^2} \right.} \\ {- \left. {{{337} \over {12}}({{\vec v}_1}\cdot{{\vec v}_2}) + {5 \over 2}v_2^2} \right)\quad \quad \quad \quad \quad \quad \quad \quad} \\ {\quad \quad \quad \quad \quad + {{{m_1}{m_2}} \over {{r_{12}}}}\left({- {5 \over {16}}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^5}({{\vec n}_{12}}\cdot{{\vec v}_2}) - {5 \over {16}}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^4}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} - {5 \over {32}}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^3}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3}} \right.} \\ {\quad + {{19} \over {16}}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^3}({{\vec n}_{12}}\cdot{{\vec v}_2})v_1^2 + {{15} \over {16}}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}v_1^2} \\ {\quad \quad \quad \quad \quad + {3 \over 4}({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3}v_1^2{{19} \over {16}}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^4}v_1^2 - {{21} \over {16}}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2})v_1^4} \\ {\quad \quad \quad \quad \quad \quad \quad \quad - 2{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}v_1^4 + {{55} \over {16}}v_1^6 - {{19} \over {16}}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^4}({{\vec v}_1}\cdot{{\vec v}_2}) - {{({{\vec n}_{12}}\cdot{{\vec v}_1})}^3}({{\vec n}_{12}}\cdot{{\vec v}_2})({{\vec v}_1}\cdot{{\vec v}_2})} \\ {\quad \quad - {{15} \over {32}}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}({{\vec v}_1}\cdot{{\vec v}_2}) + {{45} \over {16}}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}v_1^2({{\vec v}_1}\cdot{{\vec v}_2})} \\ {\quad \quad \quad \quad \quad \quad \quad + {5 \over 4}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2})v_1^2({{\vec v}_1}\cdot{{\vec v}_2}) + {{11} \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}v_1^2({{\vec v}_1}\cdot{{\vec v}_2}) - {{139} \over {16}}v_1^4({{\vec v}_1}\cdot{{\vec v}_2})} \\ {\quad \quad \quad \quad \quad \quad - {3 \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}{{({{\vec v}_1}\cdot{{\vec v}_2})}^2} + {5 \over {16}}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2}){{({{\vec v}_1}\cdot{{\vec v}_2})}^2} + {{41} \over 8}v_1^2{{({{\vec v}_1}\cdot{{\vec v}_2})}^2}} \\ {\quad \quad \quad \quad \quad \quad \quad + {1 \over {16}}{{({{\vec v}_1}\cdot{{\vec v}_2})}^3} - {{45} \over {16}}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}v_1^2v_2^2 - {{23} \over {32}}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2})v_1^2v_2^2 + {{79} \over {16}}v_1^4v_2^2} \\ {\qquad \qquad \qquad \quad \, - \left. {{{161} \over {32}}v_1^2v_2^2({{\vec v}_1}\cdot{{\vec v}_2})} \right)\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {\quad \quad \quad + {{m_1^2{m_2}} \over {r_{12}^2}}\left({- {{49} \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^4} + {{75} \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^3}({{\vec n}_{12}}\cdot{{\vec v}_2}) - {{187} \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2} + {{11} \over 2}v_1^4} \right.} \\ {\quad \quad \quad \quad + {{247}\over {24}}({{\vec n}_{12}}\cdot{{\vec v}_1}){{({{\vec n}_{12}}\cdot{{\vec v}_2})}^3} + {{49} \over 8}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}v_1^2 + {{81} \over 8}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2})v_1^2} \\ {\quad \quad \quad \quad - {{21} \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}v_1^2 - {{15} \over 2}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}({{\vec v}_1}\cdot{{\vec v}_2}) - {3 \over 2}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2})({{\vec v}_1}\cdot{{\vec v}_2})} \\ {\quad \quad \quad \quad + {{21} \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}({{\vec v}_1}\cdot{{\vec v}_2}) - 27v_1^2({{\vec v}_1}\cdot{{\vec v}_2}) + {{55} \over 2}{{({{\vec v}_1}\cdot{{\vec v}_2})}^2} + {{49} \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_1})}^2}v_2^2} \\ {\quad \quad \quad - {{27} \over 2}({{\vec n}_{12}}\cdot{{\vec v}_1})({{\vec n}_{12}}\cdot{{\vec v}_2})v_2^2 + {3 \over 4}{{({{\vec n}_{12}}\cdot{{\vec v}_2})}^2}v_2^2 + {{55} \over 4}v_1^2v_2^2 - 28({{\vec v}_1}\cdot{{\vec v}_2})v_2^2} \\ {+ \left. {\left. {{{135} \over {16}}v_2^4} \right)} \right]\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {+ (1 \leftrightarrow 2) + {\mathcal O}({\epsilon ^7}).\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ \end{array}$$
(183)

This orbital energy of the binary is computed based on that one found in Blanchet and Faye [27], the relation between their 3 PN equations of motion and our result described in Section 8.5 below, and Equation (167). (After constructing E given as in Equation (183), we have checked that our 3 PN equations of motion make E to be conserved.)

We note that Equation (171) as well gives a correct geodesic equation in the test-particle limit, is Lorentz invariant, and admits the conserved energy. These facts can be seen by the form of \(a_1^i{\vert _{{\delta _{A\;\ln}}}}\) Equation (175); it is zero when m1 → 0, is Lorentz invariant up to 3 PN order, and is the effect of the mere redefinition of the dipole moments which does not break energy conservation.

Finally, we here mention one computational detail. We have retained during our calculation \({\mathcal R}\)-dependent terms with positive powers of \({\mathcal R}\) or logarithms of \({\mathcal R}\). As stated below Equation (76), it is a good computational check to show that our equations of motion do not depend on \({\mathcal R}\) physically. In fact, we found that the \({\mathcal R}\)-dependent terms cancel each other out in the final result. There is no need to employ a gauge transformation to remove such an \({\mathcal R}\) dependence. As for terms with negative powers of \({\mathcal R}\), we simply assume that those terms cancel out the \({\mathcal R}\) dependent terms from the far zone contribution. Indeed, Pati and Will [129, 130], whose method we have adopted to compute the far zone contribution, have proved that all the \({\mathcal R}\)-dependent terms cancel out between the far zone and the near zone contributions through all post-Newtonian orders.

8.5 Comparison

By comparing Equation (181) with the Blanchet and Faye 3 PN equations of motion [27], we find the following relation:

$${m_1}\vec a_1^{\,{\rm{this}}\;{\rm{work}}} = {m_1}{(\vec a_1^{{\rm{BF}}})_{\lambda = - {{1987} \over {3080}}}} + {m_1}{\vec a_1}{\vert _{{\delta _{A\ln}}}} + \,{m_1}{\vec a_1}{\vert _{{\delta _{A,{\rm{BF}}}}}},$$
(184)

where \({m_1}\vec a_1^{{\rm{this}}\;{\rm{work}}}\) is the 3 PN acceleration given in Equation (181), \({(\vec a_1^{{\rm{BF}}})_{\lambda = - 1987/3080}}\) is the Blanchet and Faye 3 PN acceleration with λ = −1987/3080, and \({m_1}{\vec a_1}{\vert _{{\delta _{A\;\ln}}}}\) is given in Equation (175) with ϵRA replaced by \(r_A^{\prime}\) for notational consistency with the Blanchet and Faye 3 PN equations of motion shown in [27]. \({m_1}{\vec a_1}{\vert _{{\delta _{A,{\rm{BF}}}}}}\) is an acceleration induced by the following dipole moments of the stars:

$$\delta _{A,{\rm{BF}}}^i = - {{3709} \over {1260}}m_A^3a_A^i.$$
(185)

We can compute \({m_1}a_1^i{\vert _{{\delta _{A,{\rm{BF}}}}}}\) by substituting \(\delta _{A,{\rm{BF}}}^i\) instead of \(\delta _{A\Theta}^i\) Equation (166). Thus, by choosing the dipole moments,

$$D_{A\Theta, {\rm{BF}}}^i = {\epsilon ^4}\delta _{A\Theta}^i - {\epsilon ^4}\delta _{A,{\rm{BF}}}^i,$$
(186)

we have the 3 PN equations of motion in completely the same form as \({(\vec a_1^{{\rm{BF}}})_{\lambda = - 1987/3080}}\). In other words, our 3 PN equations of motion physically agree with \({(\vec a_1^{{\rm{BF}}})_{\lambda = - 1987/3080}}\) modulo the definition of the dipole moments (or equivalently, the coordinate transformation under the harmonic coordinates condition). In [93], we have shown some arguments that support this conclusion.

The value of λ that we found, λ = −1987/3080, is perfectly consistent with the relation (1) and the result of [54] (ωstatic = 0).

Finally, let us discuss the ambiguity in the 3 PN equations of motion previously derived by Blanchet and Faye in [25, 27]. In their formalism, Dirac delta distributions are used to achieve the point particle limit. The (Lorentz invariant generalized) Hadamard partie finie regularization has been extensively employed to regularize divergences caused by their use of a singular source. In fact, unless regularization is employed, divergences occur both in the evaluation of the 3 PN field where (Poisson) integrals diverge at the location of the stars and in their derivation of the equations of motion where a substitution of the metric into a geodesic equation causes divergences.

When we regularize some terms at a point, say, \(z_A^i\), where the terms are singular, using the Hadamard partie finie regularization, roughly speaking we take an angular average of the finite part of the terms in the neighborhood of the singular point. Then if there are logarithmic terms such as \(\ln (\vert \vec x - {{\vec z}_A}\vert/{r_{12}})\), we should take an angular average over some sphere centered on \(z_A^i\) with a finite radius. The radius of the sphere is arbitrary but we do not ignore it because we should ensure the argument of the logarithms to be dimensionless.

The problem here is that there is a priori no reason to expect that the radius for each star introduced to regularize the field and another radius for that star introduced to regularize the geodesic equation coincide with each other. Thus the Blanchet and Faye 3 PN equations of motion have four arbitrary constants instead of two in our equations of motion. In our framework, we can see the origin of the number if we assume that we have defined a different body zone \(B_A^{\prime}\) in the derivation of the equations of motion from BA used in the derivation of the 3 PN field. However, in reality, we have only one body zone for each star. In our formalism the field is expressed in terms of the four-momentum (and multipole moments) which are defined as volume integrals over the body zone. On the other hand our general form of the equations of motion has been derived based on the conservation law of the four-momentum, and thus we evaluate the surface integrals in the general form of the equations of motion over the boundary of the body zone.

In fact, [25, 27] have shown that two of the four arbitrary constants can be removed by using a gauge freedom remaining in the harmonic gauge condition; the two places where the singular points exist are in some sense ambiguous. The remaining two turn out to appear as the ratios \({\zeta _A} = r_A^\prime/{s_A}(A = 1,2)\) where \(r_A^\prime\) and sA are the four regularization parameters (roughly speaking, the radii of \(B_A^\prime\) and BA in the terminology in the previous paragraph). Blanchet and Faye [25, 27] then proved that assuming the equations of motion are polynomials of the two masses of the stars, those two ratios should satisfy ln ζA = σ + λ(m1 + m2)mA where σ and λ are pure numbers. Then they showed that in order for their equations of motion to admit conserved energy, then σ = 159/308, while no argument was found to fix λ.

The above argument in turn means that the Blanchet and Faye 3 PN equations of motion do not give a conserved energy unless \(B_A^\prime\) is different from BA. Damour, Jaranowski, and Schäfer [54] pointed out that there is an unsatisfactory feature in the generalized partie finie regularization which contradicts with the mathematical structure of general relativity. Indeed, by using dimensional regularization which is pointed out by them to be more satisfactory in this regard, [54] derived an unambigous ADM Hamiltonian in the ADM transeverse traceless gauge. Later, Blanchet et al. [22] used dimensional regularization and found that their new equations of motion physically agree with ours and admit a conserved energy.

8.6 Summary

To deal with strongly self-gravitating objects such as neutron stars, we have used the surface integral approach with the strong field point particle limit. The surface integral approach is achieved by using the local conservation of the energy momentum, which led us to the general form of the equations of motion that are expressed entirely in terms of surface integrals. The use of the strong field point particle limit and the surface integral approach makes our 3 PN equations of motion applicable to inspiraling compact binaries which consist of strongly self-gravitating regular stars (modulo the scalings imposed on the initial hypersurface). Our 3 PN equations of motion depend only on the masses of the stars and are independent of their internal structure such as their density profiles or radii. Thus our result supports the strong equivalence principle up to 3 PN order.

At 3 PN order, it does not seem possible to derive the field in a closed form. This is because not all the super-potentials required are available, and thus we could not evaluate all the Poisson-type N/B integrals. Some of the integrands allow us to derive super-potentials in a series form in the neighborhood of the stars. For others, we have adopted an idea that Blanchet and Faye have used in [25, 26, 27]. The idea is that while abandoning the complete derivation of the 3 PN gravitational field valid throughout N/B, one exchanges the order of integralsFootnote 15. We first evaluate the surface integrals in the evolution equation for the energy of a star and the general form of equations of motion, and then we evaluate the remaining volume integrals. Using these methods, we first derived the 3 PN mass-energy relation and the momentum-velocity relation. The 3 PN mass-energy relation admits a natural interpretation. We then evaluated the surface integrals in the general form of equations of motion, and obtained the equations of motion up to 3 PN order of accuracy.

At 3 PN order, our equations of motion contain logarithms of the body zone radii RA. We showed that we could remove the logarithmic terms by a suitable redefinition of the representative points of the stars. Thus we could transform our 3 PN equations of motion into unambiguous equations which do not contain any arbitrarily introduced free parameters.

Our so-obtained 3 PN equations of motion agree physically (modulo a definition of the representative points of the stars) with the result derived by Blanchet and Faye [27] with λ = −1987/3080, which is consistent with Equation (1) and ωstatic = 0 reported by Damour, Jaranowski, and Schäfer [54]. This result indirectly supports the validity of the dimensional regularization in the ADM canonical approach in the ADMTT gauge.

Blanchet and Faye [25, 27] introduced four arbitrary parameters. In the Hadamard partie finie regularization, one has to introduce a sphere around each singular point (representing a point mass) whose radius is a free parameter. In their framework, regularizations are employed in the evaluation of both the gravitational field having two singular points and the two equations of motion. Since, in their formalism, there is a priori no reason to expect that the spheres introduced for the evaluation of the field and the equations of motion coincide, there arise four arbitrary parameters. This is in contrast to our formalism where each body zone introduced in the evaluation of the field is inevitably the same as the body zone with which we defined the energy and the three-momentum of each star for which we derived our equations of motion.

Actually, the redefinition of the representative points in our formalism corresponds to the gauge transformation in [27], and only two of the four parameters remain in [27]. Then they have used one of the remaining two free parameters to ensure the energy conservation, and there remains only one arbitrary parameter λ which they could not fix in their formalism.

On the other hand, our 3 PN equations of motion have no ambiguous parameter, admit conservation of an orbital energy of the binary system (when we neglect the 2.5 PN radiation reaction effect), and respect Lorentz invariance in the post-Newtonian perturbative sense. We emphasize that we do not need to a posteriori adjust some parameters to make our 3 PN equations of motion to satisfy the above three physical features.

We here note that Blanchet et al. [22], who computed the 3 PN equations of motion in the harmonic gauge using the dimensional regularization, have recently obtained the same value for λ.

The gauge condition in a harmonic gauge is related to the equations of motion. One may ask if the 3 PN equations of motion that have been derived so far guarantee the harmonic gauge condition through the corresponding post-Newtonian accuracy. This has not been tested yet. Let us call the n PN accurate metric components to be the components that are needed to compute the n PN equations of motion. Then the harmonic condition for the n PN field requires that matter obeys the n − 1 PN equations of motion. Thus, we need the 4 PN field to check if our resulting 3 PN equations of motion are a necessary condition to fullfil the harmonic gauge condition. This is beyond our current knowledge.

8.7 Going further

The 4 PN templates may be required to detect gravitational wave directly with high signal to noise ratio and extract astrophysical information from the wave. We have to derive the 4 PN equations of motion to derive the 4 PN templates if we use the energy balance.

At this moment, we should say that it is difficult to derive the 4 PN equations of motion. The technical obstacles regarding the derivation of the 4 PN equations of motion are the following. First, to derive the 4 PN equations of motion, we have to derive the 4 PN gravitational field at least in the neighborhood of a star. This requires the 3 PN gravitational field valid throughout the near zone. As seen in this article, however, it seems impossible to derive the 3 PN accurate gravitational field in harmonic coordinates in a closed form completely. We have not yet found the super-potentials necessary for us to evaluate the 3 PN Poisson integrals (or, retarded integrals).

Second, the amounts of calculations required would be too large to derive the 4 PN equations of motion successfully. For example, the Newtonian field consists of only two terms. The Landau-Lifshitz pseudotensor at the Newtonian order consists of basically one term. The gravitational field at 1 PN order consists of about 10 terms, while the Landau-Lifshitz pseudotensor has about 5 terms and thus we have to handle 50 terms to derive the 1 PN equations of motion. Next, the 2 PN field has about 102 terms. The Landau-Lifshitz pseudotensor in terms of hμν consists of about 50 terms and thus ∼ 103 terms must be treated to derive the 2 PN equations of motion. The number of terms in the 3 PN field are of order 103 ∼ 104, while Landau-Lifshitz pseudotensor has about 100 terms. Thus we may encounter ∼ 105 terms to derive the 3 PN equations of motion. Then at the 4 PN order, we may expect 104 ∼ 105 terms for the field, 500 terms for the integrand, and ∼ 106 ∼ 107 terms for the equations of motion. Furthermore, for each term spatial and temporal derivative generate additional terms. The number of newly generated terms are about three or four for each term. Thus the number of terms in the intermediate expressions to be dealt with are about ten-fold the number quoted above at each orderFootnote 16. Such an enormous number of terms forces us to use algebraic computing softwares to deal with them. In this work, we have extensively used the algebraic computing software Maple [118], Mathematica [166], and MathTensor [120] to deal with tensors. However, quantity changes quality. When some equations are produced through a sequence of black-boxes (that is, calculations done by computers) we could not check the equations term by term even if the calculations follow trivial procedures such as taking a Taylor expansion of equations around a star and extracting off some coefficients from the Taylor-expanded equations. Then, how could one confirm one’s result? Some possible tests include the following: Do the resulting equations of motion respect Lorentz invariance? Do the equations of motion admit the existence of a conserved energy? Does the gravitational field (if obtainable) satisfy the harmonic condition? Do the results obtained by more than two groups agree with each other? Then do all these consist of a set of necessary and sufficient criteria for confirmation?

One way to derive the 4 PN equations of motion is to use brute force (if the first difficulty — how to derive the required super-potentials at 3 PN order — was successfully dealt with). On the other hand, one may find a scaling appropriate to the late inspiralling phase and construct an approximation scheme based on such a scaling, as the post-Newtonian approximation is based on the Newtonian scaling. In the post-Newtonian approximation, the lowest order term is just the Newtonian term, and the first order term is the 1 PN correction. In the post-Minkowskian approximation, (the lowest order is just a straight line and) the first order correction is valid for any velocity but only for a weak field. Likewise, a new approximation scheme (if any) would give (the lowest order of the conservative dynamics and) the first order correction to the radiation reaction effect. Such an approximation may produce a smaller number of terms than the post-Newtonian approximation and thus give easier-to-treat equations of motion. Also the relatively lower order equations of motion obtained by such an approximation are expected to give templates with the same accuracy as the accuracy achieved by the relatively higher order post-Newtonian equations of motion. Thus, such an approximation scheme is an attractive alternative to the post-Newtonian approximation. The construction of such an approximation remains to be done in future work.