The Post-Newtonian Approximation for Relativistic Compact Binaries

We discuss various aspects of the post-Newtonian approximation in general relativity. After presenting the foundation based on the Newtonian limit, we show a method to derive post-Newtonian equations of motion for relativistic compact binaries based on a surface integral approach and the strong field point particle limit. As an application we derive third post-Newtonian equations of motion for relativistic compact binaries which respect the Lorentz invariance in the post-Newtonian perturbative sense, admit a conserved energy, and are free from any ambiguity.


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 v/c ∼ and the deviation from the flat metric would be the Newtonian potential, say Φ. Then these are related by Φ/c 2 ∼ v 2 /c 2 ∼ 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 Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 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.

Post-Newtonian equations of motion
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 Poincaré 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 regularization 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] Therefore, it is possible to fix the λ parameter from the result of [54] as 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.

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.

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].

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 for an unknown function g, one would like to have a one-parameter (or possibly multi-parameter) family of solutions, 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 v i , i = 1, 2, 3): supplemented by an equation of state. For simplicity we have considered a perfect fluid. It can be seen that the variables {ρ(x i , t), P (x i , t), Φ(x i , t), v i (x j , t)} obeying the above equations satisfy the following scaling law: 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 v 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 (Q ij , P ij ) 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: where the functions a, b, and c i 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 R 4 being the spacetime (coordinates t, x i ). 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, and define the integral curve as the curve on which τ and x i stay constant 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.

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: where L V is the Lie derivative with respect to the tangent vector of the curves defined above, and the remainder term R µν n+1 is 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 (τ, x k ). 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 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] 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.

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 where t µν LL is the Landau-Lifshitz pseudotensor [115]. In this section we shall choose an isentropic perfect fluid for T αβ which is enough for most applications, where ρ is the rest mass density, Π the internal energy, P the pressure, and u µ the four-velocity of the fluid with normalization g αβ u α u β = −1.
The conservation of energy and momentum is expressed as Defining the gravitational field variable as where η µν is the Minkowski metric, the reduced Einstein equations (14) and the gauge condition (15) take the following form: Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 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]: where Equations (23) and (24) together imply the conservation law We shall take as our variables the set {ρ, P, v i , h αβ }, with the definition The time component of four-velocity u 0 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 : 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 v 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, where ∆ is the Laplacian in the flat space. We shall choose h ij and h ij ,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 h ij and h ij ,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, where r = |y i − x i | and C( , τ, x i ) is the past flat light cone of the event (τ, x i ) in the spacetime given by , truncated where it intersects with the initial hypersurface τ = 0. h µν H is the unique solution of the homogeneous wave equation in the flat spacetime, Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 h µν H evolves from a given initial data on the τ = 0 initial hypersurface which are subject to the constraint equations (30). The explicit form of h µν H is available via the Poisson formula (see e.g. [144]), 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 ( , τ, x i ). Accordingly one should convert all the time indices to τ indices. For example, T τ τ = 2 T tt which is of order 4 , since T tt ∼ ρ is of order 2 . Similarly the other components of stress-energy tensor T τ i = T ti and T ij 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 where we have adopted the notation and 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 √ −gg µν and η µν . The conservation law (27) also has its first nonvanishing derivatives at this order, which are 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.

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].

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 separation 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 gravity 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 = O( 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.

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 i 1 on star 1 becomes The integral region B 1 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, Note that the sphere ∂B 1 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 Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 x BA α BA Figure 1: Body zone coordinates and near zone coordinates. In the near zone coordinates (τ, x i ), 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 (τ, α  z i A (τ ) is a representative point of the star A, e.g. the center of the mass of the star A. R A , 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. dR A /dτ = 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).
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 ∂B A , 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 ∂B A is far away from the surface of star A, we can evaluate the gravitational energy momentum flux over ∂B A with the post-Newtonian gravitational field. In fact we shall express our equations of motion in terms of integrals over ∂B A and be able to evaluate them explicitly.

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 Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 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, x i ) 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 η = −2 t. 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 O( ). In other words, we assume that the star rotates slowly and is pressure supported 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 in the body zone coordinates: Here the underlined indices mean that for any tensor A i , A i = −2 A i . In [94] , we have transformed T µν N , the components of the stress-energy tensor of the matter in the near zone coordinates, to T µν A 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 (or rather Λ µν N , the source term of the relaxed Einstein equations; see Equation (63)), 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.

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: We define the mass, the dipole moment, the quadrupole moment, and the momentum of the star A as Here z i A 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, where v i A = dz i A /dτ . Using the velocity momentum relation, we calculate the time derivative of the momentum, where F i 1 is defined by Equation (43). The Newtonian potential can be expressed by the mass and multipole moment as Substituting Φ into Equation (56), we obtain the equations of motion, 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, Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 . . denotes the symmetric-tracefree operation on the indices between the brackets, and I Mp A are the Newtonian mass multipole moments of order 2p.

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].

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 (τ, x i )). Thus we define a deviation field h µν as 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, where the comma denotes the partial derivative. In the harmonic gauge, we can recast the Einstein equations into the relaxed form, where = η µν ∂ µ ∂ ν is the flat d'Alembertian and Here, T µν and t µν LL denote the stress-energy tensor of the stars and the Landau-Lifshitz pseudotensor [115]. The explicit form of t µν LL in the harmonic gauge is χ µναβ 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 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, where C(τ, x k ) means the past light cone emanating from the event (τ, 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 R/ , where 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 P N , we write the field as h µν P N (N ) is the near zone integral contribution to the near zone field, and h µν P N (F ) 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 ) 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 µν H (τ, x i ) = 0 corresponds to a no-incoming radiation condition at (Minkowskian) past null infinity (see e.g. [77]).
Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 Equation (69) can be written down as Kirchhoff's formula, Then the no-incoming radiation condition at (Minkowskian) past null-infinity, is a sufficient condition to h µν H (τ, 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, Henceforth we shall omit the subscript P N (N ) of the field h µν for notational simplicity.

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, Note that the above integral depends on the arbitrary length R in general. The cancellation between the 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 R (R −k ; k > 0). In other words, we simply let R → 0 whenever it gives a convergent result. On the other hand, we shall retain terms having positive powers of R (R k ; k > 0), and logarithmic terms (ln R) to confirm that the final result, in the end of calculation, is independent of 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 B A , and from elsewhere, N/B. Schematically we evaluate the following two types of integrals (we omit indices of the field), where r A ≡ x − z A . We shall deal with these two contributions successively.

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 , gives Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 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 ≡ | r A |.
To derive the 3 PN equations of motion, h τ τ up to O( 10 ) and h τ i as well as h ij up to O( 8 ) are required.
In the above equations we define the multipole moments as where we introduced a collective multi-index We simply call P µ A the four-momentum of the star A, P i A the momentum, and P τ A the energy 7 . Also we call D k A the dipole moment of the star A, and I kl A the quadrupole moment of the star A.
Then we transform these moments into more convenient forms. By the conservation law (67), we have where v i A ≡ż i A , an overdot denotes a time derivative with respect to τ , and y A ≡ y − z A . Using these equations and noticing that the body zone remains unchanged (in the near zone coordinates), i.e.Ṙ A = 0, we have where 7 It is just for simplicity to call P µ A a four-momemtum. P µ A is not defined (general, nor even Lorentz) covariantly nor as a volume integral of tensor density of weight −1. The integrand Λ µν is a tensor density of weight −2. In Section 6.1, we will find that P µ A equals √ −gm A u µ A in a certain sense. P τ A is an energy of the star (modulo the contribution from the integrand χ µναβ ,αβ ) if the star considered were isolated, the space-time is asymptotically flat, and finally, if we integrate over a spacelike hypersurface, rather than integrate only over the body zone of the star [115].
Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 and The operators [ ] and ( ) attached to the indices denote anti-symmetrization and symmetrization. M ij A 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 and where l is the number of indices in the multi-index K l . Now, from the above equations, especially Equation (88), we find that the body zone contributions, h µν B n=0 , are of order O( 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 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 , we have to evaluate the surface integrals appear formally at the order 2l+4 and 2l+2 . Thus where we omitted irrelevant terms and numerical coefficients. Thus one may expect that Q K l i A and R K l ij A 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 K l ij A contribute to the 3 PN equations of motion. The important thing here is that 4 . Finally, since the order of h µν Bn (n ≥ 1) is higher than that of h µν B n=0 , we conclude that h µν B = O( 4 ).

N/B contribution
As far as the N/B contribution is concerned, since the integrand Λ µν N = −gt µν LL + χ µναβ ,αβ 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( x) which satisfy ∆g( x) = f ( x), ∆ denoting the Laplacian, we have for each integral (see e.g. the n = 0 term in Equation (77)) 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πg( x)" terms 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 superpotentials 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 , and since h µν where we expanded h µν B in an series, Similarly, in the following we expand h µν in an series. From these equations we find that the deviation field in N/B, h µν , is O( 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 h µν in the body zone, we have to know the internal structure of the star.)

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 B A 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 where I ij A,NZC ≡ I ij A . I ij A,GFC are the quadrupole moments defined in the generalized Fermi normal coordinates. Note that the difference δI ij A 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 δI ij A above is a necessary term which makes the so-obtained 3 PN metric equal to the Schwarzschild metric boosted at the constant velocity 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)) The radius of a compact body A is of order of its mass m A , and the Lorentz contraction deforms the body so that the radius of the body changes by the amount 2 m A v 2 A + O( 4 ). So the apparent quadrupole moments due to Lorentz contraction is of order of 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.

