POST-NEWTONIAN HYDRODYNAMIC EQUATIONS USING THE (3+1) FORMALISM IN

BANACH CENTER PUBLICATIONS, VOLUME 41 INSTITUTE OF MATHEMATICS POLISH ACADEMY OF SCIENCES

WARSZAWA 1997

POST-NEWTONIAN HYDRODYNAMIC EQUATIONS USING THE (3+1) FORMALISM IN

GENERAL RELATIVITY

H I D E K I A S A D A

Department of Earth and Space Science

Graduate School of Science, Osaka University, Toyonaka, Osaka 560, Japan E-mail: asada@vega.ess.sci.osaka-u.ac.jp

Abstract. The post-Newtonian (PN) hydrodynamic equations are obtained in the (3+1) formalism, which include the 2.5PN order as the reaction due to the quadrupole gravitational radiation. These equations are valid in various slice conditions, while we adopt a kind of trans- verse gauge condition to determine the shift vector. In particular, we describe methods to solve the 2PN tensor potential which arises from the spatial 3-metric. Our formulaton in the PN ap- proximation using the (3+1) formalism will be useful for numerical computations providing an initial data for the final merging phase of coalescing binary neutron stars which can be treated only by fully general relativistic simulations.

1. Introduction. Aiming at the direct detection of gravitational waves from rela- tivistic astrophysical objects, kilometer-size interferometric gravitational wave detectors, such as LIGO [1] and VIRGO [2], are now under construction. Coalescing binary neutron stars are the most promising sources of gravitational waves for such detectors. After a long time emission of gravitational waves, the orbital separation of binary becomes com- parable to the radius of the neutron star. Then, each of binary neutron stars begins to behave as a hydrodynamic object, not as a point particle, because they are tidally coupled with each other. Recently, Lai, Rasio and Shapiro [3] have pointed out that such a tidal coupling of binary neutron stars is very important for their evolution in the final merging phase because the tidal effect causes the instability of their circular motion. The general relativistic gravity is also important because in such a phase, the orbital separation is less than ten times the Schwarzschild radius of the system. Thus, we need not only a hydrodynamic treatment, but also a general relativistic one to study the final phase of binary neutron stars.

1991 Mathematics Subject Classification: 83C25, 83C35.

The paper is in final form and no version of it will be published elsewhere.

[85]

Surely fully general relativistic simulation must be the best method, but it is also one of the most difficult ones. Although much effort has been focused and much progress can be expected there [4], it will take a long time until numerical relativistic calculations become reliable. One of the main reasons is that we do not know the behavior of the geometric variables in the strong field around coalescing binary neutron stars. According to this, we do not know what sort of gauge condition is useful and how to give an appropriate general relativistic initial condition for coalescing binary neutron stars.

The other reason is a technical one: In numerical relativistic simulations, gravitational waves are generated when we need to cover a region L > λ ∼ the wavelength in order to perform accurate simulations. This is in contrast with the case of Newtonian and/or PN simulations, in which we only need to cover a region (λ >)L > R ∼ the orbital separation.

At present, we had better search other methods to prepare an initial condition for binary neutron stars.

The PN approximation in the fluid was pioneered by Chandrasekhar et al. [5] who ob- tained the equation of motion up to the 2.5PN order. However, their expressions include potentials for non-compact sources which should be solved rather carefully in practical cases. On the other hand, using the ADM gauge, Blanchet, Damour and Sch¨ afer obtained the formula, including the 1PN effects and the quadrupole radiation reaction (2.5PN ef- fect) [6], which only consists of potentials for compact sources. Their formula actually works well in simulations of the coalescing binary neutron stars [7]. Our aim is to establish the formulation up to the 2.5PN hydrodynamic equation for a fairly wide class of gauge conditions and, in particular, to provide methods to solve the 2PN tensor potential which does not appear at the 1PN order [6].

We use the units of c = G = 1 in this paper. Greek and Latin indices take 0, 1, 2, 3 and 1, 2, 3, respectively.

