Regularization of spherically symmetric evolution codes in numerical relativity

0.17 MB·English

by  Miguel Alcubierre

Regularization of spherically symmetric evolution codes in numerical relativity
Miguel Alcubierre∗ and Jos´e A. Gonz´alez†
Instituto de Ciencias Nucleares, Universidad Nacional Aut´onoma de M´exico, A.P. 70-543, M´exico D.F. 04510, M´exico.
(Dated: February 6, 2008)
The lack of regularity of geometric variables at the origin is often a source of serious problem
for spherically symmetric evolution codes in numerical relativity. One usually deals with this by
restricting the gauge and solving the hamiltonian constraint for the metric. Here we present a
generic algorithm for dealing with the regularization of the origin that can be used directly on the
evolution equations and that allows very general gauge choices. Our approach is similar in spirit
to the one introduced by Arbona and Bona for the particular case of theBona-Masso formulation.
However,our algorithm is more general and can be used with a wide variety of evolution systems.
PACSnumbers: 04.20.Ex,04.25.Dm,95.30.Sf,
4
0 I. INTRODUCTION The main drawback of the standard approach is that
0
the gaugechoicehas beencompletely exhausted. Inpar-
2
ticular,thepolar/arealgaugecannotpenetrateapparent
When developing spherically symmetric codes in nu-
n horizons,sinceinsideanapparenthorizonitisimpossible
mericalrelativity,thecoordinatesingularityattheorigin
a
to keep the areas of spheres fixed without a non-trivial
J can be a source of serious problems caused by the lack
shift vector. The polar/areal gauge has, nevertheless,
8 of regularity of the geometric variables there. The prob-
been used successfully even in the study of critical col-
2 lem arises because of the presence of terms in the evolu-
lapsetoablackhole,wherethepresenceoftheblackhole
tionequations that goas 1/r near the origin. Regularity
1 of the metric (essentially local flatness) guarantees the is identified by the familiar “collapse of the lapse” even
v if no apparent horizon can be found [3].
exact cancellations of such terms at the origin, thus en-
3 Still, one would like to have a way of dealing with the
suring well-behaved solutions. This exact cancellation,
1 regularityissuethatallowsmoregenericgaugechoicesto
however, though certainly true for analytical solutions,
1
be made,eitherbecause oneisinterestedinstudying the
1 usually fails to hold for numerical solutions. One then
region inside an apparent horizon, or because one wants
0 finds thatthe1/r termsdonotcancelandthe numerical
4 solutionbecomes ill-behavednear r =0: itnot only fails to test interesting gauge conditions in the simple case of
0 to convergethere,but caneasilyturn outto be violently sphericalsymmetry. Becauseofthiswehavedevelopeda
/ generalregularizationtechniquethatcanbeuseddirectly
c unstable in just a few timesteps.
q on the Einstein evolution equations.
The usual way to deal with this problem is to use the
- Ourregularizationmethodissimilarinspirit,ifnotin
r so-called areal (or radial) gauge, where the radial coor-
g detail, to the one presented by Arbona and Bona in [4].
dinate r is chosen in such a way that the proper area of
v: spheres of constant r is always 4πr2. If, moreover, one The main difference being that the approach of Arbona
andBonawastiedtotheuseoftheBona-Massoevolution
i also choses a vanishing shift vector one ends up in the
X system [5, 6, 7, 8, 9], while our algorithm is much more
standard polar/areal gauge [1, 2], for which the lapse is
r general.
forcedtosatisfyacertainordinarydifferentialequationin
a
A final point deserves notice. Spherically symmetric
r. Thenamepolarcomesfromthefactthatforthisgauge
evolution codes that involve eternal black holes usually
choicethereisonlyonenon-zerocomponentoftheextrin-
ignore the regularity problem. For example, one can ex-
sic curvature tensor, namely K [1]. In the polar/areal
rr
cise the black hole interior and eliminate r =0 from the
gauge the problem of achieving the exact cancellation of
numerical grid. But even if one does not use excision,
the 1/r terms is reduced to imposing the boundary con-
for a black hole r=0 is not a regular point but rather a
dition g = 1 at r = 0, which can be easily done if one
rr
compactification of the asymptotic infinity on the other
solves for g from the hamiltonian constraint (which in
rr
side of the Einstein-Rosen bridge. This compactifica-
this caseisanordinarydifferentialequationinr)andig-
tion introduces geometric factors that compensate the
nores its evolution equation. If one does this in vacuum,
1/r terms making the equations regular even at r = 0.
one ends up inevitably with Minkowski spacetime in the
Thismeansthat,contrarytowhatonewouldnaivelyex-
usualcoordinates(one canalsorecoverSchwarzschildby
pect, in spherical symmetry it is easier to evolve eternal
working in isotropic coordinates and factoring out the
black holes (at least for some time) than it is to evolve
conformalfactoranalytically). Ofcourse,inthepresence
regular spacetimes.
of matter, one can still have truly non-trivial dynamics.
This paper is organized as follows. In Sec. II we dis-
cuss the regularity conditions that the metric functions
must satisfy at the origin of spherical coordinates, and
∗Electronicaddress: [email protected] we show which terms need to be regularized in the Ein-
†Electronicaddress: [email protected] stein equations. Section III describes our regularization
2
algorithm in a generic way. In Sec. IV we present exam- ∂ D = −2α[K D +∂ K ], (2.12)
t B B α r B
ple ofhowtoregularizesomespecific formulationsofthe α D D
2 α A
Einstein equations, and show some numerical examples. ∂tKA = −A ∂r(Dα+DB)+Dα− 2
We conclude in section V. (cid:20)
D2 D D
+ B − A B −AK (K +2K )
A A B
2 2
II. REGULARITY CONDITIONS 1
− (D −2D ) , (2.13)
A B
r
(cid:21)
We start by writing the general form of the spatial α D D
2 A B
∂ K = − ∂ D +D D +D −
metric in spherical symmetry as t B 2A r B α B B 2
(cid:20)
2 2 2 2 1 2(A−B)
dl =A(r,t)dr +r B(r,t)dΩ , (2.1) − (D −2D −4D )−
r A α B r2B
(cid:21)
2
withAandB positivemetricfunctionsanddΩ thesolid + αK (K +2K ) (2.14)
angle element: dΩ2 = dθ2 +sin2θdφ2. Notice that we B A B
2
have already factored out the r dependency of the an- where α is the lapse function and D := ∂ lnα. The
α r
gularmetricfunctions. Thishastheadvantageofmaking hamiltonian and momentum constraints take the form
explicit the dependency onr ofgeometric quantities and
1
makes the regularizationprocedure easier. ∂ D = (A−B)+AK (2K +K )
As we will deal with the Einstein equations in first r B r2B B A B
2
order form, we will introduce the auxiliary quantities: + 1(D −3D )+ DADB − 3DB ,(2.15)
A B
r 2 4
DA :=∂rlnA, DB :=∂rlnB . (2.2) 1 DB
∂ K = (K −K ) + , (2.16)
r B A B
r 2
We will also work with the mixed components of the ex- (cid:20) (cid:21)
trinsic curvature: KA :=Krr,KB :=Kθθ =Kφφ. Since {Dα,DA,DB} go as r near the origin, terms of
There are in fact two different types of regularitycon- the type D /r are regular and represent no prob-
{α,A,B}
ditions that the variables {A,B,DA,DB,KA,KB} must lem. However, we see that both in the hamiltonian con-
satisfy at r = 0. The first type of conditions are simple straint and in the evolution equation for K there is a
B
thoseimposedbytherequirementthatthedifferentvari- termoftheform(A−B)/r2,whileinthemomentumcon-
ables should be well defined at the origin, and imply the straint there is a term of the form (K −K )/r. Given
A B
following behavior for small r: the behavior of these variables near r = 0 these terms
would seem to blow up at the origin. The reason why
0 2
A ∼ A +O(r ), (2.3) this does not in fact happen is that, near the origin, we
0 2
B ∼ B +O(r ), (2.4) must also ask for the extra regularity conditions
D ∼ O(r), (2.5)
A 2 2
A−B ∼O(r ), K −K ∼O(r ), (2.17)
A B
D ∼ O(r), (2.6)
B
0 2
K ∼ K +O(r ), (2.7) that is
A A
0 2
K ∼ K +O(r ), (2.8)
B B 0 0 0 0
A =B , K =K . (2.18)
A B
with {A0,B0,K0,K0} perhaps functions of time, but
A B Itisnotdifficultto understandwherethese conditions
not of r. These regularity conditions are in fact quite
come from. They are just a consequence of the fact that
easy to implement numerically. For example, one can
space must remain locally flat at r = 0. This local flat-
useafinitedifferencinggridthatstaggerstheorigin,and
ness conditionimplies that, near r=0, itmust be possi-
then obtain data on the fictitious point at r = −∆r/2
ble to write the metric as
by demanding for {A,B,K ,K } to be even functions
A B
atIrt=is0thaendsefcoorn{dDtAy,pDeBo}f rteogbuelaoridtdy.conditions that is dl2R∼0 =dR2+R2dΩ2 , (2.19)
more troublesome. In order to see the problem, we
withRaradialcoordinatethatmeasuresproperdistance
will first write the Arnowitt-Deser-Misner (ADM) equa-
from the origin. A local transformation of coordinates
tions [10, 11] in first order form for the case of spherical
from R to r then takes the metric into the form
symmetry. The evolution equations are (in the case of
zero shift) 2
dR
2 2 2 2
dl r∼0 = dr +r dΩ , (2.20)
dr
∂tA = −2αAKA , (2.9) (cid:18) (cid:19)r=0
(cid:0) (cid:1)
∂ B = −2αBK , (2.10) 0 0
t B whichimpliesthatA =B and,sincethismustholdfor
∂ D = −2α[K D +∂ K ], (2.11) all time, also that K0 =K0.
t A A α r A A B
3
Itturns outthat itis nottrivialtoimplement numeri- the problem of having over-determined boundary condi-
callyboththesymmetryregularityconditions(2.3)-(2.8) tions: the extra boundary condition will be imposed on
and the local flatness regularity conditions (2.17) at the the auxiliary variable. We will then define the variable,
same time. The reason for this is that at r = 0 we now
have three boundary conditions for just two variables: 1 A
λ:= 1− . (3.1)
boththederivativesofAandB mustvanish,plus Aand r B
(cid:18) (cid:19)
B mustbe equaltoeachother(andthe samethingmust
Notice that, if the local flatness regularity condi-
happen for K and K ). The boundary conditions for
A B
tions (2.17) are satisfied, then the variable λ has the
the exact equations are, of course, also over-determined,
following behavior at the origin
but in that case the consistency of the equations implies
that if they are satisfied initially they remain satisfied
λ∼O(r), (3.2)
for all time. In the numerical case, however, this is not
trueduetotruncationerrors,andveryrapidly(typically
which,asmentionedabove,caneasilybeimposednumer-
within one or two time steps) one of the three bound-
icallyusingagridthatstaggerstheorigin,andaskingfor
ary conditions fails to hold. It is easy to convince one-
λ to be odd across r =0.
selfthat simply ignoringone conditionandimposing the
The hamiltonian constraint now becomes
other two does not work. If we impose the zero deriva-
tive condition and ignore the A=B condition, then the λ
∂ D = +AK (2K +K )
2 r B B A B
(A−B)/r term in the evolution equations rapidly be- r
comes singular. On the other hand, if we impose the 1 D D 3D2
+ (D −3D )+ A B − B , (3.3)
A=B condition and ignore one of the symmetry condi- r A B 2 4
tions, then we introduce an inconsistency with the finite
difference versionof the evolution equations, since for fi- and the evolution equation for KB becomes
nite∆ritisverydifficulttoguaranteethatthedifference
α D D
between ∂tKA and ∂tKB approaches zero at the origin. ∂tKB = −2A ∂rDB+DαDB+DB2 − A2 B
This inconsistency then very rapidly causes large (and (cid:20)
non-convergent)gradients to develop near the origin. In 1
− (D −2D −4D −2λ)
the followingsectionwe willintroduceanalgorithmthat r A α B
(cid:21)
successfullyregularizesthenumericalevolutionequations
+ αK (K +2K ). (3.4)
B A B
near r =0.
As a final comment, from the above equations we can As mentioned, the problem terms have now been trans-
also very easily see why the polar/arealgauge has no se- formed into λ/r, which will be well behaved as long as
rious regularity problem. In that gauge we have B = 1 λ is odd at r = 0. The momentum constraint still has
by construction. If we now impose the boundary condi- a term (K −K )/r, but this should cause no trouble
A B
tion A(r =0)=1, and solve for A(r) by integrating the since it does not feed back into the evolution equations
hamiltonian constraint with B = 1 and DB = 0 (ignor- (onecanalwaysmultiply the momentumconstraintwith
ing the evolution equations), then the (A−B)/r2 term r before evaluating it). Of course, multiples of the mo-
causes no trouble. mentum constraint are typically added to the evolution
equations in order to build hyperbolic formulations, and
wewilldiscuss howtodealwiththatterminsuchacase
III. REGULARIZATION ALGORITHM below.
There is one other ingredient that needs to be added:
an evolution equation for λ. This can be obtained di-
In Ref. [4], Arbona and Bona developed a regulariza-
rectly from its definition:
tion technique for the spherically symmetric version of
the Bona-Masso (BM) evolution system [5, 6, 7, 8, 9].
2αA K −K
A B
Their technique is based on redefining the auxiliary dy- ∂ λ= . (3.5)
t
B r
namical variable Vr that is part of the standard BM (cid:18) (cid:19)
formulation in a way that allows them to absorb the
The last evolution equation clearly has the dangerous
(A−B)/r2 terms and reduces the regularization prob-
(K −K )/r term, but this term can be removed with
A B
lem to applying the correct boundary condition to V at
r the help of the momentum constraint (2.16) to find
r =0.
Sinceweareinterestedindevelopingagenericregular- 2αA D
B
∂ λ= ∂ K − (K −K ) , (3.6)
izationtechnique, we willstartfromthe ADM equations t r B A B
B 2
fromtheprevioussection. However,wewilltaketheidea (cid:20) (cid:21)
fromtheregularizationtechniqueofArbonaandBonaof which is now regular at the origin.
introducinganauxiliaryvariablethatwillallowustoab- The regularized first order ADM evolution equations
sorbtheproblematicterms. Addinganauxiliaryvariable are then (2.9)-(2.13), with (2.14) replaced by (3.4), plus
seemstoustobethemorestraightforwardwayofsolving the evolution equation for λ given by (3.6).
4
Having regularized the standard ADM equations, the The finalstep is to substitute the spatialderivative ofv
j
question arises of how to regularize alternative formula- for that of v′ to find
j
tions where multiples of the constraints can be added to
the evolution equations in a number of ways. Adding ∂ v′ = Mj(u)∂ v′ +p′(u,v)−λFt(u,v)
t i i r j i i
multiples of the hamiltonian constraint represents no
f (u)B
i
problem, as the introduction of the variable λ already + ∂ λ
r
2αA
regularized this constraint, as seen in (3.3). The mo- (cid:18) (cid:19)
mentum constraint, however,is notregularizedas it still = Mj(u)∂ v′ +p′(u,v)+λ Fr(u,v)−Ft(u,v)
i r j i i i
includes the term (K −K )/r. One could try to play
A B f (u)B
i (cid:0) (cid:1)
the same game as before and introduce yetanother vari- + ∂rλ, (3.14)
2αA
abletoabsorbthisterm. However,wewillnowshowthat (cid:18) (cid:19)
this is not really necessary. with Fr(u,v) = ∂ (f (u)B/2αA). Using now the fact
i r i
Letusthenconsidersomearbitraryfirstorderformula- that
tionoftheEinsteinevolutionequationsinsphericalsym-
metry that has the generic form 1 A
∂ λ=− λ+ (D −D ) , (3.15)
r A B
r B
(cid:20) (cid:21)
∂ u = q (u,v), (3.7)
t i i
we finally find
∂ v = Mj(u)∂ v +p (u,v), (3.8)
t i i r j i
∂ v′ = Mj(u)∂ v′ +p′(u,v)+λ Fr(u,v)−Ft(u,v)
t i i r j i i i
where u = (A,B,λ) and v = (D ,D ,K ,K ). The
A B A B
f (u)B A
source terms q and p are assumed not to depend on − i λ+ (D −D (cid:0)) . (3.16(cid:1))
A B
derivativesofanyofthefields. Theformulationmightbe 2αAr B
(cid:20) (cid:21)
hyperbolic or not, depending onthe characteristicstruc-
This last equation is now regular, and has precisely
ture of the matrix M. We will assume that one has ar-
the same characteristic structure as the original system.
rived at such a formulation by adding multiples of the
What we have done is transform the original evolution
hamiltonianandmomentum constraintsto the evolution
equationsforthev variablesintoevolutionequationsfor
equations for the v’s. This means that one can expect i
thenewv′variablesforwhichtheprincipalparttermsare
that the source terms p will in general contain terms i
i the same and the source terms are regular. Notice that
proportional to (K −K )/r. We will then rewrite the
A B typically only some of the f (u) will be different from
evolution equations for the v as i
i zero, so one does not need to transform all variables.
In the following section we will consider examples of
f (u)
∂tvi =Mij(u)∂rvj+p′i(u,v)+ ir (KA−KB). (3.9) howtoregularizesomespecificsystemsofevolutionequa-
tions.
Here we are assuming that the coefficient f (u) of the
i
(K −K )/rtermsdependsoftheu’s,butnotonthev’s,
A B
whichwilltypicallybethecase. Usingnowequation(3.5) IV. EXAMPLES
we find
The regularized first order ADM equations where de-
f (u)B
∂ v =Mj(u)∂ v +p′(u,v)+ i ∂ λ, (3.10) rivedinthelastsection. Herewewillexamplesofnumer-
t i i r j i 2αA t ical evolutions using two different regularized systems.
Inthenumericalsimulationswewilltakeasinitialdata
which implies
Minkowski spacetime in the usual coordinates, so that:
f (u)B
∂ v − i λ = Mj(u)∂ v +p′(u,v) A = B =1, (4.1)
t i 2αA i r j i
(cid:18) (cid:19) D = D =0, (4.2)
A B
f (u)B
i
− λ∂t 2αA . (3.11) KA = KB =0. (4.3)
(cid:18) (cid:19)
In orderto havea non-trivialevolution,we will chosean
If we now define
initial lapse profile of the form:
f (u)B
vi′ :=vi− i2αA λ, (3.12) α(t=0)=1+r2Cˆe−(r−σˆ2rˆ)2 , (4.4)
we can transform the last equation into thatis,weaddasmallgaussiancontributiontotheinitial
Minkowski lapse. We will then evolve the lapse using a
∂tvi′ =Mij(u)∂rvj +p′i(u,v)−λFit(u,v), (3.13) Bona-Masso (BM) slicing condition [8], so the evolution
equation for the lapse will be:
with Ft(u,v) = ∂ (f (u)B/2αA). Notice that Ft(u,v)
i t i i
2
sodefinedwillinvolvenospatialderivativesofu’sorv’s. ∂ α=−α f(α)(K +2K ). (4.5)
t A B
5
FIG. 1: System I not regularized. The plots show the evo- FIG. 2: System I not regularized. The plots show the evolu-
lution of the metric function A at different times. Note the tionofthelapsefunctionαatdifferenttimes. Notethespike
spikeat r=0 for t=1. at r=0 for t=1.
In the simulations shown below we have taken the ∂ A = −2αAK , (4.10)
t A
gaussian parameters to be: Cˆ = 0.001,rˆ= 5.0, σˆ = 1.0.
∂ B = −2αBK , (4.11)
t B
We have also restricted ourselves to harmonic slicing,
f2Bλ
that is, f(α)=1.
∂ U = −αf ∂ K +α(K +2K )
t α r A A B
A
(cid:20)
A. System I − Uα(f +αf′) +αf(KB−KA)
(cid:21)
2Bλ
As a first example we will build a hyperbolic system D − , (4.12)
B
A
starting from the ADM equations and the BM slicing (cid:18) (cid:19)
condition. We constructthis systembyusing the Hamil- fBλ
∂ D = −2α ∂ K +K U − , (4.13)
tonian and Momentum constraints to remove the terms t A r A A α A
(cid:20) (cid:18) (cid:19)(cid:21)
proportionalto∂ D and∂ K fromtheevolutionequa-
r B r B fBλ
tions of KA and Dα respectively. The resulting system ∂tDB = −2α ∂rKB+KB Uα− A , (4.14)
can be easily shown to be strongly hyperbolic. (cid:20) (cid:18) (cid:19)(cid:21)
Figures1and2showtheevolutionoftheradialmetric α α fBλ DA
∂ K = − ∂ U − U − U −
t A r α α α
A andlapse function α using the systemdesribed above, A A A 2
(cid:26)(cid:18) (cid:19)(cid:20)
withnoregularization. Notethatbothplotsshowaspike Bλ fBλ D2
at r =0 for times t≈1. − (f +αf′) + (DA−DB)− B
A A 4
Next we look at the regularized case. As described (cid:21)
1
above, we first introduce the auxilliary variable + A K2 −K2 − λ+D
B A r B
(cid:20)
λ= 1 1− A . (4.6) (cid:0) Bλ (cid:1)
r B − f +DA−DB , (4.15)
(cid:18) (cid:19) A
(cid:18) (cid:19)(cid:21)(cid:27)
Also, since we have used the momentum constraint to α α fBλ
2
∂ K = − ∂ D + − U − D −D
modify the evolution equation for Dα, we will need to t B 2A r B 2A α A B B
replace this variable with (cid:26) (cid:18) (cid:19)
D D 1
A B
+ +2AK (K +2K )+ D −2U
Bλ 2 B A B r A α
U :=D + . (4.7) (cid:20)
α α
A 2fBλ
+ −4D −2λ , (4.16)
B
The set of variables to be evolved is then: A
(cid:21)(cid:27)
2αA αAD
B
{α,A,B,Uα,DA,DB,KA,KB,λ}. (4.8) ∂tλ = B ∂rKB+ B (KB −KA). (4.17)
Theevolutionequationsforthesevariablestaketheform:
with f′ :=df/dα.
2
∂ α = −α f(K +2K ), (4.9) Figures 3 and 4 shownagain the evolution of A and α
t A B
6
the variable λ. Since in this case we have used the mo-
mentum constraint to modify the evolution equation for
D, we need to replace this variable with
4Bλ
e U :=D− . (4.18)
A
The final set of dynaemicalevariables is then
{α,A,B,D ,U,D ,K,K ,λ}, (4.19)
α B B
and their evolution equations are:
e
2
∂ α = −α fK , (4.20)
t
∂ A = 2αA(2K −K), (4.21)
t B
∂ B = −2αBK , (4.22)
t B
∂ D = −∂ (αfK), (4.23)
t α r
FIG. 3: System I regularized. The plots show the evolution ∂tU = −2∂r(αK)+4αDB(K −3KB)
of the metric function A at different times.
Bλ
+ 8α D K + (3K −K) , (4.24)
e α B B
A
(cid:20) (cid:21)
∂ D = −2∂ (αK ), (4.25)
t B r B
2D
2 α 2
∂ K = α −4KK +6K − +K
t B B Ar
(cid:20)
D 4λB D2 ∂ D
+ α U + − α − r α ,(4.26)
2A A A A
(cid:18) (cid:19) (cid:21)
α Ue 2λB
∂ K = + −D −λ−D
t B B α
Ar 2 A
" #
e
α D D ∂ D
α B r B
+ − −
A 2 2
(cid:20)
D 4λB
B
+ U + +AKK , (4.27)
B
4 A
(cid:18) (cid:19) (cid:21)
2αA D
e B
∂ λ = ∂ K − (K−3K ) . (4.28)
t r B B
B 2
(cid:20) (cid:21)
FIG. 4: System I regularized. The plots show the evolution
In this case there is no need to show any plots, as the
of the lapse function α at different times.
numericalevolutionbehavesexactlyinthe wayitdidfor
System I, with the origin remaining well behaved.
for the regularizedcase. The system can now be evolved
for very long times and the origin remains well behaved.
V. DISCUSSION
The plots show the evolution up to t=12.
Lack of regularity of geometric variables at the origin
is often a problem for spherically symmetric evolution
B. System II
codes in numerical relativity. We have shown that the
problem can be traced to the existence of two types of
Oursecondexampleisagainahyperbolicsystembuilt regularity conditions at the origin. In the first place,
from the ADM equations and the BM slicing condi- there are regularity conditions that guarantee that the
tion, but now instead of using DA and KA as funda- variablesarewelldefined atthe origin. These conditions
mental variables, we have used D = DA − 2DB and can be written down as a series of symmetry conditions
K ≡trK:=KA+2KB. In order to make the system hy- atthe originforthe differentvariables,andcanbe easily
perbolic we remove the terms propeortional to ∂rDB and enforced in numerical simulations. However, there also
∂rKB fromtheevolutionequationsofK andDusingthe exist regularity conditions related to the condition that
constraints. We then obtain a new strongly hyperbolic spacetime must be locally flat at the origin. Together,
system different from System I above. e all these regularity conditions imply that we have more
For the regularization procedure we introduce again conditions to satisfy at r = 0 than dynamical variables,
7
whichmeansthatnumericallysomeregularityconditions if not in detail, to an algorithm presented by Arbona
will inevitably be violated. We have presented a generic and Bona for the particular case of the Bona-Masso for-
regularization algorithm that is based on the introduc- mulation [4]. We have also shown the effectiveness of
tionofanauxiliaryvariablethatabsorbstheproblematic ouralgorithmonacoupleofdifferentformulationsofthe
terms and on which we can impose the extra boundary evolution equations.
conditions at r = 0. Our algorithm is similar in spirit,
[1] J.M.BardeenandT.Piran,Phys.Rep.196,205(1983). Lett. 75, 600 (1995), gr-qc/9412071.
[2] M. W. Choptuik,Phys. Rev.D 44, 3124 (1991). [9] C. Bona,J.Mass´o, E.Seidel,andJ.Stela,Phys.Rev.D
[3] M. W. Choptuik,Phys. Rev.Lett. 70, 9 (1993). 56, 3405 (1997), gr-qc/9709016.
[4] A.Arbonaand C. Bona, Comput. Phys. Commun. 118, [10] R.Arnowitt,S.Deser,andC.W.Misner,inGravitation:
229 (1999), gr-qc/9805084. AnIntroductiontoCurrentResearch,editedbyL.Witten
[5] C. Bona and J. Mass´o, Phys.Rev. D 40, 1022 (1989). (John Wiley, New York,1962), pp.227–265.
[6] C.BonaandJ.Mass´o,Phys.Rev.Lett.68,1097(1992). [11] J.York,inSources of Gravitational Radiation,editedby
[7] C. Bona and J. Mass´o, International Journal of Modern L.Smarr(CambridgeUniversityPress,Cambridge,Eng-
Physics C: Physics and Computers 4, 88 (1993). land, 1979).
[8] C. Bona, J. Mass´o, E. Seidel, and J. Stela, Phys. Rev.