Harmonic coordinate method for simulating generic singularities

6 Pages·0.13 MB·English

by  David Garfinkle

Harmonic coordinate method for simulating generic singularities
David Garfinkle
Department of Physics, Oakland University, Rochester, Michigan 48309
This paper presents both a numerical method for general relativity and an application of that
method. Themethodinvolvestheuseofharmoniccoordinates ina3+1codetoevolvetheEinstein
equations with scalar field matter. In such coordinates, the terms in Einstein’s equations with the
highest number of derivatives take a form similar to that of the wave equation. The application is
an exploration of the generic approach to the singularity for this type of matter. The preliminary
results indicate that the dynamics as one approaches the singularity is locally the dynamics of the
Kasner spacetimes.
04.25.Dm, 04.20.Dw
2
0
0 I. INTRODUCTION coordinate time, the whole spacetime up to the singu-
2 larity is covered and the behavior of the metric as the
n It has long been knownthat harmonic coordinatesare singularityisapproachedsimplybecomesasymptoticbe-
a useful for mathematical relativity. In particular, these havior in the limit of large time coordinate. For Gowdy
J coordinates were used [1] to prove local existence of so- spacetimes there is a natural choice of such a time co-

9 lutions ofthe vacuumEinsteinequation. This usefulness ordinate: these spacetimes are foliated by T2s invariant
of harmonic coordinates stems from their putting Ein- under the symmetry group. The area of the symmetry
2 stein’s equation into a form that is similar to the curved T2sgoestozeroasthe singularityisapproached. There-
v
spacetime wave equation. Since many numerical tech- fore minus the logarithm of this area is a natural time
3
1 niquesworkwellonthewaveequation,onemightexpect coordinate that goes to infinity as the singularity is ap-
0 that harmonic coordinates would be used extensively in proached. While this method works well for the Gowdy
0 numerical relativity, and it is somewhat surprising that spacetimes, since it depends crucially on the symmetry
1 they are not. (However see [2] for some recent mathe- of the Gowdy spacetimes, the method does not seem to
1 matical and numerical work on the linearized case. In generalize to the case of the generic singularity with no
0 addition, harmonic time slices have been advocated and symmetries. A natural generalization comes when one
/
c used in numerical relativity [3–7]). This is perhaps due notices that this time coordinate (minus the logarithm
q tothefollowingdrawbackofharmoniccoordinates: these ofthe areaofthe symmetryT2s)isalsoharmonic. Since
r- coordinates are solutions of the wave equation, and such in some sense one expects the wave equation to become
g solutions need not have a timelike gradient at all points singular as the spacetime singularity is approached, one
v: ofspacetime,eveniftheystartoutwithatimelikegradi- mightalsoexpectasolutionofthewaveequationtoblow
i entonaninitialdatasurface. Whatthismeansisthatin up as the singularity is approached and thus one might
X
harmoniccoordinates,thetimecoordinatewillingeneral want to use such a solution as the coordinate time.
ar not remain timelike and this is likely to cause numerical Whatisexpectedtobethegenericbehaviorofaspace-
problems.[8,9]Aswewillseelater,thereisawayaround time as the singularity is approached? Based on studies
this difficulty. of spacetimes with T2 symmetry [10–12] and spacetimes
One areaofnumericalrelativity where harmoniccoor- with U(1) symmetry [13], the expected answer [14] is
dinates have been used (though somewhat unintention- the following: the singularity is expected to be spacelike
ally)isinthe studyofthe approachtothe singularity,in and as it is approached each spatial point is expected to
particular in the Gowdy spacetimes. [10]Numerical sim- “decouple”fromthe othersandundergoadynamicscor-
ulation of approach to a spacetime singularity presents responding to that of a homogeneous spacetime (though
problems of its own. Since it is expected that various adifferenthomogeneousspacetimeateachspatialpoint).
quantities become infinite atthe singularity,the numeri- Which types of homogeneous spacetime is this dynamics
calsimulationwillgenerallystopafterafinitecoordinate expected to correspond to? For vacuum spacetimes, it
time, the time when a surface of constant time first en- is thoughtthatthe dynamics willbe oscillatory,possibly
counters the singularity. Since this first encounter gen- correspondingtotheMixmasterspacetimeasconjectured
erally occurs at one spatial point, the information about in[15]. Formanyothertypesofmatteritisexpectedthat
the behavior of the singularity at other spatial points as the singularity is approached the matter terms in the
will be unavailable from the numerical simulation. The Einstein equations become negligible and the dynamics
solution to this difficulty is to choose a time coordinate approaches that of a vacuum spacetime.
that tends to infinity as the singularity is approached. However for so called “stiff matter” (i.e. a scalar field
In this way, the simulation is not forced to end at finite or a perfect fluid with equation of state P = ρ) it is ex-
1
pected that the dynamics will notbe oscillatoryandwill well on Einstein’s equations in this form. The reasonfor
correspondto thatofa Kasnerspacetime. This expecta- considering nonzero source terms Hµ in Eq. (5) is that
tion is greatly bolstered by a theorem due to Andersson these terms allow us to change the behavior of the time
and Rendall [16] which shows local existence in a neigh- coordinateandthusmayallowustoeliminate(oratleast
borhood of the singularity of solutions to the Einstein- postpone)thebehaviorwherethetime coordinateceases
scalar equations, with the expected asymptotic behavior to be timelike. In the simulations done in this paper, no
and with enough degrees of freedom to be the generic source terms were needed. I know of no systematic way
solutions. to find appropriate source terms and expect that in the
Sincethe approachto thesingularityisexpectedtobe cases where they are needed, one will have to resortto a
simpler for stiff matter than for vacuum, the stiff matter trialanderrormethodtofindappropriateHµ. Notethat
caseshouldbe easiserto treatnumerically. Therefore,in theuseofspatialharmoniccoordinatescouldalsoleadto
this paper we confine ourselves to a numerical study of coordinate problems. This would occur if the gradients
the approach to the generic singularity in the Einstein- ofthe four coordinatesfailto be linearlyindependent. If
scalar system. Section 2 presents the equations and nu- this sort of problem occurs, one might expect to be able
mericalmethodsused. Theresultsaregiveninsection3. to postpone or eliminate it by using appropriate source
Section4containsadiscussionoftheresultsandofother terms.
possibleapplicationsoftheharmoniccoordinatemethod. Thenumericalmethodusedrequiresequationsthatare
first order in time. To put the equations in such a form,
we define quantities P and P given by
αβ φ
II. EQUATIONS AND NUMERICAL METHODS
P =∂ g (7)
αβ t αβ
The equations that we wish to evolve numerically are
the Einstein-scalar equations P =∂ φ (8)
φ t
R =8π φ φ (1) Then the Einstein-scalar equation becomes
αβ α β
∇ ∇
Here,weusetheconventionsof[17]includingunitswhere g00∂tPαβ =2g0k∂kPαβ +gik∂i∂kgαβ
c=G=1. As a consequence of Eq. (1) and the Bianchi − +16π∂αφ∂βφ+2∂(αHβ)−2ΓγαβHγ
identities, the scalarfieldmustsatisfy the waveequation C µνC C µνC +2Γγ Γν (9)
− α µνβ − β µνα να γβ
αφ=0 (2) The wave equation for φ becomes
α
∇ ∇
The Riccitensorisgivenintermsofthe Christoffelsym- −g00∂tPφ =2g0k∂kPφ+gik∂i∂kφ−gαβΓγαβ∂γφ (10)
bols by
Thefullsetofequationsthatareevolvedinthecomputer
R =∂ Γγ ∂ Γγ +Γγ Γν Γγ Γν (3) code are Eqs. (7-10).
αβ γ αβ − α γβ αβ νγ − αν βγ We now turn to the numericalmethods used to evolve
while the Christoffel symbols are given in terms of the theseequations. Spatialderivativesareapproximatedby
metric by centereddifferences. Thevariablesareevolvedintimeus-
ingathreestepiteratedCrank-Nicholson(ICN)method.
Γγαβ = 12gγδ(∂αgβδ+∂βgαδ−∂δgαβ) (4) [18–20]This works as follows: evolutionequations of the
form∂ S =W(S)forsomesetofvariablesS areapprox-
t
Harmonic coordinates are solutions of the wave equa-
imated as
tion. As a generalization of harmonic coordinates, con-
sider coordinates that satisfy the wave equation with Sn+1 =Sn+ ∆t W(Sn)+W(Sn+1) (11)
source 2 (cid:2) (cid:3)
where Sn is the value of S at time step n and ∆t is the
α xµ =Hµ (5)
∇ ∇α time step. Then, using Sn as an initial guess for Sn+1,
where Hµ are specified from the beginning. Then using Eq. (11) is iterated three times.
Eqs. (3,4,5) we find that the Ricci tensor is given by The spacetimes we consider have topology T3 R.
Each spatial slice has topology T3. In terms of the×spa-
R = 1gγσ∂ ∂ g + 1C µνC tial coordinates, this means that 0 x 2π with 0 and
αβ −2 γ σ αβ 2 β µνα ≤ ≤
+1C µνC Γγ Γν ∂ H +Γγ H (6) 2π identified (and correspondingly for y and z). This
2 α µνβ − να γβ− (α β) αβ γ topology is implemented numerically as follows: a spa-
where C ∂ g . Note that the second derivative tial coordinate x has N grid points. The variables on
αµν α µν
≡
terms appear only in the wave operator. Therefore, one points from 2 to N 1 are evolved using the evolution
−
might expect that Einstein’s equations in this form be- equations. The variables on point 1 are set to the values
have similarly to the wave equation and that numerical atpointN 1whileatpointN theyaresettothevalues
−
methodsthatworkwellonthewaveequationmightwork at point 2.
2
III. RESULTS code, 3 gridpoints were used in the x direction, 3 grid-
points in the y direction and 500 gridpoints in the z
All runs were done in double precision on Compaq direction. For the comparison the initial data used is
XP1000 workstations and on the NCSA Origin 2000. P = 0, ∂τP = 5cosz, Q = cosz, ∂τQ = 0. These data
The time step was ∆t = ∆x/2. While the source terms are evolved until τ = π and the results for the compar-
Hµ shouldbe helpful in keeping the coordinateswellbe- ison are given in Fig. 1. Here, the solid line represents
haved, we did not need them for the cases studied here the1+1evolutionoftheGowdyequations,whilethedots
and all runs have Hµ =0. representthefull3+1evolutionusingtheharmoniccode.
There is clearly agreement between the two.
P
C
6
0.0035
0.003
4
0.0025
0.002
2
0.0015
z 0.001
0.0005
-2 t
0.5 1 1.5 2 2.5 3 3.5
FIG. 1. comparison of the quantity P as found using a
Gowdy code and the3+1 harmonic code
FIG.2. convergence test involving theconstraint
Before exploring the generic singularity,we would like Before presenting the results of another code test, we
to test the code. One way of doing this is to run the turn to the question of finding initial data without sym-
code on cases with symmetry where results are already metries. On an initial data surface, the intrinsic metric
known. Such a case is the Gowdy spacetimes. These h andthe extrinsic curvatureK mustsatisfy the con-
ij ij
have the form [10] straint equations
ds2 =e(τ−λ)/2 e−2τdτ2+dz2 DiKij DjK =8πφ˙Djφ (15)
+e−τ ePdx2+2ePQdxdy+ e(cid:2)P−Q2+e−P dy2(cid:3) (12) −
(cid:2) (cid:0) (cid:1) (cid:3)
(3)R+K2 KijK =8π φ˙2+DiφD φ (16)
Here,P,Qandλarefunctionsofτ andz. Itfollowsfrom ij h i i
−
Eq. (12) that the coordinates (τ,x,y,z) are all solutions
ofthewaveequation. Therefore,theyareharmoniccoor- (Note that we use the convention of [17] for the sign of
dinates. The vacuum Einstein equations for the Gowdy Kij which is opposite to the convention usually used in
spacetimes are [11] numerical relativity). Here, an overdot denotes deriva-
tive along the normal to the surface, and D and (3)R
i
∂ ∂ P e−2τ∂ ∂ P e2P (∂ Q)2 e−2τ(∂ Q)2 =0 are respectively the covariant derivative and scalar cur-
τ τ − z z − h τ − z i vature of h . Given a solution of Eqs. (15) and (16),
ij
(13) we produce initial data for evolution in harmonic coor-
dinates as follows: For spatialdirections i and j we have
g = h , g = g = 0, g = 1, P = 2K . The
ij ij i0 0i 00 ij ij
∂ ∂ Q e−2τ∂ ∂ Q+2 ∂ P∂ Q e−2τ∂ P∂ Q =0 remaining components of P are−solved for using Eq.
τ τ z z τ τ z z αβ
− (cid:0) − (cid:1)
(5).
(14)
WewanttofindasolutionofEqs. (15)and(16)thatis
simplebuthasnosymmetriesandhassomefreeparame-
plus constraint equations that determine λ once P and
ters. Wechooseφ=0andh equaltotheflatEuclidean
Q are known. ij
metric in the usual coordinates. For the extrinsic curva-
To test the 3+1 code, we evolve Eqs. (13) and (14)
ture we choose
using a 1+1 code. Then we evolve the same initial
data with the 3+1 harmonic code and compare the re-
K =(b +a cosy+a cosz)/2
xx 1 2 3
sults. To do the comparison, note that P = τ +lng
xx K =(b +a cosx a cosz)/2
yy 2 1 3
so that it is straightforward to compare the values of −
K =(b a cosx a cosy)/2
zz 3 1 2
P produced by the two codes. For these simulations − −
K =K =(c cosz)/2
xy yx 1
500 gridpoints were used in the 1+1 code. In the 3+1
K =K =(c cosy)/2
xz zx 2
3
K =K =(c cosx)/2 (17)
yz zy 3
2
Here, the quantities a , b and c are constants that are f
i i i
free parameters. It is straightforward to show that the 0
ln g_xx
extrinsic curvature of Eq. (17) with our choice of initial
h and φ satisfies Eq. (15). Equation (16) then be- -2
ij
comes an agebraic equation for P which can be solved
φ
ln g_yy
provided that the left hand side of the equation is pos- -4
itive. Note that if all the a and c are zero and all the ln g_zz
i i
bi are equal, then the initial data evolve to the spatially -6
flatRobertson-Walkerspacetimewithscalarfieldmatter.
Thus,thisfamilyofdatacanbethoughtofasRobertson-
-8
Walker with large gravitationaland scalar waves. ln (-g_tt)
We now consider a convergence test involving a con-
-10
straintthatcomesfromtheuseofharmoniccoordinates. 0 2 4 6 8
Define the quantities Cµ by t
FIG.3. behavior of metric components and scalar field as
thesingularity is approached at thespatial point (0,0,0)
Cµ =gαβΓµ (18)
αβ
(this is the appropriate form of the constraint for the
casewherethesourcetermHµ vanishes. Forthe general 2
case, the constraint would be given by Cµ = gαβΓµ + f
αβ
Hµ). Then from Eq. (5) it follows that Cµ = 0. Since 0
ln g_xx
we are solvingthe evolutionequations by approximating
them by finite difference equations, the quantities Cµ as -2
evaluatedby the computer code will notbe zerobecause
of errors due to the finite grid spacing ∆x. Define the -4 ln g_yy
quantity C by ln g_zz
-6
√ggαβCαCβdxdydz 1/2
C =(cid:12)R (cid:12) (19)
(cid:12) √gdxdydz (cid:12) -8
(cid:12) R (cid:12) ln (-g_tt)
(cid:12) (cid:12)
This quantity is a type of measure of the average size
-10
of the constraint. Figure 2 shows a plot of C vs. 0 2 4 6 8
t
time. The parameters used are a = (0.1,0.1,0.2), b = FIG.4. behavior of metric components and scalar field as
i i
(0, 0.5, 0.5) and c = (0,0,0). Here the curve corre- thesingularityisapproachedatthespatialpoint(0,π/4,π/2)
i
− −
sponds to a run done with 20 gridpoints in each spatial
direction, while the dots correspond to a run with 38
gridpoints in each spatial direction (which gives half the
is a linear function of time. Thus, if the behavior of a
grid spacing) and with C multiplied by 4. The results
genericsolutionnearthesingularityislocalKasner,then
show second order convergence.
we should expect metric components that are exponen-
We now consider the approach to the singularity. To
tial functions of time, with the exponent depending on
see what is expected, it is helpful to consider the Kasner
space. Weshouldalsoexpectascalarfieldthatisalinear
spacetime in harmonic coordinates. This is given by
function of time with the slope depending on space.
ds2 = e(q1+q2+q3)τ dτ2+eq1τdx2+eq2τdy2+eq3τdz2 Results on the approach to the singularity are shown
− inFigs. 3and4. Therunwasdonewith34gridpointsin
(20) eachspatial directions. The parametersa , b and c are
i i i
thesameasfortheconvergencetest. Herethescalarfield
where the q are constants This metric is generally a so-
i and the logarithms of the diagonal metric components
lution of the Einstein-scalar equations; but for the case
are plotted as functions of time. Figure 3 corresponds
where ( q )2 = q2 it is a vacuum spacetime. (In
P i P i to the spatial point (x,y,z) = (0,0,0) while for Fig. 4
the usual treatment of vacuum Kasner spacetimes, one
the point is (0,π/4,π/2). Note that as the singularity is
defines the quantities p = q /( q ) and then obtains
i i m approached these quantities all become linear functions
the condition p = p2 = 1P). Note that in the vac-
P i P i of time. The differences between Figs. 3 and 4 show
uum case one cannot have all three directions contract-
thatthere isa spatialdependence ofthe approachto the
ing but that this is possible for the Einstein-scalar case.
singularity. Note from Fig. 3 that though the x direc-
Note also that the metric components are exponential
tionis initially expanding,eventuallyallthree directions
functions of time. It also turns out that the scalar field
4
contract. Thisiswhatonewouldexpectifthemetricsof ds2 = α2dt2+h (dxi+βidt)(dxj +βjdt) (21)
ij
−
[16] represent the generic behavior near the singularity
Einstein’sequationsarewrittenasanevolutionequation
since these metrics have all three directions contracting
for the spatial metric h (and the extrinsic curvature
in a neighborhood of the singularity. ij
K ) while the gauge choice results in equations for the
ij
lapse α and the shift βi. This framework is sufficiently
general to accomodate the use of harmonic coordinates,
IV. DISCUSSION
whichcorrespondtothefollowingequationsforlapseand
shift
While this study is somewhat preliminary, it indicates
thatharmoniccoordinatescanbeausefultoolinnumer- ∂ α βi∂ α=Kα2 (22)
t i
−
ical relativity. Though, in principle coordinate problems
could occur, this did not happen in the cases studied ∂ βi βm∂ βi =hmnΓ¯i α2 (23)
here, even though they involved very strong fields. Fur- t − m mn
thermore, the use of the source terms Hµ may cure such (NotethatthesignofK inequation(22)comesfromthe
problems if they arise. convention of [17]). Here Γ¯imn is the Christoffel symbol
Asforthebehaviorofgenericsingularities,thenumer- associated with hij. (Equations (22) and (23) hold for
ical results indicate that solutions of the form proved in the case of vanishing source term Hµ. Similar equations
reference [16] to exist in a neighborhood of the singular- hold in the case of nonzero Hµ). The use of equations
ity also exist globally. Thus such solutions are likely to (22) and (23) in an ADM code is not precisely equiva-
describe the generic singularity in the stiff matter case. lent to the approach of this paper. The reason for this
Thereareseveralprojectsforwhichthemethodsofthis is that equations (22) and (23) directly solve the con-
papercouldbeused. Oneistodoamoreextensivestudy straint gαβΓµαβ = 0, while in our approach, this con-
ofthesingularityintheEinstein-scalarcase,withamore straint is used to change the form of the evolution equa-
thorough exploration of the evolution corresponding to tions. Nonetheless, it would be interesting to use equa-
variousvaluesoftheparametersintheinitialdataofthe tions (22) and (23) in an ADM code to see how that
previous section. It would also be helpful to evolve for a compares with other choices of lapse and shift. In par-
longer time. ticular, one might expect better stability properties and
Anotherprojectistoremovethescalarfieldandstudy compatibilitywithasimpleoutgoingwaveboundarycon-
the approach to the singularity of the generic vacuum dition. (Though note that such improvements can also
spacetime. Here, the behavior is expected to be more be obtained using the BSSN approach [23–25]).
complicated and a treatment will probably require more Anotherdesiredfeatureofanumericalcodeistheabil-
spatial resolution to resolve the expected sharp features, itytotreatblackholes,andthissometimesrequiresblack
as well as longer evolution in time to see the expected holeexcision. Whileharmoniccoordinatesaresingularity
oscillatory behavior. avoiding [3] they are just barely so and come arbitrarily
Yetanotherprojectistostudy thebehaviorofasymp- close to a singularity. Thus the need for excision in an
totically flat spacetimes rather than closed cosmologies. approach that uses harmonic coordinates should be at
Here,the closedcosmologieswerestudiedpartly forsim- least as great as in the standard approach. Nonetheless,
plicity. The periodic boundary conditions are simple one might hope that excision itself would be easier to
to implement and completely consistent with Einstein’s implement using the approach of this paper. This is be-
equations. In contrast, for an outer boundary at a finite causenoellipticequationsareinvolvedandthelightcone
distance in an asymptotically flat spacetime, one needs of the wave operator is the same as that of the physical
to put some sort of outgoing wave boundary condition. metric.
SuchconditionsareusuallynotconsistentwithEinstein’s All these projects are work in progress and prelimi-
equation(it isknownhowto haveaconsistentboundary nary results from them are promising. Thus, I expect
condition [21] but this condition is quite complicated). that harmonic coordinates will become a useful tool in
These inconsistent conditions may lead to numerical in- numerical relativity.
stability. Since harmonic coordinates make Einstein’s
equation look like the wave equation, simple outgoing
V. ACKNOWLEDGEMENTS
wave boundary conditions that work numerically with
the waveequationmightbe expectedto “work”(atleast
inthesenseofnotcausingnumericalinstability)forEin- IwouldliketothankHelmutFriedrich,BeverlyBerger,
stein’s equation. Matt Choptuik, G. Comer Duncan, Lars Andersson and
Since much work has been done on numerical simula- Alan Rendall for helpful discussions. I would also like
tionsofasymptoticallyflatspacetimesusingthestandard to thank the Albert Einstein Institute for hospitality.
Arnowitt-Deser-Misner (ADM) [22] approach, it is help- This work was partially supported by NSF grant PHY-
fultomakecomparisonswiththisapproach. IntheADM 9988790 to Oakland University. Some of the computa-
approach, the spacetime metric is written as tions were performed at the National Center for Super-
computing Applications (Illinois).
5
[1] Y.Foures-Bruhat, Acta Mathematica, 88, 141 (1952)
[2] B. Szilagyi, B. Schmidt and J. Winicour, gr-qc/0106026
[3] C. Bona and J. Masso, Phys.Rev. D38, 2419 (1988)
[4] C. Bona and J. Masso, Phys.Rev. D40, 1022 (1989)
[5] C. Bona and J. Masso, Phys.Rev. D68, 1097 (1992)
[6] A. Abrahams, A. Anderson, Y. Choquet-Bruhat and J.
York,Phys.Rev.Lett. 75, 3337 (1995)
[7] P. Anninos, K. Camarda, J. Masso, E. Seidel, W. Suen
and J. Towns, Phys. Rev. D52, 2059 (1995)
[8] M. Alcubierre, Phys. Rev. D55, 5981 (1997)
[9] M. Alcubierre and J. Masso, Phys. Rev. D57, 4511
(1998)
[10] B.BergerandV.Moncrief,Phys.Rev.D48,4676(1993)
[11] B.BergerandD.Garfinkle,Phys.Rev.D57,4767(1998)
[12] B. Berger, J. Isenberg and M. Weaver, Phys.Rev. D64,
084006 (2001)
[13] B. Berger and V. Moncrief, Phys. Rev. D58, 064023
(1998)
[14] B.Berger, D.Garfinkle,J.Isenberg,V.MoncriefandM.
Weaver,Mod. Phys.Lett. A13, 1565 (1998)
[15] V. Belinskii, I. Khalatnikov and E. Lifshitz, Sov. Phys.
Usp.13, 745 (1971)
[16] L. Andersson and A. Rendall, Commun. Math. Phys.
218, 479 (2001)
[17] R.M. Wald, General Relativity (University of Chicago
Press, Chicago, 1984)
[18] M.Choptuik,inDeterministicChaosinGeneralRelativ-
ity D. Hobill, A. Burd and A. Coley eds. (Plenum, New
York,1994) pp.155-175
[19] S.Teukolsky,Phys. Rev. D61, 087501 (2000)
[20] M. Alcubierre et al, Phys.Rev.D62, 044034 (2000)
[21] H. Friedrich and G. Nagy: Commun. Math. Phys. 201,
619 (1999)
[22] R.Arnowitt, S.DeserandC. Misner, in Gravitation: an
introduction to Current Research, L. Witten ed. (John
Wiley,New York,1962) pp. 254-280
[23] T. Nakamura, K. Oohara and Y. Kojima, Prog. Theor.
Phys.Suppl.90, 1 (1987)
[24] M. Shibata and T. Nakamura, Phys. Rev. D52, 5428
(1995)
[25] T. Baumgarte and S. Shapiro, Phys. Rev. D59, 024007
(1999)
6