2. (3+1) Formalism for post-Newtonian approximation (3+1) Formalism. In the (3+1) formalism, the metric is split as

g µν = γ µν − ˆ n µ n ˆ ν , (2.1) where ˆ n µ = (−α, 0). Then the line element is written as

ds ² = −(α ² − β i β ⁱ )dt ² + 2β i dtdx ⁱ + γ ij dx ⁱ dx ^j . (2.2) In order to distinguish the wave part from the non-wave part (for example, Newtonian potential) in the metric, we use ˜ γ ij = ψ ⁻⁴ γ ij instead of γ ij where ψ = det(γ ij ) ^1/12 so that det(˜ γ ij ) = 1 is satisfied. Using the extrinsic curvature, K ij , we define ˜ A ij as

A ˜ _ij ≡ ψ ⁻⁴ A _ij ≡ ψ ⁻⁴  K _ij − 1

3 γ _ij K 

. (2.3)

The Einstein equation is split into the constraint equations and the evolution equations.

The former are called Hamiltonian and momentum constraints which respectively become

∆ψ = ˜ 1

8 tr ˜ Rψ − 2πρ H ψ ⁵ − ψ ⁵ 8

 ˜ A ij A ˜ ^ij − 2 3 K ² 

, (2.4)

D ˜ j (ψ ⁶ A ˜ ^j _i ) − 2

3 ψ ⁶ D ˜ i K = 8πψ ⁶ J i , (2.5)

where we introduce ρ H = T µν n ^µ n ^ν and J j = −T µν n ^µ γ _j ^ν , and ˜ D i , ˜ ∆ and tr ˜ R are the covariant derivative, Laplacian and the scalar curvature with respect to ˜ γ ij . Evolution equations for the spatial metric and extrinsic curvature are respectively

∂

∂n ˜ γ ij = − 2α ˜ A ij + ˜ γ il

∂β ^l

∂x ^j + ˜ γ jl

∂β ^l

∂x ⁱ − 2 3 γ ˜ ij

∂β ^l

∂x ^l , (2.6)

∂

∂n

A ˜ ij = 1 ψ ⁴

h α 

R ij − 1

3 γ ij trR 

−  ˜ D i D ˜ j α − 1

3 ˜ γ ij ∆α ˜ 

− 2 ψ

 ψ _,i α _,j + ψ _,j α _,i − 2

3 ˜ γ _ij ˜ γ ^kl ψ _,k α _,l i + α(K ˜ A ij − 2 ˜ A il A ˜ ^l _j ) + ∂β ^m

∂x ⁱ

A ˜ mj + ∂β ^m

∂x ^j

A ˜ mi − 2 3

∂β ^m

∂x ^m A ˜ ij

− 8π α ψ ⁴

 S ij − 1

3 γ ij S ^l _l 

, (2.7)

∂

∂n ψ = ψ 6

 −αK + ∂β ⁱ

∂x ⁱ



, (2.8)

∂

∂n K = α  ˜ A _ij A ˜ ^ij + 1 3 K ² 

− 1 ψ ⁴

∆α − ˜ 2

ψ ⁵ γ ˜ ^kl ψ _,k α _,l + 4πα(S ⁱ _i + ρ _H ), (2.9) where R _ij , γ, S _ij and _∂n ^∂ denote, respectively, the Ricci tensor with respect of γ _ij , deter- minant of γ ij , T kl γ ^k _i γ ^l _j and _∂t ^∂ − β ^{i ∂} _∂x

i

.

In order to clarify the wave property of ˜ γ ij , we impose a kind of transverse gauge* to h _ij as

h ij,j = 0. (2.10)

Finally, we consider the matter as the perfect fluid, T ^µν = (ρ + ρε + P )u ^µ u ^ν + P g ^µν , where u ^µ , ρ, ε and P are the four velocity, the mass density, the specific internal energy and the pressure.

Post-Newtonian approximation. Each metric variable, which is relevant to the present paper, is expanded as