General form of the equations of motion
From the definition of the four-momentum, and the conservation law (67), we have the evolution equations for the four-momentum: 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.
Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 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, 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 i A freely to determine the representative point z i A (τ ) of star A. Up to 2.5 PN order we take D i A = 0 and simply call z i A 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 i A (see [94] and Appendix B.1). At 3 PN order, yet another choice of the value of the dipole moment D i A shall be examined (see Sections 7 and 8.2). In Equation (111), P τ A rather than the mass of star A appears. Hence we have to derive a relation between the mass and P τ A . 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 i A to be specified, we can derive the equations of motion for a strongly self-gravitating star using the post-Newtonian approximation.

On the arbitrary constant R A
Since we have introduced the body zone by hand, the arbitrary body zone size R A seems to appear in the metric, the multipole moments of the stars, and the equations of motion. More specifically R A 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.

R A dependence of the field
B and N/B contributions to the field depend on the body zone boundary R A . But h µν itself does not depend on R A . Thus it is natural to expect that there are renormalized multipole moments which are independent of R A since we use nonsingular matter sources. This renormalization would absorb the R A 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( R A ) 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 i A 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 " for the renormalized P µ A .

R A 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, R A . Actually this is not the case.
Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 In the derivation of Equation (111), if we did not use the conservation law (67) until the final step, we have Now the conservation law is satisfied for whatever value we take for R A , then the right hand side of the above equation is zero for any R A . Hence the equations of motion (111) do not depend on R A (a similar argument can be found in [73]). Along the same line, the momentum-velocity relation (86) does not depend on R A . In Section 4.8 we shall explicitly show the irrelevance of the field and the equations of motion to R A by checking the cancellation among the R A dependent terms up to 0.5 PN order.

