2008/11/28 @ 大阪市立大学

Formulation for Numerical Relativity
both for Einstein / Gauss-Bonnet
真貝寿明

大阪工業大学情報科学部

1. Introduction
定式化問題？
2. The Standard Approach to Numerical Relativity
ADM 形式，BSSN 形式，Hyperbolic 形式
3. Robust system for Constraint Violation
Adjusted systems, Adjusted ADM, Adjusted BSSN
4. 高次元数値相対論に向けて
ref. arXiv:0805.0068 and 0810.1790

http://www.is.oit.ac.jp/~shinkai/「応用数理」日本応用数理学会学会誌， vol.15, No.1 (2005)，p2-15

In the last 5 years, ...
Binary BH-BH coalescencesimulations are available!!

INSPIRAL PHASE
Newtonian / Post-Newtonian

Coalescence / Merger

・Pretorius (2005)
・Univ. Texas Brownsville (2006)
・NASA-Goddard (2006)

Innermost Stable Circular Orbit
Post-Newtonian / GR

Breakthrough suddenly occurs.

Black Hole Formation
Quasinormal Ringing

Black HoleIn the last 5 years, ...
Binary BH-BH coalescencesimulations are available!!
・Pretorius (2005)--> Princeton Univ.
・Univ. Texas Brownsville (2006) --> Rochester Univ.
・NASA-Goddard (2006)
・Louisiana State Univ.
・Jena Univ.
・Pennsylvania State Univ.
"Gold-Rush of parameter searches" (B. Bruegmann, July 2007 @GRG)
But ..... Why it works?

Goals of the Talk
何を指標にして発展方程式を選択すれば良いのか？
どうして多くのグループが BSSN
形式を使っているのか？
BSSN 形式に代わる formulation
はあるか？

The Einstein equation
geodesics

spacetime curvature

matter distribution1Rµν − gµν R + Λgµν = 8πG Tµν
2
Einstein tensor

Energy-Momentum tensor

cosmological constant

Solve for metricgµν (t,x,y,z)(10 components)flat
at spacetime (Mink(Minkowskii
wskii spacetime):
ds2 = −dt
dt2 + dx2 + dy 2 + dz 2
= −dt
dt2 + dr2 + r2(dθ2 + sin2 θdϕ2)

ds2 =

!

µ,ν

gµν dxµdxν := gµν dxµdxν


gtt



gtx gty gtz 





g
g
g


xx
xy
xz


gµν = gyy gyz 




sym.
sy
m.
gzz

First Question: How to foliate space-time?

Cauchy approachor ADM 3+1 formulation

Characteristic approach(if null, dual-null 2+2 formulation)

time direction

S: Initial 3-dimensional Surface

ingoing
direction

outgoing
direction

S: Initial 2-dimensional Surface

3+1 versus 2+2

pioneers
variables

foliation
initial data
evolution
singularity
disadvantages

Cauchy (3+1) evolution
ADM (1961), York-Smarr (1978)

Characteristic (2+2) evolution
Bondi et al (1962), Sachs (1962),
Penrose (1963)
easy to understand the concept of has geometrical meanings
time evolution
1 complex function related to 2 GW
polarization modes
has Hamilton structure
allows implementation of Penrose’s
space-time compactificationneed to solve constraints
no constraints
PDEs
ODEs with consistent conditions
need to avoid constraint violation propagation eqs along the light raysneed to avoid by some method
can truncate the gridcan not cover space-time globally difficulty in treating caustics
hard to treat matterProcedure of the Standard Numerical Relativity■ 3+1 (ADM) formulation■ Preparation of the Initial Data
◆ Assume the background metric◆ Solve the constraint equations■ Time Evolution

do time=1, time_end◆ Specify the slicing condition◆ Evolve the variables◆ Check the accuracy◆ Extract physical quantities
end doThe 3+1 decomposition of space-time:

The ADM formulation[1 ] R. Arnowitt, S. Deser and C.W. Misner, in Gravitation: An Introduction to Current Research,ed. by L.Witten, (Wiley, New York, 1962).
[2 ] J.W. York, Jr. in Sources of Gravitational Radiation, (Cambridge, 1979)
Dynamics of Space-time = Foliation of Hypersurface• Evolution of t =const. hypersurface Σ(t).
ds2 = gµν dxµdxν ,on Σ(t)... dℓ2 = γij dxi dxj ,(µ, ν = 0, 1, 2, 3)(i, j = 1, 2, 3)

time direction• The unit normal vector of the slices, nµ .
nµ = (−α, 0, 0, 0)
nµ = g µν nν = (1/α, −β i/α)
• The lapse function, α. The shift vector, β i.

Σ: Initial 3-dimensional Surface

ds2 = −α2dt2 + γij (dxi + β idt)(dxj + β j dt)

The decomposed metric:
ds2 = −α2dt2 + γij (dxi + β idt)(dxj + β j dt)= (−α2 + βl β l )dt2 + 2βidtdxi + γij dxidxj


2

l



−α + βl β βj 
gµν =
,
βi
γij






g

µν

j

2

βj ≡ g0j .

• The unit normal vector of the slices, nµ.

• The lapse function, α.



−1/α
β /α


=
i
2
ij
i j
2
β /α γ − β β /α




√ 00
where α and βj are defined as α ≡ 1/ −g ,nµ = (−α, 0, 0, 0)
nµ = g µν nν = (1/α, −β i/α)2shift vector, βisurface normal line

lapse function, αβi dtcoordinate constant line

α dt

• The shift vector, β i.

t = constant hypersurface

Projection of the Einstein equation:• Projection operator (or intrinsic 3-metric) to Σ(t),γµν = gµν + nµnνγνµ = δνµ + nµnν ≡ ⊥µν• Define the extrinsic curvature Kij ,Kij ≡ −⊥µi⊥νj nµ;ν= −(δiµ + nµni)(δjν + nν nj )nµ;ν
= −ni;j
1 α= Γij nα = · · · =−∂tγij + βi|j + βj|i .
2α• Projection of the Einstein equation:Gµν nµ nν = 8πG Tµν nµ nν ≡ 8πρHGµν nµ ⊥νi = 8πG Tµν nµ ⊥νi ≡ −8πJi
Gµν ⊥µi ⊥νj = 8πG Tµν ⊥µi ⊥νj ≡ 8πSij⇒ the Hamiltonian constraint eq.
⇒ the momentum constraint eqs.
⇒ the evolution eqs.

The Standard ADM formulation (aka York 1978):The fundamental dynamical variables are (γij , Kij ), the three-metric and extrinsic
curvature. The three-hypersurface Σ is foliated with gauge functions, (α, β i), thelapse and shift vector.
• The evolution equations:∂tγij = −2αKij + Diβj + Dj βi,∂tKij = α (3)Rij + αKKij − 2αKik K k j − DiDj α+(Diβ k )Kkj + (Dj β k )Kki + β k Dk Kij−8πGα{Sij + (1/2)γij (ρH − trS)},

where K = K ii, and (3)Rij and Di denote three-dimensional Ricci curvature,
and a covariant derivative on the three-surface, respectively.
• Constraint equations:
Hamiltonian constr.
momentum constr.
where (3)R =(3) Rii.

HADM := (3)R + K 2 − Kij K ij ≈ 0,
MADM:= Dj K j i − DiK ≈ 0,istrategy 0 The standard approach :: Arnowitt-Deser-Misner (ADM) formulation (1962)
shift vector, βi
3+1 decomposition of the spacetime.
surface normal line
Evolve 12 variables (γij , Kij )coordinate constant line
βi dtwith a choice of gauge condition.

lapse function, αα dt

t = constant hypersurface

constraints

evolution eqs.

Maxwell eqs.
div E = 4πρ
div B = 0

ADM Einstein eq.
(3)
R + (trK)2 − Kij K ij = 2κρH + 2ΛDj K ji − DitrK = κJi

14π∂tE = rot B − jcc∂tγij = −2N Kij + Dj Ni + DiNj ,∂tKij = N ( (3)Rij + trKKij ) − 2N Kil K l j − DiDj N+ (Dj N m )Kmi + (DiN m )Kmj + N m DmKij − N γij Λ− κα{Sij + 12 γij (ρH − trS)}

1∂tB = −rot EcProcedure of the Standard Numerical Relativity■ 3+1 (ADM) formulation■ Preparation of the Initial Data
◆ Assume the background metric◆ Solve the constraint equations■ Time Evolution

do time=1, time_end◆ Specify the slicing condition◆ Evolve the variables◆ Check the accuracy◆ Extract physical quantities
end doNeed to solve elliptic PDEs-- Conformal approach-- Thin-Sandwich approachProcedure of the Standard Numerical Relativity■ 3+1 (ADM) formulation■ Preparation of the Initial Data
◆ Assume the background metric

Need to solve elliptic PDEs-- Conformal approach-- Thin-Sandwich approach◆ Solve the constraint equations■ Time Evolution

do time=1, time_end◆ Specify the slicing condition◆ Evolve the variables◆ Check the accuracy◆ Extract physical quantities
end dosingularity avoidance,
simplify the system,
GW extraction, ...

Procedure of the Standard Numerical Relativity■ 3+1 (ADM) formulation■ Preparation of the Initial Data
◆ Assume the background metric

Need to solve elliptic PDEs-- Conformal approach-- Thin-Sandwich approach◆ Solve the constraint equations■ Time Evolution

do time=1, time_end◆ Specify the slicing condition◆ Evolve the variables◆ Check the accuracy◆ Extract physical quantities
end dosingularity avoidance,
simplify the system,
GW extraction, ...

Robust formulation ?
-- modified ADM / BSSN-- hyperbolization-- asymptotically constrained

Formulation Problem

strategy 0 The standard approach :: Arnowitt-Deser-Misner (ADM) formulation (1962)
shift vector, βi
3+1 decomposition of the spacetime.
surface normal line
Evolve 12 variables (γij , Kij )coordinate constant line
βi dtwith a choice of gauge condition.
lapse function, αα dt

t = constant hypersurface

constraints

evolution eqs.

Maxwell eqs.
div E = 4πρ
div B = 0

ADM Einstein eq.
(3)
R + (trK)2 − Kij K ij = 2κρH + 2ΛDj K ji − Di trK = κJi

14π∂tE = rot B − jcc∂t γij = −2N Kij + Dj Ni + DiNj ,∂t Kij = N ( (3)Rij + trKKij ) − 2N Kil K l j − Di Dj N
+ (Dj N m )Kmi + (DiN m )Kmj + N m DmKij − N γij Λ− κα{Sij + 21 γij (ρH − trS)}

1∂tB = −rot Ec

S. Frittelli, Phys. Rev. D55, 5992 (1997)
HS and G. Yoneda, Class. Quant. Grav. 19, 1027 (2002)

The Constraint Propagations of the Standard ADM:
∂tH = β j (∂j H) + 2αKH − 2αγ ij (∂iMj )+α(∂l γmk )(2γ ml γ kj − γ mk γ lj )Mj − 4γ ij (∂j α)Mi,∂tMi = −(1/2)α(∂iH) − (∂iα)H + β j (∂j Mi)+αKMi − β k γ jl (∂iγlk )Mj + (∂iβk )γ kj Mj .

From these equations, we know that
if the constraints are satisfied on the initial slice Σ,then the constraints are satisfied throughout evolution (in principle).

Primary / Secondary constraint
First-class / Second-class constraint• Primary Constraints• Secondary Constraints= when propagation of constraints require additional constraints• First-Class Constraints=

Numerical Relativity in the 20th century
1960s Hahn-Lindquist
May-White
1970s ÓMurchadha-York
Smarr
Smarr-Cades-DeWitt-Eppley
Smarr-York
ed. by L.Smarr
1980s Nakamura-Maeda-Miyama-Sasaki
MiyamaBardeen-Piran
Stark-Piran1990 Shapiro-Teukolsky
Oohara-Nakamura
Seidel-Suen
Choptuik
NCSA group
Cook et al
Shibata-Nakao-Nakamura
Price-Pullin
1995 NCSA group
NCSA group
Anninos et al
Scheel-Shapiro-Teukolsky
Shibata-Nakamura
Gunnersen-Shinkai-Maeda
Wilson-Mathews
Pittsburgh group
Brandt-Brügmann
Illinois group
Shibata-Baumgarte-Shapiro
BH Grand Challenge Alliance
Baumgarte-Shapiro
Brady-Creighton-Thorne
Meudon group
Shibata

2 BH head-on collision
spherical grav. collapse
conformal approach to initial data
3+1 formulation
2 BH head-on collision
gauge conditions“Sources of Grav. Radiation”
axisym. grav. collapse
axisym. GW collapse
axisym. grav. collapse
axisym. grav. collapse
naked singularity formation
3D post-Newtonian NS coalesence
BH excision techniquecritical behaviour
axisym. 2 BH head-on collision
2 BH initial data
BransDicke GW collapse
close limit approach
event horizon finder
hyperbolic formulation
close limit vs full numerical
BransDicke grav. collapse
3D grav. wave collapse
ADM to NP
NS binary inspiral, prior collapse?
Cauchy-characteristic approach
BH puncture datasynchronized NS binary initial data
2 NS inspiral, PN to GR
characteristic matching
Shibata-Nakamura formulationintermediate binary BH
irrotational NS binary initial data
2 NS inspiral coalesence

AnaPhys29(1964)304
PR141(1966)1232
PRD10(1974)428
PhD thesis (1975)
PRD14(1976)2443
PRD17(1978)2529
Cambridge(1979)
PTP63(1980)1229
PTP65(1981)894
PhysRep96(1983)205
unpublished
PRL66(1991)994
PTP88(1992)307
PRL69(1992)1845
PRL70(1993)9
PRL71(1993)2851
PRD47(1993)1471
PRD50(1994)7304
PRL72(1994)3297
PRL74(1995)630
PRL75(1995)600
PRD52(1995)4462
PRD51(1995)4208
PRD52(1995)5428
CQG12(1995)133
PRL75(1995)4161
PRD54(1996)6153
PRL78(1997)3606
PRL79(1997)1182
PRD58(1998)023002
PRL80(1998)3915
PRD59(1998)024007
PRD58(1998)061501
PRL82(1999)892
PRD60(1999)104052

Critical Phenomena in Gravitational Collapse

Choptuik, Phys. Rev. Lett. 70 (1993) 9

Spherical Sym., Massless Scalar Field(1) scaling(2) echoing(3) universality

Discrite
Self-Similarity

Head-on Collision of 2 Black-Holes (Misner initial data)
NCSA group 1995

S. Frittelli, Phys. Rev. D55, 5992 (1997)
HS and G. Yoneda, Class. Quant. Grav. 19, 1027 (2002)

The Constraint Propagations of the Standard ADM:
∂tH = β j (∂j H) + 2αKH − 2αγ ij (∂iMj )+α(∂l γmk )(2γ ml γ kj − γ mk γ lj )Mj − 4γ ij (∂j α)Mi,∂tMi = −(1/2)α(∂iH) − (∂iα)H + β j (∂j Mi)+αKMi − β k γ jl (∂iγlk )Mj + (∂iβk )γ kj Mj .

From these equations, we know that
if the constraints are satisfied on the initial slice Σ,then the constraints are satisfied throughout evolution (in principle).

S. Frittelli, Phys. Rev. D55, 5992 (1997)
HS and G. Yoneda, Class. Quant. Grav. 19, 1027 (2002)

The Constraint Propagations of the Standard ADM:
∂tH = β j (∂j H) + 2αKH − 2αγ ij (∂iMj )+α(∂l γmk )(2γ ml γ kj − γ mk γ lj )Mj − 4γ ij (∂j α)Mi,∂tMi = −(1/2)α(∂iH) − (∂iα)H + β j (∂j Mi)+αKMi − β k γ jl (∂iγlk )Mj + (∂iβk )γ kj Mj .

From these equations, we know that
if the constraints are satisfied on the initial slice Σ,then the constraints are satisfied throughout evolution (in principle).
But this is NOT TRUE in NUMERICS....

• By the period of 1990s, NR had provided a lot of physics:

Gravitational Collapse, Critical Behavior, Naked Singularity, Event Horizons,Head-on Collision of BH-BH and Gravitational Wavve, Cosmology, · · ·• However, for the BH-BH/NS-NS inspiral coalescence problem, · · · why ???

Many (too many) trials and errors, hard to find a definit recipe.

error

Blow upt=0

time
Constrained Surface(satisfies
Einstein's constraints)

Best formulation of the Einstein eqs. for long-term stable & accurate simulation?

“Convergence”= higher resolution runs approach to the continuum limit.
(All numerical codes must have this property.)• When the code has 2nd order finite difference scheme,
then the error should be scaled with.
• “Consistency”, Choptuik, PRD 44 (1991) 3124,

“Accuracy”= The numerical results represent the actual solutions.
(All numerical codes must have this property.)• Check the code with known results.

Gauge wave test in BSSN;
Kiuchi, HS, PRD (2008)

Best formulation of the Einstein eqs. for long-term stable & accurate simulation?
• Many (too many) trials and errors, hard to find a definit recipe.
Blow up

Blow up

error

ADM
BSSN

time

Mathematically equivalent formulations, but differ in its stability!
strategy 0: Arnowitt-Deser-Misner (ADM) formulation
strategy 1: Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation
strategy 2: Hyperbolic formulations
strategy 3: “Asymptotically constrained” against a violation of constraints
By adding constraints in RHS, we can kill error-growing modes⇒ How can we understand the features systematically?

80s

90s

2000s

87

95

99

Nakamura-Oohara

Shibata-Nakamura

Baumgarte-Shapiro

62

ADM

97

92

Alcubierre

Bona-Masso
95-97

99

ChoquetBruhat-York Anderson-York01

Kidder-Scheel-Teukolsky

96

Frittelli-Reula
99

lambda-system
97
86

Ashtekar

Iriondo-Leguizamon-Reula

99

Yoneda-Shinkai

strategy 1

Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation
T. Nakamura, K. Oohara and Y. Kojima, Prog. Theor. Phys. Suppl. 90, 1 (1987)
M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 (1995)
T.W. Baumgarte and S.L. Shapiro, Phys. Rev. D 59, 024007 (1999)

The popular approach. Nakamura’s idea in 1980s.
BSSN is a tricky nickname. BS (1999) introduced a paper of SN (1995).
• define new set of variables (φ, γ̃ij ,K,Ãij ,Γ̃i), instead of the ADM’s (γij ,Kij ) whereγ̃ij ≡ e−4φ γij ,Ãij ≡ e−4φ(Kij − (1/3)γij K),Γ̃i ≡ Γ̃ijk γ̃ jk ,and impose detγ̃ij = 1 during the evolutions.
• The set of evolution equations become(∂t − Lβ )φ = −(1/6)αK,(∂t − Lβ )γ̃ij = −2αÃij ,(∂t − Lβ )K = αÃij Ãij + (1/3)αK 2 − γ ij (∇i∇j α),(3)(∂t − Lβ )Ãij = −e−4φ(∇i∇j α)T F + e−4φαRij − e−4φ α(1/3)γij R(3) + α(K Ãij − 2Ãik Ãk j )∂tΓ̃i = −2(∂j α)Ãij − (4/3)α(∂j K)γ̃ ij + 12αÃji(∂j φ) − 2αÃk j (∂j γ̃ ik ) − 2αΓ̃k lj Ãj k γ̃ il


k
ij
kj
i
ki
j
ij
k
−∂j β ∂k γ̃ − γ̃ (∂k β ) − γ̃ (∂k β ) + (2/3)γ̃ (∂k β )Momentum constraint was used in Γi-eq.

• Calculate Riemann tensor as
φkm k
Rij = ∂k Γkij − ∂i Γkkj + Γm
Γ
−
Γ
Γ
=:
R̃
+
R
ij
ij
ij mk
kj mi

φ
Rij
= −2D̃iD̃j φ − 2g̃ij D̃ l D̃l φ + 4(D̃iφ)(D̃j φ) − 4g̃ij (D̃ l φ)(D̃l φ)R̃ij = −(1/2)g̃ lm∂lm g̃ij + g̃k(i ∂j)Γ̃k + Γ̃k Γ̃(ij)k + 2g̃ lm Γ̃kl(iΓ̃j)km + g̃ lm Γ̃kim Γ̃klj• Constraints are H, Mi.
But thre are additional ones, G i , A, S.
Hamiltonian and the momentum constraint equations
HBSSN = RBSSN + K 2 − Kij K ij ,MBSSN
= MADM
,
i
i

(1)
(2)

Additionally, we regard the following three as the constraints:
G i = Γ̃i − γ̃ jk Γ̃ijk ,A = Ãij γ̃ ij ,S = γ̃ − 1,

Why BSSN better than ADM?
Is the BSSN best? Are there any alternatives?

(3)
(4)
(5)

Some known fact (technical):• Trace-out Aij at every time step helps the stability.
Alcubierre, et al, [PRD 62 (2000) 044034]• “The essential improvement is in the process of replacing terms by the momentum
constraints”,
Alcubierre, et al, [PRD 62 (2000) 124011]• Γ̃i is replaced by −∂j γ̃ ij where it is not differentiated,
Campanelli, et al, [PRL96 (2006) 111101; PRD 73 (2006) 061501R]• Γ̃i-equation has been modified as suggested in Yo-Baumgarte-Shapiro [PRD 66(2002) 084026]
Baker et al, [PRL96 (2006) 111102; PRD73 (2006) 104002]
Some guesses:• BSSN has a wider range of parameters that give us stable evolutions in von
Neumann’s stability analysis.
Miller, [gr-qc/0008017]• The eigenvalues of BSSN evolution equations has fewer “zero eigenvalues” thanthose of ADM, and they conjectured that the instability can be caused by “zeroeigenvalues” that violate “gauge mode”.
M. Alcubierre, et al, [PRD 62 (2000) 124011]

80s

90s

2000s

87

95

99

Nakamura-Oohara

Shibata-Nakamura

Baumgarte-Shapiro

62

ADM

97

92

Alcubierre

Bona-Masso
95-97

99

ChoquetBruhat-York Anderson-York01

Kidder-Scheel-Teukolsky

96

Frittelli-Reula
99

lambda-system
97
86

Ashtekar

Iriondo-Leguizamon-Reula

99

Yoneda-Shinkai

80s

90s

Nakamura-Oohara

2000s

Shibata

87

95

99

Nakamura-Oohara

Shibata-Nakamura

Baumgarte-Shapiro

62

G-codeH-codeADM

NCSA

BSSN-code

AEI

97

92

PennState

Alcubierre

Bona-Masso
95-97

99

UWash

ChoquetBruhat-York Anderson-York
01

Cornell-Illinois
96

Kidder-Scheel-Teukolsky
CaltechFrittelli-Reula
Hern

99

LSU

lambda-system
Shinkai-Yoneda

97
86

Ashtekar

Iriondo-Leguizamon-Reula

99

Yoneda-Shinkai

strategy 2Hyperbolic formulation

Construct a formulation which reveals a hyperbolicity explicitly.
For a first order partial differential equations on a vector u,












u1   u1 












∂ u2  =  A  ∂x  u2  . ..  
.
|{z}characteristic part




t 






u1 

+
B u2 
.. 
|
{z
}
lower order part








Hyperbolic Formulation(1) Definition

Hyperbolic Formulation(2) Expectations• if strongly/symmetric hyperbolic ==> well-posed system– Given initial data + source terms -> a unique solution exists– The solution depends continuously on the data
– Exists an upper bound on (unphysical) energy norm

• Better boundary treatments<== existence of characteristic field• Known numerical techniques in
Newtonian hydro-dynamics

strategy 2Hyperbolic formulation

Construct a formulation which reveals a hyperbolicity explicitly.
For a first order partial differential equations on a vector u,












u1   u1 












∂ u2  =  A  ∂x  u2  . ..  
.
|{z}characteristic part




t 






u1 

+
B u2 
.. 
|
{z
}
lower order part








However,• ADM is not hyperbolic.

• BSSN is not hyperbolic.

• Many many hyperbolic formulations are presented.
One might ask ...
Are they actually helpful?
Which level of hyperbolicity is necessary?

Why many? ⇒ Exercise.

Exercise 1 of hyperbolic formulation

Wave equation(∂t∂t − c2∂x∂x)u = 0

Exercise 1 of hyperbolic formulation

Wave equation(∂t∂t − c2∂x∂x)u = 0[1a] use u as one of the fundamental variables.


∂t 





2





u  0 c   u  =∂x  
1 0
v
v

(6)

Eigenvalues = ±c. Not a symmetric hyperbolic, but a kind of strongly hyperbolic.
[1b]













u  0 c  u 
 =

∂t 
∂  
c 0 x v
v

(7)

Eigenvalues = ±c. Symmetric hyperbolic.
[2a] Let U = u̇, V = u′,













0 c2   U  U ∂t   =∂x  
V
1 0
V

(8)

Eigenvalues = ±c. Not a symmetric hyperbolic, but a kind of strongly hyperbolic.
[2b] Let U = u̇, V = cu′ ,













U  0 c  U 

 =
∂t 
∂  
V
c 0 x V
Eigenvalues = ±c. Symmetric hyperbolic.

(9)Exercise 1 of hyperbolic formulation[3a] Let v = u̇, w = v ′,





t



∂



(∂t∂t − c2∂x∂x)u = 0

Wave equation













u  0 0 0   u   v 









2
∂ 
 = 




0
0
c
+
v 
v
0

 x











0 1 0
w
w
0

(10)

Eigenvalues = 0, ±c. Not a symmetric hyperbolic, nor a strongly hyperbolic.
[3b] Let v = u̇, w = cv ′ ,

















0 0 0  u   v 
 u 






















∂t  v  = 00c ∂x  v  +  0 









w
0 c 0
w
0

(11)

Eigenvalues = 0, ±c. Not a symmetric hyperbolic, nor a strongly hyperbolic.
[4] Let f = u̇ − cu′ , g = u̇ + cu′ ,



∂t 











f   −c 0   f 
 =
∂  
0 c x gg

Eigenvalues = ±c. Symmetric hyperbolic, de-coupled.

(12)

Exercise 2 of hyperbolic formulation

Maxwell equations

Consider the Maxwell equations in the vacuum space,
div E = 0,
div B = 0,
1 ∂E
rot B −
= 0,
c ∂t
1 ∂B
rot E +
= 0.
c ∂t

(1a)
(1b)
(1c)
(1d)

Exercise 2 of hyperbolic formulation

Maxwell equations(cont.)• Take a pair of variables as ui = (E1, E2, E3, B1, B2, B3)T , and write (1c) and(1d) in the matrix form






l j
A
 Ei  ∼ ∂t   =  l ij
Bi
Ci
|





l j
B i   Ej .
l j  ∂l 
Bj
Di

{z

(2)

}

Hermitian?
• In the Maxwell case, we see immediately0ǫilm  ∂l um∂tui = c 
lm
−ǫi
0

or with the actual components











t









∂





E1 





E2 
0





E3 
 = c
l
l 


0
δ
−δ
B1 

3
2





l
l



−δ
0
δ

B2 
3
1 



l
l
δ2 −δ10
B3That is, symmetric hyperbolic system.









0
δ3l
−δ2l

−δ3l
0
δ1l
0

δ2l
−δ
0




l 


1


















E1 
E2 

E3 
∂.
B1 

B2 

B3










l









Exercise 2 of hyperbolic formulation

Maxwell equations(cont.)• The eigen-equation of the characteristic matrix becomes


 
l
l
l 
0
−δ
δ
−λ
0
0


3
2 
 



 




l
l
l
 

 δ

0
−λ
0
c
0
−δ







3
1
 





l j
j
l j
l
l
l
l


 A

− λ δi
Bi


0
0
−λ
−δ
δ
0
i


2
1

 = det  
det 
 = 0
l
l
l

l j
l j


l j
0
δ
−δ
−λ
0
0

Ci
D i − λ δi


 

3
2



 




 

l
l
l


 

−δ
0
δ
0
−λ
0
c


 

3
1




 
l
l
l
δ2 −δ100−λ0
We therefore obtain the eigenvalues ass

0 (2 multi),±c (δ1l )2 + (δ2l )2 + (δ3l )2 ≡ ±c (2 each)

Exercise 3 of hyperbolic formulation

Adjusted Maxwell equations

By adding constraints (1a) and (1b) in the RHS of equations, and see what will behappend.












−ǫilm 
 z 
 0
x
 ∂k Bk , ∂l um + c  ∂k Ek + c ∂tui = c  lmǫi0yw
where x, y, z, w are parameters.

(3)

Exercise 3 of hyperbolic formulation

Adjusted Maxwell equations(cont.)

By adding constraints (1a) and (1b) in the RHS of equations, and see what will behappend.












−ǫilm 
 0
x
 z 
 ∂k Ek + c  ∂l um + c ∂ B ,∂tui = c  lm
y
ǫi
0
w k k
where x, y, z, w are parameters.
• The actual components are





E1



E 

 2




E 

 3
 = c
∂t 
  l
B 
δ

 1
  1



l
B 
y
 δ1
 2
δ1l
B3

 l

δ1 δ2l δ3l
 l
l
l 

x
 δ1 δ2 δ3 
δ1l δ2l δ3l
δ2l δ3l
0
δ3l


l

δ2l δ3l 0 +  −δ3
l
l
l
δ2 δ3
δ2 −δ1l



−δ2l

δ1l 

0

 l
δ
 1l
z
 δ1
δ1l

δ2l
δ2l
δ2l





δ3l
0
−δ3l

 l

δ3l 
0
 +  δ3ll
δ3
−δ2  δ1l
δ l δ2l δ3l
 1l
l
l 

w
 δ1 δ2 δ3 
δ1l δ2l δ3l

(3)



δ2l

−δ1l 



0 






E1


E 
 2


E 
 3
.
 ∂l 
B 

1









 B2 
B3

We see that adding constraint terms break the symmetricity of the characteristic
matrix.
• The eigenvalues will be changed as
!
r
c
2
2
x + w ± x − 2xw + w + 4yz (δ1l + δ2l + δ3l ) (1 each),±c (2 each).
2
The zero eigenvalues disappear by adding constraints, and they can be also |c| if
the parameters have the relation (yz − xw − 1)2 = (x + w)2.

80s

90s

Nakamura-Oohara

2000s

Shibata

87

95

99

Nakamura-Oohara

Shibata-Nakamura

Baumgarte-Shapiro

62

G-codeH-codeADM

NCSA

BSSN-code

AEI

97

92

PennState

Alcubierre

Bona-Masso
95-97

99

UWash

ChoquetBruhat-York Anderson-York
01

Cornell-Illinois
96

Kidder-Scheel-Teukolsky
CaltechFrittelli-Reula
Hern

99

LSU

lambda-system
Shinkai-Yoneda

97
86

Ashtekar

Iriondo-Leguizamon-Reula

99

Yoneda-Shinkai

Kidder-Scheel-Teukolsky hyperbolic formulation (Anderson-York + Frittelli-Reula)Phys. Rev. D. 64 (2001) 064017• Construct a First-order form using variables (Kij , gij , dkij ) where dkij ≡ ∂k gij
Constraints are (H, Mi , Ckij , Cklij ) where Ckij ≡ dkij − ∂k gij , and Cklij ≡ ∂[k dl]ij• Densitize the lapse, Q = log(N g −σ )• Adjust equations with constraints∂ˆ0gij = −2N Kij∂ˆ0Kij = (· · ·) + γN gij H + ζN g ab Ca(ij)b∂ˆ0dkij = (· · ·) + ηN gk(i Mj) + χN gij Mk• Re-deining the variables (Pij , gij , Mkij )Pij ≡ Kij + ẑgij K,ˆ j))],
Mkij ≡ (1/2)[k̂dkij + êd(ij)k + gij (âdk + b̂bk ) + gk(i (ĉdj) + db

dk = g ab dkab , bk = g ab dabk

The redefinition parameters– do not change the eigenvalues of evolution eqs.
– do not effect on the principal part of the constraint evolution eqs.
– do affect the eigenvectors of evolution system.
– do affect nonlinear terms of evolution eqs/constraint evolution eqs.

Parabolized ADM (Kidder-Scheel-Teukolsky variables + parabolized)
Paschalidis et al, PRD75(2007)024026 , PRD78(2008)024002 , PRD78(2008)064048• Construct a First-order form using variables (Kij , gij , dkij ) where dkij ≡ ∂k gij
Constraints are (H, Mi , Ckij , Cklij ) where Ckij ≡ dkij − ∂k gij , and Cklij ≡ ∂[k dl]ij• Adjust equations with constraints∂ˆ0gij = (· · ·) + λg ab ∂b Caij
∂ˆ0Kij = (· · ·) + φgij g ab ∂aMb + θ∂(i Mj),
∂ˆ0dkij = (· · ·) + ǫg ab ∂aCbkij + ξgij ∂k H + ζCkij• 6 parameters satisfies the Petrovskii condition for well-posedness if
λ > 0, ǫ > 0, η + φ < 0, θ > 0.
• “parabolic” character ≃ constraint violating mode ∂t gij ≃ λg ab ∂a∂b gij .
• Parabolized-ADM is named since the evolution eqs form a parabolic and mixed hyperbolic.
















2
0
R12(∂x)   u 
0  u   P2 (∂x )
 u 


+

∂t   = v
R21(∂x)
0v
0
P1(∂x )v

• Better performance than KST (one parameter set) and ADM in gauge-wave and Gowdy-wave
tests.

FIG. 7. L 2 norms of the errors for the metric.

FIG. 8. Fourier components of the numerical metric for )FIG. 11. Amplification factor associated with the difference!0,4,8. Some of the components grow exponentially from the very
scheme &12' approximating the ill posed equation v t !i v x .
beginning.

Numerical experiments of KST hyperbolic formulation

lutions the coarsest one gives smaller errors after a while,and the time at which this occurs decreases as one increases
resolution. This indicates that the difference scheme is notconvergent. Note that one could be easily misled to think thatthe code is convergent if one did not evolve the system forlong enough time, or without high enough resolution. For
example, if one performs two runs, with 120 and 240 gridpoints, one has to evolve until, roughly, 150 crossing times,
in order to notice the lack of convergence. To put these numbers in context, suppose one had a similar situation in a 3D
black hole evolution. To give some usual numbers, suppose
the singularity is excised, with the inner boundary at, say, r!M , and the outer boundary is at 20M !which is quite amodest value if one wants to extract waveforms". In this case
120 and 240 gridpoints correspond to grid spacings of, approximately, M /5 and M /10, respectively !usual values aswell in some simulations". If one had to evolve up to 150
crossing times in order to notice the lack of convergence,that would correspond to t!3000M , which is several times
more than what present 3D evolutions last. Of course, thesituation presented in this simple example need not appear inexactly
the same way
in anSARBACH,
evolution AND
of a different
CALABRESE,
PULLIN,TIGLIO spacetime, or with a different discretization; in fact, in the next
subsection we show an example where the instability becomes obvious sooner. Also, there are some ways of noticing
in advance that the code is not converging. Namely, it seems
that the numerical solution has the expected power law
growth that the continuum linearized analysis predicts until
all of a sudden an exponential growth appears. But if one
looks at the Fourier components of the numerical solution,
one finds that there are nonzero components growing exponentially from the very beginning, starting at the order oftruncation error !see Fig. 8".
One might expect that, since for WH systems thefrequency-dependent growth at the continuum is a power law
one, it is possible to get convergence by adjusting the dissipation. In #2$ we show that even though certain amount ofdissipation might help, the code is never convergent and,indeed, adding too much dissipation violates the von Neumann condition, which leads to a much more severe numeriL 2 normsof the errors
for themetric. excal instability.FIG.
We 7.havesystematicallydone
numerical

Weak wave on flat spacetime.
-> No non-principal part.
-> We can observe thefeatures of hyperbolicity.

-> Using constraints in RHS
may improve the blow-up.

periments changing the value of %̃ without being able to reaches the value ˜! !0.32 the instability is even worse thanstabilize the simulations !more details are given below", veri- adding less dissipation, and the same thing happens if one
fying, thus, the discrete predictions.
PHYSICAL REVIEWkeeps
D 66,on064011!2002"˜! beyond 0.32. So we next narrow theincreasing
interval in which the dissipation is fine tuned, we start at ˜!
3. The CIP case (!ÄÀ1Õ2)
Stability properties of a formulationof increase
Einstein’s
equations!0.24, andat intervalsof 0.01. We find the same
Figure 9 shows the error in the metric, for different reso- result; at ˜! "0.25 there is already too much dissipation and†lutions. As in the WH case,Gioel
the errors
originatemostlyfrom* JorgeCalabrese,
Pullin,Olivier Sarbach,‡ and Manuel Tiglio§
the situation is worse. Fine tuning even more, we change ˜!
the nonzero
frequencies!i.e.
ones thatLouisiana
typicallyState
growin
Departmentof Physics
andtheAstronomy,
University,
202 Nicholson Hall, Baton Rouge, Louisiana 70803-4001
in intervalsof 0.001 around 0.250, but it is also found thatan unstable numerical scheme". But!Received
now they
morepublished
27 grow
May 2002;
19 September 2002"!
"0.250the effect of more dissipation is adverse, as
forthan 10 orders of magnitude in much less than one crossing
We study the stability properties of the Kidder-Scheel-Teukolsky
!KST"
also
shown
in
Fig.
14.many-parameter formulation oftime and it is quite obvious that the code is not converging.
Einstein’s equations for weak gravitational waves on flat space-timefromacontinuum and numerical point ofThis is so because in the CIP case the instability grows ex- The fact that beyond ˜! !0.250 the situation becomesview. At the continuum, performing a linearized analysisoftheequationsaroundflat space-time, it turns out
ponentially with the number of gridpoints !see Fig. 10". This worse is in perfect agreement with the discrete analysis ofthat they have, essentially, no non-principal terms. Asaconsequence,
in
the
weak
field limit the stability
can be seen performing a discrete analysis for the single ill "2#. There we show that a necessary condition for the vonproperties of this formulation depend only on the level of hyperbolicity of the system. At the discretelevel weposed equation in 1D, v t !i v x . One gets that the symbol Neumann condition to be satisfied is ˜! $%1/8. Here the uppresent some simple one-dimensional simulations using the KST family. The goal is to analyze the type of& ( ' (x) is real and cannot be bounded by 1 in magnitude, per limit of 1/8 corresponds to, precisely, ˜! !1/4. Exceeding
instabilities that appear as one changes parameter values in the formulation. Lessons learned in this analysis
making the difference scheme unstable !independently of this value results in a violation of the von Neumann condican be applied in other formulations with similar properties.
resolution". If one changed to characteristic variables exactly tion; as explained in "2#, when this happens there is a nuthis equationDOI:
would
appear in 1D as a subset of the system merical instability that grows
with the number
10.1103/PhysRevD.66.064011
PACSexponentiallynumber!s": 04.25.Dm
REVIEW
D 66, is,064011that we are considering,PHYSICAL
so this model
equationin the!2002"of gridpoints &i.e. as in the CIP case', much faster than when

FIG. 13.
stopped oncethe von Ne
growth goe
Finally,
tried with soften used
!0.25, wit
gent simula

Finally h
but now w
thus, the sy
exponential
merical sch

Perform
method yie
shows, forthe lapse, awith two it
I. INTRODUCTION
nonconvergent, but it might appear to converge for a signifi- the minim
cant range of resolutions. We will see that the addition of scheme fordissipation does not fix these problems, but actually can ex- unstable fo
Numerical simulations of the Einstein equations for situacerbate them. It is often the case in numerical relativity that factor andations of interest in the binary black hole problem do not run
discretization schemes that are convergent for strongly hy- have also t
forever. The codes either stop due to the development of
perbolic equations are applied to weakly hyperbolic formu- dissipation
floating point overflows, or even if they do not crash, theylations. The examples of this section will teach us how dan- to confirm
produce answers that after a while are clearly incorrect. It is
gerous such a practice is and confirm the analytic results of WH or CIP
usually very difficult to pinpoint a clear reason why a code
Ref. #2$. This part of the paper is also instructive in that the tions rewrit
stops working. Recently, Kidder, Scheel and Teukolsky
KST system has only been evolved with pseudospectral convergenc!KST" #1$ introduced a twelve-parameter family of evolution
methods. We use ordinary integration via the method of ten as secoequations for general relativity. Performing an empirical pasame initia
lines.
rameter study within a certain two-parameter subfamily, they
The plan of this paper is as follows. In the next section wewere able to evolve single black hole space-times for overwill discuss several notions of stability that are present in the1000 M, where M is the mass of the black hole, something
literature, mostly to clarify the notation. In Sec. III we disthat had been very difficult to achieve in the past. It is of
We havcuss the stability of the KST equations in the continuum
interest to try to understand what makes some of the paramequations ilinearized perturbations. In Sec. IV we discuss the nueter choices
betterothers,in ofparticulargiven metricthat afor )under FIG.
FIG.
8.9.Fouriercomponents
thethenumericalFIG.
L 2 thannormof the errors
formetric.
coordinates
12. L 2 norm of the errors for the metric.
mericalsimulations.
twelve!0,4,8.
dimensionalparameter
study
appears
prohibitiveatSome of the components grow exponentially from the very
lutions the coarsest one gives smaller errors after a 064011-8
while,present.
The intention of this paper is to take some steps inbeginning.
064011-10
and the time at which this occurs decreases as one increasesthis direction. We will first perform a linearized analysis of
II. DIFFERENT DEFINITIONS OF STABILITY
resolution. This indicates that the difference scheme is notthe KST
equationsin the thecontinuum,
small
without being
able toperimentschanging
value ofby%̃ considering
convergent. Note that one could be easily misled to think thatperturbations
aroundflat
space-time.
We
will
observethatthe veri- The term stability is used in numerical relativity in several
stabilize the simulations !more details are given below",
the code is convergent if one did not evolve the system forstabilityof thus,flat space-timeentirely characterized by thedifferent ways. We therefore wanted to make the notation
fying,
the discreteispredictions.
long enough time, or without high enough resolution. Forlevel of hyperbolicity of the system. Since the latter is conclear at least in what refers to this paper. Sometimes theClass. Quantum Grav. 17 (2000) 4799–4822. Printed in the UK