ψ = 1 + ₍₂₎ ψ + ₍₄₎ ψ + ₍₆₎ ψ + ₍₇₎ ψ + . . . , α = 1 − U +  U ²

2 + X 

+ ₍₆₎ α + ₍₇₎ α + . . . , β ⁱ = ₍₃₎ β _i + ₍₅₎ β _i + ₍₆₎ β _i + ₍₇₎ β _i + ₍₈₎ β _i + . . . , h ij = ₍₄₎ h ij + ₍₅₎ h ij + . . . ,

A ˜ ij = (3) A ˜ ij + (5) A ˜ ij + (6) A ˜ ij + . . . , K = ₍₃₎ K + ₍₅₎ K + ₍₆₎ K + . . . ,

(2.11)

where subscripts denote the PN order (c ⁻ⁿ ) and U is the Newtonian potential satisfying

∆ _{f lat} U = −4πρ. (2.12)

* Hereafter, we call this condition merely the transverse gauge.

From Eqs.(2.6) and (2.7), the wave equation for h ij is written as

 − ∂ ²

∂t ² + ∆ f lat

 h ij

=  1 − α ²

ψ ⁴

 ∆ _{f lat} h _ij +  ∂ ²

∂n ² − ∂ ²

∂t ²

 h _ij

+ 2α ψ ⁴

h − 2α ψ

 ˜ D _i D ˜ _j − 1 3 ˜ γ _ij ∆ ˜ 

ψ + 6α ψ ²

 ˜ D _i ψ ˜ D _j ψ − 1

3 γ ˜ _ij D ˜ _k ψ ˜ D ^k ψ 

−  ˜ D _i D ˜ _j − 1 3 ˜ γ _ij ∆ ˜ 

α − 2 ψ

 ˜ D _i ψ ˜ D _j α + ˜ D _j ψ ˜ D _i α − 2

3 γ ˜ _ij D ˜ ^k ψ ˜ D _k α i + 2α ² 

K ˜ A ij − 2 ˜ A il A ˜ ^l _j  + 2α 

β _{, i} ^m A ˜ mj + β ^m _,j A ˜ mi − 2

3 β ^m _,m A ˜ ij



− 16π α ² ψ ⁴

 S ij − 1

3 γ ij S ^l _l 

− ∂

∂n



β ^m _,i γ ˜ mj + β ^m _,j γ ˜ mi − 2 3 β ^m _,m γ ˜ ij

 + 2 ∂α

∂n A ˜ ij

≡τ ij .

(2.13)

Thus, ₍₄₎ h _ij in the near zone is determined from

∆ _{f lat(4)} h _ij = ₍₄₎ τ _ij , (2.14)

while ₍₅₎ h ij in the near zone is derived from

(5) h ij (t) = 1 4π

∂

∂t Z

(4) τ ij (t, y)d ³ y. (2.15) This integral is evaluated so that we obtain [6,8]

(5) h _ij (t) = − 4 5

d ³ dt ³

 I _ij − 1

3 δ _ij I _kk 

, (2.16)

where I ij = R ρx ⁱ x ^j d ³ x.

Up to the 2.5PN order, the hydrodynamic equations become

∂S _i

∂t + ∂(S _i v ^j )

∂x ^j = − 

1 + 2U + 5

4 U ² + 6 ₍₄₎ ψ + X  P _,i + ρ _∗ h

U ,i

n

1 + ε + P ρ + 3

2 v ² − U + 5

8 v ⁴ + 4v ² U +  3

2 v ² − U 

ε + P ρ

 + 3 ₍₃₎ β _j v ^j o

− X ,i



1 + ε + P ρ + v ²

2



+ 2v ² ₍₄₎ ψ ,i − ₍₆₎ α ,i − ₍₇₎ α ,i

+ v ^j n

(3) β j,i



1 + ε + P ρ + v ²

2 + 3U 

+ ₍₅₎ β j,i + ₍₆₎ β j,i

o

+ ₍₃₎ β _{j (3)} β j,i

+ 1