Newtonian equations of motion
Then we define the mass of star A as the integrating constant of this equation, here m A 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 m A is constant. Then we obtain the lowest order h τ τ : Second, from Equation (91) with Equations (102) and (103) we obtain Q i A = 0 at the lowest order. Thus we have the Newtonian momentum-velocity relation . 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: where y i A is the integral variable and y A = y − z A . We defined n 1 as the spatial unit vector 9 emanating from z 1 , r 12 ≡ z 1 − z 2 , and n i 12 ≡ r i 12 /r 12 . We used r 1 = R 1 n 1 and r 2 = r 12 + R 1 n 1 . For details see Figure 3).
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.

First post-Newtonian equations of motion
Next, at 1 PN order, we need 6 h τ τ and 4 h µν . 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 m A , i.e. 5 h τ τ = 0. Now we obtain 4 h τ i from the Newtonian momentum-velocity relation: We evaluate the surface integrals in the evolution equation for P τ A at 1 PN order. The result for star 1 is dP τ where we used the Newtonian equations of motion. From this equation we have the mass-energy relation at 1 PN order, Then we have to calculate the 1 PN order  (91) and (92). We then find that they depend on R A , hence we ignore them and obtain To derive 6 h τ τ and 4 h ij , we have to evaluate non-compact support integrals for 6 h τ τ N/B and 4 h ij N/B , and the n = 2 term in Equation (76) for h τ τ : 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: (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, Then it is relatively easy to find the super-potentials for these scalars. The results are ∆ ln r 1 = 1/r 2 1 and ∆ ln S = 1/(r 1 r 2 ) where S = r 1 +r 2 +r 12 [77]. Equation (101) for each integrand then becomes The second last term of Equation (128) and the terms abbreviated as O( R A ) in the above two equations arise from the surface integrals in Equation (101) (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 where n A ≡ r A /r A . Evaluating the surface integrals in Equation (111) as in the Newtonian case, we obtain the 1 PN equations of motion, where we defined the relative velocity as V ≡ v 1 − 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 . Then we solve it functionally and obtain the mass-energy relation at n + 1 PN order. Next we calculate Q i A at n + 1 PN order and derive the momentum-velocity relation at n + 1 PN order. Then we calculate Q K l i A and R K l ij A . With the n + 1 PN mass-energy relation, the n + 1 PN momentum-velocity relation, Q K l i A , and R K l ij A , 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 2 dv i 1 /dτ 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 − 2 m 2 r i 12 /r 3 12 . Basically we shall derive the 3 PN equations of motion with the procedure as described above.

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 R A in the 0.5 PN field. Returning back to the derivation of the 1 PN h τ τ with care for the R A dependence (see Equation (128)), we get Now we split P τ A as P τ A =P τ A +P τ A whereP τ A does not depend on R A whileP τ A does. In order to evaluate the R A dependent part of P τ A ,P τ A , we use the fact that the integral of Λ τ τ N over the near zone does not depend on the size of the body zone. (Notice that P τ A is defined as the volume integral of Λ τ τ N over B A .) In fact, from the relevant expression of the pseudotensor, we obtain at the lowest order Hence we findP τ A = 2 7m 2 A /(2 R A ) + O( 2 ). The above equation shows us that the 0.5 PN field is independent of R A and fully expressed by the R A independent energyP τ A , or mass up to 0.5 PN order.
In the similar manner we split P i A as P i , we find that the 0.5 PN momentum-velocity relation does not depend on R A : . 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 R A 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.

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.
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/r 1 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πg( x)" term and the surface integral in Equation (101).

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 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, 10 A super-potential here is a particular solution of a Poisson equation whose source term is non-compact.
Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 Note that the integration region is N/B and therefore g( x) is nonsingular in N/B. For this kind of source term, we have given a method in [91] to find a field F (m,n) [A,c] in the neighborhood of star A in the following sense: 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.

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 for star 1. Here the operator disc R A means to discard all the R A dependent terms other than logarithms of R A [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 R A in front of the Poisson integral as in Equation (140) above. We first exchange the order of integrals, where we defined a sphere B 1 whose center is z 1 and that has a radius r 1 which is a constant slightly larger than R 1 for any (small) ( R 1 < r 1 r 12 ). The reason we introduced r 1 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 R A , however we reasonably discard such R A -dependent terms (other than logarithmic dependence) as explained in Section 4.5. Using the so-obtained R A -independent field, we evaluate the surface integrals in the general form of the 3 PN equations of motion by discarding the R A dependence emerging from the surface integrals, and obtain the equations of motion. Thus the "discarding-R A " 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 was introduced to distinguish the two kinds of R A dependence and to discard the R A dependence in the right order. We show here a simple example. Let us consider the following integral: Using ∆ ln r 1 = 1/r 2 1 , we can integrate the Poisson integral and obtain the "field", Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 Since the "body zone contribution" must have an R 1 dependence hidden in the "moments" as 4π R 1 /r 1 [+O(( 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", On the other hand, we can derive the "equations of motion" by first evaluating the surface integral over ∂B 1 , Thus, if we take ∂B 1 as the integral region instead of ∂B 1 in the first equality in Equation (145), or if we take lim r 1 → R1 without employing disc R A beforehand, we will obtain an incorrect result, 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, 8 h τ i N/B n=0 and 10 h τ τ N/B n=0 + 8 h l N/B n=0 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 directintegration method, the super-potential method, and the super-potential-in-series method.

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 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 as dP τ (−4( n 12 · v 2 ) + 6( n 12 · v 1 )) + m 2 r 12 (−10( n 12 · v 1 ) + 11( n 12 · v 2 )) Remarkably, we can integrate Equation (147) functionally: Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 Equation (148) together with Equations (149,150,151) gives the 3 PN order mass-energy relation.

Meaning of P τ AΘ
In this section we explain the meaning of P τ AΘ . First of all, we expand in an series the four-velocity of star A normalized as g µν u µ The field should be evaluated somehow at z A . This is a formal series since the metric derived via the point particle description diverges at 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 ( x) which can be expanded around z 1 in the form Then the Hadamard partie finie at z 1 of the function f ( x) is defined by Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 For example, by this procedure h τ τ becomes (see Equation (130))  (149,150,151), we find at least up to 3 PN order: 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Θ 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 distributions 11 .
Finally, we note that up to 3 PN order is satisfied if we use the Hadamard partie finie regularization explained above. 11 In fact, we need the 2.5 PN field and a part of the 3 PN field ( 8 h τ [i,j] ) to evaluate the 3 PN mass-energy relation. This does not mean that the Hadamard partie finie regularization is sufficient to derive the 3 PN equations of motion. Indeed, when a Dirac delta distribution is adopted to achieve the point particle limit instead of the strong field point particle limit, one needs another regularization scheme such as the dimensional regularization to derive the 3 PN equations of motion in an unambiguous form [22].

Third Post-Newtonian Momentum-Velocity Relation
We now derive the 3 PN momentum-velocity relation by calculating the Q i A integral at 3 PN order. From the definition of the Q i A integral, Equation (91), we find that the calculation required is almost the same as that in the equation for dP τ A /dτ . Namely, it turns out that we do not need to use the direct-integration method to compute the Q i A integral. Therefore it is straightforward to evaluate the surface integrals in the definition of Q i A . Here we split Q i A into Q i AΘ and Q i Aχ as with and we show only Q i AΘ : where ≤n f is the quantity f up to order n inclusively. Here it should be understood that a i A in the last expression is evaluated with the Newtonian acceleration.
Q i A of O( 6 ) appears at the 4 PN or higher order field. Thus up to 3 PN order, 6 Q i A affects the equations of motion only through the 3 PN momentum-velocity relation. For this reason, not Q i Aχ but Q i AΘ is necessary to derive the 3 PN equations of motion. The explicit expression for Q i Aχ is given in [91]. Now with Q i AΘ in hand, we obtain the momentum-velocity relation. It turns out that the χ part of the momentum velocity relation is a trivial identity 12 . Thus, defining the Θ parts of P µ A and D i A in the same way as for Q i A , we obtain As explained in the previous Sections 3.4 and 4.4, we define the representative point z i A of star A by choosing the value of D i A . In other words, one can freely choose D i A in principle 13 . One may set D i A 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, 12 Defining the Θ parts of P µ A and D i A in the same way as for Q i A , it turns out that P i Aχ /dτ 2 is a trivial identity [91]. Note that P µ Aχ and D i Aχ can be converted into surface integral forms and thus can be evaluated explicitly. 13 It is natural to expect that z A is in star A unless the star is, say, crescent-shaped. One may not choose the representative point of the sun to be in the earth.
Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 Henceforth, we shall define z i A by this equation. Finally, it is important to realize that the nonzero dipole moment D i A 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 δ i AΘ appears only at 10 h τ τ as Then the corresponding acceleration becomes The last term compensates the Q i A 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, 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 8 h τ τ , the field ≤6 h µν , 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 K l i A and R K l ij A integrals appearing in the 3 PN field in h µν B contribute to the 3 PN equations of motion. 6 Q i AΘ 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 i Aχ appear independently of D i AΘ in the field, and hence 4 D i Aχ affects the 3 PN equations of motion. In summary, the ≤4 Q K l i A , ≤4 R K l ij A , 6 Q i AΘ , δ i AΘ , and 4 D i Aχ contributions to the 3 PN field can be written as where ". . . " are other contributions. On the other hand, 6 Q i AΘ and δ i AΘ affect the equations of motion through the momentum-velocity relation in Equation (111), but they cancel each other out, since we choose Equation (164). Then these contributions to a 3 PN acceleration can be summarized into the contribution to m 1 a i 1 from Collecting these contributions mentioned above, we obtain the 3 PN equations of motion. However, we found that logarithmic terms having the arbitrary constants R A in their arguments survive, where the acceleration through 2.5 PN order, (dv i 1 /dτ ) ≤2.5PN , 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 R A , 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 R A from our equations of motion. There, we show the explicit form of the 3 PN equations of motion we obtained.

Arbitrary constant R A
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 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 3 1 ( r 1 · a 1 )/r 5 1 as an integrand. Then we find Actually, we could remove the R A 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 R A dependence in Equation (171): 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)  Comparing the above equations with Equation (171), we easily conclude that the representative point z i A of star A defined by 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 i A is defined by the value of D i A , and in turn we have a freedom to assign to D i A any value as we like (though it may be natural to set the value of D i A such that z i A resides inside star A). We also note that we define z i A 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 radii 14 . It is well-known that choosing different values of dipole moments corresponds to a coordinate transformation.

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: 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 , the momentum-velocity relation, and the equations of motion (and relations among higher multipole moments, hidden in ". . . "). Thus if the logarithmic dependence of R A arises from the second term of Equation (178) There is yet another fact which supports our interpretation. Let us retain D i A = 0 for a while. We then find that the near zone dipole moment D i N defined by a volume integral of Λ τ τ N y i becomes Then if we take the old choice of D i A , the volume integral becomes where terms denoted by ". . . " have no logarithmic dependence. Notice that the near zone dipole moment can be freely determined, say, D i N = 0, because we can define the origin of the near zone freely. By taking temporal derivatives twice of D i N , we see that D i A,new gives a natural definition of the center of the mass in terms of which the 3 PN equations of motion are independent of R A .

Third post-Newtonian equations of motion
By adding m 1 a i 1 | δ 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), Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 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 i 1 | δ A ln , Equation (175); it is zero when m 1 → 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 Rdependent terms with positive powers of R or logarithms of R. As stated below Equation (76), it is a good computational check to show that our equations of motion do not depend on R physically. In fact, we found that the R-dependent terms cancel each other out in the final result. There is no need to employ a gauge transformation to remove such an R dependence. As for terms with negative powers of R, we simply assume that those terms cancel out the 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 R-dependent terms cancel out between the far zone and the near zone contributions through all post-Newtonian orders.

Comparison
By comparing Equation (181) with the Blanchet and Faye 3 PN equations of motion [27], we find the following relation: where m 1 a this work 1 is the 3 PN acceleration given in Equation (181), ( a BF 1 ) λ=−1987/3080 is the Blanchet and Faye 3 PN acceleration with λ = −1987/3080, and m 1 a 1 | δ A ln is given in Equation (175) with R A replaced by r A for notational consistency with the Blanchet and Faye 3 PN Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 equations of motion shown in [27]. m 1 a 1 | δ A,BF is an acceleration induced by the following dipole moments of the stars: We can compute m 1 a i 1 | δ A,BF by substituting δ i A,BF instead of δ i AΘ into Equation (166). Thus, by choosing the dipole moments, we have the 3 PN equations of motion in completely the same form as ( a BF 1 ) λ=−1987/3080 . In other words, our 3 PN equations of motion physically agree with ( a BF 1 ) λ=−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.
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 i A , 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(| x − z A |/r 12 ), we should take an angular average over some sphere centered on z i A 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 in the derivation of the equations of motion from B A 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 ζ A = r A /s A (A = 1, 2) where r A and s A are the four regularization parameters (roughly speaking, the radii of B A and B A 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 = σ + λ(m 1 + m 2 )/m A 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 is different from B A . Damour, Jaranowski, and Schäfer [54] pointed out that there is an unsatisfactory feature in the generalized partie finie regularization Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 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.

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 Poissontype 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 integrals 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 R A . 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.

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 10 2 terms. The Landau-Lifshitz pseudotensor in terms of h µν consists of about 50 terms and thus ∼ 10 3 terms must be treated to derive the 2 PN equations of motion. The number of terms in the 3 PN field are of order 10 3 ∼ 10 4 , while Landau-Lifshitz pseudotensor has about 100 terms. Thus we may encounter ∼ 10 5 terms to derive the 3 PN equations of motion. Then at the 4 PN order, we may expect 10 4 ∼ 10 5 terms for the field, 500 terms for the integrand, and ∼ 10 6 ∼ 10 7 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 order 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.

Acknowledgements
The authors are greatly indebted to Bernard F. Schutz for his helpful comments and valuable discussions while they visited the Max-Planck-Institut für Gravitationsphysik (Albert-Einstein-Institut), Germany. When the authors wrote this article, Y.I. stayed at Tokyo University, for which he particularly thanks Masaru Shibata for his hospitality. The authors would like to thank Takashi Fukumoto for fruitful discussion.

A Far Zone Contribution
A formal retardation expansion gives a divergent integral when we evaluate the gravitational field. Schematically the field h(τ, x i ) at the field point (τ, x i ) is given by (we omit the indices of the field for simplicity) This integral is divergent if we set the upper bound of the integral to infinity. In our formalism we take the upper bound of the integral as ∼ R/ and keep the R dependence in the field. In the last step of the derivation, we let R go to infinity if and only if this procedure gives finite result at least up to the relevant post-Newtonian order. At a high order of the retardation expansion approximation, however, some terms must have a ln R, and/or a R n (n > 0) dependence. In this case, we can not let R go infinite.
To solve this problem, it is important to realize that the field at (τ, x i ) consists of two contributions: the retarded integral inside the near zone (the near zone contribution) and the retarded integral outside the near zone (the far zone contribution).
Since R is introduced artificially, we expect that the far zone contribution to the near zone field must have some R dependent terms which cancel out completely the R dependent terms in the near zone. This expectation is the case and was proved to any post-Newtonian order by Pati and Will [129] within their formalism. They have shown that the total field (which is the sum of the near zone contribution plus the far zone contribution) is finite and independent of R.
In this section we calculate the far zone contribution to the 3 PN equations of motion to make this article self-contained. We entirely follow [129,162,165], and the result is the same as theirs: The far zone field does not have any influence on the equations of motion up to 3 PN order inclusively.
For the near zone field, we write the field as where h µν N is the field in the near zone, h µν N (N ) is the near zone integral contribution to the near zone field, and h µν N (F ) is the far zone integral contribution to the near zone field. Our task here is to evaluate h µν N (F ) to the 3 PN order. In turn, this means that we should derive h µν F to lower order because the integrand Λ µν in the retarded integral for h µν N (F ) consists of h µν F . Only in this section, we do not use our bookkeeping parameter while it should be understood that v is of order , and m of order 2 . Now we evaluate the near zone contribution to the far zone field using a multipole expansion: Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 where K l is a collective multi-index, r ≡ | x| is the distance from the near zone center to the field point, and u = t − r. The near zone multipole moments M K l µν are defined as Next, the far zone contribution to the far zone field point can be evaluated as follows. The key idea would be an introduction of a new time u and may be seen from the following transformation of the retarded integral where the radial integral is transformed into a temporal integral from the past infinity to u = t − r: where We then make STF decomposition of the integrand Λ µν in the retarded integral as Λ µν ∼ f B,L r −B n L , where n = r/r. The integrand in Equation (193) becomes a summation of terms, each of which consists of (some algebraic combinations of) near zone multipole moments (defined by Equation (192)) multiplied by terms which explicitly depend on t, u , and Ω . Roughly speaking, the idea is that we do integral by parts many times, each time increasing the number of the time derivatives of the multipole moments, and assume that the system is sufficiently stationary in the past so that the contributions from the past infinity disappear. We then have for the far zone contribution to the far zone field point: with z = R/r 17 , and for B = 2 and for B = 2. Other quantities in the above equations are given by 17 Notice that if we recover the expansion parameter , z < for the field point in the far zone.
Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 Here P l and Q l are the Legendre function of the first and the second kind. α represents an angular defect due to the fact that the far zone integral does not cover the whole spacetime due to the near zone. The function k(t) and the retarded time multi-derivatives of the STF coefficient f B,L comes from the recursive integrals by parts. Then combining the near zone and the far zone contributions to the far zone field point, we have to the required order: with for B = 2 and for B = 2.
Evaluating the coefficients E q B,L , we finally have Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 where we reintroduce our bookkeeping parameter , and Z ij N is written with the near zone quadrupole moment as Z ij N = (1/2)d 2 I ij N /dτ 2 . We note that the fields h τ τ of O( 10 ) and h µi of O( 8 ) are the 3 PN fields in our formalism. Finally assuming that in the distant past the binary was sufficiently stationary, we find that the far zone contribution becomes a function of time only at the 3 PN order. It turns out that only the spatial derivative of those 3 PN fields contributes to the 3 PN equations of motion. Thus the far zone contribution does not affect the equations of motion up to 3 PN order inclusively.
In this section, we have given a highly rough sketch about the method developed by Pati and Will [129] which is based on their previous work [165]. Readers may consult [129,130,162,163,165] for more details.

B Effects of Extendedness of Stars
In this section, we derive the spin-orbit coupling force, the quadrupole-orbit coupling force, the spin-spin coupling force, and the spin geodesic precession equation to the lowest order.
Our order-counting of the multipole couplings may need an explanation. The magnitude of a mass multipole moment and a current multipole moment of order l of a star are roughly (mass) × (radius of the star) l and (mass)×(radius of the star) l ×(velocity of steller internal motion). Since we assume slow stellar rotation where (velocity of steller internal motion) = O( ) and the strong field point particle limit where (radius of the star) = O( 2 ), the mass multipole moment and the current multipole moment are of order O( 2l+2 ) and O( 2l+3 ), respectively. For example, the spinorbit coupling force which takes a form of (mass) × (orbital velocity) × (spin) appears at O( 4 ), that is, 2 PN order, not at the usual 1.5 PN order where rapid stellar rotation is assumed.
In summary the structure of the equations of orbital motion is written schematically as Here F i Newton , F i 1 PN , F i 2 PN , F i RR , and F i 3 PN , respectively, are the Newtonian force, the 1 PN force, the 2 PN force, the 2.5 PN radiation reaction force, and the 3 PN force. F i SO and F i 1 PNSO are the spin-orbit coupling force and its 1 PN correction, while F i QO and F i 1 PN QO are the quadrupoleorbit coupling force and its 1 PN correction. F i OO , F i SS , and F i TO are the octupole-orbit coupling force, spin-spin coupling force, and tidal-orbit coupling force, respectively 18 .
Though formally the effects of the 1 PN spin-orbit coupling, the 1 PN quadrupole-orbit, the octupole-orbit coupling, and the tidal-orbit coupling appear up to 3 PN order in our ordering, we focus our attention onto the lowest order spin-orbit coupling, the spin-spin coupling, and the quadrupole-orbit coupling forces. The 1 PN spin-orbit coupling force was derived by Tagoshi, Ohashi, and Owen [148].
Before investigating the multipole-orbit coupling forces, it is worth noticing that our definition of multipole moments is operational and the relation between these multipole moments and the intrinsic multipole moments of a star has not been given. Let us discuss briefly about this problem.
First, for an appropriate frame for the definition of multipole moments it is natural to define the multipole moments in a frame attached to the star, and nonrotating with respect to an asymptotic inertial frame (see [37] for the case of an earth-satellite system in the solar system). If we do not define the multipole moments in an appropriate frame, for example, an apparent quadrupole moment would be produced by Lorentz contraction caused by the orbital motion of the star and an apparent spin would be produced by the Thomas precession. One realization of an appropriate frame are the (generalized) Fermi normal coordinates [13]. To derive the spin-orbit coupling force in the same form as in previous works [46,108,152], it is sufficient to assume the coordinate transformation of Λ τ τ N in the near zone coordinates to the Λ µ ν A in the (generalized) Fermi normal coordinates in the following implicit form ( [13,37,57,58,59,60] (see also Appendix B.4): with Γ τ A τ = 1 + O( 2 ) and Γ τ A i = v i A + O( 2 ) (think of a Lorentz transformation). The 2 in front of the second term arises from the body zone coordinates rescaling (x i A = 2 α i A ). An explicit expression of Γ µ A ν is not required for our purpose. To 2 PN order, which is the sufficient order to derive the lowest order spin-orbit, spin-spin, and quadrupole-orbit coupling forces, the transformation changes only 8 h τ τ : where ". . . " denotes irrelevant terms and where M ij A is given by Equation (90). Second, we discuss the χ part of our multipole moments. We have used Λ µν N = Θ µν N +χ µναβ N ,αβ in the definition of our multipole moments since we could not evaluate χ parts of multipole moments separately except for some low order moments. However, we have to take into account carefully the fact that χ parts of our higher multipole moments can affect the equations of motion for two point masses. An obvious example can be found from the definition of the energy and the mass. It is natural to define the mass as a volume integral of Θ τ τ N . In fact the χ part of the energy, P τ Aχ , appears at 2 PN order. As another example, the χ part of our dipole moment can be evaluated directly as Thus if we define the center of mass of the star not by D i A = 0 but by D i AΘ = 0, the form of the equations of motion for the two point masses would change 19 .
Finally, we list the relevant field and Q i up to the required order to derive the lowest order spin-orbit, the spin-spin, and the quadrupole-orbit coupling forces and the spin geodesic precession equation.

B.1 Spin-orbit coupling force
It is well known that the definition of a dipole moment of the star, which we equate to zero to determine the center of mass of the star, affects the appearance of the spin-orbit coupling force 19 We could evaluate the χ part of the spin by virtue of its antisymmetric nature, and the result is M ij Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 (see e.g. [108]). If one chooses d i A = 0 as in [94], the corresponding spin-orbit coupling force takes the usual form [46,108,152] 6( s 2 × n 12 ) · V n i 12 + 4 s 2 × V − 6 s 2 × n 12 ( n 12 · V ) + 4 m 2 r 3 12 6( s 1 × n 12 ) · V n i 12 + 3 s 1 × V − 3 s 1 × n 12 ( n 12 · V ) , where ∆ ij ≡ δ ij − 3n i 12 n j 12 and × in this section denotes the outer product for the usual Euclidean spatial three-vectors. The spin vector s i A is defined by where ijk is the totally antisymmetric symbol.

B.2 Spin-spin coupling force
It is straightforward to derive the spin-spin coupling force. The definition of the center of mass does not change the form of the spin-spin coupling force as expected.

B.3 Quadrupole-orbit coupling force
The quadrupole-orbit coupling force can be derived by evaluating the following integral with 8 [(−g)t ij LL ] = (δ i k δ j l + δ j k δ i l + −δ ij δ kl ) 4 h τ τ,k 8 h τ τ,l /4 + . . . . Then we have the quadrupole-orbit coupling force in the same form as the result in [152], This result agrees formally with the Newtonian quadrupole-orbit coupling force, see Equation (58).
Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 Now defining the intrinsic center of mass by setting D µ A = Mîτ A = 0, we obtain Notice that this relation provides the spin-orbit coupling force in the previous form (see Appendix B.1; thus we get d i A = D i A in the present treatment). Using Equation (228), we finally obtain the spin geodesic equation in the usual form (we omit the quadrupole and Z ijkl 1 term for simplicity), where S 1 is the spatial part of the intrinsic spin four-vector S 1î = ijk Mĵk 1 /2. Equation (239) is the geodesic precession equation, or called the de Sitter-Fokker precession (see e.g. Section 40.7 in [122] and [36,145]).

B.5 Remarks
We have derived the spin-orbit coupling force, the spin-spin coupling force, and the quadrupoleorbit coupling force to the lowest order. The spin geodesic equation could also be derived in a cruder way. Our results agree with the previous results modulo the definition of the center of mass. It should be noted, however, that we have defined our multipole moments M ij A and I ij A as a volume integral of Λ µν N which manifestly includes the effect of strong internal gravity of the star. Thus, our results support the applicability of the spin-orbit, spin-spin, and quadrupole-orbit coupling forces to a relativistic compact binary where the component stars have a strong internal gravitational field. On the effect of these multipole moments in the orbital evolution and/or the gravitational waveform (see e.g. [6,7,8,9,15,44,108,109,110,112,128,132,148]).

C A Generalized Equivalence Principle Including The Emission of Gravitational Wave
In this section, we give a simple and divergence free derivation of the equations of motion for a small compact object with mass m around any sort of massive body, following the paper by Fukumoto and co-authors [78]. The object moves along the geodesic determined by the smooth part of the geometry around the object up to order m. We note that the smooth part includes the gravitational waves emitted by the orbital motion of the object. Thus it generalizes the equivalence principle for such a compact object to the emission of gravitational waves. Also our derivation of the equations of motion provides a useful method to explicitly calculate the radiation reaction in the fast motion of a compact body.

C.1 Introduction
The idea that an apple and the Moon fall with the same acceleration led Newton to the discovery of the law of gravity. This is not trivial because the Moon is self-gravitating and an apple is not.
That any test particle moves on a geodesic of an external gravitational field is called the weak equivalence principle (WEP). Including self-gravitating objects in this statement is leads to the strong equivalence principle (SEP). The SEP is now experimentally verified with an accuracy of better than 1.5 × 10 −13 [4]. The verification of WEP and SEP is regarded as one of the most important experiments in physics because it plays a fundamental role in any theory of gravity. In fact, Einstein regarded the WEP as the starting point of his theory of gravity. However it is not obvious that a theory constructed that way respects the SEP, and there have been investigations to confirm this. It is then natural to ask how far one can generalize the principle within the framework of general relativity and other theories of gravity. This is not only of academic interest but also has practical importance. There are already several detectors of gravitational waves around the world and we are expecting to directly detect gravitational waves from astronomical sources in the near future and hopefully to open a new window to the universe by gravitational waves. In order to use gravitational waves as a practical tool in astronomy, it is definitely necessary to have a good understanding of the equations of motion for systems with more general situations such as small compact objects like a neutron star/black hole moving at an arbitrary speed in an arbitrary external field. In such a situation the perturbation of the external field including gravitational waves generated by the orbital motion is not negligible. This is exactly the situation we have in mind here and for which we would like to generalize the equivalence principle. In this respect it should be mentioned that Mino et al. and others derived the equations of motion for a point particle with mass m which is represented by a Dirac delta distribution source in an arbitrary background. The equation is interpreted as the geodesic equation on the geometry determined by the external field and the so-called tail part of the self-field of the particle in the first order in m [67,121,133]. Furthermore, Mino et al. used another approach, the matched asymptotic expansion, to obtain the equations of motion without employing the concept of a point particle and thus avoiding divergences in their derivation.
We avoid using a singular source and make use of the point particle limit to derive directly the geodesic equations on the smooth part of the geometry around the object. Here the point particle limit is the strong field point particle limit [81]. The smooth part includes the gravitational waves emitted by the orbital motion of the object, and thus the equivalence principle is generalized to including the emission of gravitational waves. We believe that our approach simplifies the proof that the Mino-Sasaki-Tanaka equations of motion are applicable to a nonsingular source where Mino et al. used the matched asymptotic expansion for their proof.

C.2 Equations of motion
Let us start explaining again our situation and discuss the strong field point particle limit [81]. A small spherical compact object with mass m moves at an arbitrary speed around a massive body with mass M . We would like to find the equations of motion for the object including radiation reaction. We assume that the object is stationary except for higher order tidal effects, so that we can safely neglect the emission of gravitational waves from the object itself, but of course we cannot neglect the gravitational waves emitted by the orbital motion of the object. We denote the world line of the center of mass of the object as z µ (τ ) and define the body zone of the object as follows. We imagine a spherical region around z µ (τ ) and a radius that scales as . At the same time we scale the linear dimension of the object as 2 so that the boundary of the body zone is located at the far zone of the object. Namely we are able to have a multipole expansion of the field generated by the object at the surface of the body zone. We also implicitly assume that the mass of the object scales as so that the compactness of the object remains constant in the point particle limit as → 0. This is why we call this limit the strong field point particle limit. Then we calculate the metric perturbation induced by the small object in this limit. The smallness parameter has a dimension of length which characterizes the smallness of the object. One may regard it as the ratio between the physical scale of the object and the characteristic scale of the background curvature. We assume that the background metric g µν satisfies the Einstein equations in vacuum. So the Ricci tensor of the background vanishes. Since we have assumed that the mass scale of the small object is much smaller than the scale of the gravitational field of the background geometry, we approximate the metric perturbation by the linear perturbation of the small particle h µν .
We will work in the harmonic gauge,h µν ν = 0, where the semicolon means that the covariant derivative with respect to the background metric and the trace-reversed variable is defined as usual, Then the linearized Einstein equations take the following form: This can be solved formally as follows, where we have used the retarded tensor Green's function defined by G µναβ (x, y) = 1 4π θ(Σ(x), y) u µναβ (x, y)δ(σ(x, y)) + v µναβ (x, y)θ(−σ(x, y)) .
For the general tensor Green's function and the definition of σ(x, y),ḡ µα (x, y), u µναβ (x, y), v µναβ (x, y), and Σ(x), please refer to Mino et al. [121] and DeWitt and Brehme [68]. Now we take the point particle limit. In this limit the above metric perturbation contains terms of different dependence which makes the calculation of the equations of motion simple. For example, terms with negative powers of appear from the Dirac delta distribution part of the Green's function: sh µν (x) = 2 d 4 y √ −g u µν αβ (x, y)δ(σ(x, y))T αβ (y).
Living Reviews in Relativity http://www.livingreviews.org/lrr-2007-2 to m 2 a µ / and can be renormalized to the mass of the small object, vanish in the limit or the angular integration. Remembering that the LL tensor is bilinear in the Christoffel symbols and the Christoffel symbols are the derivatives of the metric tensor, one may realize that the only remaining terms come from the combination of the 0th order of the smooth part of the metric and 1/ part of the self-field which is given by Equation (252). The remaining part of the self-field is the so-called tail part that is regular at the object which is the only relevant point in our calculation. The field point (x, τ x ) is now on the body zone boundary, which is defined by 2σ(x, z(τ x )) = and [σ ;α (x, z(τ ))ż α (τ )] τ =τx = 0.
Using this expression we have the following result in the point particle limit, where the Christoffel symbols here are calculated in terms of the smooth part of the metric g s . Then the ADM mass is related to the four-momentum as follows, which is supported by the higher order post-Newtonian approximation [91,95,91]: Finally we have which is the geodesic equation on the geometry determined by the smooth part of the metric around the compact object. In fact, the spin effect on the equations of motion can be derived in a similar way and the standard result [152] can be obtained.
We have proved that a small compact object moves on the geodesic determined by only the smooth part of the geometry around the object. Thus the equations of motion are automatically obtained by determining the geometry around the object which is of course an implicit functional of the world line of the object. The smooth part contains the gravitational waves emitted by the orbital motion so that this equation includes the damping force due to radiation reaction. Our method avoids using a singular source in the first place by making use of the strong field point particle limit. All the quantities should be evaluated at the surface of the body zone boundary and thus we only need the dependence of the distance from the center of the object, namely the dependence of the field. In this way we are able to avoid using any divergent quantities in any part of our calculation. This strongly suggests that our method may be used to get unique equations of fast motion with radiation reaction. This will be investigated in future publications.
In this section, we have assumed spherical symmetry of the compact object except for the tidal effect. It is straightforward to generalize the case to multipole moments in our formalism. This will also be studied in future publications.