PII: S0264-9381(00)50711-4

Hyperbolic formulations and numerical relativity

4811

Hyperbolic formulations and numerical relativity:
experiments using Ashtekar’s connection variables
Hisa-aki Shinkai† and Gen Yoneda‡† Centre for Gravitational Physics and Geometry, 104 Davey Laboratory, Department of Physics,
The Pennsylvania State University, University Park, PA 16802-6300, USA‡ Department of Mathematical Sciences, Waseda University, Shinjuku, Tokyo, 169-8555, Japan

Hyperbolic
formulations andnumerical relativity
E-mail: shinkai@gravity.phys.psu.edu
and yoneda@mn.waseda.ac.jp

4803

Received 3 May 2000, in final form 13 September 2000

N i , Aa0 , which we call the densitized lapse function, shift vector and the triad lapse function.
Abstract. In order to perform accurate and stable long-time numerical integration of the Einstein
system has threeconstraint equations,
equation, several hyperbolic systems have been proposed. Here we present a numerical comparison
between weakly hyperbolic, strongly hyperbolic and symmetric hyperbolic systems based onASHabi primarycadvantage for using this connection formulation in
Ashtekar’s
connection variables.
The
CH
:=
(i/2)!
≈ 0, variables for all levels of hyperbolicity.
ẼẼtheFsamecaij dynamical
this experiment is that we can keep using b(7)

and we compare thethe simulation by monitoring the violation of the constraints.
ASH
a accuracy
CMi
:= −F
≈of 0,
Ẽaj obtained
By comparing
withfrom the weakly hyperbolic system, we observe that theij results(8)j

Our numerical code demonstrates gravitational wave propagation in plane-symmetric spacetimes,

strongly and symmetric hyperbolic system show better numerical performance (yield less constraint
violation), but not so much difference between the latter two. Rather, we find that the symmetricASHihyperbolic system is not always the best in terms of numerical performance.
i a first to present full numerical simulations using Ashtekar’s variables. We
GaThis study is thealso describe our procedures in detail.

C:= D Ẽ ≈ 0,

(9)

ch are called the
Hamiltonian, momentum and Gauss constraint equations, respectively.
(Some figures in this article are in colour only in the electronic version; see www.iop.org)
Figure 2. Images of gravitational wave propagation and comparisons of dynamical behaviour ofdynamical equations
for a set of (Ẽai , Aai ) are
Ashtekar’s variables and ADM variables. We applied the same initial data of two +-mode pulse
PACS numbers: 0420C, 0425, 0425D∂t Ẽai = −iDj (! cb a N Ẽcj Ẽbi ) + 2Dj (N [j Ẽai] ) + iAb0 !ab c Ẽci ,1. Introductiona∂t Ai = −i!

ab

˜

j c
j a
a
c N Ẽb Fij + N Fj i + Di A0 ,˜

Numerical relativity—solving the Einstein equation numerically—is now an essential field inagravity research. aAs is well aknown, criticalc collapse in gravity systems was first discovered by
Abi A[i Aj ] [1]. Thebc
ijnumerical simulationjcurrent
mainstreamof numerical relativity is to demonstrate thefinal phase of compact binary objects related to gravitational wave observations†, and these
efforts are now again shedding light on the mathematical structure of the Einstein equations.
Up to a couple of years ago, the standard Arnowitt–Deser–Misner (ADM) decompositionof the Einstein equation was taken as the standard formulation for numerical relativists.
Difficulties in accurate/stable long-term evolutions were supposed to be overcome by choosing
proper gauge conditions and boundary conditions. Recently, however, several numericalwaves (a = 0.2, b = 2.0, c = ±2.5 in equation (21) and K0 = −0.025), and the same slicing
condition, the standard geodesic slicing condition (N = 1). (a) Image of the 3-metric component
gyy of a function of proper time τ and coordinate x. This behaviour can be seen identically bothin ADM and Ashtekar evolutions, and both with the Brailovskaya and Crank–Nicholson timeintegration scheme. Part (b) explains this fact by comparing the snapshot of gyy at the same proper
time slice (τ = 10), where four lines at τ = 10 are looked at identically. Parts (c) and (d) are of theyreal part of the densitized triad Ẽ2 , and the real part of the connection A2y , respectively, obtained
from the evolution of the Ashtekar variables.

(10)

(11)

ere F := 2∂− i!
is the curvature 2-form.
We have to consider the reality conditions when we use this formalism
to describe theWhen the pulses collide, then the amplitude seems simply to double, as they are superposed,
and will
the pulses
keep travellingsical Lorentzian spacetime. As we review in appendix A.2, the metricremainon itsin their original propagation direction. That is, we observe
something like solitonic wave pulse propagation.
-valued constraint surface during time evolution automatically if we prepare
initialdata3.2, we have to assume our background not to be flat, therefore
As we mentioned
in section
there theare notriad
exact solutions.
The reader might think that if we set | tr K| to be small and pulse
ch satisfy the reality condition. More practically, we also require thatbe real-4816
Class. Quantum Grav. 17 (2000) 4799–4822. Printed in the UK

H Shinkai and G Yoneda

PII: S0264-9381(00)50711-4