2 v ^j v ^k ( ₍₄₎ h jk,i + ₍₅₎ h jk,i ) + O(c ⁻⁸ ) i

, (2.17)

∂H

∂t + ∂(Hv ^j )

∂x ^j = − P h

v ^j _,j + ∂

∂t

 1

2 v ² + 3U  + ∂

∂x ^j n 1

2 v ² + 3U  v ^j o

+ O(c ⁻⁵ ) i

, (2.18)

where S i = αψ ⁶ (ρ + ρε + P )u ⁰ u i , S ⁰ = αψ ⁶ (ρ + ρε + P )(u ⁰ ) ² and H = αψ ⁶ ρεu ⁰ .

3. Strategy to obtain 2PN tensor potential. We describe methods to solve the equation for the 2PN tensor potential ₍₄₎ h ij . Although Eq.(2.14) is formally solved as

(4) h ij (t, x) = − 1 4π

Z (4) τ _ij (t, y)

|x − y| d ³ y, (3.1)

it seems difficult to estimate this integral in practice since (4) τ ij → O(r ⁻³ ) for r → ∞ and the integral is taken all over the space. Thus it is desirable to replace this equation by some tractable forms in numerical evaluation. In the following, we show two approaches:

In section 3.1, we change Eq. (3.1) into the form in which the integration is performed only over the matter distribution like the Newtonian potential. In section 3.2, we propose a method to solve Eq. (2.14) as a boundary value problem.

Direct integration method. The explicit form of ₍₄₎ τ _ij is

(4) τ ij = − 2 ˆ ∂ ij (X + 2 ₍₄₎ ψ) + U ˆ ∂ ij U − 3U ,i U ,j + δ ij U ,k U ,k − 16π 

ρv ⁱ v ^j − 1

3 δ ij ρv ² 

− 

(3) β ˙ i,j + ₍₃₎ β ˙ j,i − 2

3 δ _{ij (3)} β ˙ k,k



, (3.2)

where ˆ ∂ _ij ≡ _∂x ^∂

i

∂x

^{2 j}

− ¹ ₃ δ _ij ∆ _{f lat} . Since ₍₃₎ β ˙ _i is written as

(3) β ˙ _i = ˙ p _i − (X + 2 ₍₄₎ ψ) _,i − n Z



ρv ² + 3P − ρU/2 

|x − y| d ³ y o

,i

, (3.3)

the source term, ₍₄₎ τ ij , is split into five parts ₍₄₎ τ ij = ₍₄₎ τ _ij ^(S) + ₍₄₎ τ _ij ^{(U )} + ₍₄₎ τ _ij ^(C) + ₍₄₎ τ _ij ^(ρ) +

(4) τ _ij ^{(V )} , where we introduce

(4) τ _ij ^(S) = − 16π 

ρv ⁱ v ^j − 1

3 δ ij ρv ²  ,

(4) τ _ij ^{(U )} =U U ,ij − 1

3 δ ij U ∆ f lat U − 3U ,i U ,j + δ ij U ,k U ,k ,

(4) τ _ij ^(C) =4 ∂

∂x ^j

Z (ρv ⁱ ) ^·

|x − y| d ³ y + 4 ∂

∂x ⁱ

Z (ρv ^j ) ^·

|x − y| d ³ y − 8 3 δ ij

∂

∂x ^k

Z (ρv ^k ) ^·

|x − y| d ³ y,

(4) τ _ij ^(ρ) = ˆ ∂ ij

Z

¨

ρ|x − y|d ³ y,

(4) τ _ij ^{(V )} =2 ˆ ∂ ij

Z



ρv ² + 3P − ρU/2 

|x − y| d ³ y.

(3.4)

Thus it becomes clear that ₍₄₎ h ij and ₍₅₎ h ij are expressed in terms of only matter vari- ables, without metric, through ₍₄₎ τ _ij and thus they do not depend on slicing conditions.

Then we define ₍₄₎ h ^(∗) _ij = ∆ ⁻¹ _{f lat(4)} τ _ij ^(∗) where the symbol ‘∗’ denotes S, U , C, ρ and V . First, since ₍₄₎ τ _ij ^(S) is a compact source, we immediately obtain

