{'fraction_non_alphanumeric': 0.03556553784426176, 'fraction_numerical': 0.015569462918947207, 'mean_word_length': 4.112027644632706, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Third-and-a-half order post-Newtonian equations of motion for relativistic compact binaries using the strong field point particle limit

We report our rederivation of the equations of motion for relativistic compact binaries through the third-and-a-half post-Newtonian (3.5 PN) order approximation to general relativity using the strong field point particle limit to describe self-gravitating stars instead of the Dirac delta functional. The computation is done in harmonic coordinates. Our equations of motion describe the orbital motion of the binary consisting of spherically symmetric non-rotating stars. The resulting equations of motion fully agree with the 3.5 PN equations of motion derived in the previous works. We also show that the locally defined energy of the star has a simple relation with its mass up to the 3.5 PN order.


I. INTRODUCTION
A relativistic compact binary (e.g., a neutron star binary) loses its orbital angular momentum by emitting gravitational waves and coalesces in the end. Such a system is a promising source for the gravitational wave detectors such as CLIO, GEO600, LCGT, LIGO, TAMA300, and VIRGO [1]. However, even with those advanced detectors, a direct detection is not easy. Indeed, because the amplitude of gravitational wave from such a source is expected to be tiny at the Earth compared to the detectors' noises, an efficient detection method is required. One method widely used is matched filtering. When using this technique, it is known that the more accurately we know the shape of the signal, the larger the signal to noise ratio becomes. This in turn means that it is desirable to know the dynamics of the binary accurately when one hopes to increase the number of the detectable events. This purpose can be achieved by using higher order post-Newtonian (PN) equations of motion for two point particles, because up to the last several orbits before the coalescence the binary components have moderately slow orbital velocities and are affected negligibly by tidal effects [2,3]. In fact, it is suggested that at least a third order post-Newtonian correction in the equations of motion are necessary for extraction of astronomical information of the sources [4]. To obtain higher order post-Newtonian corrections, mainly three methods are employed in the literature.
In the first method, the Arnowitt-Deser-Misner (ADM) Hamiltonian is derived by means of a direct PN iteration of the Einstein equations [5,6,7,8,9,10,11,12,13,14,15]. Indeed, the ADM Hamiltonian in the ADM Transverse-Tarceless (TT) gauge is completed up to the 3.5 PN order inclusively [10,15]. In the second method, one assumes that the energy and the angular momentum fluxes at infinity are balanced by the corresponding loses of those in the binary orbital motion. Here the known PN expressions for the energy and angular momentum fluxes at infinity are used. With this second method, the n PN (n = 2.5, 3.5, and 4.5) order corrections to the Newtonian equations of motion are derived [16,17,18,19].
Finally in the third method, PN equations of motion are derived by a direct PN iteration of the Einstein equations [20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40]. (see also [3,4,41,42,43,44,45] for reviews). With this method the equations of motion are completed up to the 3.5 PN order inclusively [34,38]. This paper reports a rederivation of the 3.5 PN correction to the Newtonian acceleration for a spherically symmetric non-rotating self-gravitating star using the same method as that in our previous papers [30,33,35,36,37] which derived the corrections up to and including the 3 PN order (see [45] for a review). We use the third method mentioned in the previous paragraph. However our method is different in the several aspects from those of the previous works [34,38] which also used the third method to derive the 3.5 PN order corrections. Here we briefly mention the two most important points.
The first point is regarding the way to describe the binary components as point particles. Among the previous works, Damour, Blanchet, Schäfer and their collaborators have used Dirac delta functionals to describe binary stars as point particles. However, divergent integrals appear when using Dirac delta functionals in general relativity. Then to regularize those divergences, for example, the work [28] successfully used the Hadamard Partie Finie (HPF) regularization to derive the 2.5 PN equations of motion in harmonic coordinates. Blanchet and Faye have developed the generalized HPF regularization [46] in a Lorentz invariant manner [47] and derived the 3 PN correction in the same gauge except for one and only one numerical coefficient (denoted by λ) which could not be determined within their method [29,32]. Interestingly, two works [48,49] have shown that this coefficient λ actually corresponds to one of the undetermined coefficients (ω static ) reported previously by Jaranowski and Schäfer in their derivation of the 3 PN ADM Hamiltonian in the ADM Transverse-Traceless (TT) gauge [11,12,13]. Indeed these latter works also used Dirac delta functionals and the HPF. It has become clear that the HPF is not an appropriate method for regularization of the divergent integrals at the 3 PN order. Damour, Jaranowski and Schäfer [14] then used the dimensional regularization to derive the 3 PN ADM Hamiltonian in the ADMTT gauge and they finally completed the 3 PN correction that contains no undetermined coefficient. The 3 PN equations of motion in harmonic coordinates were later derived using the combination of the HPF and the dimensional regularization [50]. Their result physically agree with the result of [14]. Incidentally, there found no trouble when using the HPF at the 3.5 PN order and the works [10,15,38] derived the corrections at that order using the HPF. See [4] for a review. Now what is our method to achieve a point particle limit? Futamase [27] proposed the strong field point particle limit. In this limit the binary star is first described as an extended object and then its radius is taken to zero in a specific manner which will be explained in the Section II below. By this limit we obtain a point particle with strong internal gravity that is appropriate as a model of a relativistic compact star. With this limit and the method mentioned below we have successfully derived the 3 PN equations of motion that contains no undetermined coefficient [35,36] and confirmed the earlier result [14].
The second aspect of our method that is different from the others is the way to derive equations of motion. In our method, using the local energy momentum conservation law, we compute the force acting on the star by computing the gravitational energy momentum flux going through a suitably defined surface around the star. This idea was used by Einstein, Infeld and Hoffmann in their derivation of the 1 PN equations of motion [20,21]. On the other hand, it was assumed in some works ( [28] at the 2.5 PN order, [29,32] at the 3 PN order, and [34,38] at the 3.5 PN order) that the binary star follows a "geodesic" of spacetime.
However for an equal mass binary, the metric components diverge at the position of the star when it is described as a point particle. In those works [28,29,32,38], the divergence was regularized by the HPF (at the 2.5 and 3.5 PN order), or the combination of the HPF and the dimensional regularization (at the 3 PN order). Our results up to the 3 PN order inclusively have fully agreed with these results. This agreement has confirmed that the star follows the regularized geodesics at least through the 3 PN order. One of the aims of the current paper is to extend it to the 3.5 PN order. (In this respect, Fukumoto et al [51] showed that the star should follows the geodesic of the "smooth part of the full metric" where the smoothing is done by means of the surface integral approach. However, this work showed it for the case of an extreme mass ratio binary.) This paper is organized as follows. The next section explains our method of deriving the equations of motion, sketching each steps in our method with emphases on the strong field point particle limit, the surface integral approach and the scalings of the stress energy tensor of the matter assumed in our approach. Sec. III explains how to derive the PN gravitational field and the PN equations of motion. After showing the formal structure of the 3.5 PN equations of motion in Sec. IV, the 3.5 PN gravitational field in harmonic coordinates is derived in Sec. V. Since we shall not use a Dirac delta functional nor assume any specific form of a momentum velocity relation, we have to derive the relation between the star's mass and its energy and the one between the star's three momentum and its three velocity.
Those relations at the 3.5 PN order are reported in the sections VI and VII. We show the resulting 3.5 PN equations of motion in Sec. VIII.
Throughout this paper, we will use harmonic coordinates and the unit of c = 1 = G where c is the velocity of light and G is the Newton's gravitational constant. A tensor with alphabetical indexes such as x i , denotes a Euclidean three-vector. We raise or lower its indexes with a Kronecker delta. For an object with Greek indexes such as x µ , its indexes are raised or lowered with a flat Minkowskian metric. We shall call our previous papers [30,33,36] Paper I, II, and III, respectively.

II. KEY IDEAS IN OUR DERIVATION
In this section we briefly review our method to obtain post-Newtonian equations of motion, listing the procedures in our method with emphases on the strong field point particle limit, the surface integral approach and the associated scalings of the stress energy tensor of the matter. See [27,30,33,42,44,52,53,54] for more details.
We first introduce an adimensional parameter ǫ which represents the slowness of the star's where z i is the star's representative point and we assume v i ≡ dx i /dτ of order unity. The time coordinate τ is called the Newtonian dynamical time [27,44]. The post-Newtonian scaling implies that the square modulus of the orbital velocity is of order of the inter-body gravitym/r wherem and r are the star's mass and the orbital radius, respectively. In terms of ǫ, this means that the mass of the starm is of order of ǫ 2 [30,33,36]. In what follows, m (without a tilde) denotes the star's mass that is of order unity. The smallness of the orbital velocity and the inter-body gravity are represented in terms of ǫ. Thus, ǫ is a post-Newtonian expansion parameter. Henceforth we call the (τ, x i ) coordinate the near zone coordinate.
The typical magnitude of the velocity of the system far away from the binary is the velocity of the gravitational wave emitted by the binary, or the velocity of light c. Hence there must be different coordinates appropriate in the far zone [55,56]. The near zone coordinate is useful within, say, one wavelength of the gravitational waves (of order ∼ c× (orbital period of the binary)) from the binary system center of the mass. Now, the next subsection sketches the principal steps of our derivation of post-Newtonian equations of motion.

A. Principal steps of our derivation
Suppose we have the n PN metric and the n PN equations of motion. To derive the n + 1 PN equations of motion, we first solve the Einstein equations to the n + 1 PN order and then compute the gravitational force acting on the binary star. When solving the Einstein equations, we use the following techniques.
• Gauge: We use harmonic coordinates to solve the Einstein equations.
• Einstein equations: We solve the harmonically relaxed Einstein equations where we use the Minkowskian wave operator for computational simplicity and where we rewrite the Einstein equations as a set of retarded integrals in the same manner as in [57].
• Boundary condition at infinity: We take the no-incoming radiation condition at the Minkowskian past null-infinity [41,58].
Other previous works [34,38] take essentially the same approaches on all of the above points.
Now it is well-known that the post-Newtonian approximation breaks down far away from the binary at high post-Newtonian orders (e.g., [59]). To obtain the metric, it is thus convenient to divide the spacetime into the near zone and the far zone. The near zone is a time-like world tube surrounding the binary stars. The center of the near zone is defined (if necessary) as the center of the mass of the binary. The intersection between the τ =constant spatial hypersurface and the world tube becomes a spatial 3-sphere. The radius of the sphere is of order of the largest wavelength of the gravitational wave emitted by the binary due to its orbital motion. Its mathematical definition shall be given later in Sec. III A. Then the far zone is the outside of the near zone. The domain of the integration of the retarded integrals is now divided into that in the near zone and that in the far zone. We compute two contributions separately.
In the near zone we can use the post-Newtonian approximation safely. We use the following method to obtain the near zone contribution.
• Near zone contribution: To obtain the near zone metric, we first expand the retarded integrals in ǫ and change the domain of the integration from the Minkowskian null-cone to the τ = constant spatial hypersurface. This hypersurface is denoted as NZ.
• Body zone: The region NZ is divided into three regions. Two of the three are spheres called the body zones each of which surrounds each of the binary star. The remaining is the region outside of the two body zones.
• Multipole expansion: In each body zone, the integrals are evaluated by a multipole expansion.
In the multipole expansion, the multipole moments are defined as the volume integrals on the τ = constant spatial hypersurface whose integrands include the matter's stress energy tensor plus the gravitational stress energy tensor. In the usual post-Newtonian approximation, one assumes that the gravitational field is everywhere weak. Since the gravitational wave detectors search for relativistic compact binaries that have strong internal gravity, we would like to have a method where we use the post-Newtonian approximation only outside of the stars. One procedure which we adopt here to solve this problem is proposed by Futamase [27].
• Strong field point particle limit [27]: With this limit the star's internal gravity can be assumed to be strong. Namely if the companion star were absent, the star's mass would become the ADM mass. We shall explain this limit in more detail in the next subsection.
As in the Newtonian dynamics, the multipole moments are defined with respect to some reference point [60].
• Star's representative point: The representative point of the star (e.g., the center of the mass of the star) is defined by setting the star's dipole moment appropriately in the same manner as in the Newtonian dynamics. The star's multipole moments are defined with respect to this representative point.
As is mentioned in the introduction, this paper studies an intrinsically spherically symmetric star. As the multipole moments above are defined as volume integrals on the τ = constant spatial hypersurface, we cannot obtain a spherically symmetric star by simply making the multipole moments vanish. This is because the Lorentz contraction makes such a star acquire apparent multipole moments when an observer sees the star moving with respect to the observer. We define the star's intrinsic multipole moments in the generalized Fermi normal coordinate [61] and set those zero to obtain an intrinsically spherically symmetric star.
Now, in the far zone, it is well-known that a simple post-Newtonian expansion gives divergent integrals (e.g., [59]). To solve this problem, we use the following method.
• Far zone contribution: we use the Direct Integration of the Relaxed Einstein equations (DIRE) method to evaluate the far zone contribution.
This method was proposed by Will and his collaborators [31,34,62]. This method assumes that the binary is sufficiently stationary in the long past. On the other hand, the work [38] used the multipolar post-Minkowskian (MPM) approach [63] to obtain the far zone metric.
The asymptotic matching is done between the PN inner metric and the MPM outer metric [64].
Having the metric components to the n + 1 PN order, we now derive the n + 1 PN equations of motion for relativistic compact binaries. The following steps are used only in our method.
• Surface integral approach: Using the local conservation law of the stress energy, we obtain expressions for the time derivative of the four momentum in terms of the surface integrals over the star's body zone surface. The integrands shall be the Landau-Lifshitz pseudo-tensor. Thus, the surface integrals give the net flux going through the star's body zone surface and amount to the gravitational force acting on the star.
• Four momentum -mass and velocity relation: As is written in the previous point, we compute the time derivative of the four momentum. To obtain equations of motion, we need to compute the time derivative of the coordinate velocity. In our paper the star's four momentum is defined as the volume integral of the sum of the matter's and gravitational stress energy tensors. It is non-trivial that the four momentum is proportional to the four velocity [60]. We develop a method to derive the expression for the four momentum in terms of the stars' masses and their coordinate velocities.
• Gravitational energy momentum flux: Using the metric components (the solution of the Einstein equations to the n+1 PN order), we evaluate the Landau-Lifshitz pseudotensor to the n + 1 PN order.
• Surface integrals: We evaluate the surface integrals to the n + 1 PN order and obtain the n + 1 PN equations of motion.
The previous works [31,34] use the volume integral approach where they assume a geodesic, multiply it by the conserved baryon density, and integrate it over the star on the time constant hypersurface. Nissanke and Blanchet [38] also assume a geodesic, but in their method they substitute the metric components into the geodesic and regularize divergence due to their use of Dirac delta functionals.
In the following subsections, we shall explain in more detail the strong field point particle limit, the surface integral approach, and the associated scalings of the stress energy tensor components. Other features listed above shall be explained in Sec. III. We apply the procedures listed in this section to the derivation of the 3.5 PN order equations of motion starting from Sec. IV.

B. Strong Field Point Particle Limit
One would describe a star as a point particle by making the radius of the star zero.
However, by this procedure, one would obtain a black hole rather than a point particle.
The radius of the star cannot be made smaller than of order of its mass (in the unit of G = 1 = c). Futamase [27] proposed that we may obtain a point particle model for the star by taking a limit where both the radius of the star and its mass shrink at the same rate.
Since the post-Newtonian scaling implies that the star's massm is O(ǫ 2 ), we assume that its radius is also O(ǫ 2 ). The point particle limit is now achieved in the limit where ǫ goes to zero.
The usual post-Newtonian approximation assumes the gravitational field is everywhere weak even inside the stars. On the other hand, with Futamase's procedures we obtain a point particle with finite internal gravity, since a typical magnitude of the self-gravitational field (the mass over the radius) of the point particle achieved by this procedure is finite irrespective of ǫ. For this reason, this limit is called the strong field point particle limit.
Note that the inter-body gravity (the mass over the orbital separation) is O(ǫ 2 ) and the post-Newtonian approximation applies for the aorbital motion. With this limit, we can derive post-Newtonian equations of motion for a binary star whose component stars have strong internal gravity.

C. Surface Integral Approach and Body Zone
One way to derive equations of motion of a star is a surface integral approach where we compute the total gravitational energy momentum flux going through a surface around the star. When we let the radius of the surface shrink to zero and when the radius of the star goes to zero faster than that of the surface in the point particle limit, we would then obtain equations of motion for the star as a point particle. For this purpose, we introduce two body zones B A for the stars A = 1, 2 as is a representative point of the star A, e.g., the center of the mass of the star A. We shall later fix z i A by specifying the star's dipole moment as in the Newtonian dynamics. The body zone radius ǫR A are much smaller than the orbital separation but is larger than the radius of the star for any ǫ (Recall that the radius of the star decreases proportionally to ǫ 2 in the strong field point particle limit while the body zone radius does so proportionally to ǫ). Note that the two body zones do not overlap each other. Finally, R A are constant, i.e., dR A /dτ = 0.
Other than these conditions, R A are arbitrary. Now the above scalings of the masses and the radii of the stars motivate us to introduce a body zone coordinate for the star A as (τ, ). The scalings of the body zones and the body zone coordinates give us a situation where in the body zone coordinate A the star A does not shrink as ǫ → 0 while the boundary of the body zone expands to infinity. Thus, it is appropriate to define the star's characteristic quantities such as its mass using the body zone coordinate. Moreover, since the body zone boundary ∂B A is far away from the surface of the star A (in its body zone coordinate), we can evaluate explicitly the gravitational energy momentum flux on ∂B A using the post-Newtonian gravitational field. After evaluating the surface integrals, we make the body zone shrink to derive the equations of motion for the compact star.
Possible effects of the internal structures of the compact stars are coded in the multipole moments of the stars. These moments in turn appear in the gravitational energy momentum flux and would affect the orbital motion. However, in this paper we shall concentrate on spherically symmetric stars and ignore those multipole moments.

D. Scalings of the Matter Stress Energy Tensor
The scalings of the radii and the masses of the stars indicate that the matter density is O(ǫ −4 ) in the (t, x i ) coordinate (or ǫ −2 in the near zone coordinate (τ, x i )). We further assume that the internal time scale of the star is comparable to that of the binary orbital motion and is O(ǫ). Namely, in this paper we assume that the star is pressure supported and non-rotating, although an extension to rapidly rotating stars is straightforward [27].
In terms of the stress energy tensor of the matter T µν (or the source terms of the harmonically relaxed Einstein equations Λ µν in Eq. (2) below), these scalings imply T τ τ = O(ǫ −2 ), where the underlined indexes mean that for any tensor S i , S i = ǫ −2 S i and reminds the scaling of the body zone spatial coordinate [27]. Because a star moves in a gravitational field, the tensor components in the body zone coordinate are different from those in the near zone coordinate. Let us denote the matter's stress energy tensor in the body zone coordinate (τ, α i A ) by T µν A and that in the near zone coordinate (τ, Hence the stress energy tensor components of the matter in the near zone coordinate varies with respect to ǫ in the same way as those in the body zone coordinate, or in short,

III. MATHEMATICAL FORMULATION
Based on the idea explained in the previous section, this section formulates our method to derive the equations of motion. We first explain how to solve the Einstein equations, then show how we achieve the surface integral approach.

A. Field Equation
In the surface integral approach, we need to compute the gravitational field near the body zone boundary where the field is well described by the post-Newtonian approximation and slightly deviates from the flat metric. We define a deviation field h µν as where η µν = diag(−ǫ 2 , 1, 1, 1) is the flat metric in the near zone coordinate (τ, x i ) and g is the determinant of the metric. The indexes of h µν are raised or lowered by the flat metric.
Now we impose a harmonic coordinate condition h µν ,ν = 0 where the comma denotes a partial derivative. In the harmonic gauge, we can recast the Einstein equations into a relaxed form, where = η µν ∂ µ ∂ ν is the flat spacetime d'Alembertian. The source term Λ µν of the relaxed Einstein equations consists of two pseudo-tensors. The first is the sum of the stress energy tensor of the stars denoted by T µν and the Landau-Lifshitz pseudo-tensor t µν LL [65]. The second arises due to our use of the flat spacetime d'Alembertian instead of the curved spacetime one. The explicit expressions are The harmonic condition on h µν implies a local energy momentum conservation law; Note that χ µναβ ,αβ itself is divergence free, i.e., χ µναβ ,αβν = 0. We can formally rewrite the relaxed Einstein equations as retarded integrals; where C(τ, x k ) means the past light cone emanating from the event (τ, x k ). We have assumed no homogeneous solution of the relaxed Einstein equations. It is well-known that this condition can be deduced from the Minkowskian no-incoming radiation condition (See, e.g., the section 92 of [58] or the section 6 of [41]).
We solve the Einstein equations as follows. First we split the domain of the integration into two zones: the near zone and the far zone. The near zone is the neighborhood of the gravitational wave source where the wave character of the gravitational radiation is not manifest. In this paper, as in our previous papers, we define the near zone as a time-like world tube surrounding the binary stars. The center of the near zone is defined (if necessary) as the center of the mass of the binary. The intersection between the τ =constant spatial hypersurface and the world tube becomes a spatial 3-sphere. Mathematically, denoting the harmonic coordinate distance from the center by r, the near zone is defined as r < R/ǫ where R/ǫ is of order of the largest wavelength of the gravitational wave emitted by the binary due to its orbital motion. We assume R/ǫ sufficiently large so that the near zone covers the binary stars. Finally, R is constant in time, i.e., dR/dτ = 0. Otherwise R is arbitrary. The ǫ −1 scaling of the near zone radius is derived from the ǫ dependence of the wavelength of the gravitational wave emitted by the binary. The outside of the near zone is the far zone where the retardation effect of the field is crucial.
For the retarded integrals in the far zone, we evaluate it using the Direct Integration of the Relaxed Einstein equations (DIRE) method. This method was proposed by Will and his collaborators [31,34,62]. DIRE directly and nicely fits into our formalism since it utilizes the relaxed Einstein equations in the harmonic gauge. Using DIRE, Pati and Will [31] showed that the far zone contribution to the near zone field affects (physically) the orbital motion starting at the 4 PN order In fact, this result was first obtained by Blanchet and Damour [66]. Although we do not show our explicit computation in this paper, we have followed the DIRE method and checked that the far zone contribution does not affect the equation of motion through the 3.5 PN order. We shall thus focus our attention on the near zone contribution to the near zone field and neglect all the far zone contributions to the near zone field.

B. Near Zone Contribution
For the near zone contribution, we first make retardation expansion and change the domain of the integration to a τ = constant spatial hypersurface where NZ denotes the near zone and we attach the subscript NZ to Λ µν to clarify that they are quantities in the near zone. 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 through all the post-Newtonian order was shown by Pati and Will [31]. This paper uses their method, and hence we safely neglect all the R dependent terms other than those where R appears in arguments of logarithms. This is to make the arguments of all possible logarithmic terms adimensional.
Second we split the integral into two parts; contribution from the body zone B = B 1 ∪B 2 , and from elsewhere, NZ/B. We thus evaluate the following two types of integrals where r A ≡ x − z A . We shall deal with these two contributions successively in the followings.

Body Zone Contribution
As for the body zone contribution, we make multipole expansion using the scaling of the integrand (Λ µν N Z ) in the body zone. For example, the n = 0 part in Eq. (10), h µν Bn=0 , gives where r A ≡ | r A |. The quantity with <> denotes a symmetric and tracefree (STF) operation on the indexes between the brackets. To derive the 3.5 PN equations of motion, we need h τ τ up to O(ǫ 11 ) and h µi up to O(ǫ 9 ).
In the above equations we defined the multipole moments of the star A as where the capital index denotes a set of collective indexes, I l ≡ i 1 i 2 ···i l and α We simply call P µ A the four momentum of the star A, P i A the (three) momentum, and P τ A the energy. Also we call D i A the dipole moment and I ij A the quadrupole moment.
Then we transform these moments into more convenient forms using the conservation law Eq. (6). In the following, v i A ≡ż i A , an overdot denotes a τ time derivative, and y A ≡ y − z A . Noticing that the body zone radii are constant , i.e.,Ṙ A = 0, we have where and The symbol [ ] (or ( )) attached to some of the indexes denotes anti-symmetrization (or symmetrization) on the indexes between the brackets. M ij A is the spin of the star A which is set to be zero for the purpose of the current paper. Eq. (18) gives a momentum-velocity relation.
Thus our momentum-velocity relation is a direct analogue of the Newtonian momentumvelocity relation [60]. In general, we have and Because the energy of the star, P τ A , is not the mass of the star, we shall find the relation between these two in Sec. VI. We may set the dipole moment D i A to be zero to define the center of the mass of the star z i A , but we will later choose a different value for D i A for convenience. As for the multipole moments , we may equate those to zero, as it is our aim to derive equations of motion for a spherically symmetric star.
However, there is a subtlety here, which we will explain in the next section.

A Spherically Symmetric Star and the Star's Multipole Moments
As mentioned in the introduction, this paper studies a binary consisting of two spherically symmetric compact stars. In other words, all the multipole moments of the star defined in an appropriate reference coordinate where effects of its orbital motion and the companion star are removed (modulo, namely, the tidal effect) vanish. We adopt the generalized Fermi normal coordinate (GFC) [61] as the reference coordinate.
We have defined the multipole moments of the star A as Eqs. (15), (16), and (17). Those are defined on the τ = constant three-surface and differ from ones defined in the GFC. Then a question specific to our formalism is if the differences between the multipole moments defined in Eqs. (15), (16), and (17) and the those in the GFC give purely monopole terms.
Roughly speaking, a body that is spherical in its GFC may acquire apparent multipole moments due to, for example, the Lorentz contraction.
At the 3 PN order, this problem was addressed in the appendix C of Paper III and the differences were mainly attributed to the shape of the body zone. The body zone B A which is spherical in the near zone coordinate (NZC) is not spherical in the GFC because of a kinematic effect (Lorentz contraction). In fact, the difference between the GFC quadrupole moment and the NZC one was found to contain monopole terms.
At the 3.5 PN order, no NZC multipole moment was found to contain purely monopole terms. We shall be back to this point later in Sec. V A 2. We equate all those NZC multipoles (l ≥ 1) to zero from now on in this paper.

N Z/B Contribution
About the NZ/B contribution, since the integrand Λ µν N Z = −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 a (super-)potential g( x) which satisfies ∆g( x) = f ( x), ∆ = ∂ i ∂ i denoting Laplacian, we have for each integral (, e.g., n = 0 term in Eq. (11)) Eq. (30) can be proved without using a Dirac delta functional (see Appendix B of Paper III). For n ≥ 1 terms in Eq. (11), we use appropriate (super-)potentials many times to convert all the volume integrals into surface integrals and bulk terms ("−4πg( x)") [67].
Finding the super-potentials is one of the most formidable task especially when we proceed to a high post-Newtonian order. Fortunately, at the 3.5 PN order, all the required superpotentials are available (See [34,38] and Sec. V B below).

C. General Form of the Equations of Motion
From the definition of the four momentum (Eq. (15) with l = 0) and the conservation law Eq. (6) we obtain an evolution equation for the four momentum; We have defined the four momentum by a volume integral of Λ τ ν N Z . Because χ τ ναβ ,αβ = χ τ ναi ,αi , we can transform the volume integral of χ τ ναβ ,αβ into a surface integral form and can evaluate it explicitly in terms of the mass, velocity and the orbital separation. As a result, it is straightforward to see that χ part of Eq. (31) is a trivial identity. The identity to the 3 PN order was shown in Appendix E of Paper III by an explicit calculation. This observation implies that equations of motion can be derived from Θ part of Eq. (31). We thus define the Θ part of the four momentum, dipole moment, and Q i A integral as Correspondingly, we split the momentum-velocity relation Eq. (18)  The Θ part of (31) is Substituting the Θ part of the momentum-velocity relation into the spatial components of Eq. (35), we obtain the general form of the equations of motion for the star A; All the right hand side terms in Eq. (36) except for 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 the star A. In Eq. (36), P τ AΘ rather than the mass of the star A appears. Hence we have to derive a relation between the mass and P τ AΘ . We shall derive the relation by solving the temporal component of the evolution equation (35) functionally. In fact, at the lowest order, we have shown in Paper II that Then we define the mass of the star A as the integrating constant of this equation; m A is the ADM mass that the star A had if the star A were isolated. We took ǫ zero limit in Eq. (38) to ensure that the mass defined above does not include the effect of the companion star and the orbital motion of the star itself. Some subtleties about this definition were discussed in Paper II. By definition m A is constant. The procedure that we solve the evolution equation of P τ AΘ and obtain the mass energy relation is achieved up to the 3.5 PN order successfully and the result will be shown in Sec. VI.

D. On the Arbitrary Constant R A
Our final remark in this section is on the two arbitrary constants R A . Since we introduce the body zones by hand, the arbitrary body zone radii R A seem to appear in the metric, the multipole moments of the stars, and the equations of motion. Paper II has proved that the surface integrals in the general equations of motion Eq. (36) do not depend on R A through any order of the post-Newtonian iteration. The appendix D of Paper III explained that the field and the multipole moments are independent of ǫR A .
Practically, those two accounts justify that we safely discard all the ǫR A dependent terms except for logarithms of ǫR A that appear in the course of computations. We keep ln ǫR A dependent terms to make the arguments of the logarithms adimensional.
We here emphasize that we discard the ǫR A dependent terms in the field first and then evaluate the surface integrals in the general form of the equations of motion using the field independent of ǫR A . We then discard the ǫR A dependent terms arising in the computation of the surface integrals. The details of this procedure were explained in Paper III.

IV. STRUCTURE OF THE 3.5 PN EQUATIONS OF MOTION
In the following sections, we shall derive an acceleration for two spherical compact stars through the third and a half post-Newtonian accuracy. For this purpose, we evaluate the surface integrals in Eqs. (35) and (36) to the appropriate order. As the mass of the star is O(ǫ 2 ), the Newtonian force appears as ǫ 2 correction to the lowest order equations of motion (m A dv i A /dτ = 0). The 3.5 PN order correction, or (v/c) 7 correction to the Newtonian force appears at O(ǫ 7 ) and hence the equations that we have to evaluate to derive an evolution equation for the energy and the equations of motion are See the appendix B for the complete expressions. The field components up to the 2.5 PN order inclusively, ≤9 h τ τ , ≤7 h τ i , and ≤7 h ij , are listed in Paper II. In addition to those, to derive the 3.5 PN mass-energy relation and the 3.5 PN momentum-velocity relation, we have to derive 9 h τ i . To derive the 3 PN equations of motion, we further need 11 h τ τ + 9 h k k . We do not need to use the 3 PN field 8 h τ i , 8 h ij , and 10 h τ τ for the purpose of this paper.
Up to the 2.5 PN order inclusively, the super-potentials required to compute the field could be found [11,28]. At the 3 PN order, it is quite difficult to complete the required superpotentials. We took another method to overcome this problem in Paper III. Fortunately, we can find the super-potentials at the 3.5 PN order and hence derive the 3.5 PN field in a closed form. The next section shows our method of the derivation of the 3.5 PN field.

V. 3.5 PN GRAVITATIONAL FIELD IN HARMONIC COORDINATES
As written in the previous section, our derivation of the 3.5 PN equations of motion requires 9 h τ i and 11 h τ τ + 9 h k k . We start with the body zone contribution.

A. 3.5 PN Body Zone Contribution for a Spherically Symmetric Star
This paper studies equations of motion for a spherically symmetric star and hence discard all the GFC multipole moments of the stars. As stated in Sec. III B 2 and will be shown in Sec. V A 2, up to the 3.5 PN order we can safely discard the NZC multipole moments except for the NZC quadrupole moments I ij A . Using Eqs. (18) to (29), the body zone contribution for a spherically symmetric star becomes In the above equations, D i A may be specified freely to define the star's representative point z i A . Although we explicitly wrote Q K l i A and R K l ij A integrals above, the next section shall show those do not affect the 3.5 PN acceleration.
The surface integrals Q K l i A and R K l ij A did contribute to the field and the equations of motion at the 3 PN order (Paper III).
At the 3.5 PN order, we compute Q K l i A and R K l ij A up to O(ǫ 5 ) to derive the 3.5 PN field. Then the integrands of the surface integrals in Eqs. (23) and (24) are 9 Λ µν N Z . Now suppose those integrands behaves near the star A as 1/r p A with p some positive integer. Then Inspection shows that p = 6 for 9 Λ τ µ N Z and p = 7 for 9 Λ iµ N Z . For those degree l of the collective indexes K l , by explicitly computing those integrals in Eqs. (23) and (24), we found that Q K l i A and R K l ij A do not contribute to the field nor equations of motion at the 3.5 PN order.

A Spherical Body and Multipole Moments
As promised, this section explains the reason we can safely discard the NZC multipole moments except for the NZC quadrupole moments through the 3.5 PN order. Denoting the NZC (GZC) quadrupole moment as I ij A,NZC (I ij A,GFC ), Eq. (C14) of Paper III gave an expression of the difference between these two moments; for a star which is spherically symmetric in the GFC. The coefficients A ij (τ ) and B i (τ ) can be read off from [61] (see also Eq. (C16) and (C17) of Paper III) and are for the star A = 1 and r 12 = z 1 − z 2 = r 2 − r 1 . Note that there is no ǫ 3 term in A ij (τ ) nor B i (τ ) because there is no 0.5 PN field in the harmonic gauge (i.e., 5 h τ τ = 0) as was shown in the appendix E of Paper II. The surface integral in Eq. (47) can be expanded in ǫ as The first integral in the right hand side of the equality is found to be non-zero and contribute to the 3 PN metric 10 h τ τ and hence the 3 PN equations of motion. The second integral could in principle affect the 3.5 PN field, as the coefficient ǫ in front of it indicates. So let us evaluate the second integral. Its integrand is 9 Λ τ τ N Z = − 4 h τ τ ,ij 5 h ij which behaves near the star A = 1 as 9 Λ τ τ N Z ∼ 3r i 1 r j 1 /r 5 1 − δ ij /r 3 1 . The surface integral then results in and thus there is no monopole term in δI ij A at the 3.5 PN order. Similarly, using Eq. (C14) of Paper III, we can check that none of the NZC multipole moments hide any monopole terms at the 3.5 PN order for a star spherically symmetric in the GFC.

B. 3.5 PN N/B Field
This section explains how we derive the 3.5 PN N/B field components 9 h τ i N Z/B and 11 h τ τ N Z/B + 9 h k N Z/B k . The integrals for the latter which we evaluate may be written as Because only a spatial derivative of 11 h τ τ N Z/B + 9 h k N Z/B k appears in 11 [−gt ij LL ] as shown in Eq. (43), the second and the third term in the right hand of the equality in the above equation does not contribute to the 3.5 PN equations of motion. Similarly, 9 h τ i N Z/B may be written as Our task now is to find super-potentials f that satisfies, e.g., 11 Λ τ τ N Z + 9 Λ k N Z k = ∆f . With the super-potentials, we use Eq. (30) and find explicit expressions of the 3.5 PN field components in closed forms. For the higher order retarded expansion terms, for example for the second retarded expansion term of Eq. (52) (the n = 2 term of the first term in Eq. (52)) , we use super-super-potential f ( y) satisfying 9 Λ τ τ N Z + 7 Λ k N Z k = ∆ 2 f . Then the volume integral may be converted into the bulk term and the surface integral terms as The source terms for which we need to find particular solutions of Poisson equations have the following spatial coordinate dependence.
and 1 ↔ 2, i.e., like 1/r 2 2 for 1/r 2 1 . (Specific members of the list depends on how one simplifies the expressions of the source terms of the Einstein equations. For instance, one can always erase r i 2 by using r i 2 = r i 12 + r i 1 .) Our method to derive the super-potentials are heuristic; there are few guidelines available to find the required super-potentials. We proceed as follows. First, we convert all the tensorial sources into scalars with spatial derivatives. For example, 11 Λ τ τ N Z + 9 Λ k N Z k include the following term Paper III gives the explicit expression of (dP τ 1Θ /dτ ) ≤3P N . We can integrate Eq. (57) functionally as The k Γ A up to 3 PN order are given in Paper III. The new result we show in this paper is 12 3 , where V = v 1 − v 2 . The mass-energy relation for the χ part up to the 3.5 PN order is given in the appendix A. Eqs. (58) and (A2) give the 3.5 PN order mass-energy relation in our formalism.
From the definition of P τ AΘ , we expect it equal to √ −gm A u τ A where u τ A is the time component of the 4-velocity of the star A normalized as g µν u µ √ −gu τ A at 3.5 PN order can be written in terms of the deviation field as This expression should be somehow evaluated at the star A. However, because the deviation field h µν diverges at the star in the point particle description, the above expectation is not trivial. As in Paper I, II, and III, we checked that up to the 3.5 PN order inclusively a relation holds where [f ] ext A means that we regularize the quantity f at the star A by the Hadamard's Partie Finie (see e.g. [46]) or whatever regularization which gives the same result. We here emphasize that we have never assumed this "natural" relation in advance. The relation Eq.
The computation is straightforward and we found 7 Q i AΘ = 0. From the 3 PN accurate Q i AΘ calculation, we have a momentum-velocity relation.
where a i 1 is the acceleration of the star 1 and should be replaced by the 3 PN order expression of the acceleration. As in Paper III, we define the representative points of the stars z i A by choosing The first term in Eq. (64) makes the three momentum proportional to v i A at the 3 PN order. The second term gauges away the logarithmic terms from the 3 PN equations of motion [36].
In any case, our choice of D i AΘ affects only the 3 PN order correction and P i 1Θ = P τ 1Θ v i 1 holds as long as we are concerned with the 3.5 PN order correction to the equations of motion.
As the second term in Eq. (64) depends on the arbitrary parameter ǫR A , one may suspect that our equations of motion lose their predictive power on the binary dynamics. This is not the case. When we define the star's representative point by D i AΘ = 0 instead of Eq. (64), logarithmic terms that depend on ǫR A appear in equations of motion [36]. But those terms are mere gauge. Indeed, it can be shown that an observable such as the orbital energy does not depend on ǫR A when that observable is written in terms of gauge independent variables (such as the gravitational wave frequency) [29,35].
The explicit expressions for the 2 PN and 3 PN terms, 4 P τ 1χ , and 6 P τ 1χ , are given in Paper III. The 3.5 PN χ part of the star's energy 7 P τ 1χ does contribute the 3.5 PN field 11 h τ τ (See Eq. (44)) and affects the 3.5 PN equations of motion.