Hyperbolic formulations and numerical relativity:
experiments using Ashtekar’s connection variables
Hisa-aki Shinkai† and Gen Yoneda‡† Centre for Gravitational Physics and Geometry, 104 Davey Laboratory, Department of Physics,
The Pennsylvania State University, University Park, PA 16802-6300, USA‡ Department of Mathematical Sciences, Waseda University, Shinjuku, Tokyo, 169-8555, Japan
E-mail: shinkai@gravity.phys.psu.edu and yoneda@mn.waseda.ac.jp

Hyperbolic formulations and numerical relativity

Received 3 May 2000, in final form 13 September 2000

4815

Abstract. In order to perform accurate and stable long-time numerical integration of the Einstein
5. Experiments
2: another way to control the accuracy/stability
equation, several hyperbolic systems have been proposed. Here we present a numerical comparison
between weakly hyperbolic, strongly hyperbolic and symmetric hyperbolic systems based on
Ashtekar’s connection variables. The primary advantage for using this connection formulation inthis experiment is that we can keep using the same dynamical variables for all levels of hyperbolicity.
Our numerical code demonstrates gravitational wave propagation in plane-symmetric spacetimes,
and we compare the accuracy of the simulation by monitoring the violation of the constraints.
By comparing with results obtained from the weakly hyperbolic system, we observe that thestrongly and symmetric hyperbolic system show better numerical performance (yield less constraint
violation), but not so much difference between the latter two. Rather, we find that the symmetric
hyperbolic system is not always the best in terms of numerical performance.
This study is the first to present full numerical simulations using Ashtekar’s variables. We
also describe our procedures in detail.

The results we have presented in the previous section indicate that both strongly and symmetric
hyperbolic systems show better performance than the original weakly hyperbolic system.
These systems are obtained by adding constraint terms (or ‘adjusted’ terms) to the right-handside of the original equations, (10) and (11). In this section, we report on simple experimentsin changing the magnitude of the multipliers of such adjusted terms.
We consider the following system, where the equations of motion are adjusted in the same(Somefiguresin thisare in colour onlyin the electronic
version; seeκ:www.iop.org)
way as before,butwith
a article
real-valued
constant
multiplier
PACS numbers: 0420C, 0425, 0425D

ASH b
Figure(31)6. Comparisons of the ‘adjusted’ system with the different multiplier, κ, in equations (31)∂t Ẽai = −iDj (# cb a N Ẽcj Ẽbi ) + 2Dj (N [j Ẽai] ) + iAb0 #ab c Ẽci + κP i ab CG,
and (32). The model uses +-mode pulse waves (a = 0.1, b = 2.0, c = ±2.5) in equation (21) in a
˜background K0 = −0.025. Plots are of the L2 norm of the Hamiltonian and momentum constraint
ASH
ASH
where
P i ab ≡ N i δab + iN #ab c Ẽci ,equations, CH
and CM((a) and (b), respectively). We see some κ produce a better performancethanthesymmetric
hyperbolic
system.
1. Introduction
˜
j c
a
ab
j a
a
a ASH
j a ASH∂tNumerical
Ai = −i#
N Fj i +equation
Di A0 numerically—is+ κQi CH now+ κR
CMj field, in(32)
c N Ẽb Fij +relativity—solving
the Einsteinaniessential
˜
gravity research. As is well known, criticalcollapse inagravity systems
was first−2
discovered
byb j Our numerical code demonstrates gravitational wave propagation in plane-symmetricN # ac b Ẽ
where
Qaimainstream≡ e−2 NofẼnumerical
Ri j a ≡.
numerical simulation[1]. The current
relativity
is toiedemonstratethei Ẽcspacetime,i,and we compare the ‘accuracy’ and/or ‘stability’ by monitoring the violation of˜
˜
final phase of compact binary objects related to gravitational wave observations†, and thesetheconstraints.
The
setareofnow
equations(31)light
andon(32)
becomes structure
the original
weaklyhyperbolic system if κActually,= 0, our experiments in section 4 were the comparisons of accuracy
efforts
again shedding
the mathematicalof the Einsteinequations.
in evolutions, while in section 5 we observed cases of unstable evolution. By comparing
Up to thea coupleof years ago,the standard system
Arnowitt–Deser–Misner
becomessymmetric
hyperbolicif κ = 1 and(ADM)
N =decomposition
constant, andwith
remainsa strongly
the results
obtained from the weakly hyperbolic system, we observe that the stronglyof the Einstein equation was taken as the standard formulation for
1 numerical relativists.
hyperbolic
systemforotherchoicesofκ
except
κ
=
which
only
forms
a
weaklyhyperbolic
andsymmetric
hyperbolic
system show better properties with little differences between them.
2 overcome by choosing
Difficulties in accurate/stable long-term evolutions were supposed to be
Therefore, we
may conclude that higher levels of hyperbolic formulations help the numerics
system.
Weconditions
again remarkthat conditions.
the coefficients
constructing
the symmetric
hyperbolic
proper gauge
and boundary
Recently, forhowever,
several numericalNo drastic differences in stability
between 3 levels of hyperbolicity.

BSSN Pros:• With Bona-Masso-type α (1+log), and frozon β (∂tΓi ∼ 0), BSSN plus auxiliary
variables form a 1st-order symmetric hyperbolic system,
Heyer-Sarbach, [PRD 70 (2004) 104004]• If we define 2nd order symmetric hyperbolic form, principal part of BSSN can beone of them,
Gundlach-MartinGarcia, [PRD 70 (2004) 044031, PRD 74 (2006) 024016]

BSSN Cons:• Existence of an ill-posed solution in BSSN (as well in ADM)
Frittelli-Gomez [JMP 41 (2000) 5535]• Gauge shocks in Bona-Masso slicing is inevitable. Current 3D BH simulation islack of resolution.
Garfinke-Gundlach-Hilditch [arXiv:0707.0726]

strategy 2Hyperbolic formulation (cont.)

Are they actually helpful?
“YES” group“Well-posed!”, ||u(t)|| ≤ eκt||u(0)||Mathematically Rigorous Proofs
IBVP in future

Initial Boundary Value Problem
Consistent treatment is available
only for symmetric hyperbolic
systems.
GR-IBVP
Stewart, CQG15 (98) 2865

Tetrad formalism
Friedrich & Nagy, CMP201 (99) 619

Linearized Bianchi eq.
Buchman & Sarbach, CQG 23 (06) 6709

Constraint-preserving BC
Kreiss, Reula, Sarbach & Winicour, CQG 24 (07) 5973

Higher-order absorbing BC
Ruiz, Rinne & Sarbach, CQG 24 (07) 6349

strategy 2Hyperbolic formulation (cont.)

Are they actually helpful?
“YES” group“Really?” group“Well-posed!”, ||u(t)|| ≤ eκt||u(0)||“not converging”, still blow-up

Mathematically Rigorous Proofs

Proofs are only simple eqs.
Discuss only characteristic part.
Ignore non-principal part.
...

IBVP in futurestrategy 2Hyperbolic formulation (cont.)

Are they actually helpful?
“YES” group“Really?” group“Well-posed!”, ||u(t)|| ≤ eκt||u(0)||“not converging”, still blow-up

Mathematically Rigorous Proofs

Proofs are only simple eqs.
Discuss only characteristic part.
Ignore non-principal part.
...

IBVP in future
Which level of hyperbolicity is necessary?

symmetric hyperbolic ⊂ strongly hyperbolic ⊂ weakly hyperbolic systems,
Advantages in Numerics (90s)
Advantages in sym. hyp.
– KST formulation by LSU

strategy 2Hyperbolic formulation (cont.)

Are they actually helpful?
“YES” group“Really?” group“Well-posed!”, ||u(t)|| ≤ eκt||u(0)||“not converging”, still blow-up

Mathematically Rigorous Proofs

Proofs are only simple eqs.
Discuss only characteristic part.
Ignore non-principal part.
...

IBVP in future
Which level of hyperbolicity is necessary?

symmetric hyperbolic ⊂ strongly hyperbolic ⊂ weakly hyperbolic systems,
Advantages in Numerics (90s)

These were vs. ADM

Advantages in sym. hyp.
– KST formulation by LSU

Not much differences in hyperbolic 3 levels– FR formulation, by Hern
– Ashtekar formulation, by HS-Yoneda
sym. hyp. is not always the best

80s

90s

Nakamura-Oohara

2000s

Shibata

87

95

99

Nakamura-Oohara

Shibata-Nakamura

Baumgarte-Shapiro

62

G-codeH-codeADM

NCSA

BSSN-code

AEI

97

92

PennState

Alcubierre

Bona-Masso
95-97

99

UWash

ChoquetBruhat-York Anderson-York
01

Cornell-Illinois
96

Kidder-Scheel-Teukolsky
CaltechFrittelli-Reula
Hern

99

LSU

lambda-system
Shinkai-Yoneda

97
86

Ashtekar

Iriondo-Leguizamon-Reula

99

Yoneda-Shinkai

90sNakamura-Oohara

2000s

Shibata

87

95

99

Nakamura-Oohara

Shibata-Nakamura

Baumgarte-Shapiro

62

G-codeH-codeADM

NCSA

BSSN-code

AEI

97

92

PennState

Alcubierre

01

Bona-Masso

Illinois

95-97

99

adjusted-system

80s

UWash

ChoquetBruhat-York Anderson-York
01

Cornell-Illinois
96

Kidder-Scheel-Teukolsky
CaltechFrittelli-Reula
Hern

99

LSU

lambda-system
Shinkai-Yoneda

97
86

Ashtekar

Iriondo-Leguizamon-Reula

99

Yoneda-Shinkaiここまでのまとめ (前半)
[Keyword 1]
Formulation Problem
数学的に等価な式であっても，数値的には安定性が異なる．
[Keyword 2]
ADM formulation
歴史的に ADM 形式からスタートした．90 年代までは成功していた．
BH-BH/NS-NS 連星シミュレーションで問題発生[Keyword 3]
BSSN formulation
ADM に新変数を導入してゲージも一部 ﬁx．
どうして安定性が改善されたのかは謎のまま．Technical tips がいろいろ．
[Keyword 4]
hyperbolic formulations
PDE の適切性(wellposedness)という根拠はあるが限定的．
Einstein 方程式の双曲形式はいくらでも作れる．
BSSN の利点を双曲形式で説明しようというグループもあるが，．．
But are they really helpful in numerics?

2008/11/28 @ 大阪市立大学

Formulation for Numerical Relativity
both for Einstein / Gauss-Bonnet
真貝寿明

大阪工業大学情報科学部

1. Introduction
定式化問題？
2. The Standard Approach to Numerical Relativity
ADM 形式，BSSN 形式，Hyperbolic 形式
3. Robust system for Constraint Violation
Adjusted systems .... better than lambda system!
Adjusted ADM
Adjusted BSSN
4. 高次元数値相対論に向けてstrategy 3 “Asymptotically Constrained” system /“Constraint Damping” system
Formulate a system which is “asymptotically constrained” against a violation of constraints
Constraint Surface as an Attractor
method 1: λ-system (Brodbeck et al, 2000)• Add aritificial force to reduce the violation of
constraints• To be guaranteed if we apply the idea to a symmetric hyperbolic system.
t=0

method 2: Adjusted system (Yoneda HS, 2000,
2001)

Constrained Surface(satisfies
Einstein's constraints)• We can control the violation of constraints byadjusting constraints to EoM.
• Eigenvalue analysis of constraint propagation
equations may prodict the violation of error.

Blow uperror?
Stabilize?

time• This idea is applicable even if the system is not
symmetric hyperbolic. ⇒for the ADM/BSSN formulation, too!!

Idea of λ-system
Brodbeck, Frittelli, Hübner and Reula, JMP40(99)909
We expect a system that is robust for controlling the violation of constraints
Recipe1. Prepare a symmetric hyperbolic evolution system∂t u = J∂iu + K
2.

Introduce λ as an indicator of violation of constraint
which obeys dissipative eqs. of motion

3.

Take a set of (u, λ) as dynamical variables

4.

Modify evolution eqs so as to forma symmetric hyperbolic system∂t λ = αC − βλ(α 6= 0, β > 0)






u
A
0
u
∂ 

∂t   ≃ 
i
F 0
λ
λ












u
A F̄   u ∂t   = 
∂i
λ
F 0
λ

Remarks• BFHR used a sym. hyp. formulation by Frittelli-Reula [PRL76(96)4667]• The version for the Ashtekar formulation by HS-Yoneda [PRD60(99)101502]
for controlling the constraints or reality conditions or both.
• Succeeded in evolution of GW in planar spacetime using Ashtekar vars. [CQG18(2001)441]• Do the recovered solutions represent true evolution? by Siebel-Hübner [PRD64(2001)024021]• The version for Z4 hyperbolic system by Gundlach-Calabrese-Hinder-MartinGarcia [CQG22(05)3767]⇒ Pretorius noticed the idea of ”constraint damping” [PRL95(05)121101]

INSTITUTE OF PHYSICS PUBLISHING

CLASSICAL AND QUANTUM GRAVITYClass. Quantum Grav. 18 (2001) 441–462

www.iop.org/Journals/cq

446

G Yoneda and H Shinkai

PII: S0264-9381(01)15923-X

Hyperbolic formulations and numerical relativity: II.
yperbolic formulations
and numerical relativity:II
445
asymptoticallyconstrained
systems of Einsteinequations
ith the initial data λ = λ = 0 and take (E, B, λ , λ ) as a set of variables to evolve:
E

B

E

B

 j 
i l
E Gen Yoneda
0 1 and−c#
0 02
E0j ShinkaiHisa-aki B i   c# i l Bj   0 
0 Sciences,
0 Waseda0 1  ofjMathematical Shinjuku, JapanUniversity,
Tokyo, 169-8555,∂t  2 Department=+∂for Gravitational
 Lab.,
Department
.The
104l DaveyCentrePhysics and Geometry,of Physics,l λE Pennsylvania
 α1State




−β
δ j University,
0 University
0 Park,
0
λ
λ
E
1
E
PA 16802-6300, USA


i 





Hyperbolic
λB formulations0and numerical
α2 δ l j relativity:
0 0 II

λB
E-mail: yoneda@mn.waseda.ac.jp and shinkai@gravity.phys.psu.edu(2.14)
445−β2 λBWe immediately
obtainan
expectedsymmetric
form
as
with
the initial
data
λB in
= final
0 and
take(E, B, λE2000, λB ) as a set of variables to evolve:
Received
27λJuly2000,
form
13 December
E = i    j il    j  
i l  i li
αE
0 0
0
−c#
−c#
E
0
E
E j 0 00j1δAbstract
 lic#i l  systems B i   c# i j l We0study numericalB j 
j  for
B0i 
0
0 
j  B0  0asymptotically α2δ constrainedintegration of the0 

  


∂t equations,(2.14)
Einsteinwhich
perturbative
errors (2.15)
 =
 are
.
 ∂l . +
∂l  tobe+robust = againstlintended
 λevolution
 −β1 λE δ initial00 −β0 1examine
0E  
0ofαthej λ1 λE   α1 δjl
λEλEthepreviouslyfor0the freedata.
weE  First,
λB

Maxwell-lambda
systemworks
proposed ‘λλsystem’, which0 introduces
artificial
flows
surfacesαδ00 λto constraintλ−β λ0 hyperbolic λformulation.−β0α δ on the0symmetricbased
We show that this system
asWeexpected.
immediately
anfor
expected
as in the Maxwell system
works asobtain
expected
the wavesymmetricpropagation form
problem
l
2 j

B

l
2B j

2 B

B

2 B

relativity
connectionand ingeneral  Second,  formulation.
2.2. Analysis ofeigenvalues.
Nowthei j evolution
Cwel using Ashtekar’s
E and CB
α1control
δ li equations
0 stability,for−c#
E i propose a0new
E jtheweconstraints
mechanism
to
the
which
call the0‘adjusted
ecome
i l is simply obtained by addingli constraint
 i system’.
 This j terms 
in the dynamical00α2 δ   B   0 
 B   c# jandadjusting their multipliers. We
whyparticular choice
=('λ
∂∂t C
∂explain =αequations ), +a 
. (2.16)ll)t EλE  1  α
E1 δt CB =0α2 ('λ0∂the0B from λE  orpure-imaginary
 non-positive−β1 λE j reducesof multipliernumerical errors(2.15)

eigenvalues
constraintThis
here ' = ∂i ∂ i . We
the Fourier
and2‘adjusted
λB takeλequations.
λλs,δjl
0 for constraints
0propagation Cs
0 of theαadjusted
B (2.16) −β
B (2.12) and2integrals
system’ is also tested in
the Maxwell system and in the Ashtekar system. This
2.13), in the form of (2.7),
to obtain
mechanism affects more than the system’s symmetric hyperbolicity.

of eigenvalues. Now  equations

2.2.2.ĈAnalysis
thefor the constraints CE and CB
0 evolution
0 0 −α1 k 2
ĈE
E
PACS numbers: 0420C, 0425, 0425D
become Ĉ   0 0
  Ĉ −αin2 kthe2 electronicB  see www.iop.org)
 B (Some figures in this article are0in colour only  version;∂t (2.17) (2.16)
=,)∂t CE = α1 ('λE ),∂t CB=α2 ('λ λ̂E   α1 0
−β10   λ̂E B

Figure 1. Demonstration of the λ system in the Maxwell equation. (a) Constraint violation (L2where ' = ∂i ∂ i . We take the Fourier integrals for constraints Cs (2.16) and λs, (2.12) and
0 α2
λ̂B
0
−β2
λ̂B
norm of CE ) versus time with constant β (= 2.0) but changing α. Here α = 0 means no λ system.
(2.13), in the form of (2.7), to obtain(b) The same plot with constant α (= 0.5) but changing β. We see better performance for β > 0,
which is the case of negative eigenvalues of the constraint propagation equation. The constants in 
 expression.
here k 2 = ki k i . We find the matrixNote
2 that thisisan exact(2.18) were chosen as A = 200 and B = 1.
ĈE to be 0constant.
0
0 −α
ĈE
1k
1. Introductionince the eigenvalues are Ĉ   0 0

2 0−α2 k   ĈB B( ∂ant '= solving,(2.17)
Numerical
relativity,
approachthe Einstein equations numerically,)to  issupposed λ̂2 2  α20−β0   λ̂ where αi "= 0 (possibly complex) and βi > 0 (real) are constants.
If we take yα := (Ẽai , Aai , λ, λi , λa ) as a set of dynamical variables, then the principal
Hyperbolic formulations and numerical relativity: II451
G
YonedaH Shinkai
CLASSICALAND
UANTUMGRAVITY can be written as
partof Qand(2.26)–(2.28)

448 INSTITUTE OF PHYSICS PUBLISHING

448

Class.
Quantum Grav.
(2001) 441–462
www.iop.org/Journals/cq
PII: S0264-9381(01)15923-X
bcd j lWe introduce
new18variables(λ, λi , λa ), obeying
the dissipative evolution
equations, ∂ λ ∼Ẽc Ẽd (∂l Abj ),t = −iα1 %(2.29)G Yoneda andH Shinkai(2.26)∂t λ = α1 CH − β1 λ,jj∂t λi ∼
(2.30)= α2 [−eδil Ẽb + eδi Ẽbl ](∂l Abj ),α2 C˜Mi −(λ,β2 λλii,, λa ), obeying the dissipative evolution equations,(2.27)t λi =
We introduce ∂newvariablesl∼∂II.
(2.31)t λa = α3 ∂l Ẽa .
β 3 λa ,(2.28)
Hyperbolic
andnumerical
relativity:t λaα=Cαformulations
3C
Gaβ− λ,∂t λ∂=−(2.26)
1 H
1asymptoticallyconstrained
The characteristic matrix of the system uα is not Hermitian. However, if we modify
where
αi "= 0 (possibly ˜complex)
and βi > 0 (real)systems
are constants.of Einstein∂λ=α(2.27)
Ci −aβ2 λi ,
t
i
2
Mi
If
we
take
y
:=
(
Ẽ
,
A
,
λ,
λ
,
λ
)
as
a
set
of
dynamical
variables,then
the principal
theright-handsideof the evolution equation of (Ẽai , Aai ), then the set becomes a symmetric
α
i
a
a
i
equationspart of (2.26)–(2.28) can be written as
hyperbolic system.
This is done by adding ᾱ 3 γ il (∂l λa ) to the equation of ∂t Ẽai , and by adding∂t λa = α3 CGa − β3 λa ,(2.28)
c
l
m la
a
lm a
iᾱ 1 % a c d Ẽi Ẽ(2.29)∂t λ ∼= −iα1 % bcd Ẽ1 cj Ẽdl (∂l Abj ),2d (∂l λ)+ ᾱ 2 (−eγ Ẽi +eδi Ẽ )(∂l λm ) to the equation of ∂t Ai . The final principal
Gen
YonedaandHisa-aki
Shinkai
where αi "= 0 (possibly complex) and βi > 0 (real) are constants.
part is then written asjb∼
If we take yα ∂:=(Department
Ẽαai2,[−eδAai ,ilofλ,
λ , λj ) lasSciences,al Asetof dynamical variables,
principal(2.30)
ẼMathematicalt λi1 =j ), Waseda University, Shinjuku, Tokyo,169-8555,
Japanb +i eδiaẼb ](∂ thenthe
2 CentreẼai of Physics, Thepart of (2.26)–(2.28) can
be forwrittenas Physics and Geometry, 104 Davey Lab., Department
Gravitational∼∂t λaPennsylvania. University, University Park, PA 16802-6300, USA  a  (2.31)= α3 ∂l Ẽal State∼ Ai 
∂t λ = −iα1 % bcd Ẽcj Ẽdl (∂l Abj ), λ we modify(2.29)∂
E-mail:yoneda@mn.waseda.ac.jp
shinkai@gravity.phys.psu.edu
t  if
The characteristic
matrixof the system uandnot Hermitian. However,α is λa symmetricj lthe right-hand sideof the evolution
equationof b(Ẽai , Aai ), then the set becomesl j∼i∂t λi =Received
α2 [−eδi27+ eδ2000,),(2.30)
ẼbJulylA
i Ẽb ](∂
infinal
il j form 13 December 2000 i
hyperbolic system. This is done by adding ᾱ 3 γ (∂l λa ) to the equation of ∂t Ẽa , andλa by adding
l lm Ẽia +eδim Ẽ la )(∂l λm ) to the equation of ∂t Aai . The final∼ αᾱ 2∂(−eγiᾱ 1 % a c d Ẽic Ẽ∂dl λ(∂l λ)+
principal
t a =Abstract
3 l Ẽa .
Ml a bi(2.31)000ᾱ 3 γ il δa bjpart is then written as
We study asymptotically constrained systems for numerical integration of thej characteristicThematrixequations,of the systemuαintended
is not Hermitian.
However,
if we0errors
modifyN l ai bj
iᾱ 1 % a c d Ẽic Ẽdl ᾱ 2 e(δi Ẽ la − γ lj Ẽia )0
Einsteinwhich areto be robust against
perturbative


Ẽai
a
i


the right-handthetheevolution
equationof initial(Ẽa , Adata.
thewesetexamine
becomes
apreviously
symmetric Aa  side offor
free evolutionof theFirst,thej l∼i ), then
cd


 i
=
Ẽ
Ẽ
0
−iα
%
0
0
0
c
1
b
ilintroduces artificial flows to constrainti  surfaces
d
proposed‘λbysystem’,which
hyperbolic
system.
This
is
done
adding
ᾱ
γ
(∂
λ
)
to
the
equation
of
∂
Ẽ
,
and
by
adding
∂t 
3
l
a
t
λ

a 


j l
j
l


l
a
m
a
a d c lmlabased on thesymmetric
hyperbolic formulation. We show that this system0α2 e(δi Ẽb − δi Ẽb )000iᾱ 1 % c ẼiλẼ
Ẽi +eδi d (∂l λ)+ ᾱ 2 (−eγi Ẽ )(∂l λm ) to the equation of ∂t Ai . The final principal
systemλa written asworks as expected for the wave propagation problem in the Maxwell
part is then
α3 δab δjlwe0000
and in general relativity using Ashtekar’s connection formulation. Second,l bi
il b 
 ij
M a j propose a new
0 mechanism to 0control the stability,0 which we callᾱ 3 γtheδ‘adjusted
a
Ẽa
Ẽb
j la


c
l
a
l
a
j
a
d
lj
system’.
Thisissimplyobtainedbyaddingconstraint
terms
in
the
dynamical
 Aa   0
N ibiᾱ 1 % c Ẽi Ẽd ᾱ 2 e(δi Ẽ − γ Ẽi )0  Ab  i ∼   j
equations andadjustingtheir multipliers. We explain why a particular choice
j l


cd


∂t  λ  =  0−iα1 %b Ẽc Ẽd000×∂l  λ 
(2.32)
,ofmultiplierreducesthenumerical
errorsfrom
non-positiveorpure-imaginary λ  
j l
j


l

i
λj
0 eigenvalues
α2 e(δi Ẽb of− δthe)00 equations. This0 ‘adjusted
constraint
propagationi Ẽbadjusted
λab l system’ is also tested in the Maxwell system and in the Ashtekar system. This λα3 δ a δ j
0
0
0
0
b
 l bimechanism
affects more than the system’s symmetric hyperbolicity.

il b 
j
Ma j
00 whereᾱ 3 γ δa
Ẽb 0


j
b
 0
c l
la j
a d0425D
lj a
PACS
numbers:
0420C,iᾱ0425,
Figure 3. Demonstration of the λ system in the Ashtekar equation. We plot the violation of theᾱ 2 e(δi Ẽ la − γlabij
Ẽi ) abc0 l  Aj Nij
1 % c Ẽi ẼdibM= i%(2.32)
NẼc γ
+ N l γ ij δ ab , constraint (the L2 norm of the Hamiltonian constraint equation, CH ) for the cases of(2.33)
plane-wave


×∂
,
l
λ
∼
 figures
 cd in this(Some
are in colour only in the electronic version; see www.iop.org)
j article∼l
propagation under the periodic boundary. To see the effect more clearly, we added an artificial error
=


Ẽ
Ẽ
0
−iα
%
0
0
0
c d


λj 1 bat t = 6. Part (a) shows how the system goes bad depending on the amplitude of artificial error.
labij
abc j li
abc l j Thei error−2l are
j b of the evolution of Ashtekar’s
j lformẼ
Aj2y →lA2y (1 + −2
AllẼ
thei curves
l j
λ
b
 0

= iN
Ẽ
Ẽ
(%
γ
−
%
γ
− e wasẼofiathe
% bcd
− (b)
e error).
% acd
α2 e(δi Ẽb − δi Ẽb )
0
0N
0
c
c
c ẼdPartd Ẽofc Ẽ
∼original equation (no λ system).
shows the effectthe λ system. All the curves have

Ashtekar-lambda system works
as expected, as well.

where α3 δ b δ l

a j0Ẽ l γjij+ N l γ ij δ ab ,Mlabij = i% abc N
c bẼ
∼
1. Introduction

0020% error amplitude, but show the difference of the evolution equations. The full curve is for0
original
equation (the same as in (a)), the dotted curve is for the strongly hyperbolic−2 acd i j lb Ashtekar’s+e % Ẽd Ẽc Ẽ ) +Ashtekar
N l δ abequation.
γ ij . Other curves are of λ systems, which produce a better performance(2.34)
than that(2.33)of the strongly hyperbolic system.

aaii Ab 
abc
Ẽ j γj li− % abc Ẽ l γ j i − e−2 Ẽ ia % bcd Ẽ j Ẽ l − e−2 % acd Ẽ i Ẽ l Ẽ j b Clearly, the solution (Ẽa , Ai , λ, λi , λa ) = (Ẽa , Ai , 0, 0, 0) represents the original
N labij
=
iN
(%
c
c
c
d
d
c


Numerical
supposed∼×∂ relativity, ,solution ofisthe
Ashtekar(2.32) system.
l  λ  an approach to solving the Einstein equations numerically,

Idea of “Adjusted system” and Our Conjecture
CQG18 (2001) 441, PRD 63 (2001) 120419, CQG 19 (2002) 1027
General Procedure1.

prepare a set of evolution eqs.

∂t ua = f (ua , ∂bua , · · ·)2.

add constraints in RHS

a
a
(C
,
∂
C
, · · ·)}
∂t ua = f (ua , ∂bua , · · ·) +F
b
|
{z

3.

choose appropriate F (C a, ∂bC a, · · ·)to make the system stable evolution

How to specify F (C a, ∂b C a, · · ·) ?
4.

prepare constraint propagation eqs.

∂t C a = g(C a, ∂b C a, · · ·)5.

and its adjusted versionaa∂t C a = g(C a, ∂b C a, · · ·) +G(C,∂C
, · · ·)}
b
|
{z

6.

a
Fourier transform and evaluate eigenvalues ∂tĈ k = A(Ĉ) Ĉ k| {z }

Conjecture: Evaluate eigenvalues of (Fourier-transformed) constraint propagation eqs.
If their (1) real part is non-positive, or (2) imaginary part is non-zero, then the system is more stable.

Example: the Maxwell equations
Yoneda HS, CQG 18 (2001) 441
Maxwell evolution equations.
jk∂tEi = c"i ∂j Bk + Pi CE + Qi CB ,∂tBi = −c"ijk ∂j Ek + Ri CE + Si CB ,CE = ∂iE i ≈ 0, CB = ∂iB i ≈ 0,

Constraint propagation equations

CAFs?











sym. hyp⇔ Pi = Qi = Ri = Si = 0,strongly hyp ⇔ (Pi − Si)2 + 4RiQi > 0,





 weakly hyp⇔ (Pi − Si)2 + 4RiQi ≥ 0∂tCE = (∂iP i)CE + P i(∂iCE ) + (∂iQi)CB + Qi(∂iCB ),(∂iRi)CE + Ri(∂iCE ) + (∂iS i)CB + S i(∂iCB ),∂tCB = sym. hyp⇔ Qi = Ri,
strongly hyp ⇔ (Pi − Si)2 + 4RiQi > 0,





 weakly hyp⇔ (Pi − Si)2 + 4RiQi ≥ 0Ĉ∂t  E  =
ĈB



















ĈE 
∂ P i + P iki ∂iQi + Qiki   ĈE   P iki Qiki   ĈE 
 i

∂
≈
=:
T
l
ĈB
ĈB
ĈB
∂iRi + Riki ∂iS i + S iki
R i k i S i ki
ii)

⇒ CAFs = (P ki + S ki ± (P iki + S iki)2 + 4(QikiRj kj − P ikiS j kj ))/2

Therefore CAFs become negative-real when
P iki + S iki < 0,and

QikiRj kj − P ikiS j kj < 0

INSTITUTE OF PHYSICS PUBLISHINGClass. Quantum Grav. 18 (2001) 441–462

CLASSICAL AND QUANTUM GRAVITY
www.iop.org/Journals/cq

PII: S0264-9381(01)15923-X

Hyperbolic formulations and numerical relativity: II.
asymptotically constrained systems of EinsteinequationsHyperbolic formulations and numerical relativity: II

455

Gen Yoneda1 and Hisa-aki Shinkai2
1 Department of Mathematical Sciences, Waseda University, Shinjuku, Tokyo, 169-8555, Japan
2 Centre for Gravitational Physics and Geometry, 104 Davey Lab., Department of Physics, The

454

Pennsylvania State University, University Park, PA 16802-6300, USA

G Yoneda and H Shinkai

E-mail: yoneda@mn.waseda.ac.jp
and shinkai@gravity.phys.psu.edu
a symmetric hyperbolic
system (e.g. [2–5, 10]).
We believe, however, that the above guidelines
yield the essential mechanism for our purpose of constructing a stable dynamical system.
Received 27 July 2000, in final form 13 December 2000
In the following subsections and appendix A, we demonstrate that this adjusted system
actually worksAbstract
as desired in the Maxwell system and in the Ashtekar system of the Einsteinequations, in which
the above
two guidelines
are applied,
respectively.
We study
asymptoticallyconstrained
systems
for numerical integration of theEinstein equations, which are intended to be robust against perturbative errors
Adjusted-Maxwell
systemas well.
3.2. Example for
1: Maxwell
the free equations
evolution of the initial
data. First, weworks
examine the previously
proposed
system’,whichconsider
introduces
to constraint
surfaces
3.2.1. Adjusted
system.‘λ Herewe againtheartificial
Maxwellflows
equations(2.9)–(2.11).
We
based
on
the
symmetric
hyperbolic
formulation.
We
show
that
this
systemstart from the adjusted dynamical equations
works as expected for the wave propagation problem in the Maxwell systemj∂t Eini =c"i j k ∂relativity
Qi CB + q jformulation.
j Bk + P i C
E + pAshtekar’si (∂j CE ) +connectioni (∂j CB ),andgeneralusing
Second,(3.7)
we
j kmechanism to controljjproposeanew
thestability,
which
we
call
the
‘adjusted(3.8)∂t Bi = −c"i ∂j Ek + Ri CE + r i (∂j CE ) + Si CB + s i (∂j CB ),system’. This is simply obtained by adding constraint terms in the dynamical
where P , Q, equations
R, S, p, q,andr ands are theirmultipliers.
Thesedynamical
adjusttheadjustingmultipliers.
We explain
whyequations
a particular
choiceFigure 4. Demonstrations of the adjusted system in the Maxwell equation. We perform the same
constraint propagation
equationsasexperiments with section 2.2.3 (figure 1). Constraint violation (L2 norm of CE ) versus time areof multiplier reduces the numerical errors from non-positive or pure-imaginaryiiplotted for various κ (= p j i = s j i ). We see that κ > 0 gives a better performance (i.e. negative∂t CE = (∂i P ieigenvalues)CE + P i (∂i C) + (∂
Q
)C
+
Q
(∂
C
)
ofEtheadjusted
constraint
propagation
equations.
This‘adjustediB
i B
real part eigenvalues for the constraint propagation equation), while excessively large positive κsystem’testedjin
and in the
i the Maxwell systemji
j i Ashtekar system. This
makes the system divergent again.
)(∂also(3.9)+(∂i p j iis
j CE ) + p (∂i ∂j CE ) + (∂i q )(∂j CB ) + q (∂i ∂j CB ),mechanism affects more than the system’s symmetric hyperbolicity.
∂t CB = (∂i R i )CE + R i (∂i CE ) + (∂i S i )CB + S i (∂i CB )PACS numbers:
0420C, 0425, 0425D(3.10) 3.3. Example 2: Einstein equations (Ashtekar equations)+(∂i r j i )(∂j CE ) + r j i (∂i ∂j CE ) + (∂i s j i )(∂j CB ) + s j i (∂i ∂j CB ).
(Some figures in this article are in colour only in the electronic version; see www.iop.org)

3.3.1. Adjusted system for controlling constraint violations. Here we only consider theThis will be expressed using Fourier components by
!
" !
i
i
ji
ji
i
i
ji
ji "
adjusted system which controls the departures from the constraint surface. In the appendix,ĈE∂i P + iP ki + ikj (∂i p ) − ki kj p
∂i Q + iQ ki + ikj (∂i q ) − ki kj q∂t=
we present an advanced system which controls the violation of the reality condition together
ĈB∂i R i + iR i ki + ikj (∂i r j i ) − ki kj r j i∂i S i + iS i ki + ikj (∂i s j i ) − ki kj s j iwith a numerical demonstration.
"
!
"
!
ĈE
ĈE
Even if we restrict ourselves to adjusted equations of motion for (Ẽai , Aai ) with constraint=: T.
(3.11)
×
1. Introduction
terms (no adjustment with derivatives of constraints), generally, we could adjust them as
ĈBĈB
Assuming the multipliers are constants or functions of E and B, we can truncate the principalicbẼcj Ẽbi ) + 2Dj (N [j Ẽai] ) + iAb0 #ab c Ẽci + Xai CH + Yaij CMj + Paib CGb ,Numerical relativity, an approach to solving the Einstein equations numerically, is supposed ∂t Ẽa = −iDj (# a N∼
matrix as