(4) h ^(S) _ij = 4 Z

 ρv ⁱ v ^j − ¹ ₃ δ _ij ρv ² 

|x − y| d ³ y. (3.5)

Second, we consider the following equation

∆ _{f lat} G(x, y ₁ , y ₂ ) = 1

|x − y 1 ||x − y 2 | . (3.6) It is possible to write ₍₄₎ h ^{(U )} _ij as integrals over the matter using the function G. Eq. (3.6) has solutions [9], ln(r 1 + r 2 ± r 12 ), where r 1 = |x − y 1 |, r 2 = |x − y 2 |, r 12 = |y 1 − y 2 |. Note that ln(r 1 +r 2 −r 12 ) is not regular on the interval between y 1 and y 2 , while ln(r 1 +r 2 +r 12 ) is regular on the matter. Thus we use ln(r ₁ + r ₂ + r ₁₂ ) as a Green function. Using this function, U U ,ij and U ,i U ,j are rewritten as

U U ,ij =∆ f lat

Z

d ³ y 1 d ³ y 2 ρ(y 1 )ρ(y 2 ) ∂ ²

∂y ₁ ⁱ ∂y ₁ ^j ln(r 1 + r 2 + r 12 ), U _,i U _,j =∆ _{f lat}

Z

d ³ y ₁ d ³ y ₂ ρ(y ₁ )ρ(y ₂ ) ∂ ²

∂y ₁ ⁱ ∂y ₂ ^j ln(r ₁ + r ₂ + r ₁₂ ).

(3.7)

Thus we can express (4) h ^{(U )} _ij using the integral over the matter as

(4) h ^{(U )} _ij = Z

d ³ y ₁ d ³ y ₂ ρ(y ₁ )ρ(y ₂ ) h ∂ ²

∂y ₁ ⁱ ∂y ^j ₁ − 1 3 δ _ij 4 1

 −3  ∂ ²

∂y ₁ ⁱ ∂y ₂ ^j − 1 3 δ _ij 4 12

i ln(r ₁ +r ₂ +r ₁₂ ), (3.8)

where we introduce 4 ₁ = _∂y ^∂

_k²

1

∂y

₁^k

and 4 ₁₂ = _∂y ^∂

_k²

1

∂y

₂^k

. Using relations ∆ _{f lat} |x − y| = 2/|x − y| and ∆ f lat |x − y| ³ = 12|x − y|, ₍₄₎ h ^(C) _ij , ₍₄₎ h ^(ρ) _ij and ₍₄₎ h ^{(V )} _ij are solved as

(4) h ^(C) _ij = 2 ∂

∂x ⁱ Z

(ρv ^j ) ^· |x − y|d ³ y + 2 ∂

∂x ^j Z

(ρv ⁱ ) ^· |x − y|d ³ y + 4 3 δ ij

Z

¨

ρ|x − y|d ³ y, (3.9)

(4) h ^(ρ) _ij = 1 12

∂ ²

∂x ⁱ ∂x ^j Z

¨

ρ|x − y| ³ d ³ y − 1 3 δ ij

Z

¨

ρ|x − y|d ³ y,

(3.10)

(4) h ^{(V )} _ij = ∂ ²

∂x ⁱ ∂x ^j Z 

ρv ² +3P − ρU 2

 |x−y|d ³ y− 2 3 δ ij

Z



ρv ² + 3P − ρU/2 

|x − y| d ³ y. (3.11) In this way, we obtain ₍₄₎ h ij = ₍₄₎ h ^(S) _ij + ₍₄₎ h ^{(U )} _ij + ₍₄₎ h ^(C) _ij + ₍₄₎ h ^(ρ) _ij + ₍₄₎ h ^{(V )} _ij .

