Table Of ContentIllustrating Stability Properties of Numerical Relativity in Electrodynamics
A. M. Knapp1, E. J. Walker1 and T. W. Baumgarte1,2
1 Department of Physics and Astronomy, Bowdoin College, Brunswick, ME 04011 and
2 Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL, 61801
WeshowthatareformulationoftheADMequationsingeneralrelativity,whichhasdramatically
improved the stability properties of numerical implementations, has a direct analogue in classical
electrodynamics. We numerically integrate both the original and the revised versions of Maxwell’s
equations,andshowthattheirdistinctnumericalbehaviorreflectsthepropertiesfoundinlinearized
general relativity. Our results shed further light on the stability properties of general relativity,
illustrate them in a very transparent context, and may provide a useful framework for further
improvement of numerical schemes.
PACSnumbers: 04.25.Dm,02.60.Lj,95.30.Sf
2
0
Motivated by the prospect of gravitational wave de- and the two evolution equations
0
2 tectionsandtheaccompanyingneedfortheoreticalgrav-
itational wave templates, much effort has recently gone dtγij =−2αKij (3)
n
into the development of numerical relativity algorithms
a
and
J that arecapable ofmodeling the most promisingsources
of gravitational radiation, in particular the inspiral and
5 d K =−D D α+α(R −2K Kl +KK −M ).
1 coalescence of binary black holes and neutron stars. In t ij i j ij il j ij ij
(4)
thepast,progresshasbeenhamperedbynumericalinsta-
1 bilities that arise in straight-forwardimplementations of
v the traditional Arnowitt-Deser-Misner (ADM, [1]) 3+1 Hereγij isthespatialmetric,Kij theextrinsiccurvature,
1 decompositionofEinstein’sequations(e.g.[2,3]). These K = γijKij its trace, α and βi are the lapse function
5 andtheshiftvector,andρ,S andM aremattersource
instabilities have been associatedwith the mathematical i ij
0 terms. The time derivative is defined as d = ∂ −L ,
structureofthe ADMequations,andasacureanumber t t β
1 and R
0 ofhyperbolicformulationshavebeensuggested(e.g.[2,4] ij
2 as well as [5] and references therein). Alternatively, Shi-
1
0 bataandNakamura[6]andlaterBaumgarteandShapiro Rij = − γkl γij,kl+γkl,ij −γkj,il −γil,kj (5)
2
c/ [7]suggestedareformulationofthe ADMequationsthat (cid:16) (cid:17)
q has been demonstrated to dramatically improve the sta- +γkl ΓmilΓmkj −ΓmijΓmkl ,
r- bility properties of numerical implementations (e.g. [8]). (cid:16) (cid:17)
g Whiletheexactreasonforthisimprovementisstillsome- D andΓi aretheRiccitensor,thecovariantderivative,
i jk
: whatmysterious(see[9],hereafterAABSS,and[10]),this and the connection coefficients associated with γ . Fi-
v ij
i newformulation,nowoftencalledtheBSSNformulation, nally, R is the scalar curvature R = γijRij. Equations
X is quite widely used. (1) through (4) are commonly refered to as the ADM
r InthisBriefReportweshowthattheBSSNreformula- equations [1].
a
tionofthe ADM equationshas a directanalogueinclas- The first term in the Ricci tensor (5) is an elliptic op-
sicalelectrodynamics(E&M).Wenumericallyimplement erator acting on the components of the spatial metric
both the original and the revised versions of Maxwell’s γ . If the Ricci tensor contained only this term, the
ij
equations, and find that their distinct numerical prop- two evolution equations (3) and (4) could be combined
erties reflect those found in general relativity (GR). As to forma waveequationfor γ . This property is spoiled
ij
suggested by AABSS, these properties can be identified by the appearance of the three other second derivative
with the propagation of constraint violating modes. We termsin(5),suggestingthattheseterms maybe respon-
presentourfindingsinthehopethat,inadditiontobeing sible forthe appearanceofinstabilities inmany straight-
a very transparent illustration of the stability properties forward,three-dimensionalimplementationsoftheADM
ofGR,they mayproveusefulforthefuture development equations. Thisproblemcanbeavoidedbyeitherusinga
and improvement of numerical algorithms. hyperbolicformulationofGR(e.g.[5]),orbyeliminating
In a 3+1 decomposition of GR, Einstein’s equations themixedsecondderivativesasintheBSSNformulation.
split into the two constraint equations In this formulation,the conformally related metric γ¯
ij
is defined as γ¯ = e−4φγ , where the conformal factor
ij ij
R−KijKij +K2 =2ρ (1) eφ is chosen so that the determinant of γ¯ij is unity, γ¯ =
1. The conformal exponent φ as well as the trace of
and the extrinsic curvature, K, are evolved as independent
variables. Their evolution equations can be found from
D Kj −D K =S (2)
j i i i the traces of equations (3) and (4), while the trace-free
2
parts of those equations form evolution equations for γ¯ In terms of a vector potential A , Maxwell’s equations
ij i
and the trace-free part of the extrinsic curvature, A¯ . can be written as the evolution equations
ij
The latter equation still contains the Ricci tensor R¯
ij
associated with γ¯ij, which contains all the mixed second ∂tAi = −Ei−Diψ (10)
derivativesof(5). Thecrucialstepistorealizethatthese ∂ E = −DjD A +DjD A −4πj , (11)
t i j i i j i
second derivatives can be absorbed in a first derivative
and the constraint equation
of the “conformalconnection functions”
D Ei =4πρ . (12)
Γ¯i ≡γ¯jkΓ¯i =−γ¯ij , (6) i e
jk ,j
Here Ei is the electrical field, ρ the charge density, j
e i
where the last equality holds because γ¯ = 1. Here and
the flux, and ψ the scalar gauge potential. Identifying
in the following all barred quantities are associated with A with γ , E with K , and ψ with βi, we see that
γ¯ . In terms of Γ¯i, the Ricci tensor can be written i ij i ij
ij the structure of the evolution equations (10) and (11) is
very similar to that of equations (3) and (4), and that
1
R¯ = − γ¯lmγ¯ +γ¯ ∂ Γ¯k+Γ¯kΓ¯ + (7) the constraint equation (12) can be similarly identified
ij 2 ij,lm k(i j) (ij)k
withthemomentumconstraint(2). InanalogywithGR,
γ¯lm 2Γ¯k Γ¯ +Γ¯k Γ¯ we refer to equations (10) through (12) as System I.
l(i j)km im klj
In further analogy,we willeliminate the mixed second
(cid:16) (cid:17)
derivative in (11) by introducing a new variable Γ
(compare [11]). Evidently, the only remaining second
derivativetermistheellipticoperatorγ¯lmγ¯ ,iftheΓ¯i
ij,lm Γ≡D Ai (13)
are considered as independent functions. For that pur- i
pose, an evolution equation is derived by permuting a (compare (6)), in terms of which (11) becomes
time and space derivative in (6)
∂ E =−DjD A +D Γ−4πj . (14)
t i j i i i
2
∂tΓ¯i =−∂j 2αA¯ij −2γ¯m(jβi),m+ 3γ¯ijβl,l+βlγ¯ij,l . As in (7), the mixed second derivatives have been ab-
(cid:16) (cid:17)(8) sorbed in a first derivative of the new variable. An evo-
lution equationfor Γ can againbe derived by permuting
The divergence of the extrinsic curvature can now be a time and space derivative in the definition of the new
eliminated with the help of the momentum constraint variable (13)
(2), which yields the evolution equation
∂ Γ=D ∂ Ai =−D Ei−DiD ψ =−4πρ −DiD ψ.
t i t i i e i
∂ Γ¯i =−2A¯ij∂ α (9) (15)
t j
2
+2α Γ¯i A¯kj − γ¯ij∂ K−γ¯ijS +6A¯ij∂ φ Here we have used the constraint (12), which, as in GR,
jk 3 j j j will turn out to be crucial (see (21) below). Equations
+βj∂(cid:16)¯ Γi−Γ¯j∂ βi+ 2Γ¯i∂ βj + 1γ¯liβj +γ¯(cid:17)ljβi . (10), (14), (15) form the evolution equations of what we
j j 3 j 3 ,jl ,lj call System II. The definition (13) together with (12)
form the constraint equations of System II.
As in[7]we willreferto the ADM equations(1)through
For vanishing scalar potential ψ = 0, an analytical
(4)asSystemI,andtothenewBSSNsystemofequations
vaccuum solution to Maxwell’s equations is given by the
as System II. A complete listing and derivation of the
purely toriodal dipole field
latter can be found in [7].
Inanefforttobetter understandthe improvednumer- e−λv2 −e−λu2 ve−λv2 +ue−λu2
ical behavior of System II, AABSS linearized Systems I Aφˆ=Asinθ −2λ .
r2 r
!
andII andidentifiedtheir characteristicstructure. They
(16)
found that SystemI has constraintviolating modes with
a characteristic speed of zero. In System II, the charac-
Here A is the amplitude, λ parameterizes the size of the
teristicspeedofthesemodeschangestothespeedoflight.
wavepackage,anduandv aretheretardedandadvanced
AABSSfurther demonstratedthatina non-linearmodel time u≡t+r and v ≡t−r. According to (10), Eφˆ can
problem the existence of non-propagating,constraint vi-
befoundbytakingatimederivativeof(16). Sinceρ =0
olating modes may lead to numerical instabilities, as en- e
and D Ai = 0 in this solution, ψ = 0 is consistent with
counteredinimplementationsofSystemI.Theiranalysis i
both the Coulomb gauge and the Lorentz gauge
alsodemonstratedthatthe usageofthe momentumcon-
straint in the derivation of (9) is crucial for the charac- ∂ ψ =−D Ai−4πρ =−Γ−4πρ (17)
t i e e
teristicspeedfortheconstraint-violatingmodetochange
to a non-zero value, and hence for the stability of the (where the second equality applies for System II), if ap-
system. In the following we will show that very similar propriateboundaryandinitialconditionsarechosen. All
properties can be found in E&M. results shown below were obtained with Lorentz gauge.
3
FIG. 1: The integrated constraint violation kCk≡( C2dV)1/2 for evolutions usingLorentz gauge with resolution 323 and 643
with outer boundaries at 4 (OB4) and 8 (OB8).
R
Asinitialdataforourdynamicalsimulationsweadopt time t=9.375, C is larger in System II than in System
the analytical solution (16) at t=0 I, as one expects from Fig. 1. At a much later time
(t=200), however, System I has settled down to a con-
Aφˆ=0, Eφˆ =8Arsinθλ2e−λr2 (18) stant shape, while C in System II has almost completely
dissipated.
with A = λ = 1 and transformed into cartesian coordi- These numerical results demonstrate that, as in GR,
nates. the constraintviolationC behaves very differently in the
We numerically implement System I and II following two systems. In E&M, this different behavior can be
the algorithm of [7] as closely as possible. In particular, understood very easily from an analytic argument.
wewrotethecodeinthreespatialdimensionsusingcarte- For System I, it is easy to show that a time derivative
sian coordinates. We use an iterative Crank-Nicholson of the constraint violation (19) vanishes
scheme [12] to update the evolution equations, and im-
pose outgoing wave boundary conditions on Ei, Ai and ∂tC =Di∂tEi−4π∂tρe =−4π(Diji+∂tρe)=0, (20)
Γ as in [7]. We verified that the numerical solution con-
wherewehaveusedthecontinuityequationD ji+∂ ρ =
verges to the analytical solution to second order as long i t e
0. ThisexplainswhyatlatetimestheprofileofCremains
as the solution is not affected by the outer boundaries
unchanged in System I. This property is the analogue of
(which are first order accurate).
AABSS’sfindingthatthelinearizedADMequationshave
We compare the performance of System I and II by
non-propagating,constraint violating modes.
monitoring the constraint violation
ForSystemII,ontheotherhand,itcanbeshownthat
C ≡D Ei−4πρ . (19) C satisfies a wave equation
i e
∂2C = ∂ Di∂ E −4π∂2ρ
In Fig. 1, we show integrated values kCk ≡ ( C2dV)1/2 t t t i t e
for Systems I and II for two different gridsizes (323 and = ∂tDi(−DjDjAi+DiΓ−4πji)−4π∂t2ρe
643) and two different locations of the outerRboundary = −Di(D Dj∂ A −D ∂ Γ)−4π∂ (Dij +∂ ρ )
j t i i t t i t e
(at 4 (OB4) and 8 (OB8)).
= Di D Dj(E +D ψ)−D (4πρ +DjD ψ)
At early times, System I violates the constraint (12) j i i i e j
to a lesser degree than System II. After about a light- = D (cid:16)Dj(DiE −4πρ )=D DjC, ((cid:17)21)
j i e j
crossing time, when the electro-magnetic wave has left
thenumericalgrid,kCksettlesdowntoanearlyconstant which explains why the constraint violations propagate
value,whichprimarilydependsonthegrid-resolution. In awayinSystemII.ThisistheanalogueofAABSS’sfind-
System II, kCkis also largestafter about a light-crossing ing that in the relativistic System II (the BSSN equa-
time, but after that kCk decreases exponentially. As one tions)the constraintviolatingmodespropagatewith the
might expect, the decay time of this exponential fall-off speed of light. Moreover, we now realize why the usage
scales with the location of the outer boundary. of the constraint (12) in (9) was crucial – without hav-
To further compare these Systems, we compare snap- ing made this substitution the terms on the last line of
shotsoftheconstraintviolationC forthe323 OB4evolu- (21) would have cancelled, ∂2C = 0, leading to a non-
t
tion in Fig. 2. With identical initial data, both Systems propagatingconstraint violation as in System I. The ad-
have identical values for C at t = 0. At an intermediate dition of constraint equations to the evolution equations
4
0.03 0.03
0.02 0.02
0.01 0.01
0 0
−0.01 −0.01
−0.02 −0.02
−0.03 −0.03
4 4
3 4 3 4
2 3 2 3
2 2
1 1 1 1
0 0 0 0
0.03 0.03
0.02 0.02
0.01 0.01
0 0
−0.01 −0.01
−0.02 −0.02
−0.03 −0.03
4 4
3 4 3 4
2 3 2 3
2 2
1 1 1 1
0 0 0 0
0.03 0.03
0.02 0.02
0.01 0.01
0 0
−0.01 −0.01
−0.02 −0.02
−0.03 −0.03
4 4
3 4 3 4
2 3 2 3
2 2
1 1 1 1
0 0 0 0
FIG. 2: C in the z = .0625 plane (the first of our grid points) for System I (left column) and System II (right column) for a
323 OB4 evolution in Lorentz gauge. We show results for t=0 (top row), t=9.375 (middle row) and t=100 (bottom row).
C settles down toa constant profile in System I,but dissipates away in System II.
hasbeenfoundtobecrucialinseveralotherformulations and the boundary. Similar problems are likely to occur
of Einstein’s equations as well (e.g. [4, 13]). in GR, but might be easier to analyze in the simpler
To summarize, we see that the numerical stability framework of E&M.
properties of two formulations of GR are beautifully re- Acknowledgments
flected in similar formulations of E&M. Maxwell’s equa-
tions therefore provide a verytransparentframeworkfor
analyzing these properties, which may be useful for fu- This paper was supported in part by NSF Grant
ture algorithm development. For example, Fig. 2 shows PHY 99-02833 to the University of Illinois at Urbana-
that the outer boundaries produce much more noise in Champaign and Bowdoin College as a subrecipient.
System I than in System II, which points to an incon- AMK gratefully acknowledges support from the Mellon
sistency between the treatment of the interior equations Foundation.
[1] R. Arnowitt, S. Deser & C. W. Misner, in Gravitation: 75, 600 (1995).
AnIntroductiontoCurrentResearch,editedbyL.Witten [3] A.M.Abrahamset.al.,Phys.Rev.Lett.80,1812(1998).
(Wiley,New York,1962). [4] A.Anderson&J.W.York,Jr.,Phys.Rev.Lett.82,4384
[2] C.Bona,J.Mass´o, E.Seidel&J.Stela,Phys.Rev.Lett. (1999).
5
[5] O.Reula, Living Rev.Rel. 1, 3 (1998). M. Suen,Phys.Rev. D 62, 124011 (2000) (AABSS).
[6] M. Shibata & T. Nakamura, Phys. Rev. D 52, 5428 [10] S. Frittelli & O. Reula, J. Math. Phys. 40, 5143 (1999);
(1995). Miller, M., gr-qc/0008017 (2000).
[7] T. W. Baumgarte & S. L. Shapiro, Phys. Rev. D 59, [11] T. De Donder, La gravifique einsteinienne (Gauthier-
024007 (1999). Villars,Paris,1921);C.Lanczos,Phys.Z.23,537(1922);
[8] M.Alcubierreet.al.,Phys.Rev.D62,044034(2000);L. [12] S. A. Teukolsky,Phys. Rev.D 61, 087501 (2000).
Lehner,M.Huq,&D.Garrison,Phys.Rev.D62,084016 [13] B. Kelly et. al., Phys. Rev. D 64, 084013 (2001); G.
(2000);M.Alcubierre&B.Bru¨gmann,Phys.Rev.D63, Yoneda & H.Shinkai, Phys.Rev.D 63 124019 (2001).
104006 (2001).
[9] M.Alcubierre,G.Allen,B. Bru¨gmann, E.Seidel,&W.-