Example: the Ashtekar equations
HS Yoneda, CQG 17 (2000) 4799
Adjusted dynamical equations:
j i
[j i]
b
c i
i
ij
ib
)
+
2D
(N
)
+
iA
"
+X
C
+
Y
C
+
P
CGb
Ẽ
Ẽ
Ẽ
Ẽ
∂tẼai = −iDj ("cb aN
j
H
M
j
c
b
a
0
ab
c
a
a
a
∼
!

"#

adjustẼbj Fijc + N j Fjia + DiAa0 + ΛNẼia +QaiCH + Riaj CM j + ZiabCGb∂tAai = −i"abcN∼∼!

Adjusted and linearized:"#

adjust$$Ẽia), Raj i = κ3(−ie−2N"acdẼidẼcj )X = Y = Z = 0, Pbia = κ1(iN iδba), Qai = κ2(e−2N∼∼
Fourier transform and extract 0th order of the characteristic matrix:ˆH  
ˆH 
C
C
0
i(1
+
2κ
)k
0
3
j










kj
ˆ
ˆ






i(1
−
2κ
C
C
)k
κ
"
k
0
=
∂t 




Mi 
2 i
3
i k
Mj 




0
2κ3δ j
CˆGa
CˆGb
0
a

Eigenvalues:+

0, 0, 0, ±κ3,,−kx2 − ky 2 − kz 2, ±(−1 + 2κ2)(1 + 2κ3)(kx2 + ky 2 + kz 2)

In order to obtain non-positive real eigenvalues:(−1 + 2κ2)(1 + 2κ3) < 0-

INSTITUTE OF PHYSICS PUBLISHINGClass. Quantum Grav. 18 (2001) 441–462

CLASSICAL AND QUANTUM GRAVITY
www.iop.org/Journals/cq

456

G Yoneda and H Shinkai

PII: S0264-9381(01)15923-X

Hyperbolic formulations and numerical relativity: II.
asymptotically constrained systems of Einsteinequations
Gen Yoneda1 and Hisa-aki Shinkai2
1 Department of Mathematical Sciences, Waseda University, Shinjuku, Tokyo, 169-8555, Japan
2 Centre for Gravitational Physics and Geometry, 104 Davey Lab., Department of Physics, The

Figure 4. Demonstrationsof theadjusted Park,
systemthe MaxwellUSA
equation. We perform the same
Pennsylvania
State University,
University
PAin16802-6300,
experiments with section 2.2.3 (figure 1). Constraint violation (L2 norm of CE ) versus time areplotted forvarious κ (= p j i = s j iand). We
see that κ > 0 gives a better performance (i.e. negative
E-mail:yoneda@mn.waseda.ac.jp
shinkai@gravity.phys.psu.edu
real part eigenvalues for the constraint propagation equation), while excessively large positive κ
makes the system divergent again.

Received 27 July 2000, in final form 13 December 2000

Abstract
Adjusted-Ashtekarsystem works as well.
3.3. Example 2: Einstein equations (Ashtekar equations)
We study asymptotically constrained systems for numerical integration of theEinsteinwhich constraint
are intended
to be robustHere
againsterrors
3.3.1. Adjusted
systemequations,for controlling
violations.
we perturbative
only consider
thethe freeevolutionof the initial
First, wesurface.
examineIntheadjusted systemforwhich
controls
the departures
fromdata.
the constraint
thepreviously
appendix,
proposed system‘λ system’,introduces
artificial
flows
to constraint
we present an advanced
whichwhich
controls
the violationof therealitycondition surfaces
togetherbasedon the symmetric hyperbolic formulation. We show that this systemwith a numericaldemonstration.
expectedtofor
the wavepropagation
probleminai thesystemEven if weworks
restrictasourselves
adjusted
equationsof motionfor (Ẽ, AaiMaxwell) with constraint
andingeneralrelativity
using
Ashtekar’s
connection
formulation.
Second,terms (no adjustment with derivatives of constraints), generally, we could adjust them
as wepropose a new mechanism to control the stability, which we call the ‘adjusted
cb
b
c i
ij
ib
Ẽcj ThisẼbi ) +is2simply
Dj (N [jobtainedẼai] ) + iAby
Ẽc +constraint
Xai CH + Yterms+ Pthe∂t Ẽai = −iDj (#system’.
Gb ,addingaN
0 #ab
a CMj ina Cdynamical∼
equations and adjusting their multipliers. We explain why a particular choice(3.14)
reducesthe numericalor pure-imaginaryj cjaaab of multiplierj aaaerrorsa from non-positiveab∂t Ai = −i# ceigenvaluesẼb Fij + NofFthe
Ẽi + QpropagationN
Di A0 + $N
, ‘adjusted(3.15)
Mj + Zi CGb
j i +adjusted
i CH + Ri Cequations.
∼∼constraintThis
is alsotested inij system’
aj the Maxwell system and in the Ashtekar system. Thisawhere Xai , Ya ,mechanism
Ziab , Paib , Qaffects
Ri areorder to simplify thei and morethanmultipliers.
the system’s However,symmetricin
hyperbolicity.
discussion, we restrict multipliers so as to reproduce the symmetric hyperbolic equationsof motion [10, 11],i.e.numbers: 0420C, 0425, 0425D
PACS(Some figures in this article are in colour only in the electronic version; see www.iop.org)

X = Y = Z = 0,

Paib = κ1 (N i δab + iN
#a bc Ẽci ),∼Qai = κ2 (e−2 NẼia ),∼
1. IntroductionRi

ja(3.16)= κ3 (ie−2 N# ac b Ẽib Ẽcj ).
∼Numerical relativity, an approach to solving the Einstein equations numerically, is supposed

Figure 5. Demonstration of the adjusted system in the Ashtekar equation. We plot the violation of
the constraint for the same model as figure 3(b). An artificial error term was added at t = 6, in theform of A2y → A2y (1 + error), where error is + 20% as before. (a), (b) L2 norm of the Hamiltonian
constraint equation, CH , and momentum constraint equation, CMx , respectively. The full curve isthe case of κ = 0, that is the case of ‘no adjusted’ original Ashtekar equation (weakly hyperbolic
system). The dotted curve is for κ = 1, equivalent to the symmetric hyperbolic system. We seethat the other curve (κ = 2.0) shows better performance than the symmetric hyperbolic case.

The Adjusted system (essentials):
Purpose:Control the violation of constraints by reformulating the system so as to have aconstrained surface an attractor.

Procedure:

Add a particular combination of constraints to the evolution equations, and adjust
its multipliers.

Theoretical support: Eigenvalue analysis of the constraint propagation equations.
Advantages:

Available even if the base system is not a symmetric hyperbolic.

Advantages:

Keep the number of the variable same with the original system.

Conjecture on Constraint Amplification Factors (CAFs):(A) If CAF has a negative real-part (the constraints are forced to be diminished), then we see morestable evolution than a system which has positive CAF.
(B) If CAF has a non-zero imaginary-part (the constraints are propagating away), then we see morestable evolution than a system which has zero CAF.

2001
so-called BSSN
62

Shibata

87,95,99

ADM

AEI

BSSN

PennState
Caltechhyperbolic formulation
01
92
92

Bona-Masso

04

Kidder-Scheel-Teukolsky

04

Z4 (Bona et.al.)Z4-lambda 05(Gundlach-Calabrese)harmonic
99

Nagy-Ortiz-Reulalambda system

Shinkai-Yoneda

asymptotically constrained /constraint damping
01
87

Detweiler02

adjusted-system
02

Illinois

adjusted BSSN (Yoneda-Shinkai)

adjusted ADM (Shinkai-Yoneda)

2001
so-called BSSN
62

Shibata

87,95,99

ADM

2005

AEI

BSSN

PennState
Caltechhyperbolic formulation
01
92
92

Bona-Masso

04

Kidder-Scheel-Teukolsky

04

LSU

Z4 (Bona et.al.)Z4-lambda 05(Gundlach-Calabrese)harmonic
99

Nagy-Ortiz-Reulalambda system

Shinkai-Yoneda

asymptotically constrained /constraint damping
01
87

Detweiler02

adjusted-system
02

Pretorius
Illinois

adjusted BSSN (Yoneda-Shinkai)

adjusted ADM (Shinkai-Yoneda)

General-covariant Z4 system
Bona, Palenzuela et al, PRD66(2002)084013, PRD69(2004)104003, arXiv:0710.4425, 0811.1691• Introduce 4-dim vector field Zµ as
1Rµν + ∇µZν + ∇ν Zµ = 8π(Tµν − T gµν ) .
2• Construct a general-covariant first-order form using variables (Kij , gij , Θ ≡ αZ 0 , Z i)(∂t − Lβ ) γij = −2αKij

1(∂t − Lβ ) Kij = −∇iαj + α[(3)Rij + ∇iZj + ∇j Zi − 2Kij2 + (K − 2Θ)Kij − Sij + (S − τ )γij ]2α (3)(∂t − Lβ ) Θ = [ R + 2∇k Z k + (K − 2Θ)K − K 2 − 2Z k αk /α − 2τ ]2(∂t − Lβ ) Zi = α [∇j (Kij − δi j K) + ∂iΘ − 2Kij Zj − Θαi /α − Si ],where τ ≡ 8πα2 T 00, Si ≡ 8πα T 0i , Sij ≡ 8π Tij , together with dynamical lapse and shift:(∂t − β i ∂i)α = . . . ,(∂t − β i ∂i )β i = . . . ,• Recover BSSN and KST. (!?)• Constraints Propagation? Adjustments in eqns?
• long-term single BH evolution without excision. ∼ Bona-Massó in 90s.

2008/11/28 @ 大阪市立大学

Formulation for Numerical Relativity
both for Einstein / Gauss-Bonnet
真貝寿明

大阪工業大学情報科学部

1. Introduction
定式化問題？
2. The Standard Approach to Numerical Relativity
ADM 形式，BSSN 形式，Hyperbolic 形式
3. Robust system for Constraint Violation
Adjusted systems .... better than lambda system!
Adjusted ADM......... why ADM blows up?
Adjusted BSSN
4. 高次元数値相対論に向けてstrategy 0 The standard approach :: Arnowitt-Deser-Misner (ADM) formulation (1962)
shift vector, βi
3+1 decomposition of the spacetime.
surface normal line
Evolve 12 variables (γij , Kij )coordinate constant line
βi dtwith a choice of gauge condition.
lapse function, αα dt

t = constant hypersurface

constraints

evolution eqs.

Maxwell eqs.
div E = 4πρ
div B = 0

ADM Einstein eq.
(3)
R + (trK)2 − Kij K ij = 2κρH + 2ΛDj K ji − Di trK = κJi

14π∂tE = rot B − jcc∂t γij = −2N Kij + Dj Ni + DiNj ,∂t Kij = N ( (3)Rij + trKKij ) − 2N Kil K l j − Di Dj N
+ (Dj N m )Kmi + (DiN m )Kmj + N m DmKij − N γij Λ− κα{Sij + 12 γij (ρH − trS)}

1∂tB = −rot E
cThe Standard ADM formulation (aka York 1978):The fundamental dynamical variables are (γij , Kij ), the three-metric and extrinsic
curvature. The three-hypersurface Σ is foliated with gauge functions, (α, β i), thelapse and shift vector.
• The evolution equations:∂tγij = −2αKij + Diβj + Dj βi,∂tKij = α (3)Rij + αKKij − 2αKik K k j − DiDj α+(Diβ k )Kkj + (Dj β k )Kki + β k Dk Kij−8πGα{Sij + (1/2)γij (ρH − trS)},

where K = K ii, and (3)Rij and Di denote three-dimensional Ricci curvature,
and a covariant derivative on the three-surface, respectively.
• Constraint equations:
Hamiltonian constr.
momentum constr.
where (3)R =(3) Rii.

HADM := (3)R + K 2 − Kij K ij ≈ 0,
MADM:= Dj K j i − DiK ≈ 0,i

S. Frittelli, Phys. Rev. D55, 5992 (1997)
HS and G. Yoneda, Class. Quant. Grav. 19, 1027 (2002)

The Constraint Propagations of the Standard ADM:
∂tH = β j (∂j H) + 2αKH − 2αγ ij (∂iMj )+α(∂l γmk )(2γ ml γ kj − γ mk γ lj )Mj − 4γ ij (∂j α)Mi,∂tMi = −(1/2)α(∂iH) − (∂iα)H + β j (∂j Mi)+αKMi − β k γ jl (∂iγlk )Mj + (∂iβk )γ kj Mj .

From these equations, we know that
if the constraints are satisfied on the initial slice Σ,then the constraints are satisfied throughout evolution (in principle).
But this is NOT TRUE in NUMERICS....

Original ADM The original construction by ADM uses the pair of (hij , π ij ).
√√1
L = −gR = hN [(3)R − K 2 + Kij K ij ], where Kij = £n hij
2√∂Lthen π ij == h(K ij − Khij ),∂ ḣijThe Hamiltonian density gives us constraints and evolution eqs.
"
√ !
ij
j
−1/2
ij
Nj π ) ,H = π ḣij − L = h N H(h, π) − 2Nj M (h, π) + 2Di(h















δH
N
1
√
=
2
−
(π
hij π) + 2D(iNj),
ij
δπ ij
2
h
√
δH
1 (3) ij
1 N ij
1 2
N in j 1 ij
(3) ij
mn
ij

√
√
π
=
−
=
−
hN
(
R
−
Rh
)
+
(π
π
−
)
−
2
h
(π πn − ππ )
∂
π

t
mn



δh
2
2
2
2

h√
h
ij

√





+ h(DiDj N − hij DmDmN ) + hDm(h−1/2N mπ ij ) − 2π m(iDmN j)∂thij =

Standard ADM (by York) NRists refer ADM as the one by York with a pair of (hij , Kij ).
∂thij = −2N Kij + Dj Ni + DiNj ,

(3)
 ∂ K
Rij + KKij ) − 2N Kil K l j − DiDj N + (Dj N m)Kmi + (DiN m)Kmj + N mDmKij
t ij = N (In the process of converting, H was used, i.e. the standard ADM has already adjusted.

Adjusted ADM systems
PRD 63 (2001) 120419, CQG 19 (2002) 1027

We adjust the standard ADM system using constraints as:∂t γij = −2αKij + ∇iβj + ∇j βi ,+Pij H + Qk ij Mk + pk ij (∇k H) + q kl ij (∇k Ml ),

(1)
(2)

(3)

∂t Kij = αRij + αKKij − 2αKik K k j − ∇i∇j α + (∇iβ k )Kkj + (∇j β k )Kki + β k ∇k Kij ,(3)+Rij H + S k ij Mk + rk ij (∇k H) + skl ij (∇k Ml ),(4)with constraint equations
H := R(3) + K 2 − Kij K ij ,Mi := ∇j K j i − ∇iK.
We can write the adjusted constraint propagation equations as(7)∂tH = (original terms) + H1mn [(2)] + H2imn ∂i [(2)] + H3ijmn ∂i ∂j [(2)] + H4mn [(4)],∂t Mi = (original terms) + M1i mn[(2)] + M2i jmn∂j [(2)] + M3i mn[(4)] + M4i jmn∂j [(4)].(8)

(5)
(6)

3

Constraint propagation of ADM systems

3.1

Original ADM vs Standard ADM

Try the adjustment Rij = κ1αγij and other multiplier zero, where κ1 = 0
the standard ADM−1/4 the original ADM• The constraint propagation eqs keep the first-order form (cf Frittelli, PRD55(97)5992):





l
jl 
H
β
−2αγ
H
 " 
∂ 
.
∂t 
(5)
l
l
l
l
l j
Mi
Mj
−(1/2)αδi + R i − δi R
β δi
The eigenvalues of the characteristic matrix:
l









l

l

l

)

λ = (β , β , β ± α2γ ll (1 + 4κ1))

symmetric hyperbolic when κ1 = 3/2
The hyperbolicity of (5): strongly hyperbolic when α2γ ll (1 + 4κ1) > 0 weakly hyperbolic
when α2γ ll (1 + 4κ1) ≥ 0• On the Minkowskii background metric, the linear order terms of the Fourier-transformed
constraint propagation equations gives the eigenvalues
l





That is, )Λ = (0, 0, ± −k 2(1 + 4κ1)).
(two 0s, two pure imaginary) for the standard ADM(four 0s)
for the original ADM

BETTER STABILITY

Comparisons of Adjusted ADM systems (Teukolsky wave)
3-dim, harmonic slice, periodic BC
HS original Cactus/GR code

L2 norm of Hamiltonian constraint

100
10-1
Original ADM

10-2Standard ADM

10

-3

10

-4

Adjusted ADM (SimpDet)

Adjusted ADM (Det)

10-5
10-6
0

50

100

150

200

250

300

time
Figure 1: Violation of Hamiltonian constraints versus time: Adjusted ADM systems applied for Teukolsky wave initial data evolution
with harmonic slicing, and with periodic boundary condition. Cactus/GR/evolveADMeq code was used. Grid = 243 , ∆x = 0.25, iterative
Crank-Nicholson method.

4Constraint propagations in spherically symmetric spacetime

4.1

The procedure

The discussion becomes clear if we expand the constraint Cµ := (H, Mi)T using vector harmonics.
Cµ =

! "

l,m

lm

lm

lm

lm

#

A (t, r)alm(θ, ϕ) + B blm + C clm + D dlm ,(1)

where we choose
the basis of thevector harmonicsas






0 
0
 Ylm 
 0 



















r
r
0
0
 0 
 Ylm 











.
*
*
alm =  , blm =  , clm = , dlm =


 0 
 0 
∂ Y

− 1 ∂ Y

ϕ
lm




 θ lm 


l(l
+
1)
l(l
+
1)
sin θ








0
0
∂ϕYlm
sin θ ∂θ Ylm
The basis are normalized so that they satisfy
+ 2π

+ π

#Cµ, Cν $ = 0 dϕ 0 Cµ∗Cρ η νρ sin θdθ,

where η νρ is Minkowskii metric and the asterisk denotes the complex conjugate. Therefore
Alm = #alm(ν) , Cν $,∂tAlm = #alm(ν) , ∂t Cν $,etc.

We also express these evolution equations using the Fourier expansion on the radial coordinate,
Alm =

!

k

ikr
Âlm
etc.
(k) (t) e

(2)

lm
T
So that we will be able to obtain the RHS of the evolution equations for (Âlm(k) (t), · · · , D̂(k) (t))
in a homogeneous form.

4.2

Constraint propagations in Schwarzschild spacetime

1. the standard Schwarzschild coordinate
2M 2
dr2
ds = −(1 −)dt ++ r2dΩ2,r
1 − 2M/r2(the standard expression)

2. the isotropic coordinate, which is given by, r = (1 + M/2riso)2riso:
ds2 = −(

1 − M/2riso 2 2
M 4 22) dt + (1 +) [driso + risodΩ2],
1 + M/2riso
2riso(the isotropic expression)

3. the ingoing Eddington-Finkelstein (iEF) coordinate, by tiEF = t + 2M log(r − 2M ) :ds2 = −(1 −

2M 2
4M
2M
)dtiEF +dtiEF dr + (1 +)dr2 + r2dΩ2rrr(the iEF expression)

4. the Painlevé-Gullstrand (PG) coordinates,




,
.

2M
2M 2
dtP G dr + dr2 + r2dΩ2 ,dtP G + 2
(the PG expression)
r
r
*
*
√
which is given by tP G = t + 8M r − 2M log{( r/2M + 1)/( r/2M − 1)}
ds2 = − 1 −

Example 1: standard ADM vs original ADM (in Schwarzschild coordinate)
(b)

(a)

original ADM (! = - 1/4)
F

Real / Imaginary parts of Eigenvalues (AF)Real / Imaginary parts of Eigenvalues (AF)no adjustments (standard ADM)
1

0.5

0

-0.5

-1
0

5

10

r

15

sch

20

0.5

0

-0.5
0

5

10

r

15

20

sch

Figure 1: Amplification factors (AFs, eigenvalues of homogenized constraint propagation equations) are shown for the standard
Schwarzschild coordinate, with (a) no adjustments, i.e., standard ADM, (b) original ADM (κF = −1/4). The solid lines andthe dotted lines with circles are real parts and imaginary parts, respectively. They are four lines each, but actually the two
eigenvalues are zero for all cases. Plotting range is 2 < r ≤ 20 using Schwarzschild radial coordinate. We set k = 1, l = 2, andm = 2 throughout the article.

∂tγij = −2αKij + ∇iβj + ∇j βi,(3)∂tKij = αRij + αKKij − 2αKik K k j − ∇i∇j α + (∇iβ k )Kkj + (∇j β k )Kki + β k ∇k Kij + κF αγij H,

Example 2: Detweiler-type adjusted (in Schwarzschild coord.)(b)

Detweiler type, ! = + 1/2(c)

Detweiler type, ! = - 1/2

1

L

Real / Imaginary parts of Eigenvalues (AF)Real / Imaginary parts of Eigenvalues (AF)

L

0.5

0

-0.5

-1
0

5

10

r

sch

15

20

1

0.5

0

-0.5

-1
0

5

10

r

15

20

sch

Figure 2: Amplification factors of the standard Schwarzschild coordinate, with Detweiler type adjustments. Multipliers used in
the plot are (b) κL = +1/2, and (c) κL = −1/2.

∂tγij = (original terms) + Pij H,∂tKij = (original terms) + Rij H + S k ij Mk + skl ij (∇k Ml ),where Pij = −κLα3γij , Rij = κLα3(Kij − (1/3)Kγij ),kl− (∂l α)γij γ kl ], skl ij = κLα3[δ(ik δj)− (1/3)γij γ kl ],S k ij = κLα2[3(∂(iα)δj)Detweiler’s criteria vs Our criteria• Detweiler calculated the L2 norm of the constraints, Cα , over the 3-hypersurface and imposed
its negative definiteness of its evolution,
Detweiler’s criteria ⇔ ∂t

! "αCα2 dV < 0,This is rewritten by supposing the constraint propagation to be ∂tĈα = Aα β Ĉβ in the Fouriercomponents,
! "
! "
¯
¯
¯
⇔ ∂t
Ĉα Ĉ α dV =Aα β Ĉβ Ĉ α + Ĉα Āα β Ĉ β dV < 0, ∀ non zero Ĉααα†⇔ eigenvalues of (A + A ) are all negative for ∀k.

• Our criteria is that the eigenvalues of A are all negative. Therefore,
Our criteria # Detweiler’s criteria• We remark that Detweiler’s truncations on higher order terms in C-norm corresponds our
perturbative analysis, both based on the idea that the deviations from constraint surface (theerrors expressed non-zero constraint value) are initially small.

Constraint propagation of ADM systems(2) Detweiler’s system
Detweiler’s modification to ADM [PRD35(87)1095] can be realized in our notation as:
Pij = −Lα3γij ,Rij = Lα3(Kij − (1/3)Kγij ),k− (∂l α)γij γ kl ],Sijk = Lα2[3(∂(iα)δj)
3
k l
kl
skl
ij = Lα [2δ(i δj) − (1/3)γij γ ],and else zero, where L is a constant.

• This adjustment does not make constraint propagation equation in the first order form, sothat we can not discuss the hyperbolicity nor the characteristic speed of the constraints.
• For the Minkowskii background spacetime, the adjusted constraint propagation equations
with above choice of multiplier become∂tH = −2(∂j Mj ) + 4L(∂j ∂j H),∂tMi = −(1/2)(∂iH) + (L/2)(∂k ∂k Mi) + (L/6)(∂i∂k Mk ).
Constraint Amp. Factors (the eigenvalues of their Fourier expression) are
l

2

2

!

Λ = (−(L/2)k (multiplicity 2), −(7L/3)k ± (1/3) k 2(−9 + 25L2k 2).)This indicates negative real eigenvalues if we chose small positive L.

Example 3: standard ADM (in isotropic/iEF coord.)(b)
iEF coordinate, no adjustments (standard ADM)

Real / Imaginary parts of Eigenvalues (AF)isotropic coordinate, no adjustments (standard ADM)
1

0.5

0

-0.5

-1
0

5

10

r

iso
￿

15

20

Real / Imaginary parts of Eigenvalues (AF)

(a)
2

1.5

1

0.5

0

-0.5

-1
0

5

10

r

15

20

sch

Figure 3: Comparison of amplification factors between different coordinate expressions for the standard ADM formulation (i.e.
no adjustments). Fig. (a) is for the isotropic coordinate (1), and the plotting range is 1/2 ≤ riso . Fig. (b) is for the iEF
coordinate (1) and we plot lines on the t = 0 slice for each expression. The solid four lines and the dotted four lines with circles
are real parts and imaginary parts, respectively.

Example 4: Detweiler-type adjusted (in iEF/PG coord.)
(b)

(c)

iEF coordinate, Detweiler type ! =+0.5

PG coordinate, Detweiler type ! =+0.5
L2

Real / Imaginary parts of Eigenvalues (AF)Real / Imaginary parts of Eigenvalues (AF)

L

1.5

1

0.5

0

-0.5

-1
0

5

10

r

sch

15

20

2

1.5

1

0.5

0

-0.5

-1
0

5

10

r

sch

Figure 4: Similar comparison for Detweiler adjustments. κL = +1/2 for all plots.

15

20

“Einstein equations” are time-reversal invariant. So ...
Why all negative amplification factors (AFs) are available?
Explanation by the time-reversal invariance (TRI)• the adjustment of the system I,
α γij !"#$
H
adjust term to !"#$∂t Kij = κ1 !"#$! "# $(−) (−)!"#$(+) (+) (+)

preserves TRI. ... so the AFs remain zero (unchange).
• the adjustment by (a part of) Detweiler
α γij !"#$
H
adjust term to !"#$
∂t γij = −L !"#$!"#$(−) (+)

violates TRI. ... so the AFs can become negative.

!"#$(+) (+) (+)

Therefore
We can break the time-reversal invariant feature of the “ADM equations”.

Adjusted ADM systems
PRD 63 (2001) 120419, CQG 19 (2002) 1027

We adjust the standard ADM system using constraints as:∂t γij = −2αKij + ∇iβj + ∇j βi ,+Pij H + Qk ij Mk + pk ij (∇k H) + q kl ij (∇k Ml ),

(1)
(2)

(3)

∂t Kij = αRij + αKKij − 2αKik K k j − ∇i∇j α + (∇iβ k )Kkj + (∇j β k )Kki + β k ∇k Kij ,(3)+Rij H + S k ij Mk + rk ij (∇k H) + skl ij (∇k Ml ),(4)with constraint equations
H := R(3) + K 2 − Kij K ij ,Mi := ∇j K j i − ∇iK.
We can write the adjusted constraint propagation equations as(7)∂tH = (original terms) + H1mn [(2)] + H2imn ∂i [(2)] + H3ijmn ∂i ∂j [(2)] + H4mn [(4)],∂t Mi = (original terms) + M1i mn[(2)] + M2i jmn∂j [(2)] + M3i mn[(4)] + M4i jmn∂j [(4)].(8)

(5)
(6)

No intable 1

0
P-1
P-2
P-3
P-4
P-5
Q-1
Q-2
Q-3
Q-4
R-1
R-2
R-3
S-1
S-2
p-1
p-2
p-3
q-1
q-2
r-1
r-2
r-3
s-1
s-2
s-3

0
2-P
3
–
–
–
–
–
–
–
1
4
–
2-S
–
–
–
–
–
–
–
–
–
2-s
–
–

Schwarzschild/isotropic coordinates
Adjustment
–
Pij
Pij
Pij
Pij
Pij
Qk ij
Qk ij
Qk ij
Qk ij
Rij
Rij
Rij
S k ijS k ijp k ijp k ijp k ijq kl ij
q kl ijr k ijr k ijr k ijs kl ijs kl ijs kl ij

no adjustments−κL α 3 γij−κL αγij
Prr = −κ or Prr = −κα
−κγij
−κγrr
καβ k γij
Qr rr = κQr ij = κγij or Qr ij = καγij
Qr rr = κγrr
κF αγij
Rrr = −κµ α or Rrr = −κµRrr = −κγrrk
κL α 2 [3(∂(i α)δj)− (∂l α)γij γ kl ]lkκαγ (∂l γij )p r ij = −καγij
p r rr = καp r rr = καγrr
q rr ij = καγij
q rr rr = −καγrr
r r ij = καγij
r r rr = −καr r rr = −καγrrl
κL α 3 [δ(ik δj)− (1/3)γij γ kl ]rrs ij = −καγij
s rr rr = −καγrr

1st?

TRS

Real

Imaginary

yes
no
no
no
no
no
no
no
no
no
yes
yes
yes
yes
yes
no
no
no
no
no
no
no
no
no
no
no

–
no
no
no
no
no
no
yes
yes
yes
yes
no
no
no
no
no
no
no
no
yes
yes
yes
yes
no
no
no

–
–
makes 2 Neg.
not apparent
makes 2 Neg.
not apparent
slightly enl.Neg.
not apparent
makes 2 Neg.
not apparent
red. Pos./enl.Neg.
not apparent
N/A
N/A
red. abs vals.
not apparent
red. abs vals.
not apparent
red. abs vals.
not apparent
κF = −1/4 min. abs vals.
not apparent
not apparent
enl. vals.
not apparent
not apparent
not apparent
makes 2 Neg.
not apparent
red. Pos.
red. vals.
red. Pos.
red. vals.
makes 2 Neg.
enl. vals.
κ = 1/2 min. vals.
red. vals.
red. abs vals.
not apparent
not apparent
not apparent
red. abs vals.
enl. vals.
red. abs vals.
enl. vals.
makes 4 Neg.
not apparent
makes 2 Neg.
red. vals.
makes 2 Neg.
red. vals.

iEF/PG coordinates
RealImaginary––
makes 2 Neg.
not apparent
makes 2 Neg.
not apparent
slightly enl.Neg.
not apparent
makes 2 Neg.
not apparent
red.Pos./enl.Neg.
not apparentκ ∼ 1.35 min. vals.
not apparent
red. abs vals.
not apparent
enl.Neg.
enl. vals.
red. abs vals.
not apparent
κF = −1/4 min. vals.
red.Pos./enl.Neg.
enl. vals.
red.Pos./enl.Neg.
enl. vals.
not apparent
not apparent
makes 2 Neg.
not apparent
red. Pos.
enl. vals.
red.Pos/enl.Neg.
enl. vals.
red. Pos. vals.
red. vals.
not apparent
enl. vals.
not apparent
not apparent
not apparent
enl. vals.
red. abs vals.
enl. vals.
red. abs vals.
enl. vals.
makes 4 Neg.
not apparent
makes 2 Neg.
red. vals.
makes 2 Neg.
red. vals.

H Shinkai and G Yoneda

No

1044

Table 3. List of adjustments we tested in the Schwarzschild spacetime. The column of adjustments are nonzero multipliers in terms of (13) and (14). The column ‘1st?’ and ‘TRS’ arethe same as in table 1. The effects to amplification factors (when κ > 0) are commented for each coordinate system and for real/imaginary parts of AFs, respectively. The ‘N/A’ meansthat there is no effect due to the coordinate properties; ‘not apparent’ means the adjustment does not change the AFs effectively according to our conjecture; ‘enl./red./min.’ means
enlarge/reduce/minimize, and ‘Pos./Neg.’ means positive/negative, respectively. These judgements are made at the r ∼ O(10M) region on their t = 0 slice.

Numerical Tests (method)• Cactus-based original “GR” code
http://www.cactuscode.org/[CactusBase+CactusPUGH+GR]• 3+1dim，linear wave evolution(Teukolsky wave)• harmonic slice• periodic boundary, [-3,+3]• iterative Crank-Nicholson method• 12^3, 24^3, 48^3, 96^3
Towards standard testbeds for numerical relativity
Mexico Numerical Relativity Workshop 2002 Participants
CQG 21 (2004) 589-613

Numerical Tests (Detweiler-type)
PRD35(1987)1095

24^3

ADM (standard)
Detweiler k=0.025
Detweiler k=0.01
Detweiler k=0.005

10-1

error (norm of Hamiltonian constraint)

error (norm of Hamiltonian constraint)

10-1

10-2

10-3

10-4

10-5

ADM (standard)
Detweiler 48^3 k=0.02
Detweiler 48^3 k=0.01
Detweiler 48^3 k=0.005
Detweiler 48^3 k=0.0025

48^3

10-2

10-3

10-4

10-5

10-6

10-6
0

50

100

150

200

time

250

300

350

0

50

100

150

time

200

250

Numerical Tests (Simpliﬁed Detweiler)

ADM (standard) 24^3
simp. Det. 24^3 k=0.1
simp. Det. 24^3 k=0.05
simp. Det. 24^3 k=0.01

ADM (standard) 48^3
simp. Det. 48^3 k=0.01
simp. Det. 48^3 k=0.0025

10-1

error (norm of Hamiltonian constraint)

error (norm of Hamiltonian constraint)

100

10-1

10-2

10-3

10-4

10-5

10-2

10-3

10-4

10-5
0

50

100

150

time

200

250

300

0

50

100

150

time

200

250

300

Numerical Tests (Detweiler, k-adjust)

24^3

10-1

error (norm of Hamiltonian constraint)

ADM (standard)
Detweiler k=0.05
Detweiler k=0.05 cut 10^{-3}
Detweiler k=0.05 cut 10^{-4}
Detweiler k=0.05 cut 10^{-5}

10-2

10-2

10-3

10-3

10-4

10-4

10-5

10-5

10-6

10-6
0

50

100

150

200

time

250

300

350

zero
all
cut3
cut4
cut5

10-1

0
1
kappa*factor

Numerical Tests (Detweiler, k-adjust)ADM (standard)
Detweiler k=0.05, slope1-4
Detweiler k=0.05, slope1-5
Detweiler k=0.05, slope1-6
Detweiler k=0.05

12^3

error (norm of Hamiltonian constraint)

10-1

10-2

10-2

10-3

10-3

10-4

10-4

10-5

zero
all
slope14slope15
slope16

10-1

10-5

10-6

10-6
0
1
kappa*factor

10-7
0

100

200

300

400

time

500

600

700

Numerical Tests (Detweiler, k-adjust)ADM (standard)

12^3

Detweiler k=0.05, off @ 10^(-3)!"#$%&'(
Detweiler k=0.05, off @ 10^(-4)

error (norm of Hamiltonian constraint)

10-1

Detweiler k=0.05, off @ 10^(-5)
Detweiler k=0.05

10-2

10-2

10-3

10-3

10-4

10-4

10-5

10-5

10-6

10-6

10-7
0

100

200

300

400

time

500

600

700

800

zero
all
offH3
offH4
offH5

10-1

0
1
kappa*factor

2008/11/28 @ 大阪市立大学

Formulation for Numerical Relativity
both for Einstein / Gauss-Bonnet
真貝寿明

大阪工業大学情報科学部

1. Introduction
定式化問題？
2. The Standard Approach to Numerical Relativity
ADM 形式，BSSN 形式，Hyperbolic 形式
3. Robust system for Constraint Violation
Adjusted systems .... better than lambda system!
Adjusted ADM......... why ADM blows up?
Adjusted BSSN......... why BSSN works well?
4. 高次元数値相対論に向けてstrategy 1

Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation
T. Nakamura, K. Oohara and Y. Kojima, Prog. Theor. Phys. Suppl. 90, 1 (1987)
M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 (1995)
T.W. Baumgarte and S.L. Shapiro, Phys. Rev. D 59, 024007 (1999)

The popular approach. Nakamura’s idea in 1980s.
BSSN is a tricky nickname. BS (1999) introduced a paper of SN (1995).
• define new set of variables (φ, γ̃ij ,K,Ãij ,Γ̃i), instead of the ADM’s (γij ,Kij ) whereγ̃ij ≡ e−4φ γij ,Ãij ≡ e−4φ(Kij − (1/3)γij K),Γ̃i ≡ Γ̃ijk γ̃ jk ,and impose detγ̃ij = 1 during the evolutions.
• The set of evolution equations become(∂t − Lβ )φ = −(1/6)αK,(∂t − Lβ )γ̃ij = −2αÃij ,(∂t − Lβ )K = αÃij Ãij + (1/3)αK 2 − γ ij (∇i∇j α),(3)(∂t − Lβ )Ãij = −e−4φ(∇i∇j α)T F + e−4φαRij − e−4φ α(1/3)γij R(3) + α(K Ãij − 2Ãik Ãk j )∂tΓ̃i = −2(∂j α)Ãij − (4/3)α(∂j K)γ̃ ij + 12αÃji(∂j φ) − 2αÃk j (∂j γ̃ ik ) − 2αΓ̃k lj Ãj k γ̃ il


k
ij
kj
i
ki
j
ij
k
−∂j β ∂k γ̃ − γ̃ (∂k β ) − γ̃ (∂k β ) + (2/3)γ̃ (∂k β )Momentum constraint was used in Γi-eq.

• Calculate Riemann tensor as
φkm k
Rij = ∂k Γkij − ∂i Γkkj + Γm
Γ
−
Γ
Γ
=:
R̃
+
R
ij
ij
ij mk
kj mi

φ
Rij
= −2D̃iD̃j φ − 2g̃ij D̃ l D̃l φ + 4(D̃iφ)(D̃j φ) − 4g̃ij (D̃ l φ)(D̃l φ)R̃ij = −(1/2)g̃ lm∂lm g̃ij + g̃k(i ∂j)Γ̃k + Γ̃k Γ̃(ij)k + 2g̃ lm Γ̃kl(iΓ̃j)km + g̃ lm Γ̃kim Γ̃klj• Constraints are H, Mi.
But thre are additional ones, G i , A, S.
Hamiltonian and the momentum constraint equations
HBSSN = RBSSN + K 2 − Kij K ij ,MBSSN
= MADM
,
i
i

(1)
(2)

Additionally, we regard the following three as the constraints:
G i = Γ̃i − γ̃ jk Γ̃ijk ,A = Ãij γ̃ ij ,S = γ̃ − 1,

Why BSSN better than ADM?
Is the BSSN best? Are there any alternatives?

(3)
(4)
(5)

Constraints in BSSN system
The normal Hamiltonian and momentum constraints
HBSSN = RBSSN + K 2 − Kij K ij ,MBSSN
= MADM
,
i
i

(1)
(2)

Additionally, we regard the following three as the constraints:
G i = Γ̃i − γ̃ jk Γ̃ijk ,A = Ãij γ̃ ij ,S = γ̃ − 1,

(3)
(4)
(5)

Adjustments in evolution equations∂tB ϕ = ∂tA ϕ + (1/6)αA − (1/12)γ̃ −1(∂j S)β j ,(6)∂tB γ̃ij = ∂tA γ̃ij − (2/3)αγ̃ij A + (1/3)γ̃ −1(∂k S)β k γ̃ij ,(7)∂tB K = ∂tA K − (2/3)αKA − αHBSSN + αe−4ϕ (D̃j G j ),(8)
∂tB Ãij = ∂tA Ãij + ((1/3)αγ̃ij K − (2/3)αÃij )A + αe−4ϕ ((1/2)(∂k γ̃ij ) − (1/6)γ̃ij γ̃ −1(∂k S))G k+αe−4ϕ γ̃k(i (∂j)G k ) − (1/3)αe−4ϕγ̃ij (∂k G k )(9)∂tB Γ̃i = ∂tA Γ̃i − ((2/3)(∂j α)γ̃ ji + (2/3)α(∂j γ̃ ji) + (1/3)αγ̃ jiγ̃ −1(∂j S) − 4αγ̃ ij (∂j ϕ))A−(2/3)αγ̃ ji(∂j A) + 2αγ̃ ij Mj − (1/2)(∂k β i)γ̃ kj γ̃ −1(∂j S) + (1/6)(∂j β k )γ̃ ij γ̃ −1(∂k S)+(1/3)(∂k β k )γ̃ ij γ̃ −1(∂j S) + (5/6)β k γ̃ −2γ̃ ij (∂k S)(∂j S) + (1/2)β k γ̃ −1(∂k γ̃ ij )(∂j S)(10)+(1/3)β k γ̃ −1(∂j γ̃ ji)(∂k S).

A

Full set of BSSN constraint propagation eqs.










HBS  A11
A12 A13
A14
A15
HBS 







 M 
 −(1/3)(∂ α) + (1/6)∂
 M 
αK A23
0
A25
i 
i
i
j 











ij
j 
 Gi  = 



0
αγ̃
0
A
A
G
∂tBS 
34
35










k





0
0
0
β
(∂
S)
−2αγ̃
 S 

 S 
k





k
A
0
0
0
0
αK + β ∂kA
A11 = +(2/3)αK + (2/3)αA + β k ∂kA12 = −4e−4ϕ α(∂k ϕ)γ̃ kj − 2e−4ϕ (∂k α)γ̃ jk
A13 = −2αe−4ϕ Ãk j ∂k − αe−4ϕ (∂j Ãkl )γ̃ kl − e−4ϕ (∂j α)A − e−4ϕ β k ∂k ∂j − (1/2)e−4ϕ β k γ̃ −1 (∂j S)∂k+(1/6)e−4ϕ γ̃ −1 (∂j β k )(∂k S) − (2/3)e−4ϕ (∂k β k )∂j
A14 = 2αe−4ϕ γ̃ −1 γ̃ lk (∂l ϕ)A∂k + (1/2)αe−4ϕ γ̃ −1 (∂l A)γ̃ lk ∂k + (1/2)e−4ϕ γ̃ −1 (∂l α)γ̃ lk A∂k + (1/2)e−4ϕ γ̃ −1 β m γ̃ lk ∂m ∂l ∂k−(5/4)e−4ϕ γ̃ −2 β m γ̃ lk (∂m S)∂l ∂k + e−4ϕ γ̃ −1 β m (∂m γ̃ lk )∂l ∂k + (1/2)e−4ϕ γ̃ −1 β i (∂j ∂i γ̃ jk )∂k+(3/4)e−4ϕ γ̃ −3 β i γ̃ jk (∂i S)(∂j S)∂k − (3/4)e−4ϕ γ̃ −2 β i (∂i γ̃ jk )(∂j S)∂k + (1/3)e−4ϕ γ̃ −1 γ̃ pj (∂j β k )∂p ∂k−(5/12)e−4ϕ γ̃ −2 γ̃ jk (∂k β i )(∂i S)∂j + (1/3)e−4ϕ γ̃ −1 (∂k γ̃ ij )(∂j β k )∂i − (1/6)e−4ϕ γ̃ −1 γ̃ mk (∂k ∂l β l )∂m
A15 = (4/9)αKA − (8/9)αK 2 + (4/3)αe−4ϕ (∂i ∂j ϕ)γ̃ ij + (8/3)αe−4ϕ (∂k ϕ)(∂l γ̃ lk ) + αe−4ϕ (∂j γ̃ jk )∂k+8αe−4ϕ γ̃ jk (∂j ϕ)∂k + αe−4ϕ γ̃ jk ∂j ∂k + 8e−4ϕ (∂l α)(∂k ϕ)γ̃ lk + e−4ϕ (∂l α)(∂k γ̃ lk ) + 2e−4ϕ (∂l α)γ̃ lk ∂k+e−4ϕ γ̃ lk (∂l ∂k α)kl
A23 = αe−4ϕ γ̃ km (∂k ϕ)(∂j γ̃mi ) − (1/2)αe−4ϕ Γ̃mkl γ̃ (∂j γ̃mi )+(1/2)αe−4ϕ γ̃ mk (∂k ∂j γ̃mi ) + (1/2)αe−4ϕ γ̃ −2 (∂i S)(∂j S) − (1/4)αe−4ϕ (∂i γ̃kl )(∂j γ̃ kl ) + αe−4ϕ γ̃ km (∂k ϕ)γ̃ji ∂m
kl−4ϕ mk+αe−4ϕ (∂j ϕ)∂i − (1/2)αe−4ϕ Γ̃mγ̃ Γ̃ijk ∂m + (1/2)αe−4ϕ γ̃ lk γ̃ji ∂k ∂lkl γ̃ γ̃ji ∂m + αe+(1/2)e−4ϕ γ̃ mk (∂j γ̃im )(∂k α) + (1/2)e−4ϕ (∂j α)∂i + (1/2)e−4ϕ γ̃ mk γ̃ji (∂k α)∂m
A25 = −Ãk i (∂k α) + (1/9)(∂i α)K + (4/9)α(∂i K) + (1/9)αK∂i − αÃk i ∂k
A34 = −(1/2)β k γ̃ il γ̃ −2 (∂l S)∂k − (1/2)(∂l β i )γ̃ lk γ̃ −1 ∂k + (1/3)(∂l β l )γ̃ ik γ̃ −1 ∂k − (1/2)β l γ̃ in (∂l γ̃mn )γ̃ mk γ̃ −1 ∂k+(1/2)β k γ̃ il γ̃ −1 ∂l ∂k
A35 = −(∂k α)γ̃ ik + 4αγ̃ ik (∂k ϕ) − αγ̃ ik ∂k

BSSN Constraint propagation analysis in flat spacetime• The set of the constraint propagation equations, ∂t (HBSSN , Mi, G i, A, S)T ?
• For the flat background metric gµν = ηµν , the first order perturbation equations of (6)-(10):∂t(1)ϕ = −(1/6)(1)K + (1/6)κϕ(1)A∂t(1)γ̃ij = −2(1)Ãij − (2/3)κγ̃ δij(1)A∂t(1)K = −(∂j ∂j(1)α) + κK1∂j(1)G j − κK2(1)HBSSN
BSSN T F∂t(1)Ãij = (1)(Rij) − (1)(D̃iD̃j α)T F + κA1δk(i(∂j)(1)G k ) − (1/3)κA2δij (∂k(1)G k )∂t(1)Γ̃i = −(4/3)(∂i(1)K) − (2/3)κΓ̃1(∂i(1)A) + 2κΓ̃2(1)Mi

(11)
(12)
(13)
(14)
(15)

We express the adjustements asκadj := (κϕ, κγ̃ , κK1, κK2, κA1, κA2, κΓ̃1, κΓ̃2).

(16)• Constraint propagation equations at the first order in the flat spacetime:∂t(1)HBSSN = (κγ̃ − (2/3)κΓ̃1 − (4/3)κϕ + 2) ∂j ∂j(1)A + 2(κΓ̃2 − 1)(∂j(1)Mj ),∂t(1)Mi = (−(2/3)κK1 + (1/2)κA1 − (1/3)κA2 + (1/2)) ∂i ∂j(1)G j+(1/2)κA1∂j ∂j(1)G i + ((2/3)κK2 − (1/2)) ∂i(1)HBSSN ,∂t(1)G i = 2κΓ̃2(1)Mi + (−(2/3)κΓ̃1 − (1/3)κγ̃ )(∂i(1)A),∂t(1)S = −2κγ̃(1)A,∂t(1)A = (κA1 − κA2)(∂j(1)G j ).

(17)
(18)
(19)
(20)
(21)

Effect of adjustments
No.

Constraints (number of components)
H (1) Mi (3) G i (3)
A (1)

S (1)

Amplification Factors (AFs)in Minkowskii background

0.
1.
2.

standard ADMuse
use(0, 0, ℑ, ℑ)
BSSN no adjustment
use
use
use
use
use
(0, 0, 0, 0, 0, 0, 0, ℑ, ℑ)
the BSSN
use+adj use+adj use+adj use+adj use+adj (0, 0, 0, ℑ, ℑ, ℑ, ℑ, ℑ, ℑ)

3.
4.
5.
6.
7.

no S adjustment
no A adjustment
no G i adjustment
no Mi adjustment
no H adjustment

use+adj
use+adj
use+adj
use+adjuse

use+adj
use+adj
use+adjuse
use+adj

8.
9.
10.
11.
12.

ignore G i , A, S
ignore G i , A
ignore G iignore A
ignore S

use+adj
use+adj
use+adj
use+adj
use+adj

use+adj(0, 0, 0, 0)
use+adj use+adj(0, ℑ, ℑ, ℑ, ℑ, ℑ, ℑ)
use+adj
use+adj use+adj (0, 0, 0, 0, 0, 0)
use+adj use+adj
use+adj (0, 0, ℑ, ℑ, ℑ, ℑ, ℑ, ℑ)
use+adj use+adj use+adj(0, 0, ℑ, ℑ, ℑ, ℑ, ℑ, ℑ)

use+adj
use+adjuse
use+adj
use+adj

use+adjuse
use+adj
use+adj
use+adjuse
use+adj
use+adj
use+adj
use+adj

no difference in flat background(0, 0, 0, ℑ, ℑ, ℑ, ℑ, ℑ, ℑ)(0, 0, 0, 0, 0, 0, 0, ℑ, ℑ)(0, 0, 0, 0, 0, 0, 0, ℜ, ℜ) Growing modes!
(0, 0, 0, ℑ, ℑ, ℑ, ℑ, ℑ, ℑ)

New Proposals :: Improved (adjusted) BSSN systems
TRS breaking adjustments
In order to break time reversal symmetry (TRS) of the evolution eqs, to adjust ∂t φ,∂t γ̃ij ,∂t Γ̃i using S, G i , or to adjust∂t K, ∂tÃij using Ã.
∂t φ = ∂tBS φ + κφH αHBS + κφG αD̃k G k + κφS1 αS + κφS2 αD̃j D̃j S∂t γ̃ij = ∂tBS γ̃ij + κγ̃H αγ̃ij HBS + κγ̃G1 αγ̃ij D̃k G k + κγ̃G2 αγ̃k(i D̃j) G k + κγ̃S1 αγ̃ij S + κγ̃S2 αD̃i D̃j S∂t K = ∂tBS K + κKMαγ̃ jk (D̃j Mk ) + κK Ã1 αÃ + κK Ã2αD̃j D̃j Ã
∂t Ãij = ∂tBS Ãij + κAM1 αγ̃ij (D̃k Mk ) + κAM2 α(D̃(i Mj) ) + κAÃ1 αγ̃ij Ã + κAÃ2 αD̃i D̃j Ã

∂t Γ̃i = ∂tBS Γ̃i + κΓ̃H αD̃i HBS + κΓ̃G1 αG i + κΓ̃G2 αD̃j D̃j G i + κΓ̃G3 αD̃i D̃j G j + κΓ̃S αD̃i HBSor in the flat background∂tADJ(1)φ =∂tADJ(1)γ̃ij =∂tADJ(1)K =∂tADJ(1)Ãij =∂tADJ(1)Γ̃i =+κφH(1)HBS + κφG ∂k(1)G k + κφS1(1)S + κφS2 ∂j ∂j(1)S+κγ̃H δij(1)HBS + κγ̃G1 δij ∂k(1)G k + (1/2)κγ̃G2 (∂j(1)G i + ∂i(1)G j ) + κγ̃S1 δij(1)S + κγ̃S2 ∂i∂j(1)S+κKM∂j(1)Mj + κK Ã1(1)Ã + κK Ã2 ∂j ∂j(1)Ã+κAM1 δij ∂k(1)Mk + (1/2)κAM2(∂iMj + ∂j Mi) + κAÃ1δij Ã + κAÃ2∂i ∂j Ã
+κΓ̃H ∂i(1)HBS + κΓ̃G1(1)G i + κΓ̃G2 ∂j ∂j(1)G i + κΓ̃G3 ∂i∂j(1)G j + κΓ̃S ∂i(1)SConstraint Amplification Factors with each adjustment

∂t φ
∂t φ
∂t γ̃ij∂t γ̃ij∂t γ̃ij∂t γ̃ij∂t γ̃ij∂t K∂t Ãij∂t Ãij∂t Ãij∂t Ãij∂t Γ̃i∂t Γ̃i∂t Γ̃i∂t Γ̃iadjustment
κφH αHκφG αD̃k G k
κSD αγ̃ij H
κγ̃G1 αγ̃ij D̃k G k

CAFs
√
2
(0, 0, ±√−k (∗3), 8κφHk 2 )(0, 0, ±√−k 2 (∗2), long expressions)(0, 0, ±√−k 2 (∗3), (3/2)κSD k 2 )(0, 0, ± −k 2 (∗2), longexpressions)q2(0,0, (1/4)k κγ̃G2 ± k 2 (−1 + k 2 κγ̃G2 /16)(∗2),
κγ̃G2 αγ̃k(i D̃j) G klong expressions)√
κγ̃S1 αγ̃ij S(0, 0, ±√−k 2 (∗3), 3κγ̃S1 )2
2
κγ̃S2 αD̃i D̃j S(0, 0, ± −k
√ (∗3), −κγ̃S2 k )(0, 0, 0, ± −k 2 (∗2), q
κKM αγ̃ jk (D̃j Mk )(1/3)κKM k 2 ± (1/3) k 2 (−9 + k 2 κ2KM ))√2
κAM1 αγ̃ij (D̃ k Mk ) (0, 0, ± −k 2 (∗3), −κq
AM1 k )
(0,0, −k 2 κAM2 /4 ± k 2 (−1 + k 2 κAM2 /16)(∗2) ,κAM2 α(D̃(i Mj))
long expressions)√(0, 0, ±√−k 2 (∗3), 3κAA1)
κAA1 αγ̃ij A
κAA2 αD̃i D̃j A(0, 0, ±√−k 2 (∗3), −κAA2 k 2 )(0, 0, ± −k 2 (∗3),q−κAA2 k 2 )κΓ̃H αD̃ i H
κΓ̃G1 αG i(0, 0, (1/2)κΓ̃G1 ± −k 2 + κ2Γ̃G1 (∗2) , long.)κΓ̃G2 αD̃ j D̃j G iiκΓ̃G3 αD̃ D̃j G

j(0, 0, −(1/2)κΓ̃G2 ±(0, 0, −(1/2)κΓ̃G3 ±q−k 2 + κ2Γ̃G2 (∗2) , long.)q−k 2 + κ2Γ̃G3 (∗2) , long.)

diag?
no
yes
yes
yes

effect of the adjustment
κφH < 0 makes 1 Neg.
κφG < 0 makes 2 Neg. 1 Pos.
κSD < 0 makes 1 Neg.
Case (B)κγ̃G1 > 0 makes 1 Neg.

yes

κγ̃G2 < 0 makes 6 Neg. 1 Pos.

nono

κγ̃S1 < 0 makes 1 Neg.
κγ̃S2 > 0 makes 1 Neg.

no

κKM < 0 makes 2 Neg.

yes

κAM1 > 0 makes 1 Neg.

yes

κAM2 > 0 makes 7 Neg

Case (D)

yes
yes
no
yes

κAA1 < 0 makes 1 Neg.
κAA2 > 0 makes 1 Neg.
κΓ̃H > 0 makes 1 Neg.
κΓ̃G1 < 0 makes 6 Neg. 1 Pos.

Case (E2)yesκΓ̃G2 > 0 makes 2 Neg. 1 Pos.

yesκΓ̃G3 > 0 makes 2 Neg. 1 Pos.

Case (E1)

Yoneda-HS, PRD66 (2002) 124003

An Evolution of Adjusted BSSN Formulation
by Yo-Baumgarte-Shapiro, PRD 66 (2002) 084026
Mass/M

Constraint residual

0.94

10

4

10

6

10

0.7
Innersurface + volume
Outersurface

0.93

0.6
0.5

Ham.
Mom.

0.92

0.4
0.91

8

0.3

0.9

!Krms

10
10

0.2

10

12

10

14

10

16

10

18

0.1

0.89

A3
A4
A5
A6
A7
A8

Angular momentum/M2

rms of !f0.96
Innersurface + volume
Outersurface

10

4

10

6

10

8

10

10

10

12

10

14

10

16

10
6000

18

"
trK

0.94
0.92
0.9

0

2000

4000

6000

t/M

0.88

0

2000

4000

6000 0

t/M

2000

4000

t/M

Kerr-Schild BH (0.9 J/M), excision with cube, 1 + log-lapse, Γ-driver shift.
22∂tΓ̃i = (· · ·) + Γ̃iβ i,j − (χ + )G iβ j ,j33∂tγ̃ij = (· · ·) − καγ̃ij H

χ = 2/3 for (A4)-(A8)κ = 0.1 ∼ 0.2 for (A5), (A6) and (A8)

2001
so-called BSSN
62

Shibata

87,95,99

ADM

2005

AEI

BSSN

PennState
Caltechhyperbolic formulation
01
92
92

Bona-Masso

04

Kidder-Scheel-Teukolsky

04

LSU

Z4 (Bona et.al.)Z4-lambda 05(Gundlach-Calabrese)harmonic
99

Nagy-Ortiz-Reulalambda system

Shinkai-Yoneda

asymptotically constrained /constraint damping
01
87

Detweiler02

adjusted-system
02

Pretorius
Illinois

adjusted BSSN (Yoneda-Shinkai)

adjusted ADM (Shinkai-Yoneda)

2001
so-called BSSN
62

UTB-Rochester

BSSN

Caltech92

Bona-Masso

04

Kidder-Scheel-Teukolsky

04

PennState

BSSN is “well-posed” ?
(Sarbach / Gundlach ...)LSU

Z4 (Bona et.al.)Z4-lambda 05(Gundlach-Calabrese)harmonic
99

Nagy-Ortiz-ReulaSouthampton

Jena

LSU

PennState

01

Parma

NASA-Goddard

AEI

hyperbolic formulation
92

AEI

Shibata

87,95,99

ADM

2005

lambda system

Shinkai-Yoneda

asymptotically constrained /constraint damping
01
87

Detweiler02

adjusted-system
02

Pretorius
Illinois

Kiuchi-Shinkai

adjusted BSSN (Yoneda-Shinkai)

adjusted ADM (Shinkai-Yoneda)

ShinkaiSome known fact (technical):• Trace-out Aij at every time step helps the stability.
Alcubierre, et al, [PRD 62 (2000) 044034]• “The essential improvement is in the process of replacing terms by the momentum
constraints”,
Alcubierre, et al, [PRD 62 (2000) 124011]• Γ̃i is replaced by −∂j γ̃ ij where it is not differentiated,
Campanelli, et al, [PRL96 (2006) 111101; PRD 73 (2006) 061501R]• Γ̃i-equation has been modified as suggested in Yo-Baumgarte-Shapiro [PRD 66(2002) 084026]
Baker et al, [PRL96 (2006) 111102; PRD73 (2006) 104002]

Some known fact (technical):• Trace-out Aij at every time step helps the stability.
Alcubierre, et al, [PRD 62 (2000) 044034]
This is because A-violation affects to all other constraint violations.
• “The essential improvement is in the process of replacing terms by the momentum
constraints”,
Alcubierre, et al, [PRD 62 (2000) 124011]
This is because M-replacement in Γi equation kills the positive real eigenvaluesof CAFs. eigenvalues• Γ̃i is replaced by −∂j γ̃ ij where it is not differentiated,
Campanelli, et al, [PRL96 (2006) 111101; PRD 73 (2006) 061501R]
This is because G-violation affects to H, Mi-violation constraint violations.

• Γ̃i-equation has been modified as suggested in Yo-Baumgarte-Shapiro [PRD 66(2002) 084026]
Baker et al, [PRL96 (2006) 111102; PRD73 (2006) 104002]
No doubt about this.

Numerical Experiments of Adjusted BSSN Systems
Kenta Kiuchi
木内 健太

Waseda University
早稲田大学 理工学部
kiuchi@gravity.phys.waseda.ac.jp
Hisa-aki Shinkai Osaka Institute of Technology
真貝寿明
大阪工業大学 情報科学部
shinkai@is.oit.ac.jp• BSSN vs adjusted BSSN Numerical tests• gauge-wave, linear wave, and Gowdy-wave tests, proposed by the Mexico workshop 2002
• 3 adjusted BSSN systems.
• Work as Expected– When the original BSSN system already shows satisfactory good evolutions (e.g., linear wave test),
the adjusted versions also coincide with those evolutions.
– For some cases (e.g., gauge-wave or Gowdy-wave tests) the simulations using the adjusted systems
last 10 times longer than the standard BSSN.
Phys. Rev. D77, 044010 (2008)

Adjusted BSSN systems; we tested
from the proposals in Yoneda & HS, Phys. Rev. D66 (2002) 124003

1. Ã-equation with the momentum constraint:∂tÃij = ∂tB Ãij + κAαD̃(iMj) ,(1)with κA > 0 (predicted from the eigenvalue analysis).
2. γ̃-equation with G constraint:∂t γ̃ij = ∂tB γ̃ij + κγ̃ αγ̃k(iD̃j)G k ,(2)∂t Γ̃i = ∂tB Γ̃i + κΓ̃αG i.

(3)with κγ̃ < 0.
3. Γ̃-equation with G constraint:with κΓ̃ < 0.

Numerical Testbed Models

A: Gauge-wave testbed

from the proposals in Mexico Workshop 2002, Class. Quant. Gravity 21 (2004) 589

The trivial Minkowski space-time, but time-dependent tilded slice.
2ds2 = −Hdt2 + Hdx2 + dy 2 + dz,
2π(x − t)
H = H(x − t) = 1 − A sin,dParameters:• Gauge-wave parameters: d = 1 and A = 10−2
• Simulation domain: x∈[−0.5, 0.5], y = z = 0• Grid: xi = −0.5 + (n − 12 )dx with n = 1, · · · 50ρ, where dx = 1/(50ρ) with ρ = 2, 4, 8
• Time step: dt = 0.25dx• Periodic boundary condition in x direction• Gauge conditions: ∂tα = −α2K,

β i = 0.

The 1D simulation is carried out for a T = 1000 crossing-time or until the code crashes, where one
crossing-time is defined by the length of the simulation domain.

Error evaluation methods
It should be emphasized that the adjustment effect has two meanings, improvement of stability and of
accuracy. Even if a simulation is stable, it does not imply that the result is accurate.
• We judge the stability of the evolution by monitoring the L2 norm of each constraint,s
1 X(C(t; x, y, z))2,||δC||2(t) ≡N x,y,zwhere N is the total number of grid points,• We judge the accuracy by the difference of the metric components gij (t; x, y, z) from the exact(exact)solution gij(t; x, y, z),
s
2
1 X
(exact).
gij − gij||δgij ||2(t) ≡N x,y,zNumerical Results

A: Gauge-wave test (1)A.1 The plain BSSN system
10000

100000
ρ=2
ρ=4
ρ=8Log10[||H||2]

100

ρ=2
ρ=4
ρ=8

10000
1000
Log10[||Mx||2]

1000
10
1
0.1
0.01

100
10
1
0.1

0.001

0.01

1e-04

0.001

1e-05

1e-04
0

20

40

60
T

80

100

120

0

20

40

60

80

100

120

T

FIG. 1: The one-dimensional gauge-wave test with the plain BSSN system. The L2 norm of H and Mx , rescaled by ρ2 /4, are plotted with a function of the crossing-time. The
amplitude of the wave is A = 0.01. The loss of convergence at the early time, near the 20 crossing-time, can be seen, and it will produce the blow-ups of the calculation in theend.

• The poor performance of the plain BSSN system has been reported.
Jansen, Bruegmann, & Tichy, PRD 74 (2006) 084022.
• The 4th-order finite differencing scheme improves the results.
Zlochower, Baker, Campanelli, & Lousto, PRD 72 (2005) 024021.

Numerical Results

A: Gauge-wave test (2)A.2 Adjusted BSSN with Ã-equation
0.1

0.1

Log10[||H||2]

Log10[||Mx||2]

ρ=2
ρ=4
ρ=8

0.01

0.001

ρ=2
ρ=4
ρ=8

0.01

0.001

1e-04

1e-05

1e-04
0

200

400

600

800

1000 1200

0

T

200

400

600

800

1000 1200

T

FIG. 2: The one-dimensional gauge-wave test with the adjusted BSSN system in the Ã-equation (1). The L2 norm of H and Mx , rescaled by ρ2 /4, are plotted with a functionof the crossing-time. The wave parameter is the same as with Fig. 1, and the adjustment parameter κA is set to κA = 0.005. We see the higher resolution runs show convergencelonger, i.e., the 300 crossing-time in H and the 200 crossing-time in Mx with ρ = 4 and 8 runs. All runs can stably evolve up to the 1000 crossing-time.

• We found that the simulation continues 10 times longer.

• Convergence behaviors are apparently improved than those of the plain BSSN.
• However, growth of the error in later time at higher resolution.

2∂t Ãij = −e−4φ [Di Dj α + αRij ]TF + αK Ãij − 2αÃik Ãkj + ∂i β k Ãkj + ∂j β k Ãki − ∂k β k Ãij + β k ∂k Ãij +κA αD̃(i Mj)
3

Numerical Results

A: Gauge-wave test (4)A.4 Evaluation of Accuracy• L2 norm of the error in γxx , (4), with the function of time.
• The error is induced by distortion of the wave; the both phase and amplitude errors.
1
(a)
plain
A-adj.
Γ-adj.

0.01
0.001

γxx

Log10[||δγxx||2]

0.1

1e-04
1e-05
1e-06
1e-07
0.1

1

10
Log10[T]

100

1000

1.05
1
0.95
(b)
0.9
Exact
0.85
plain
A-adj.
0.8
0.75
0.7
0.65
0.6
0.55
0.5
-0.6 -0.4 -0.2

0
x

0.2

0.4

0.6

FIG. 4: Evaluation of the accuracy of the one-dimensional gauge-wave testbed. Lines show the plain BSSN, the adjusted BSSN with A-equation, and with Γ̃-equation. (a) The
L2 norm of the error in γxx , using (4). (b) A snapshot of the exact and numerical solution at T = 100.

Numerical Testbed Models

B: Linear wave testbed

from the proposals in Mexico Workshop 2002, Class. Quant. Gravity 21 (2004) 589

Check the ability of handling a travelling gravitational wave.
ds2 = −dt2 + dx2 + (1 + b)dy 2 + (1 − b)dz 2,2π(x − t)b = A sindParameters:• Linear wave parameters: d = 1 and A = 10−8
• Simulation domain: x∈[−0.5, 0.5], y = 0, z = 0• Grid: xi = −0.5 + (n − 12 )dx with n = 1, · · · 50ρ, where dx = 1/(50ρ) with ρ = 2, 4, 8
• Time step: dt = 0.25dx• Periodic boundary condition in x direction• Gauge conditions: α = 1 and β i = 0

The 1D simulation is carried out for a T = 1000 crossing-time or until the code crashes.

Numerical Results

B: Linear Wave Test

1e-08

2e-10
1e-10
0
γzz-γzz(exact)

γzz-1

5e-09
0
-5e-09Exact
ρ=2-1e-08ρ=4
ρ=8-1.5e-08-0.6 -0.4-1e-10-2e-10-3e-10-4e-10plain
A-adj.

-5e-10
-0.2

0

0.2

0.4

0.6

-6e-10
-0.6

-0.4

-0.2

x

Snapshots of the one-dimensional linear wave at different resolutions with the plain
BSSN system at the simulation time 500 crossing-time. We see there exists phaseerror, but they are convergent away at higher resolution runs.

0

0.2

0.4

0.6

x
Snapshot of errors with the exact solution for the Linear Wave testbed with the plain
BSSN system and the adjusted BSSN system with the Ã equation at T = 500. The
highest resolution ρ = 8 is used in both runs. The difference between the plain andthe adjusted BSSN system with the Ã equation is indistinguishable. Note that the
maximum amplitude is set to be 10−8 in this problem.

• The linear wave testbed does not produce a significant constraint violation.
• The plain BSSN and adjusted BSSN results are indistinguishable.
This is because the adjusted terms of the equations are small due to the small violations of constraints.

Numerical Testbed Models

C: Collapsing polarized Gowdy-wave testbed

from the proposals in Mexico Workshop 2002, Class. Quant. Gravity 21 (2004) 589

Check the formulation in a strong field context using the polarized Gowdy metric.
ds2 = t−1/2eλ/2(−dt2 + dz 2 ) + t(eP dx2 + e−P dy 2).
P = J0(2πt) cos(2πz),λ = −2πtJ0(2πt)J1(2πt) cos2(2πz) + 2π 2t2[J02(2πt) + J12(2πt)]1− [(2π)2[J02(2π) + J12(2π)] − 2πJ0(2π)J1(2π)],2where Jn is the Bessel function.
Parameters:• Perform the evolution in the collapsing (i.e. backward in time) direction.
• Simulation domain: z ∈ [−0.5, 0.5], x = y = 0• Grid: z = −0.5 + (n − 12 )dz with n = 1, · · · 50ρ, where dz = 1/(50ρ) with ρ = 2, 4, 8
• Time step: dt = 0.25dz• Periodic boundary condition in z-direction• Gauge conditions: the harmonic slicing ∂t α = −α2K,

β i = 0. and β i = 0• Set the initial lapse function is 1, using coordinate transformation.
The 1D simulation is carried out for a T = 1000 crossing-time or until the code crashes.

Numerical Results

C: Collapsing polarized Gowdy-wave testbed (1)C.1 The plain BSSN
1000
100

100

1
0.1
0.01
0.001

10
1
0.1
0.01

1e-04

0.001

1e-05

1e-04

1e-06
-800

-600

-400
-200
t(harmonic)

ρ=2
ρ=4
ρ=8

1000
Log10[||Mz||2]10
Log10[||H||2]

10000

ρ=2
ρ=4
ρ=8

0

200

1e-05
-800

-600

-400

-200

0

200

t(harmonic)

FIG. 5: Collapsing polarized Gowdy-wave test with the plain BSSN system. The L2 norm of H and Mz , rescaled by ρ2 /4, are plotted with a function of the crossing-time.
(Simulation proceeds backwards from t = 0.) We see almost perfect overlap for the initial 100 crossing-time, and the higher resolution runs crash earlier. This result is quite
similar to those achieved with the Cactus BSSN code, reported by [? ].

• Our result shows similar crashing time with that of Cactus BSSN code.
Alcubierre et al. CQG 21, 589 (2004)• Higher order differencing scheme with Kreiss-Oliger dissipation term improves the results.
Zlochower, Baker, Campanelli & Lousto, PRD 72, 024021 (2005)

Numerical Results

C: Collapsing polarized Gowdy-wave testbed (2)C.2 Adjusted BSSN with Ã-equation100
Log10[||H||2]

10

10000

ρ=2
ρ=4
ρ=8

1000
100
Log10[||Mz||2]

1000

1
0.1
0.01
0.001

ρ=2
ρ=4
ρ=8

10
1
0.1
0.01

1e-04

0.001

1e-05

1e-04

1e-06
-1000 -800 -600 -400 -200

1e-05
-1000 -800 -600 -400 -200

0200t(harmonic)

0200

t(harmonic)

FIG. 6: Collapsing polarized Gowdy-wave test with the adjusted BSSN system in the Ã-equation (1), with κA = −0.001. The style is the same as in Fig. 5 and note that bothconstraints are normalized by ρ2 /4. We see almost perfect overlap for the initial 1000 crossing-time in both constraint equations, H and Mz , even for the highest resolution run.

• Adjustment extends the life-time of the simulation 10 times longer.

• Almost perfect convergence upto t = 1000tcross for both H and Mz , while we find oscillations in Mz
later time.
2∂t Ãij = −e−4φ [Di Dj α + αRij ]TF + αK Ãij − 2αÃik Ãkj + ∂i β k Ãkj + ∂j β k Ãki − ∂k β k Ãij + β k ∂k Ãij +κA αD̃(i Mj)
3

Numerical Results

C: Collapsing polarized Gowdy-wave testbed (3)

10000
ρ=2
1000
ρ=4
100
ρ=8
10
1
0.1
0.01
0.001
1e-04
1e-05
1e-06
-1000 -800 -600 -400 -200
t(harmonic)

10000
1000
100
Log10[||Mz||2]

Log10[||H||2]

C.3 Adjusted BSSN with γ̃-equation
ρ=2
ρ=4
ρ=8

10
1
0.1
0.01
0.001
1e-04

0

200

1e-05
-1000 -800 -600 -400 -200

0200t(harmonic)

FIG. 7: Collapsing polarized Gowdy-wave test with the adjusted BSSN system in the γ̃-equation (2), with κγ̃ = 0.000025. The figure style is the same as Figure 5. Note thealmost perfect overlap for 200 crossing-time in the both the Hamiltonian and Momentum constraint and the ρ = 2 run can evolve stably for 1000 crossing-time.

• Almost perfect convergence up to t = 200tcross in both H and Mz .
2∂t γ̃ij = −2αÃij + γ̃ik ∂j β k + γ̃jk ∂iβ k − γ̃ij ∂k β k + β k ∂k γ̃ij + κγ̃ αγ̃k(i D̃j)G k
3

Numerical Results

C: Collapsing polarized Gowdy-wave testbed (4)C.4 Adjustment works for Accuracy
Error of γzz to the exact solution normalized by γzz .
• Accurate Evolution ⇔ Error < 1 %.
(Zlochower, et al., PRD72 (2005) 024021 )the Plain BSSN ≈ t = 200tcrossadjusted BSSN Ã-eq ≈ t = 1000tcrossadjusted BSSN γ̃-eq ≈ t = 400tcross

Log10[||δγzz||2/||γzz||2]

1
0.1
0.01
0.001
1e-04
1e-05
1e-06

plain
A-adj.
γ-adj.

1e-07
1e-08-1000 -800 -600 -400 -200
t(harmonic)

0200

Comparisons of systems in the collapsing polarized Gowdy-wave test.
The L2 norm of the error in γzz , rescaled by the L2 norm of γzz , for theplain BSSN, adjusted BSSN with Ã-equation, and with γ̃-equation are
shown. The highest resolution run, ρ = 8, is depicted for the plots. We
can conclude that the adjustments make longer accurate runs available.
Note that the evolution is backwards in time.

A

Full set of BSSN constraint propagation eqs.










HBS  A11
A12 A13
A14
A15
HBS 







 M 
 −(1/3)(∂ α) + (1/6)∂
 M 
αK A23
0
A25
i 
i
i
j 











ij
j 
 Gi  = 



0
αγ̃
0
A
A
G
∂tBS 
34
35










k





0
0
0
β
(∂
S)
−2αγ̃
 S 

 S 
k





k
A
0
0
0
0
αK + β ∂kA
A11 = +(2/3)αK + (2/3)αA + β k ∂kA12 = −4e−4ϕ α(∂k ϕ)γ̃ kj − 2e−4ϕ (∂k α)γ̃ jk
A13 = −2αe−4ϕ Ãk j ∂k − αe−4ϕ (∂j Ãkl )γ̃ kl − e−4ϕ (∂j α)A − e−4ϕ β k ∂k ∂j − (1/2)e−4ϕ β k γ̃ −1 (∂j S)∂k+(1/6)e−4ϕ γ̃ −1 (∂j β k )(∂k S) − (2/3)e−4ϕ (∂k β k )∂j
A14 = 2αe−4ϕ γ̃ −1 γ̃ lk (∂l ϕ)A∂k + (1/2)αe−4ϕ γ̃ −1 (∂l A)γ̃ lk ∂k + (1/2)e−4ϕ γ̃ −1 (∂l α)γ̃ lk A∂k + (1/2)e−4ϕ γ̃ −1 β m γ̃ lk ∂m ∂l ∂k−(5/4)e−4ϕ γ̃ −2 β m γ̃ lk (∂m S)∂l ∂k + e−4ϕ γ̃ −1 β m (∂m γ̃ lk )∂l ∂k + (1/2)e−4ϕ γ̃ −1 β i (∂j ∂i γ̃ jk )∂k+(3/4)e−4ϕ γ̃ −3 β i γ̃ jk (∂i S)(∂j S)∂k − (3/4)e−4ϕ γ̃ −2 β i (∂i γ̃ jk )(∂j S)∂k + (1/3)e−4ϕ γ̃ −1 γ̃ pj (∂j β k )∂p ∂k−(5/12)e−4ϕ γ̃ −2 γ̃ jk (∂k β i )(∂i S)∂j + (1/3)e−4ϕ γ̃ −1 (∂k γ̃ ij )(∂j β k )∂i − (1/6)e−4ϕ γ̃ −1 γ̃ mk (∂k ∂l β l )∂m
A15 = (4/9)αKA − (8/9)αK 2 + (4/3)αe−4ϕ (∂i ∂j ϕ)γ̃ ij + (8/3)αe−4ϕ (∂k ϕ)(∂l γ̃ lk ) + αe−4ϕ (∂j γ̃ jk )∂k+8αe−4ϕ γ̃ jk (∂j ϕ)∂k + αe−4ϕ γ̃ jk ∂j ∂k + 8e−4ϕ (∂l α)(∂k ϕ)γ̃ lk + e−4ϕ (∂l α)(∂k γ̃ lk ) + 2e−4ϕ (∂l α)γ̃ lk ∂k+e−4ϕ γ̃ lk (∂l ∂k α)kl
A23 = αe−4ϕ γ̃ km (∂k ϕ)(∂j γ̃mi ) − (1/2)αe−4ϕ Γ̃mkl γ̃ (∂j γ̃mi )+(1/2)αe−4ϕ γ̃ mk (∂k ∂j γ̃mi ) + (1/2)αe−4ϕ γ̃ −2 (∂i S)(∂j S) − (1/4)αe−4ϕ (∂i γ̃kl )(∂j γ̃ kl ) + αe−4ϕ γ̃ km (∂k ϕ)γ̃ji ∂m
kl−4ϕ mk+αe−4ϕ (∂j ϕ)∂i − (1/2)αe−4ϕ Γ̃mγ̃ Γ̃ijk ∂m + (1/2)αe−4ϕ γ̃ lk γ̃ji ∂k ∂lkl γ̃ γ̃ji ∂m + αe+(1/2)e−4ϕ γ̃ mk (∂j γ̃im )(∂k α) + (1/2)e−4ϕ (∂j α)∂i + (1/2)e−4ϕ γ̃ mk γ̃ji (∂k α)∂m
A25 = −Ãk i (∂k α) + (1/9)(∂i α)K + (4/9)α(∂i K) + (1/9)αK∂i − αÃk i ∂k
A34 = −(1/2)β k γ̃ il γ̃ −2 (∂l S)∂k − (1/2)(∂l β i )γ̃ lk γ̃ −1 ∂k + (1/3)(∂l β l )γ̃ ik γ̃ −1 ∂k − (1/2)β l γ̃ in (∂l γ̃mn )γ̃ mk γ̃ −1 ∂k+(1/2)β k γ̃ il γ̃ −1 ∂l ∂k
A35 = −(∂k α)γ̃ ik + 4αγ̃ ik (∂k ϕ) − αγ̃ ik ∂kWhich constraint should be monitored?
Yoneda & HS, PRD 66 (2002) 124003

Kiuchi & HS, PRD 77 (2008) 044010

Order of constraint violation?
• A and S constraints propagate independently ofthe other constraints.
• G-constraint is triggered by the violation of the
momentum constraint.
• H and M constraints are affected by all the otherconstraints.

Normalized Constraint

1e+08

S(t)/S(t=0)A(t)/A(t=0)gx(t)/gx(t=0)
H(t)/H(t=0)
Mx(t)/Mx(t=0)

1e+07
1e+06
100000
10000
1000
100
10
1
0.1

0

20

40

60
T

80

100

120

The violation of all constraints normalized with their initial values, ||δC||2 (t)/||δC||2 (0),are plotted with a function of time. The evolutions of the gauge-wave testbeds withthe plain BSSN system are shown.

By observing which constraint triggers the other constraint’s violation from the constraint propagation
equations, we may guess the mechanism by which the entire system is violating accuracy and stability.

ここまでのまとめ (後半)
[Keyword 1]
Adjusted Systems
運動方程式を拘束条件を用いて補正するのは普通に使われている手段．
これを手法として明確化し，lambda システムのように対称双曲型にこだわらない．
[Keyword 2]
Constraint Propagation Analysis -> Constraint Damping System
拘束条件の発展方程式を計算すれば，運動方程式をどう補正すればよいかがわかる．
(Step 1)
Fourier mode expression of all terms of constraint propagation eqs.
(Step 2)
Eigenvalues and Diagonalizability of constraint propagation matrix.
Eigenvalues = Constraint Ampliﬁcation Factors(Step 3)
If CAF=negatives -> Constraint surface becomes the attractor.
[Keyword 3]
Adjusted ADM systems
Standard ADM は，拘束条件を破るモードが存在していることがわかった．
いくつかの補正候補があり，実際に計算は安定化して寿命は伸びる．
[Keyword 3]
Adjusted BSSN systems
BSSN の長所は momentum constraint を用いた補正にあった．
いくつかのさらなる補正候補があり，実際に計算は安定化して寿命は伸びる．

Goals of the Talk
何を指標にして発展方程式を選択すれば良いのか？
どうして多くのグループが BSSN
形式を使っているのか？
BSSN 形式に代わる formulation
はあるか？

Goals of the Talk
何を指標にして発展方程式を選択すれば良いのか？
--- Constraint Propagation eqs.
どうして多くのグループが BSSN
形式を使っているのか？
BSSN 形式に代わる formulation
はあるか？

Goals of the Talk
何を指標にして発展方程式を選択すれば良いのか？
--- Constraint Propagation eqs.
どうして多くのグループが BSSN
形式を使っているのか？
--- Just rush, not to be late.
BSSN 形式に代わる formulation
はあるか？

Goals of the Talk
何を指標にして発展方程式を選択すれば良いのか？
--- Constraint Propagation eqs.
どうして多くのグループが BSSN
形式を使っているのか？
--- Just rush, not to be late.
BSSN 形式に代わる formulation
はあるか？
--- Yes, there are. But we do not know the best.

2008/11/28 @ 大阪市立大学

Formulation for Numerical Relativity
both for Einstein / Gauss-Bonnet
真貝寿明

大阪工業大学情報科学部

1. Introduction
定式化問題？
2. The Standard Approach to Numerical Relativity
ADM 形式，BSSN 形式，Hyperbolic 形式
3. Robust system for Constraint Violation
Adjusted systems .... better than lambda system!
Adjusted ADM......... why ADM blows up?
Adjusted BSSN......... why BSSN works well?
4. 高次元数値相対論に向けてDiscussion
Application 1 : Constraint Propagation in N + 1 dim. space-time
HS-Yoneda, GRG 36 (2004) 1931

Dynamical equation has N -dependency
Only the matter term in ∂tKij has N -dependency.
1
0 ≈ CH ≡ (Gµν − 8πTµν )nµnν = ((N) R + K 2 − K ij Kij ) − 8πρH − Λ,2µ ν
0 ≈ CMi ≡ (Gµν − 8πTµν )n ⊥i = Dj Kij − DiK − 8πJi,∂tγij = −2αKij + Dj βi + Diβj ,∂tKij = α(N)Rij + αKKij − 2αK ℓj Kiℓ − DiDj α2α1γij T  −γij Λ,+β k (Dk Kij ) + (Dj β k )Kik + (Diβ k )Kkj − 8πα Sij −
N −1
N −1

Constraint Propagations remain the same
From the Bianchi identity, ∇ν Sµν = 0 with Sµν = Xnµnν + Yµnν + Yν nµ + Zµν , we get
0 = nµ ∇ν Sµν = −Zµν (∇µnν ) − ∇µYµ + Yν nµ∇µ nν − 2Yµnν (∇ν nµ ) − X(∇µnµ ) − nµ (∇µX),
0 = hiµ ∇ν Sµν = ∇µZiµ + Yi(∇µnµ) + Yµ (∇µni) + X(∇µni)nµ + nµ(∇µ Yi).
• (Sµν , X, Yi, Zij ) = (Tµν , ρH , Ji, Sij ) with ∇µ Tµν = 0 ⇒ matter eq.
• (Sµν , X, Yi, Zij ) = (Gµν − 8πTµν , CH , CMi , κγij CH ) with ∇µ(Gµν − 8πTµν ) = 0 ⇒ CP eq.

N + 1 formalism in Einstein-Gauss-Bonnet Gravity
Takashi Torii
鳥居 隆
Hisa-aki Shinkai
真貝寿明

Osaka Institute of Technology
大阪工業大学 工学部
Osaka Institute of Technology
大阪工業大学 情報科学部• (N + 1)-dimensional space-time decomposition of Einstein-Gauss-Bonnet gravity• Due to the quasi-linear property of the Gauss-Bonnet gravity, we ﬁnd that theevolution equations can be in a treatable form in numerics.
• We also show the conformally-transformed constraint equations for constructing an initial data.
• Both for timelike and spacelike foliations.
Phys Rev D 78, 084037 (2008).

Einstein-Gauss-Bonnet action• (N + 1)-dimensional spacetime (M, gµν )∫N +1

S = Md√1X −g  2 (R − 2Λ + αGB LGB ) + Lmatter2κ(1)

LGB = R2 − 4Rµν Rµν + Rµνρσ Rµνρσ• The action gives the gravitational equation
Gµν + αGB Hµν = κ2 Tµν

(2)

where
1
Gµν = Rµν − gµν R + Λgµν ,
2
Hµν =

[

2 RRµν − 2Rµα Rαν − 2Rαβ Rµανβ + Rµαβγ Rναβγ

Tµν ≡ −2

δLmatter+ gµν Lmatter.
δg µν

]

1
− gµν LGB ,2

Projections to Hypersurface ΣN (spacelike or timelike)(1)• the projection operator,⊥µν = gµν − εnµnν ,nµ n µ = ε(3)where nµ is the unit-normal vector to Σ with nµ is timelike (if ε = −1) or spacelike (if ε = 1).
Σ is spacelike (timelike) if nµ is timelike (spacelike).
• The projections of the gravitational equation:(Gµν + αGB Hµν ) nµ nν = κ2 Tµν nµ nν =: κ2ρH ,(4)(Gµν + αGB Hµν ) nµ ⊥νρ = κ2 Tµν nµ ⊥νρ =: −κ2Jρ,(5)(Gµν + αGB Hµν ) ⊥µρ ⊥νσ = κ2 Tµν ⊥µρ ⊥νσ =: κ2Sρσ ,(6)where we deﬁned
Tµν = ρH nµnν + Jµnν + Jν nµ + Sµν , T = −ρH + S ``• Introduce the extrinsic curvature Kij
1
Kij := − £nhij = −⊥αi⊥βj ∇α nβ ,(7)2where £n denotes the Lie derivative in the n-direction and ∇ and Di is the covariant diﬀerentiation with respect to gµν and γij , respectively.

Projections to Hypersurface ΣN (spacelike or timelike)(2)• Projection of the (N + 1)-dimensional Riemann tensor onto ΣN
Gauss eq.
Codacci eq.

Rαβγδ ⊥αi ⊥βj ⊥γk ⊥δl = Rijkl − εKik Kjl + εKil Kjk ,(8)Rαβγδ ⊥αi ⊥βj ⊥γk nδ = −2D[iKj]k ,(9)Rαβγδ ⊥αi ⊥γk nβ nδ = £nKik + Ki` K `k ,(10)

• Curvature relations
Rµνρσ = Rµνρσ − ε(KµρKνσ − Kµσ Kνρ − nµDρKνσ + nµDσ Kρν + nν DρKσµ − nν Dσ Kρµ−nρDµKνσ + nρDν Kµσ + nσ DµKνρ − nσ Dν Kµρ)+nµnρKνα K ασ − nµnσ Kνα K αρ − nν nρKµα K ασ + nν nσ Kµα K αρ+nµnρ£nKνσ − nµnσ £nKνρ − nν nρ£nKµσ + nν nσ £nKµρ,Rµν =Rµν − ε[KKµν − 2Kµα K αν + nµ (Dα K αν − Dν K) + nν
+nµnν Kαβ K αβ + ε£nKµν + nµnν γ αβ £nKαβ ,R = R − ε(K 2 − 3Kαβ K αβ − 2γ αβ £nKαβ ).

(Dα K αµ − DµK

(11)
)

]
(12)
(13)

N + 1 Einstein-Gauss-Bonnet equations

Substituting (11)-(13) into (3) or (4)-(6), we ﬁnd:(a) dynamical equations for γij :1Mij − M γij − ε(−KiaK aj + γij KabK ab − £nKij + γij γ ab£nKab)2+2αGB [Hij + ε(M £nKij − 2Mi a£nKaj − 2Mja£nKai − Wij ab£nKab)] = κ2Tµν γ µiγ νj(b) Hamiltonian constraint equation:
M + αGB (M 2 − 4MabM ab + MabcdM abcd) = −2εκ2Tµν nµnν(c) momentum constraint equation:(

Ni + 2αGB M Ni − 2Mi Na + 2M Niab − MiM=

Nijk =Ni =Wijkl =

ab

cab

)

Nabc = −κ2Tµν nµγ νiHij = M Mij − 2(Mia M aj + M ab Miajb ) + Miabc Mjabc[
1aba−2ε −Kab K ab Mij − M Kia K aj + Kia K ab M bj + Kja K ab M bi + K ac Kcb Miajb
γ Miajb = Rij − ε(KKij − Kia K j )2]
γ ab Mab = R − ε(K 2 − Kab K ab )
1
a
ab
ab
+Ni Nj − N (Naij + Naji ) − Nabi N j − Niab Nj
2
Di Kjk − Dj Kikaba1γ Naib = Da Ki − Di K
− γij [M 2 − 4Mab M ab + Mabcd M abcd ]
4
kl
kl
kl
ak bl
M γij γ − 2Mij γ − 2γij M + 2Miajb γ γ−εγij [Kab K ab M − 2Mab K ac Kcb − 2Na N a + Nabc N abc ]Mijkl = Rijkl − ε(Kik Kjl − Kil Kjk )
Mij =

a

Conformal Approach to solve constraints : Eqs. for Initial Data construction• We generalized the Conformal approach by York and ÓMurchadha (1974) to N-dim & forGauss-Bonnet gravity.
• Conformal transformation
solutionγij = ψ 2mγ̂ij γ ij = ψ −2mγ̂ ij

trial metric

this gives{}

R = ψ −2m R̂ − 2(N − 1)mψ −1 (D̂a D̂a ψ) + (N − 1)[2 − (N − 2)m]mψ −2 (D̂ψ)2 ,Rij= R̂ij − mγ̂ij ψ −1 D̂a D̂a ψ − (N − 2)mψ −1 D̂i D̂j ψ + (N − 2)m(m + 1)ψ −2 D̂i ψ D̂j ψ − m[(N − 2)m − 1]ψ −2 (D̂ψ)2 γ̂ij ,{Rijkl = ψ 2m R̂ijkl + mψ −1 γ̂il [D̂j D̂k ψ − (m + 1)ψ −1 D̂j ψ D̂k ψ] − mψ −1 γ̂ik [D̂j D̂l ψ − (m + 1)ψ −1 D̂j ψ D̂l ψ]}+ mψ −1 γ̂jk [D̂i D̂l ψ − (m + 1)ψ −1 D̂i ψ D̂l ψ] − mψ −1 γ̂jl [D̂i D̂k ψ − (m + 1)ψ −1 D̂i ψ D̂k ψ] + m2 ψ −2 (D̂ψ)2 (γ̂il γ̂jk − γ̂ik γ̂jl ) .

1
γij K, and assumeN
Aij = ψ `Âij , Aij = ψ `−4mÂij ,• Decompose the extrinsic curvature Kij as Kij ≡ Aij +K = ψ τ K̂
• When matter exists, deﬁne also the conformal transformationρ = ψ −pρ̂,

J i = ψ −q Jˆi

Hamiltonian constraint2(N − 1)mD̂aD̂aψ − (N − 1)[2 − (N − 2)m]m(D̂ψ)2ψ −1
N − 1 2m+2τ +1 2
= R̂ψ −εψK̂ + εψ −2m+2`+1ÂabÂab + 2εκ2ρ̂ψ −p − 2Λ̂
(N
)
2
ab
abcd
+ αGB M − 4MabM + MabcdM
ψ 2m+1.
(Θ̂ =M 2 − 4Mab M ab + Mabcd M abcd

= (N − 3)mψ

−4m

{

(14)

)
[

]

[

]

4(N − 2)mψ −2 (D̂a D̂a ψ)2 − D̂a D̂b ψ D̂a D̂b ψ − 4ψ −1 M̂ − (N − 2)[(N − 3)m − 2]mψ −2 D̂a ψ D̂a ψ D̂a D̂a ψ][+ 8ψ −1 M̂ ab + (N − 2)m(m + 1)ψ −2 D̂a ψ D̂b ψ D̂a D̂b ψ + (N − 1)2 m2 [(N − 4)m − 4]ψ −4 (D̂a ψ D̂a ψ)2 − 2ψ −2 [(N − 4)m − 2]M̂ D̂c ψ D̂c ψ− 8(m + 1)ψ

−2

[

where

}

M̂ D̂a ψ D̂b ψ + ψ −4m (Υ̂2 − 4Υ̂ab Υ̂ab + Υ̂abcd Υ̂abcd ),
ab

]

[

]

N − 1 2m+2τ 2
N − 1 2m+2τ
N − 2 `+τψψ K̂ Âij − ψ 2`−2m Âia Âaj ,Υ̂ = R̂ − εK̂ − ψ 2`−2m Âab Âab , Υ̂ij = R̂ij − εψγ̂ij K̂ 2 +
2
N
N
N
]
[
1 2m+2τ
1 `+τ
2
2`−2m
Υ̂ijkl = R̂ijkl − εψ(γ̂ik γ̂jl − γ̂il γ̂jk )K̂ + ψ (Âik γ̂jl − Âil γ̂jk + Âjl γ̂ik − Âjk γ̂il ) + ψ(Âik Âjl − Âil Âjk ) .
N2N(A) If we specify τ = ` − 2m and m = 2/(N − 2), then (14) becomes4(N − 1)D̂a D̂a ψ = R̂ψ − εψ 2`+1−4/(N −2) (K̂ 2 − K̂ab K̂ ab ) + 2εκ2 ρ̂ψ −p − 2Λ̂ + αGB Θ̂ψ 1+4/(N −2) .
N −2

(15)

(B) If we specify τ = 0 and m = 2/(N − 2), then (14) becomes4(N − 1)N − 1 1+4/(N −2) 2D̂a D̂a ψ = R̂ψ − εψK̂ + εψ 2`+1−4/(N −2) Âab Âab + 2εκ2 ρ̂ψ −p − 2Λ̂ + αGB Θ̂ψ 1+4/(N −2) . (16)N −2N

Momentum constraint• Introduce the TT part and the longitudinal part of Âij , and its vector potential as
ij
ij
ij
D̂j Âij
T T = 0, ÂL = Â − ÂT T ,

i
j
j
i
Âij
L = D̂ W + D̂ W −2 ijγ̂ D̂k W k .
N• Conformal transformations: Dj Ai j = ψ `−2m {D̂j Âi j + ψ −1 [` + m(N − 2)]Âi j D̂j ψ}

D̂a D̂a Wi +N −2
D̂i D̂k W k + R̂ik W kN

2 abγ̂ D̂k W k )γ̂bi D̂a ψ
NN−1−ψ 2m−`
D̂i (ψ τ K̂) + ψ 2m−` 2αGB Ξ̂i = κ2 ψ 4m−`−q Jˆi (See next page for Ξ̂i .)N+ψ −1 [` + (N − 2)m](D̂a W b + D̂b W a −(17)(A) If we specify τ = ` − 2m and m = 2/(N − 2), then (17) becomes
(
)
N −2
2
D̂i D̂k W k + R̂ik W k + ψ −1 (` + 2) D̂a W b + D̂b W a − γ̂ ab D̂k W k γ̂bi D̂a ψ
N
N
[(
]
)
4
N −1
`−
(D̂i ψ)K̂ + D̂i K̂ + ψ −`+4/(N −2) 2αGB Ξ̂i = κ2 ψ 8/(N −2)−`−q Jˆi
−
N
N −2

D̂a D̂a Wi +(18)

(B) If we specify τ = 0 and m = 2/(N − 2), then (17) becomes
N −2
2
D̂i D̂k W k + R̂ik W k + ψ −1 (` + 2)[D̂a W b + D̂b W a − γ̂ ab D̂k W k ]γ̂bi D̂a ψ
N
N
N
−
1
D̂i K̂ + ψ 4/(N −2)−` 2αGB Ξ̂i = κ2 ψ 8/(N −2)−`−q Jˆi−ψ 4/(N −2)−`ND̂a D̂a Wi +

(19)

{

}

N 2 − 3N + 4 +2m+2τ 2
R̂ − 2(N − 3)mψ D̂b D̂ ψ − (N − 3)m[(N − 4)m + 2]ψ D̂b ψ D̂ ψ −Ξi = ψεψK̂ − εψ 2`−2m Âbc Âbc D̂a Âai
N2}{
2(N − 3) `+τ
`−4m
b
2`−2m
c
b
b
−1 b
−2
b
+ψ
εψ K̂ Âi − 2εψÂi Âc D̂a Âab
−2R̂ i + 2(N − 3)mψ D̂ D̂i ψ − 2(N − 3)m(m + 1)ψ D̂i ψ D̂ ψ +N}{
2(N − 3) `+τεψ K̂ Âab + 2εψ 2`−2m Âac Âcb (D̂i Âab − D̂a Âib )+ψ `−4m 2R̂ab − 2(N − 3)mψ −1 D̂b D̂a ψ − 2(N − 1)m(m + 1)ψ −2 D̂a ψ D̂b ψ −
N
3`−6m
a bc
(1)
(2)
(3)
+2εψ
Âi Â (D̂a Âbc − D̂b Âac ) + Ri + Di + A D̂i ψ + A D̂i K̂ + A D̂a ψ Âai
2(N − 2)3 `−2m+2τ
2(N − 3)−[(N − 4)m + 2` + τ ]εψ 2`−4m+τ −1 K̂ D̂b ψ Âba ÂaiεψK̂(D̂a K̂ + τ ψ −1 K̂ D̂a ψ)Âai +2
NN2(N − 3) 2`−4m+τ
+εψD̂b K̂ Âba Âai − 2[(N − 6)m + 3`]εψ 3`−6m−1 D̂c ψ Âcb Âba Âai ,
N
where
{
}
N − 3 −2m+τ
`−4m−1a−1Ri =[(N − 3)m + `]ψψ(D̂i K̂ + τ ψ K̂ D̂i ψ) R̂
Âi D̂a ψ −N}{2(N − 3) −2m+τ −1
2(N − 3) −2m+ττψK̂ D̂a ψ +ψD̂a K̂ − 2[(N − 3)m + `]ψ `−4m−1 Âab D̂b ψ R̂i a+NN`−4m−1−2(m − `)ψ(Âab D̂i ψ − Âib D̂a ψ)R̂ab + 2(m − `)ψ `−4m−1 D̂a ψ Âbc R̂i cab ,{ 2}
N − 8N + 112`−4m−2b−2m+τ −1−1Di =D̂b ψ Â i D̂a D̂a ψmψ
(D̂i K̂ + τ ψ K̂ D̂i ψ) − 2m[(N − 6N + 7)m + (N − 3)`]ψN{}
2(N − 2)3−mψ −2m+τ −1 (D̂a K̂ + τ ψ −1 K̂ D̂a ψ) − 2m[(N 2 − 4N + 5)m + (N − 2)(` − 2)]ψ `−4m−2 D̂b ψ Âba D̂a D̂i ψN+2(N − 3)m(m − `)ψ `−4m−2 (Âab D̂i ψ − Âia D̂b ψ)D̂b D̂a ψ,{
N −2
m(m + 1)ψ −2m+2τ −2 (D̂a K̂ + τ ψ −1 K̂ D̂a ψ)(D̂a K̂ + τ ψ −1 K̂ D̂a ψ)
A(1) = 2 −N
N −2(N − 2)2m(m + 1)ψ −2m+τ −2 D̂a K̂ D̂a ψ +m[(N 2 − 4N + 5)m + 2]τ ψ −2m+τ −3 K̂ D̂a ψ D̂a ψ
+
N
2N
}
N −3
− (N − 2)3 m2 (m + 1)ψ `−4m−3 D̂a ψ D̂b ψ Âab +(m − ` − τ )εψ 2`−4m+τ −1 K̂ Âab Âab − (m − `)εψ 3`−6m−1 Âab Âbc Âca ,N{}
1(N − 1)2 (N + 1) 3τ 2
(2)−2m+τ −2a2`−4m+τabA=(N − 2)3 m[(N − 3)m − 2]ψ
εψ K̂ − (N − 3)εψD̂a ψ D̂ ψ −Âab Â ,N
N2
A(3) = −m[(N − 2)2 (N − 5)m2 + (N − 2)3 (` − 2)m + (N − 1)(3` − 2)]ψ `−4m−3 D̂a ψ D̂a ψ
1− 2 [(N − 1)(N 2 − 8)m + (N 2 − N + 2)`]εψ `−2m+2τ −1 K̂ 2 + [(N − 6)m + 3`]εψ 3`−6m−1 Âab Âab .
N
`−4m

−1

b

−2

b

(A) Hamiltonian constraint
4(N − 1)D̂a D̂a ψ = R̂ψ − εψ 2`+1−4/(N −2) (K̂ 2 − K̂ab K̂ ab ) + 2εκ2 ρ̂ψ −p − 2Λ̂ + αGB Θ̂ψ 1+4/(N −2) .
N −2(A) momentum constraint
2N −2
D̂i D̂k W k + R̂ik W k + ψ −1 (` + 2)(D̂a W b + D̂b W a − γ̂ ab D̂k W k )γ̂bi D̂a ψ
N
N
[
]
N −1
4
−
(` −
)(D̂i ψ)K̂ + D̂i K̂ + ψ −`+4/(N −2) 2αGB Ξ̂i = κ2 ψ 8/(N −2)−`−q JˆiN
N −2

D̂a D̂a Wi +Procedures to construct the initial hypersurface data (γij , Kij , ρ, J i )ˆ
1. Give the initial assumption (trial values) for γ̂ij , trK, ÂTijT and ρ̂, J.
2. Solve above 2 equations for ψ and W i .
3. inverse conformal transformations,
γij = ψ 4/(N −2) γ̂ij ,ρ = ψ −p ρ̂,Kij = ψ ` [ÂTijT + (l̂W )ij ] +J i = ψ −q Jˆi

1 `−4/(N −2)ψγ̂ij trK,N

N + 1 Einstein-Gauss-Bonnet evolution equations(1 + 2αGB M )£nKij − (γij γ ab + 2αGB Wij ab )£nKab − 8αGB M(ia£nK|a|j)()
1= −ε Mij − M γij − KiaK aj + γij KabK ab + εκ2Sij − εγij Λ − 2εαGB Hij , (20)2• £nKµν terms appear only in the linear form, due to the quasi-linear property of the GaussBonnet gravity.
• Iterative scheme is necessary, but treatable in numerics.




































£nγ11 
£nγ12 









£nγ13 




...






 = 



£nK11 







£nK12 








£nK13 



...

O

O

O

Mixing






































£nγ11 
 K11 









£nγ12  K12 









 K13 £nγ13 






.
...



.



.




+





£nK11 









 Source £nK12 












£nK13 








...

• Coding is in progress, .... but .... Are the evolution eqs always invertible??

現時点でのまとめ (高次元への応用)[高次元 GR 初期値決定問題]
conformal 変換を用いて Constraint を解く方法は一般的に導出できた．
[高次元 GR 時間発展問題]
4 次元と同じく，formulation problem が存在し，4 次元と同じ解決方法が適用されてゆくであろう．
[Gauss-Bonnet 初期値問題]
conformal 変換を用いて Constraint を解く方法は一般的に導出できたが，解けるかどうかはわからない．
[Gauss-Bonnet 時間発展問題]
ADM 的な定式化はできた．
理論上，状況によっては時間発展が追えなくなる可能性も考えられる．

Discussion
Future : Construct a robust adjusted system
HS-Yoneda, in preparation(1) dynamic & automatic determination of κ under a suitable principle.
e.g.) Efforts in Multi-body Constrained Dynamics simulations∂C a∂
p i = Fi + λ a
,
∂t
∂xi

with C a(xi, t) ≈ 0• J. Baumgarte (1972, Comp. Methods in Appl. Mech. Eng.)
Replace a holonomic constraint ∂t2C = 0 as ∂t2C + α∂t C + β 2C = 0.
• Park-Chiou (1988, J. Guidance), “penalty method”
Derive “stabilization eq.” for Lagrange multiplier λ(t).
• Nagata (2002, Multibody Dyn.)
Introduce a scaled norm, J = C T SC, apply ∂tJ + w2 J = 0, and adjust λ(t).
e.g.) Efforts in Molecular Dynamics simulations• Constant pressure· · · · · · potential piston!

• Constant temperature · · · · · · potential thermostat!! (Nosé, 1991, PTP)

(2) target to control each constraint violation byadjusting multipliers.

Momentum constr.
grow

CP-eigenvectors indicate directions of constraint grow/decay, if CP-matrix is diagonalizable.

grow
decay

(3) clarify the reasons of non-linear violation in thelast stage of current test evolutions.

Hamiltonian constr.

decay(4) Alternative new ideas?
– control theories, optimization methods (convex functional theories), mathematical programming methods, or ....
(5) Numerical comparisons of formulations, links to other systems, ...
– “Comparisons of Formulations” (e.g. Mexico NR workshop, 2002-2003); more formulations to be tested, ...

Find a RECIPE for all. Avoid un-essential techniques.

Discussion
Application 2 : Constraint Propagation of Maxwell field in Curved space
HS-Yoneda, in preparation
Towards a robust GR-MHD system:• Maxwell eqs in curved space-time∂tE i = ǫijk Dj (αBk ) − 4παJ i + αKE i + £β E i∂t B i = −ǫijk Dj (αEk ) + αKB i + £β B i
CE := DiE i − 4πρe
CB := DiB i• CP of Maxwell system in curved space-time∂t CE = αKCE + β j Dj CE∂t CB = αKCB + β j Dj CB• CP of ADM+Maxwell







t




∂

• CP of ADM+Maxwell+Hydroin progress.





CE 
∗



CB 

∗
 = 


0
H 




Mi
0

∗
∗
0
0

0
0
∗
∗





0   CE 




0
  CB 






∗ H Mi∗