Treatment as a boundary value problem. The above expression for ₍₄₎ h _ij is quite inter- esting because it only consists of integrals over the matter. However, in actual numerical simulations, it will take a very long time to perform the direct integration. Therefore, we also propose other strategies where Eq. (2.14) is solved as the boundary value problem.

Here, we would like to emphasize that the boundary condition should be imposed at r(= |x|)  |y ₁ |, |y ₂ |, but r does not have to be greater than λ, where λ is a typical wave length of gravitational waves. We only need to impose r > R (a typical size of matter distribution). This means that we do not need a large amount of grid numbers compared with the case of fully general relativistic simulations, in which we require r > λ  R.

First of all, we consider the equation

∆ _{f lat} 

(4) h ^(S) _ij + ₍₄₎ h ^{(U )} _ij 

= ₍₄₎ τ _ij ^(S) + ₍₄₎ τ _ij ^{(U )} . (3.12)

Since its source term behaves as O(r ⁻⁶ ) at r → ∞, this equation can be accurately solved under the boundary condition at r > R as

(4) h ^(S) _ij + ₍₄₎ h ^{(U )} _ij = 2 r

 ¨ I _ij − 1 3 δ _ij I ¨ _kk  + 2

3r ²



n ^k I ¨ _ijk − 1

3 δ _ij n ^k I ¨ _llk + 2n ^k ( ˙ S _ikj + ˙ S _jki ) − 4

3 δ _ij n ^k S ˙ _lkl 

+ O(r ⁻³ ), (3.13)

where I ijk = R ρx ⁱ x ^j x ^k d ³ x and S ijk = R ρ(v ⁱ x ^j − v ^j x ⁱ )x ^k d ³ x.

(4) h ^(C) _ij , ₍₄₎ h ^(ρ) _ij and ₍₄₎ h ^{(V )} _ij can be rewritten as

(4) h ^(C) _ij = 2 Z

(ρv ^j ) ^· x ⁱ − y ⁱ

|x − y| d ³ y + 2 Z

(ρv ⁱ ) ^· x ^j − y ^j

|x − y| d ³ y − 4 3 δ ij

Z

(ρv ^k ) ^· x ^k − y ^k

|x − y| d ³ (3.14) y,

(4) h ^(ρ) _ij = 1 4

∂ ²

∂x ⁱ ∂x ^j Z

ρv ^k v ^l (x ^k − y ^k )(x ^l − y ^l )

|x − y| d ³ y + 1 3 δ _ij

Z

(ρv ^k ) ^· x ^k − y ^k

|x − y| d ³ y + 1

2

 ∂

∂x ⁱ Z

P ⁰ (x ^j − y ^j )

|x − y| d ³ y + ∂

∂x ^j Z

P ⁰ (x ⁱ − y ⁱ )

|x − y| d ³ y



− 1 8

 2

Z

ρ U _,j (x ⁱ − y ⁱ ) + U _,i (x ^j − y ^j )

|x − y| d ³ y + x ^k ∂

∂x ⁱ Z

ρ U ,k (x ^j − y ^j )

|x − y| d ³ y + x ^k ∂

∂x ^j Z

ρ U ,k (x ⁱ − y ⁱ )

|x − y| d ³ y

 ,(3.15)

(4) h ^{(V )} _ij = 1 2

h ∂

∂x ⁱ Z 

ρv ² + 3P − ρU 2

 x ^j − y ^j

|x − y| d ³ y + ∂

∂x ^j Z 

ρv ² + 3P − ρU 2

 x ⁱ − y ⁱ

|x − y| d ³ y i

− 2 3 δ ij

Z



ρv ² + 3P − ρU/2 

|x − y| d ³ y, (3.16)

where P ⁰ = P + ρv ² /4 + ρU ,l y ^l /4. From the above relations, ₍₄₎ h ^(C) _ij , ₍₄₎ h ^(ρ) _ij and ₍₄₎ h ^{(V )} _ij become

(4) h ^(C) _ij = 2(x ⁱ ₍₃₎ P ˙ ^j + x ^j ₍₃₎ P ˙ ⁱ − Q ij ) + 4

3 δ _ij  Q _kk

2 − x ^k (3) P ˙ ^k 

, (3.17)

(4) h ^(ρ) _ij = 1 4

∂ ²

∂x ⁱ ∂x ^j

 V _kl ^(ρv) x ^k x ^l − 2V _k ^(ρv) x ^k + V ^(ρv)  + 1

3 δ _ij 

x ^k ₍₃₎ P ˙ _k − Q _kk 2



+ 1 2

n ∂

∂x ⁱ



V ^{(P )} x ^j − V _j ^{(P )}  + ∂

∂x ^j



V ^{(P )} x ⁱ − V _i ^{(P )} o

− 1 8 n

2 

x ⁱ V _j ^{(ρU )} + x ⁱ V _i ^{(ρU )} − V _ij ^{(ρU )} − V _ji ^{(ρU )}  + x ^k ∂

∂x ⁱ



x ^j V _k ^{(ρU )} − V _lj ^{(ρU )} 

+ x ^k ∂

∂x ^j



x ⁱ V _l ^{(ρU )} − V _li ^{(ρU )} o ,

(3.18)

(4) h ^{(V )} _ij = 1 2



Q ^(I) _,j x ⁱ + Q ^(I) _,i x ^j − Q ^(I) _i,j − Q ^(I) _j,i  + 1

3 Q ^(I) δ ij , (3.19)

where

∆ _{f lat(3)} P i = −4πρv ⁱ ,

∆ f lat Q ij = −4π n

x ^j (ρv ⁱ ) ^· + x ⁱ (ρv ^j ) ^· o ,

∆ f lat Q ^(I) = −4π 

ρv ² + 3P − 1 2 ρU 

,

∆ f lat Q ^(I) _i = −4π 

ρv ² + 3P − 1 2 ρU 

x ⁱ ,

∆ f lat V _ij ^(ρv) = −4πρv ⁱ v ^j ,

∆ f lat V _i ^(ρv) = −4πρv ⁱ v ^j x ^j ,

∆ _{f lat} V ^(ρv) = −4πρ(v ^j x ^j ) ² ,

∆ f lat V ^{(P )} = −4πP ⁰ ,

∆ _{f lat} V _i ^{(P )} = −4πP ⁰ x ⁱ ,

∆ f lat V _i ^{(ρU )} = −4πρU ,i ,

∆ f lat V _ij ^{(ρU )} = −4πρU ,i x ^j .

(3.20)

Therefore, ₍₄₎ h ^(C) _ij , ₍₄₎ h ^(ρ) _ij and ₍₄₎ h ^{(V )} _ij can be derived from the above potentials which satisfy the Poisson equations with compact sources.

We note that, instead of the above procedure, we may solve the Poisson equation for

(4) h _ij carefully imposing the boundary condition for r  R as

(4) h ij = 1 r

n 1

4 I _ij ⁽²⁾ + 3 4 n ^k 

n ⁱ I _kj ⁽²⁾ + n ^j I _ki ⁽²⁾ 

− 5

8 n ⁱ n ^j I _kk ⁽²⁾ + 3

8 n ⁱ n ^j n ^k n ^l I _kl ⁽²⁾ + 1

8 δ ij I _kk ⁽²⁾ − 5

8 δ ij n ^k n ^l I _kl ⁽²⁾ o + 1

r ² hn − 5

12 n ^k I _ijk ⁽²⁾ − 1

24 (n ⁱ I _jkk ⁽²⁾ + n ^j I _ikk ⁽²⁾ ) + 5

8 n ^k n ^l (n ⁱ I _jkl ⁽²⁾ + n ^j I _ikl ⁽²⁾ )

− 7

8 n ⁱ n ^j n ^k I _kll ⁽²⁾ + 5

8 n ⁱ n ^j n ^k n ^l n ^m I _klm ⁽²⁾ + 11

24 δ ij n ^k I _kll ⁽²⁾ − 5

8 δ ij n ^k n ^l n ^m I _klm ⁽²⁾ o + n 2

3 n ^k ( ˙ S ikj + ˙ S jki ) − 4

3 (n ⁱ S ˙ jkk + n ^j S ˙ ikk ) + 2n ^k n ^l (n ⁱ S ˙ jkl + n ^j S ˙ ikl ) + 2n ⁱ n ^j n ^k S ˙ kll + 2

3 δ ij n ^k S ˙ kll

oi

+ O(r ⁻³ ). (3.21) It should be noted that ₍₄₎ h ij obtained in this way becomes meaningless at the far zone because Eq. (2.14), from which ₍₄₎ h _ij is derived, is valid only in the near zone.

4. Summary. We have formulated the hydrodynamic equation accurate up to 2.5PN order. We carefully consider methods to obtain various metric quantities, especially the 2PN tensor potential ₍₄₎ h ij for the sake of an actual numerical simulation. It is possible to solve the Poisson equations for ₍₄₎ h _ij by using standard numerical methods. Thus, the formalism developed here will be useful also in numerical calculations.

Moreover, we would like to emphasize that, from the 2PN order, the tensor part of the

3-metric, ˜ γ _ij , cannot be neglected even if we ignore gravitational waves. Recently, Wilson

and Mathews [10] presented numerical equilibrium configurations of binary neutron stars

using a semi-relativistic approximation, in which they assume the spatially conformal flat metric as the spatial 3-metric, i.e., ˜ γ ij = δ ij . Thus, in their method, a 2PN term, h ij , was completely ignored. This means that their results unavoidably have an error of the 2PN order which will become ∼ (M/R) ² ∼ 1 − 10%. Thus if we want to obtain a general relativistic eqiulibrium configuration of binary neutron stars with a better accuracy (say less than 1%), we should take into account the tensor part of the 3-metric.

Acknowledgments. The author thanks T. Futamase, T. Nakamura, M. Sasaki, M. Shibata and T. Tanaka for invaluable discussions. He is also very grateful to S. Ikeuchi and M. Sasaki for continuous encouragements. He thanks the Stefan Banach International Mathematical Center for its hospitality during the semester. This work is based on a col- laboration with M. Shibata and T. Futamase [8].

References

[1] A. A b r a m o v i c i et al., Science, 256 (1992), 325;

K. S. T h o r n e, in proceedings of the eighth Nishinomiya-Yukawa memorial symposium on Relativistic Cosmology, edited by M. Sasaki (Universal Academy Press, Tokyo, 1994), pp.67.

[2] C. B r a d a s c h i a et al., Nucl. Instrum. Method Phys. Res. Sect. A289 (1990), 518.

[3] D. L a i, F. R a s i o and S. L. S h a p i r o, Astrophys. J. 420 (1994), 811; Astrophys. J.

supplement 88 (1993), 205.

[4] T. N a k a m u r a, in the proceedings of the eighth Nishinomiya-Yukawa memorial sympo- sium on Relativistic Cosmology (ref.1), pp. 155.

[5] S. C h a n d r a s e k h a r and F. P. E s p o s i t o, Astrophys. J. 160 (1970), 153; see their series of papers cited therein.

[6] L. B l a n c h e t, T. D a m o u r and G. S c h ¨ a f e r, Mon. Not. R. Astr. Soc. 242 (1990), 289.

[7] K. O o h a r a and T. N a k a m u r a, Prog. Theor. Phys. 83 (1990), 906; 86 (1991), 73;

88 (1992), 307;

M. S h i b a t a, T. N a k a m u r a and K. O o h a r a, Prog. Theor. Phys. 88 (1992), 1079;

89 (1993), 809.

[8] H. A s a d a, M. S h i b a t a and T. F u t a m a s e, Prog. Theor. Phys. 96 (1996), 81.

[9] T. O h t a, H. O k a m u r a, T. K i m u r a and K. H i i d a, Prog. Thor. Phys. 51 (1974), 1598.

[10] J. R. W i l s o n and G. J. M a t h e w s, Phys. Rev. Lett. 75 (1995), 4161.