Introduction to Numerical Relativity
Hisa-aki Shinkai
hshinkai@postman.riken.go.jp
Computational Science Division,
The Institute of Physical and Chemical Research (RIKEN),
Wako, Saitama, 351-0198 Japan
January 14, 2003
April 10, 2003 updated

Contents
1 Subjects of and for Numerical Relativity
1.1 Why Numerical Relativity? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Overview of Numerical Relativity . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Gravitational Wave Physics (Why Blackholes/Neutron Stars?) . . . . . . . . . .

3
3
5
7

2 The Standard Cauchy Approach to Numerical Relativity
2.1 The ADM formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 How to construct initial data 1: conformal approach . . . . . . . . . . . . . . . .
2.3 How to construct initial data 2: thin sandwich approach . . . . . . . . . . . . . .
2.4 How to choose gauge conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5 How to evolve the system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9
9
12
16
20
23

3 Alternative Approaches to Numerical Relativity
35
3.1 Full numerical, but diﬀerent foliations . . . . . . . . . . . . . . . . . . . . . . . . 35
3.2 Several approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4 Unsolved problems
4.1 Gravitational Wave Physics and related problems . . . . . . . . . . . . . . . . . .
4.2 Conjecture Hunting – reported and unreported issues – . . . . . . . . . . . . . .
4.3 Concluding Remark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42
42
42
43

Lecture note for APCTP Winter School on Gravitation and Cosmology, Korea, January 17-18, 2003.
This ﬁle and viewgraphs of the lecture are available at http://atlas.riken.go.jp/˜shinkai/winterAPCTP/1

Lecture Note at APCTP winter school 2003: H.Shinkai

2

Notations:• signature (− + ++).
• Covariant derivatives, Christoﬀel symbol∇µ Aα ≡ Aα;µ ≡ Aα,µ + Γαµν Aν
∇ µ Aα ≡Γαµν=Aα;µ ≡ Aα,µ − Γναµ Aν
(1/2)g αβ (gβµ,ν + gβν,µ − gµν,β )

(0.1)
(0.2)
(0.3)

• Riemann tensor, Ricci tensor, Weyl tensor
Rabcd ≡ ∂c Γabd − ∂d Γabc + Γaec Γebd − Γaed Γebc
Rab ≡ Rµaµb

Γµab,µ − Γµaµ,b + Γµνµ Γνab − Γµνb Γνaµ≡1
Cabcd = Rabcd − ga[c Rd]b + gb[c Rd]a − Rga[c gd]b ,
3

(0.4)
(0.5)
(0.6)

• ADM decomposition, the extrinsic curvature (§2)
ds2 = gµν dxµ dxν ,on Σ(t)...

2

d

i

j

= γij dx dx ,(µ, ν = 0, 1, 2, 3)(i, j = 1, 2, 3)

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

(0.7)

1
Kij ≡ −⊥µi ⊥νj nµ;ν = − £n γij .
2

(0.8)

Lecture Note at APCTP winter school 2003: H.Shinkai

1

Subjects of and for Numerical Relativity

1.1

Why Numerical Relativity?

3

The Einstein equation:1Rµν + gµν R + Λgµν = 8πGTµν
2

(1.9)

What are the diﬃculties? (# 1)• for 10-component metric, highly nonlinear partial diﬀerential equations.
• completely free to choose coordinates, gauge conditions, and even for decomposition of thespace-time.
• mixed with 4 elliptic eqs and 6 dynamical eqs if we apply 3+1 decomposition.
• has singularity in its nature.
How to solve it?
• ﬁnd exact solutions– assume symmetry in space-time, and decomposition of space-timespherically symmetric, cylindrical symmetric, ...
– assume simple situation and mattertime-dependency, homogeneity, algebraic speciality, ...
We know many exact solutions (O(100)) by this ”Spherical Cow” approach.
• approximations– weak-ﬁeld limit, linearization, perturbation, ...
We know correct prediction in the solar-system, binary neutron stars, ...
We know post-Newtonian behavior, ﬁrst-order correction, BH stability, ...
Why don’t we solve it using computers?
• dynamical behavior• strong gravitational ﬁeld
• no symmetry in space• gravitational wave!
• higher-dimensional theories, and/or other gravitational theories, ...
The most robust way to study the strong gravitational ﬁeld. Great.

Lecture Note at APCTP winter school 2003: H.Shinkai

Numerical Relativity= Solve the Einstein equations numerically.
= Necessary for unveiling the nature of strong gravity. For example:4

Box 1.1• gravitational waves from colliding black holes, neutron stars, supernovae, ...
• relativistic phenomena like cosmology, active galactic nuclei, ...
• mathematical feedback to singularity, exact solutions, chaotic behavior, ...
• laboratory for gravitational theories, higher-dimensional models, ...

What are the diﬃculties? (# 2)• How to construct a realistic initial data?
• How to treat black-hole singularity?
• We cannot evolve the system stably in long-term evolution. Why?
More general and recent introductions to numerical relativity are available, e.g. by d’Inverno(1996) [3], Seidel (1996/98/99) [5], Brügmann (2000) [2], Lehner (2001) [4], van Putten (2001) [6], andBaumgarte-Shapiro (2002) [1].

References[1] T.W. Baumgarte and S.L. Shapiro, gr-qc/0211028.
[2] B. Brügmann, Annalen Phys. 9, 227 (2000), available as gr-qc/9912009.
[3] R. d’Inverno, in General Relativity, ed. by G.S. Hall and J.R. Pulham (Institute of Physics
Publishing, 1996)[4] L. Lehner, Class. Quant. Grav. 18, R25 (2001); Proceedings of GR16, available as gr-qc/0202055.
[5] E. Seidel, in Relativity and Scientiﬁc Computing, ed. by by G.S. Hall, R.A. Puntigam, and H.
Ruder (Springer-Verlag, 1996);
E. Seidel, in Proceedings of GR15 conference, eds. by N. Dadhich and J. Narlikar (IUCA, Pune,
1998), available as gr-qc/9806088.
E. Seidel and W-M. Suen, gr-qc/9904014.
[6] M.H.P.M. van Putten, Proceedings of Asian Paciﬁc CTP Winter School on black hole astrophysics,
Pohang, Korea, available as gr-qc/0203076.

Lecture Note at APCTP winter school 2003: H.Shinkai

1.2

5

Overview of Numerical Relativity

Several milestones of NR
New proposals, developments, physical results.
1960s
1970s

1980s

1990

1995

2000

Hahn-Lindquist
May-WhiteÓMurchadha-York
Smarr
Smarr-Cades-DeWitt-Eppley
Smarr-York
ed. by L.Smarr
Nakamura-Maeda-Miyama-Sasaki
MiyamaBardeen-Piran
Stark-Piran
Shapiro-Teukolsky
Oohara-Nakamura
Seidel-Suen
Choptuik
NCSA group
Cook et al
Shibata-Nakao-Nakamura
Price-Pullin
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
ShibataYork
Brodbeck et al
Kidder-Finn
Shinkai-Yoneda
AEI group
AEI group
Shibata-Uryu
Shinkai-Yoneda
Meudon group
PennState group

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 ﬁnder
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
conformal thin-sandwich formulation
λ-system
BH, Spectral methods
planar GW, Ashtekar variables
full numerical to close limit
2 BH grazing collision
2 NS inspiral coalesence
adjusted ADM systems
irrotational BH binary initial dataisolated horizon

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
PRL82(1999)1350
JMathPhys40(1999)909
PRD62(2000)084026
CQG17(2000)4729
CQG17(2000)L149
PRL87(2001)271103
PTP107(2002)265
CQG19(2002)1027
PRD65(2002)044020
gr-qc/0206008

Lecture Note at APCTP winter school 2003: H.Shinkai

6

Issues to consider

Numerical Relativity – open issues
0. How to foliate space-time
Cauchy (3 + 1),
Hyperboloidal (3 + 1),
characteristic (2 + 2),
or combined?

Box 1.2⇒ see e.g. [24]⇒ see e.g. [6]⇒ if the foliation is (3 + 1), then · · ·1. How to prepare the initial data

⇒ see e.g. [3]Theoretical: Proper formulation for solving constraints?
How to prepare realistic initial data?
Eﬀects of background gravitational waves?
Connection to the post-Newtonian approximation?
Numerical: Techniques for solving coupled elliptic equations?
Appropriate boundary conditions?
2. How to evolve the data
Theoretical: Free evolution or constrained evolution?
Proper formulation for the evolution equations?
Suitable slicing conditions (gauge conditions)?
Numerical: Techniques for solving the evolution equations?
Appropriate boundary treatments?
Singularity excision techniques?
Matter and shock surface treatments?
Parallelization of the code?

⇒ see e.g. [40]3. How to extract the physical information
Theoretical: Gravitational wave extraction?
Connection to other approximations?
Numerical: Identiﬁcation of black hole horizons?
Visualization of simulations?

References[1] S. Husa, in the Proceedings of the conference “The Conformal Structure of Spacetimes: Geometry,
Analysis, Numerics”, ed. by J. Frauendiener and H. Friedrich, by Springer Verlag, Lecture Notes
in Physics series), available as gr-qc/0204043;
in the Proceedings of the 2001 Spanish Relativity meeting, eds. L. Fernandez and L. Gonzalez, tobe published by Springer, Lecture Notes in Physics series), available as gr-qc/0204057.
[2] J. Winicour, Livng Rev. Relativ. 2001-3 at http://www.livingreviews.org/[3] G. Cook, Livng Rev. Relativ. 2000-5 at http://www.livingreviews.org/[4] H. Shinkai and G. Yoneda, in Progress in Astronomy and Astrophysics (Nova Science Publ) to be
published, available as gr-qc/0209111.

Lecture Note at APCTP winter school 2003: H.Shinkai

1.3
1.3.1

7

Gravitational Wave Physics (Why Blackholes/Neutron Stars?)
General References

A resource guide by Centrella [1] might be quite useful. Essential references are by e.g. Thorne(1987/1997) [2], and Abramovici et al [3]. The latest review is by Cutler and Thorne [4]. Viewgraphsof the lecture (2002) by Dr. A. J. Weinstein (http://www.ligo.caltech.edu/˜ajw/) may be also
useful.

References[1] J. Centrella, Resource Letter GrW-1: Gravitational Waves, gr-qc/0211084.
[2] K. S. Thorne, in 300 Years of Gravitation, ed. by S. Hawking and W. Israel (Cambridge, 1987).
K. S. Thorne, gr-qc/9706079.
[3] A. Abramovici et al, Science 256, 325 (1992).
[4] C. Cutler and K. S. Thorne, in Proceedings of GR16 (Durban, South Africa, 2001), available asgr-qc/0204090.
1.3.2

Laser Interferometers

Current Projects, under operations! (target frequency: 10 − 103 Hz)
TAMA
LIGO
LIGO
GEO
VIRGO

Japan
USA
USA
Germany/UK
France/Italy

300m @ Tokyo
4Km @ Hanford, WA
4Km @ Livingston, LA
600m @ Hannover
3Km @ Pisa

19972001200120012002-

http://tamago.mtk.nao.ac.jp/
http://www.ligo-wa.caltech.edu/
http://www.ligo-la.caltech.edu/
http://www.geo600.uni-hannover.de/
http://www.pi.infn.it/virgo/virgoHome.html

Future Planning
CLIO
LCGT
LIGO 2LIGO 2
ACIGA

JAPAN
JAPAN
USA
USA
Australia

100m @ Kamioka
3Km @ Kamioka
4Km @ Hanford, WA
4Km @ Livingston, LA
3Km @

2004?
2006?
2006?
2006?
2003?

http://www.ligo.caltech.edu/
http://www.anu.edu.au/Physics/ACIGA/

Future Planning as Space Satellites (target frequency: 10−4 − 10−1 Hz)
LISADECIGO

US/ESAJAPAN

5 106 Km5 105 Km2011?
??

http://lisa.jpl.nasa.gov/(not approved yet)

Lecture Note at APCTP winter school 2003: H.Shinkai
1.3.3

8

What Information can we extract from gravitational waveform of binary neutronstar coalescence?

1. INSPIRAL phase [∼ 3 mins.] (Post-Newtonian Approx.)
df
dt

⇒

’chirp mass’ Mc ≡amplitude (incleasing)⇒

Mc , distance

amplitude (h+ /h× )⇒inclination

waveform⇒

eccentricity

moduration⇒

spin, · · ·’chirps’(M1 M2 )3/5
(M1 + M2 )1/5

2. ISCO phase (Post-Newtonian & Numerical Relativity)
frequency⇒⇒

Mass-Radius relation
Equation of States

3. COALESCE phase [∼ 3 millisecs.] (Numerical Relativity)
waveform · · ·other elements?

⇒
⇒
⇒
⇒
⇒

?
BH parameters (m, a), · · ·GR testγ-ray burst ?
r-process elements?

4. BLACK HOLE formation phase [∼ 10 msecs.] (Perturbation)
Quasi-Normal Modes⇒

Black Hole formation

with optical identiﬁcation
statics⇒⇒

Hubble parametercosmological parameters

5. STATISTICS

1.3.4

Requirements for Numerical Relativity• Where to start the simulation? How to construct physically reasonable initial data?
• How can we evolve the system stably?
• How to treat black hole singularity if it appears?
• How to extract gravitational wave?
• How can we manage the large-scale simulations?

Lecture Note at APCTP winter school 2003: H.Shinkai

92

The Standard Cauchy Approach to Numerical Relativity

2.1

The ADM formulation

2.1.1

The 3+1 decomposition of space-time

The idea of space-time evolution was ﬁrst formulated by Arnowitt, Deser, and Misner (ADM) [9]. The
formulation was ﬁrst motivated by a desire to construct a canonical framework in general relativity,
but it also gave the community to the fundamental idea of time evolution of space and time: suchas foliations of 3-dimensional hypersurface (Figure 1). This scheme is often called ‘3+1 formulation’,‘the ADM formulation’, or ‘Cauchy approach’.
Let us denote the hypersurface Σ(t) which is the three-dimensional spatial space with a parametert. The evolution of spacetime is expressed as the dynamics of Σ(t). The formulation begins bydecomposing the metic as
ds2 = gµν dxµ dxν ,on Σ(t)...

2i(µ, ν = 0, 1, 2, 3)j= γij dx dx ,d(i, j = 1, 2, 3)

Let the unit normal vector of the slices be nµ , wherenµ = g µν nν = (1/α, −β i /α).

nµ = (−α, 0, 0, 0),
We then have a 3+1 decomposed metric as

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

2

l

i

i

= (−α + βl β )dt + 2βi dtdx + γij dx dx


gµν =

−α2 + βl β l βj
βiγij

where α and βj are deﬁned as





,

g

µν



α ≡ 1/ −g 00 ,

=

(2.1)
j

−1/α2
β j /α2β i /α2 γ ij − β i β j /α2βj ≡ g0j .



(2.2)

and called the lapse function and shift vector, respectively.
shift vector, βisurface normal line

lapse function, αβi dtcoordinate constant lineαt = constant hypersurface

Figure 1: Concept of time evolution of space-time: foliations of 3-dimensional hypersurface. The lapse
and shift functions are often denoted α or N , and β i or N i , respectively.

Lecture Note at APCTP winter school 2003: H.Shinkai
2.1.2

10

The Standard ADM formulation

In order to decompose the Einstein equation into 3+1, we introduce the projection operator ⊥µν normal
to nµ ,
γµν = gµν + nµ nν ,γνµ = δνµ + nµ nν ≡ ⊥µν .
(2.3)
We also call the spatial components of γij the intrinsic 3-metric gij .1
The projections of the Einstein equation can be the following three:Gµν nµ nν= 8πG Tµν nµ nν ≡ 8πρH

(2.4)

Gµν n ⊥νiGµν ⊥µi ⊥νj

8πG Tµν n ⊥νi ≡ −8πJi
8πG Tµν ⊥µi ⊥νj ≡ 8πSij

(2.5)

µ

µ

=
=

(2.6)

To express these equation, we introduce the extrinsic curvature Kij as
Kij ≡ −⊥µi ⊥νj nµ;ν = · · · =
1 1−∂t γij + βi|j + βj|i = − £n γij .
2α
2

(2.7)

Projection of the Einstein equation on to the 3-hypersurface Σ is given using the Gauss-Codacci
relation: The Gauss equation,(3)Rαβγδ = (4) Rµνρσ ⊥µα ⊥βν ⊥γρ ⊥δσ − K αγ Kβδ + K αδ Kβγ ,(2.8)

and the Codacci equation,
Dj Ki j − Di K = −(4) Rρσ nσ ⊥i ρ ,(2.9)

where K = K i i , and Dµ is the covariant diﬀerentiation with respect to γij .
The projections (2.4)-(2.6) can be derived as follows.

The Standard ADM formulation [42, 51]:
Box 2.1
The fundamental dynamical variables are (γij , Kij ), the three-metric and extrinsic curvature.
The three-hypersurface Σ is foliated with gauge functions, (α, β i ), the lapse and shift vector.
• The evolution equations:∂t γij= −2αKij + Di βj + Dj βi ,∂t Kij= α (3)Rij + αKKij − 2αKik K k j − Di Dj α+(Di β k )Kkj + (Dj β k )Kki + β k Dk Kij(2.10)(2.11)

where K = K i i , and (3)Rij and Di denote three-dimensional Ricci curvature, and acovariant derivative on the three-surface, respectively.
• Constraint equations:
HADM

:=

MADM
i

:= Dj K i − Di K ≈ 0,

R + K 2 − Kij K ij ≈ 0,(3)

j(2.12)(2.13)

where (3)R =(3) Ri i : these are called the Hamiltonian (or energy) and momentum constraint equations, respectively.

1

If nµ is space-like, then γµν = gµν − nµ nνLecture Note at APCTP winter school 2003: H.Shinkai

11

The formulation has 12 free ﬁrst-order dynamical variables (γij , Kij ), with 4 freedom of gauge choice(α, βi ) and with 4 constraint equations, (2.12) and (2.13). The rest freedom expresses 2 modes ofgravitational waves.
We should remark here the ‘original’ ADM formulation. The evolution equations in Box 2.1 isthe version by Smarr and York which is now the standard convention for numerical relativists. They
adapted Kij as a fundamental variable instead of the conjugate momentum π ij , which was in theoriginal Arnowitt-Deser-Misner’s canonical formulation. Note that there is one replacement in (2.11)
using (2.12) in the process of conversion from the original ADM to the standard ADM equations.
More detail description: The Hamiltonian density can be written as√√
HGR = π ij γ̇ij − L, where L = −gR = α γ[(3)R − K 2 + Kij K ij ],where π ij is the canonically conjugate momentum to γij ,π ij =∂L√= − γ(K ij − Kγ ij ),∂ γ̇ijomitting the boundary terms. The variation of HGR with respect to α and βi yields the constraints, and the
δHGR
ij
GR
dynamical equations are given by γ̇ij = δH
δπ ij and π̇ = − δhij .
∂t γij∂t π ijN= 2 √ (πij − (1/2)γij π) + 2D(i Nj) ,γ
NN√= − γN ((3)Rij − (1/2)(3)Rγ ij ) + (1/2) √ hij (πmn π mn − (1/2)π 2 ) − 2 √ (π in πn j − (1/2)ππ ij )
γ
γ
√
√
i j
ij m
−1/2 m ij
m(ij)+ γ(D D N − γ D Dm N ) + γDm (γN π ) − 2π Dm N

The ADM formulation is a kind of constrained system, like Maxwell equations.
constraints
evolution eqs.

Maxwell eqs.
div E = 4πρ
div B = 0∂t E = rot B − 4πj∂t B = −rot E

ADM Einstein eq.
Hamiltonian constraint (2.12)
Momentum constraints (2.13)∂t γij = · · ·(2.10)∂t Kij = · · ·(2.11)In order to see the constraints are conserved during the evolution or not, we have to check how theconstraints evolve. The constraint propagation equations, which are the time evolution equations of
the Hamiltonian constraint (2.12) and the momentum constraints (2.13), can be written as [22, 39]

The Constraint Propagations of the Standard ADM:

Box 2.2∂t H = β j (∂j H) + 2αKH − 2αγ ij (∂i Mj )+α(∂l γmk )(2γ ml γ kj − γ mk γ lj )Mj − 4γ ij (∂j α)Mi ,(2.14)∂t Mi = −(1/2)α(∂i H) − (∂i α)H + β (∂j Mi )j

+αKMi − β k γ jl (∂i γlk )Mj + (∂i βk )γ kj Mj .

(2.15)

From these equations, we know that if the constraints are satisﬁed on the initial slice Σ, thenthe constraints are satisﬁed throughout evolution (in principle).

Lecture Note at APCTP winter school 2003: H.Shinkai
2.1.3

12

Numerical Procedures

In numerical relativity, this free-evolution approach is also the standard. This is because solving theconstraints (non-linear elliptic equations) is numerically expensive, and because free evolution allows
us to monitor the accuracy of numerical evolution.
The normal numerical scheme (free evolution scheme):
1. preparation of the initial data
solve the elliptic constraints for preparing the initial data (γij , Kij ).
2. time evolution(a) specify the gauge conditions (slicing conditions) for the lapse α and shift βi .
(b) evolve (γij , Kij ) by using the evolution equations.
(c) monitor the accuracy of simulations by checking the constraints.
(d) extract physical quantities.
3. step back to 2 and repeat.

References[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] L. Smarr, J.W. York, Jr., Phys. Rev. D 17, 2529 (1978).
[3] J.W. York, Jr., in Sources of Gravitational Radiation, ed. by L.Smarr, (Cambridge, 1979).
[4] S. Frittelli, Phys. Rev. D 55, 5992 (1997).
[5] H. Shinkai and G. Yoneda, Class. Quant. Grav. 19, 1027 (2002).

2.2

How to construct initial data 1: conformal approach

Initial Data Construction Problem
Prepare all metric and matter components by solving the two constraints:

Box 2.3• The Hamiltonian constraint equation
R + (trK)2 − Kij K ij = 2κρ + 2Λ(3)

(2.16)• The momentum constraint equations
Dj (K ij − γ ij trK) = κJ i

(2.17)

Lecture Note at APCTP winter school 2003: H.Shinkai

13

Conformal Approach – York-ÓMurchadha (1974)
N.ÓMurchadha and J.W.York Jr., Phys. Rev. D 10, 428 (1974)
2.2.1

Formulation

The idea issolutionγij = ψ 4 γ̂ij

trial metric(2.18)

We introduce the decomposition of Kij ,


Kij ⇒

trK = γ ij Kij
trace part
Aij = Kij − 13 γij trK trace-free part(2.19)

Then, other conformal transformations as consistent with (2.18) are:
γij
A

= ψ

−10

ρ = ψ

−n

ij

γ ij = ψ −4 γ̂ ij ,= ψ 4 γ̂ij ,
ij

Â ,

Aij = ψ

ρ̂,

i

−2

Âij ,(2.20)(2.21)−10 ˆiJ ,(2.22)

ˆ K̂, trA = trˆ Â = 0.
trK = tr(2.23)

Γi jk = Γ̂i jk + 2ψ −1 (δ i j D̂k ψ + δ i k D̂j ψ − γ̂jk γ̂ im D̂m ψ),ˆ
R = ψ −4 R̂ − 8ψ −5 ∆ψ.

(2.24)J =ψ

and we suppose
From (2.20), we get(2.25)

ˆ = γ̂ jk D̂j D̂k and R̂ = R(γ̂), and also Dj Aij = ψ −10 D̂j Âij .
where ∆
We further decompose Âij to divergence-free (transverse-traceless, TT) part and longitudinal part:ij
Âij = Âij
T T + (l̂W ) ,(2.26)

where we supposeD̂j Âij
TT = 0andˆ ÂT T = 0.
tr(2.27)

and2(l̂W )ij = D̂i W j + D̂j W i − γ̂ ij D̂k W k .
3Using these terms, we can writeD̂j Âij

ˆ l W )i ,= D̂j (l̂W )ij ≡ (∆ˆ )i + 1 D̂i (D̂j W j ) + R̂i j W j .
= (∆W3(2.28)(2.29)

With above transformation, the two constraints, (2.16) and (2.17), can be expressed as follows.
• The Hamiltonian constraint equation
ˆ = R̂ψ − (Âij Âij )ψ −7 + [ 2 (trK)2 − 2Λ]ψ 5 − 16πGρ̂ψ 5−n8∆ψ
3(2.30)• The momentum constraint equations
ˆ i + 1 D̂i D̂k W k + R̂i k W k = 2 ψ 6 D̂i trK + 8πGJˆi∆W
33(2.31)

Lecture Note at APCTP winter school 2003: H.Shinkai
2.2.2

14

Summary

Conformal approach (York-ÓMurchadha, 1974)
Box 2.4i
One way to set up the metric and matter components (γij , Kij , ρ, J ) so as to satisfy the constraints (2.16) and (2.17) is as follows.
1. Specify metric components γ̂ij , trK, ÂTijT , and matter distribution ρ̂, Jˆ in the conformal
frame.
2. Solve the next equations for (ψ, W i )ˆ = R̂ψ − (Âij Âij )ψ −7 + [ 2 (trK)2 − 2Λ]ψ 5 − 16πGρ̂ψ 5−n8∆ψ31
ˆ i + D̂i D̂k W k + R̂i k W k = 2 ψ 6 D̂i trK + 8πGJˆi
∆W
3
3
where

2 ij
i
j
j
i
k
Âij = Âij
T T + D̂ W + D̂ W − γ̂ D̂k W .
3(2.30)(2.31)(2.32)

3. Apply the inverse conformal transformation and get the metric and matter components
γij , Kij , ρ, J i in the physical frame:γij= ψ 4 γ̂ij ,1
Kij = ψ −2 [ÂTijT + (l̂W )ij ] + ψ 4 γ̂ij trK,3
ρ = ψ −n ρ̂,J i = ψ −10 Jˆi(2.33)(2.34)(2.35)(2.36)

Comments• Using the idea of conformal rescaling, we have a way to ﬁx 12 components of (γij , Kij ) thatsatisfy 4 constraints.
• The Hamiltonian constraint, (2.30), is a non-linear elliptic equation for ψ, so that we have tosolve it by an iterative method.
• The momentum constraints, (2.31), are PDEs for W i and coupled with (2.30). If we assumetrK = 0, then two constraints are decoupled. Normally people assume trK = 0 (maximal slicing
condition) or (trK) =const. (constant mean curvature slicing) for this purpose.
• For simplicity, people assume the background metric γ̂ij is conformally ﬂat γ̂ij = δij . The
physical appropriateness of conformal ﬂatness is often debatable.
• Two freedom of ÂTijT corresponds to the one of gravitational wave. However, there have beenno systematic discussion how to specify them, except applying tensor harmonics in a linearizedsituation.

Lecture Note at APCTP winter school 2003: H.Shinkai
2.2.3

15

Numerical procedures – Several tips

Solving the Hamiltonian constraint

Two Methods:

1. Solve the non-linear equation (2.30) directly.
2. Solve the linearized equation ψ = ψ0 + δψ iteratively.
ˆ
8∆ψ= E ψ + F ψ −7 + G ψ 5 + H ψ −3 + I ψ −1
= [E − 7F ψ0−8 + 5Gψ04 − 3Hψ0−4 − 2Iψ0−2 ] ψ + [8F ψ0−7 − 4Gψ05 + 4Hψ0−3 + 2Iψ0−1 ]Under an appropriate boundary condition, such as Robin BC ψ = 1 + const./r, or Dirichlet BC
ψ = 1 + Mtotal /2r.
Solve the momentum constraints

A couple of methods:

1. Solve the non-linear equations (2.31) directly.
2. Bowen’s method for conformally ﬂat case [GRG14(1982)1183]
Under the (∇i K = 0) condition, (2.31) becomes1∆W i + ∇i ∇j W j = 8πS i .
3
By introducing a decomposition of W i into vector and gradient terms
1
W i = V i − ∇i θ,
4the equations to solve are:∆V i = 8πS i ,(2.37)∆θ = ∇i V ,(2.38)i

If the source is of ﬁnite extent, then the the asymptotic behavior of V i and θ are given by
V i = −2
θ = −
+

∞


1
Qij1 ···jl nj1 · · · njl l+1 ,r
l=0(2.39)∞∞12(l + 1)1
1 ···jlQ{ij1 ···jl−1 } ni nj1 · · · njl−1 l−1 +nj1 · · · njl l+1
Qkj
k
r
(2l
+
1)(2l
+
3)
r
l=1
l=0

∞

2l − 11
M {ij1 ···jl−1 } ni nj1 · · · njl−1 l+1
2l + 1rl=1where ni = xi r−1 in the Cartesian cordinate, the multipoles Q and M are deﬁned as
ij1 ···jlQ≡

M ij1 ···jl≡(2l − 1)!!
S i (r)x{j1 xj2 · · · xjl } dV,l!
(2l − 1)!!
r2 S i (r)x{j1 xj2 · · · xjl } dV,l!

and where brackets denote the completely symmetric trace-free partZ {ij1 ···jl } = Z (ij1 ···jl ) −lk(j ···j
Z 1 l−1 δ jl i)
2l + 1 k(2.40)

Lecture Note at APCTP winter school 2003: H.Shinkai

2.3

16

How to construct initial data 2: thin sandwich approach
J.W.York Jr., Phys. Rev. Lett. 82, 1350 (1999)

The name “sandwich” comes from the proposal that this method prepares two spatial slices at
t = 0 and t = ∆t. There may be the following beneﬁts:• The input function is more friendly (3-metric and its time derivative) than the previous conformal
approach.
• The input quantity also requires the lapse function, N . (Actually this is the inverse and densitized lapse function. See below.)• The similar conformal transformation is applied. But the relation Āij = ψ −10 Aij is derived inthis version.
However, the numerical solvability is still debatable. Partial applications are seen in constructing
quasi-equilibrium binary neutron stars/black-holes. Matter terms are inserted by H. Shinkai.
2.3.1

metric, conformal metric, weighted conformal metricFirst, I list three types of 3-metric along to the conformal transformation.
(a) The metric, g ij which satisﬁes the constraints. (That is, the solution to seek.)(b) The conformal metric gij , whereg ij = ψ 4 gij“this world”“that world”(c) The “weighted” (−2/3) conformal metricĝij = g −1/3 g ij = g −1/3 gij(2.41)

where g = det(g ij ) and g = (det gij ). For small variation, g ij δĝij = 0 is always hold, and we getg ij ∂t ĝij = g ij ∂t ĝij = ĝ ij ∂t ĝij = 0 .

(2.42)In York’s paper, he does not use the weighted conformal metric, but he imposes that the conformal
metric does have the property similar to (2.42).
2.3.2

introduction of velocity tensor

On the second slice t = δt, we write the conformal metric= gij + uij δt,
gij(2.43)

where we introduced the velocity tensor (suppose to be a given quantity)uij = ġij ,(2.44)

together with the “weighted” condition,g ij uij = 0,andg ij ġij = 0.

(2.45)

Lecture Note at APCTP winter school 2003: H.Shinkai

17

By taking the traceless part of the evolution equation,∂t g ij ≡ ġ ij = −2N K ij + (Di β j + Dj β i )(2.46)

together with (2.24) and (2.25), we get
1ġ ij − g ij g kl ġ kl ≡ uij = −2N Aij + (Lβ)ij(2.47)
3kwhere Aij ≡ K ij − (1/3)Kg ij , and (Lβ)ij ≡ Di β j + Dj β i − (2/3)g ij D β k . (2.48)
From (2.47), we obtain
uij = ψ 4 uij(2.49)

Similarly, we obtainiβ i = ψ 4 βi ,β = βi,
4ij(Lβ)ij = ψ (Lβ)ij ,2.3.3(2.50)(Lβ) = ψ−4ij(Lβ) .

(2.51)

Redeﬁnition of the lapse and its conformal transformation

We call the standard α(t, x) > 0 slicing function, and deﬁne the lapse function N as
N = g 1/2 α.

(2.52)The slicing function is now α = g −1/2 N = N, which may be called the inverse densitized lapse. Note∼that the lapse here, N , depends g, so that N is not a pure gauge quantity. Therefore we treat that
α is the freely speciﬁed function and let

α = α . In result, we obtain a new relation

N = ψ6Nfrom (2.52).
We also impose
A

ij

K=K

as before. The relation (2.47) then derives= ψ −6 (2N )−1 ψ −4 (Lβ)ij − ψ −4 uij= ψ −10 (2N )−1 (Lβ)ij − uij

2.3.4= ψ −10 Aijthat isij

A = ψ −10 Aij

Constraints to solve

By using above boxed conformal transformations, two constraints can be transformed as 2
• The Hamiltonian constraint equation (the same with before)28∆g ψ − R(g)ψ + Aij Aij ψ −7 − [ K − 2Λ]ψ 5 − 16πGρψ 5−n = 0 ,3(2.53)• The momentum constraint equations
Dj (2N )−1 (Lβ)ij
22= Dj (2N )−1 uij + ψ 6 Di K + 8πGJ i
3

LHS of (2.54) isDa (2N )−1 (Di β a + Da β i −1 ia k1
g D βk ) + (2N )−1 ∆β i + Di Dk β k + Ri k β k33(2.54)

Lecture Note at APCTP winter school 2003: H.Shinkai
2.3.5

18

SummaryThin-Sandwich approach (York, 1999)
Box 2.5ii
One way to set up the metric, gauge values and matter components (g ij , K ij , N , β , ρ, J ) so as
to satisfy the constraints (2.16) and (2.17) is as follows.
1. Specify metric components gij , uij (= ġij ), K, the lapse function N , and matter distribution
ρ, J i in the conformal frame.
2. Solve the next equations for (ψ, β i )28∆g ψ − R(g)ψ + Aij Aij ψ −7 − [ K − 2Λ]ψ 5 − 16πGρψ 5−n = 0 ,
3
2
−1
ij
−1 ij
Dj (2N ) (Lβ) = Dj (2N ) u + ψ 6 Di K + 8πGJ i ,3
where

Aij = (2N )−1 (Lβ)ij − uij .

(2.53)(2.54)(2.55)

3. Apply the inverse conformal transformation and get the metric and matter componentsii(γ ij , K ij , N , β , ρ, J ) in the physical frame:= ψ 6 N,

(2.56)i= βi,(2.57)

g ij

4(2.58)Nβ= ψ gij ,1= ψ −2 Aij + ψ 4 gij K,3
ρ = ψ −8 ρ,

K ij

J

i

= ψ

−10 i

J .

(2.59)(2.60)(2.61)

Comments• The two equations, (2.53) and (2.54), are coupled, but they will be decoupled if we assume theconstant mean curvature condition, (trK) =const. (This is the same as the conformal approach,but we have to solve the momentum constraints ﬁrst here.)• The (general) solvability of (2.54) is still debatable.

Lecture Note at APCTP winter school 2003: H.Shinkai

19Comparison between two approaches

input
functions
treatment of
gauge
functionsusage of theconstraints

counting thefreedom

conformal approach
gij , K, ATijT (components: 6, 1, 2)
GW components are separated out
lapse and shift are not appearing in
the formulation.

Hamiltonian constraint is for the conformal factor ψmomentum constraints are for thelongitudinal part of Aij .
(input 9 functions) plus (3 functions
by solving momentum constraints)= 12 = (3-metric) plus (extrinsic
curvature).

thin-sandwich approach
gij , K, uij , N (comp.: 6, 1, 5, 1)
can specify time-derivatives
lapse is given by the conformal transformation.
shift is given by solving theconstraints.
Hamiltonian constraint is for the conformal factor ψmomentum constraints are for shiftfunction β i .
(input 13 functions) plus (3 functions
by solving momentum constraints) =16 = (3-metric) plus (extrinsic curvature) plus (gauge functions).

Interpretation of quasi-equilibrium data constructioncf) G. Cook, Living Reviews in Relativity, 2000-5
The thin-sandwich formulation oﬀers the easiest way to interpret a method constructing the initial
data for binary NSs or BHs. For example, if we assume• the binary has a constant orbital angular velocity of Ω.
• the existence of the helical Killing vector, ξ µ = tµ + Ωφµ . (ξ µ is also assumed to the time vector
in the rotating frame that we are on, while tµ is in the rest frame).
• align shift vector with Killing vector⇒kψ 4 (Lβ)ij = (Lβ)ij ≡ Di β j + Dj β i − (2/3)g ij D β k = 0• gravitational radiation is negligible, i.e. conformally ﬂat background, R = 0.
• maximally slicing condition, K = ∂t K = 0.
Thin-sandwich formulation gives us following explicit relations:

1. The momentum constraints, Dj (2N )−1 (Lβ)ij = Dj (2N )−1 uij + 23 ψ 6 Di K + 8πGJ i ,⇒uij = 02. The Hamiltonian constraint, 8∆g ψ − R(g)ψ + Aij Aij ψ −7 − [ 23 K − 2Λ]ψ 5 − 16πGρψ 5−n = 0 ,⇒∆(αψ 7 ) = (αψ 7 )7
ij −8 + 2πGψ 4 (ρ + 2S)
8 Aij A ψ

Lecture Note at APCTP winter school 2003: H.Shinkai

2.4

20

How to choose gauge conditions

The standard 3 + 1 formulation allows us to choose gauge conditions (slicing conditions) for every time
step. The fundamental guidelines for ﬁxing the lapse function α and the shift vector βi :• to avoid the foliation hitting the physical and coordinate singularity in its evolution.
• to make system suitable for physical situation.
• to make the evolution system as simple as possible.
• to enable the gravitational wave extraction easy.
I list several essential slicing conditions below. The notations hereafter follows those of §2.1 (ADM
formulation).
2.4.1

Lapse conditions

geodesic slice

α=1

harmonic slice∇a ∇a xb = 0

maximal slice

K=0

maximal slice(K-driver)constant mean
curvaturepolar slicing∂t K = −c2 K

algebraic

K = const.
Kθθ + Kϕϕ = 0, or
K = Krr√α ∼ γ,α ∼ 1 + log γ

GOOD
BAD
GOOD
GOOD
BAD
GOOD
BAD
G&B
GOOD
G&B
GOOD
GOOD
BAD
GOOD
BAD

simple, easy to understandno singularity avoidance
simplify eqs.,
easy to compare analytical investigations
no singularity avoidance or coordinate
pathologies
singularity avoidance
have to solve an elliptic eq.
same with maximal slice,
easy to maintain K = 0
same with maximal slice,
suitable for cosmological situation
singularity avoidance in isotropic coord.
trouble in Schwarzschild coord.
easy to implement
not avoiding singularity[2]-[7][1],[8][15][12][16]-[18][19]-[21]Maximal slicing This is always the ﬁrst one to be mentioned as a singularity avoiding
gauge
√ 3
condition. The name of ‘maximal’ comes from the fact that the deviation of the 3-volume V =γd xalong to the normal line becomes maximal when we set K = 0. This is simply written as
K=0

on Σ(t).

(2.62)

Pioneering idea can be seen in Lichnerowicz [8], and it was extended by York [1]. This conditionis supposed to be applied in simulations that a singularity will appear during evolutions such as
gravitational collapses. The actual equation for determining the lapse function α can be obtained
from ∂t K = ∂t (Kij γ ij ) = 0. By substituting the evolution equations, we get
Di Di α = { (3)R + K 2 + 4πG(S − 3ρH ) − 3Λ}α,(2.63)

or by using the Hamiltonian constraint further,
Di Di α = {Kij K ij + 4πG(S + ρH ) − Λ}α.

(2.64)

Lecture Note at APCTP winter school 2003: H.Shinkai

21This is an elliptic equation. When the curvature is strong (i.e. close to the appearance of a singularity),
the RHS of equation become larger, hence the lapse becomes smaller. Therefore the foliation near the
singularity evolves slowly.
For Schwarzschild black-hole space-time, Estabrook et al. [10] showed that the maximal slicing
condition allows the 3-surface to reach into r = 1.5M in the limit t → ∞, that is inside of theevent horizon, r = 2M . However, it is also reported that the diﬀerence of α-evolution causes thegrid-stretching problem.
2.4.2

Shift conditions

geodesic slice

βi = 0minimaldistortion
minimal strainmin Σij Σij
min Θij Θij

simple, easy to understandtoo simple
geometrical meaning
elliptic eqs., hard to solve
same with minimal distortion

GOOD
BAD
GOOD
BAD
G&B

[1]
[1]

Minimal distortion condition, minimal strain condition Any singularity avoiding slice conditions causes the grid stretching problem. Smarr and York [1] proposed the condition which minimize
the distortion in a global sense.
Let us deﬁne the expansion tensor Θµν and the distortion tensor Σij . Let the normal direction
to the surface nµ , and the coordinate-constant congruence tµ = αnµ + βµ . By projecting tµ onto thehypersurface using the projection operator ⊥ab = δba + na nb ,1Θµν = ⊥∇(ν tµ) = −αKµν + D(µ βν)2(2.65)We then extract this traceless part and deﬁne,


Σij







1
1
1
1
= Θij − Θγij = −2α Kij − γij K +D(i βj) − Dk βk .
3323(2.66)The minimal distortion condition is to choose β i which minimize the action1δS[β] = δ{2Σij Σij d3 x} = 0.

(2.67)

This condition can be written as Dj Σij = 0, or




2
1
Dj Dj βi + Dj Di βj − Di Dj β j = Dj 2α Kij − trKγij
3
3
or





1
1
∆βi + Di (Dj βj ) + Rij βj = Dj 2α Kij − trKγij
3
3



,

(2.68)

,(2.69)

where ∆ = Di Di .
Similarly, we can deﬁne the minimal strain condition by minimizing Θij Θij .
The both requires non-linear elliptic equations and hard to solve. Several group solves “pseudo”minimal distortion condition by replacing the covariant derivatives to the partial derivatives [22]. This
simpliﬁcation also works for inspiral binary neutron star evolution.

Lecture Note at APCTP winter school 2003: H.Shinkai

22

References[1] J.W. York,Jr., “Kinematics and Dynamics of General Relativity”, in Sources of Gravitational Radiation,
ed. by L.Smarr, (Cambridge, 1979) ;L.Smarr and J.W.York,Jr., Phys. Rev. D 17, 2529 (1978)[2] C. Bona and J. Massó, Phys. Rev. D 38, 3419 (1988).
[3] C. Bona and J. Massó, Phys. Rev. Lett. 68, 1097 (1992).
[4] G.B. Cook and M.A. Scheel, Phys. Rev. D 56, 4775 (1997).
[5] M. Alcubierre, Phys. Rev. D 55, 5981 (1997).
[6] M. Alcubierre and J. Massó, Phys. Rev. D 57, R4511 (1998).
[7] M. Alcubierre, gr-qc/0210050.
[8] A. Lichnerowicz, J. Math. Pures Appl. 23, 37 (1944).
[9] D. Bernstein, D.W. Hobill, and L. Smarr, in Frontiers in Numerical Relativity, edited by C.R. Evans, L.S.
Finn, and D.W. Hobill (Cambridge Univ. Press, 1989)[10] F. Estabrook, H. Wahlquist, S. Christensen, B. DeWitt, L. Smarr, and E. Tsiang, Phys. Rev. D. 7, 2814
(1973).
[11] M.J. Duncan, Phys. Rev. D 31, 1267 (1985).
[12] J. Balakrishna, G. Daues, E. Seidel, W-M. Suen, M. Tobias and E. Wang, Class. Quant. Grav. 13, L135
(1996).
[13] A. Geyer and H. Herold , Phys. Rev. D. 52, 6182 (1995).
[14] A. Geyer and H. Herold , Gen. Rel. Grav. 29, 1257 (1997).
[15] L. I. Petrich, S.L. Shapiro, and S.A. Teukolsky, Phys. Rev. D 31, 2459 (1985).
[16] D.M. Eardley, L. Smarr, Phys. Rev. D 19, 2239 (1979)[17] D.R. Brill, J.M. Cavallo, J.A. Isenberg, J. Math. Phys. 21, 2789 (1980)[18] K. Nakao, T. Nakamura, K. Oohara, K. Maeda, Phys. Rev. D 44, 1326 (1991)[19] J.M. Bardeen and T. Piran, Phys. Rep. 96, 205 (1983).
[20] R. A. d’Inverno, Class. Quantum Grav. 12, L75 (1995).
[21] L. I. Petrich, S.L. Shapiro, and S.A. Teukolsky, Phys. Rev. D 33, 2100 (1986).
[22] M. Shibata and K. Uryu, Prog. Theor. Phys. 107, 265 (2002).

Lecture Note at APCTP winter school 2003: H.Shinkai

2.5

23

How to evolve the system

Many trials for longterm stable and accurate simulations of binary compact objects have revealed
that mathematically equivalent sets of evolution equations show diﬀerent numerical stability in freeevolution schemes. Thus, the stability problem or the formulation problem is now shedding light on
the mathematical structure of the Einstein equations. More detail review is available as [40].
2.5.1

Overview

Up to a couple of years ago, the “standard ADM” decomposition (§2.1) of the Einstein equation was
taken as the standard formulation for numerical relativists. However, numerical simulations were often
interrupted by unexplained blow-ups. This was thought due to the lack of resolution, or inappropriate
gauge choice, or the particular numerical scheme which was applied. However, after the accumulationof much experience, people have noticed the importance of the formulation of the evolution equations,since there are apparent diﬀerences in numerical stability although the equations are mathematicallyequivalent Figure 2 is a chronological map of the researches. See Column 1 for the meaning of“stability”.

90s

Nakamura-Oohara

2000s

Shibata

87

95

99

Nakamura-Oohara

Shibata-Nakamura

Baumgarte-Shapiro

62

G-codeH-codeADM

NCSA

BSSN-code

Illinois

AEI

97

92

PennState

Alcubierre

01

Bona-Masso
95-97

99

adjusted-system

80s

UWash

ChoquetBruhat-York Anderson-York
01

Cornell-Illinois

Kidder-Scheel-Teukolsky

96

CaltechFrittelli-Reula
Hern

99

LSU

lambda-system
Shinkai-Yoneda

97
86

Ashtekar

Iriondo-Leguizamon-Reula

99

Yoneda-Shinkai

Figure 2: Chronological table of formulations and their numerical tests. Boxed ones are of proposals of
formulation, circled ones are related numerical experiments. Please refer Table 1 in [40] for references.
At this moment, there are three major ways to obtain longer time evolutions: (1) modiﬁcationsof the standard Arnowitt-Deser-Misner equations initiated by the Kyoto group, (2) rewriting of theevolution equations in hyperbolic form, and (3) construction of an “asymptotically constrained” system. Of course, the ideas, procedures, and problems are mingled with each other. The purpose of thissection is to review all three approaches and to introduce our idea to view them in a uniﬁed way. The
third idea has been generalized by us as an asymptotically constrained system. The main procedure
is to adjust the evolution equations using the constraint equations [47, 48, 39]. The method is also
applied to explain why the above approach (1) works, and also to propose alternative systems basedon the ADM [48, 39] and BSSN [49] equations.

Lecture Note at APCTP winter school 2003: H.Shinkai

24

Column 1
The word stability is used quite diﬀerent ways in the community.
• We mean by numerical stability a numerical simulation which continues without any blowups and in which data remains on the constrained surface.
• Mathematical stability is deﬁned in terms of the well-posedness in the theory of partial
diﬀerential equations, such that the norm of the variables is bounded by the initial data.
See eq. (2.83) and around.
• For numerical treatments, there is also another notion of stability, the stability of ﬁnite
diﬀerencing schemes. This means that numerical errors (truncation, round-oﬀ, etc) are
not growing by evolution, and the evaluation is obtained by von Neumann’s analysis.
Lax’s equivalence theorem says that if a numerical scheme is consistent (converging tothe original equations in its continuum limit) and stable (no error growing), then thesimulation represents the right (converging) solution. See [18] for the Einstein equations.

.
2.5.2

Strategy 0: The ADM formulation

As we see in §2.1, we know that if the constraints are satisﬁed on the initial slice Σ, then the constraints
are satisﬁed throughout evolution. The normal numerical scheme is to solve the elliptic constraints for
preparing the initial data, and to apply the free evolution (solving only the evolution equations). The
constraints are used to monitor the accuracy of simulations.
The origin of the problem was that the above statement in Italics is true in principle, but is notalways true in numerical applications. A long history of trial and error began in the early 90s. Shinkai
and Yoneda showed that the standard ADM equations has a constraint violating mode in its constraint
propagation equations even for a single black-hole (Schwarzschild) spacetime [39].
2.5.3

Strategy 1: Modiﬁed ADM formulation by Nakamura et al

Up to now, the most widely used formulation for large scale numerical simulations is a modiﬁed ADM
system, which is now often cited as the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation.
This reformulation was ﬁrst introduced by Nakamura et al. [29, 30, 35]. The usefulness of thisreformulation was re-introduced by Baumgarte and Shapiro [11], then was conﬁrmed by other groups
to show a long-term stable numerical evolution [4, 6].
2.5.4

Basic variables and equations

The widely used notation[11] introduces the variables (ϕ, γ̃ij ,K,Ãij ,Γ̃i ) instead of (γij ,Kij ), whereϕ = (1/12) log(detγij ),−4ϕÃij = eγ̃ij = e−4ϕ γij ,(Kij − (1/3)γij K),Γ̃i

K = γ ij Kij ,= Γ̃ijk γ̃ jk .

(2.70)(2.71)The new variable Γ̃i was introduced in order to calculate Ricci curvature more accurately. In BSSN
ADM = ∂ Γk − ∂ Γk + Γl Γk − Γl Γk , but asformulation, Ricci curvature is not calculated as Rij
i kj
k ij
ij lk
kj li
ϕ
BSSN
= Rij + R̃ij , where the ﬁrst term includes the conformal factor ϕ while the second term does
Rij

Lecture Note at APCTP winter school 2003: H.Shinkai

25

BSSN does have wave operator apparently in the ﬂat
not. These are approximately equivalent, but Rijbackground limit, so that we can expect more natural wave propagation behavior.
Additionally, the BSSN requires us to impose the conformal factor as γ̃(:= detγ̃ij ) = 1, during
evolution. This is a kind of deﬁnition, but can also be treated as a constraint. 3

The BSSN formulation [29, 30, 35, 11]:
Box 2.3
The fundamental dynamical variables are (ϕ, γ̃ij ,K,Ãij ,Γ̃i ).
The three-hypersurface Σ is foliated with gauge functions, (α, β i ), the lapse and shift vector.
• The evolution equations:∂tB ϕ = −(1/6)αK + (1/6)β i (∂i ϕ) + (∂i β i ),∂tB γ̃ij∂tB K∂tB Ãij∂tB Γ̃i(2.72)= −2αÃij + γ̃ik (∂j β ) + γ̃jk (∂i β ) − (2/3)γ̃ij (∂k β ) + β (∂k γ̃ij ), (2.73)
k

k

k

k

= −Di Di α + αÃij Ãij + (1/3)αK 2 + β i (∂i K),
=
=

−4ϕ

(2.74)−4ϕ

BSSN T F−e(Di Dj α) + eα(Rij) + αK Ãij − 2αÃik Ãk j+(∂i β k )Ãkj + (∂j β k )Ãki − (2/3)(∂k β k )Ãij + β k (∂k Ãij ),−2(∂j α)Ãij + 2α(Γ̃ijk Ãkj − (2/3)γ̃ ij (∂j K) + 6Ãij (∂j ϕ))(2.75)−∂j (β k (∂k γ̃ ij ) − γ̃ kj (∂k β i ) − γ̃ ki (∂k β j ) + (2/3)γ̃ ij (∂k β k )).

(2.76)

TF• Constraint equations:
HBSSN
MBSSN
i
i
G

= RBSSN + K 2 − Kij K ij ,
=
=

A =

MADM
,
i
i
jk i
Γ̃ − γ̃ Γ̃jk ,Ãij γ̃ ij ,S = γ̃ − 1.

(2.77)(2.78)(2.79)(2.80)(2.81)(2.77) and (2.78) are the Hamiltonian and momentum constraints (the “kinematic” constraints), while
the latter three are “algebraic” constraints due to the requirements of BSSN formulation.
Remarks Why BSSN is better than the standard ADM? Together with numerical comparisons withthe standard ADM case[6], this question has been studied by many groups using diﬀerent approaches.
Using numerical test evolution, Alcubierre et al [4] found that the essential improvement is in theprocess of replacing terms by the momentum constraints. They also pointed out that the eigenvaluesof BSSN evolution equations have fewer “zero eigenvalues” than those of ADM, and they conjectured
that the instability might be caused by these “zero eigenvalues”. An eﬀort was made to understandthe advantage of BSSN from the point of hyperbolization of the equations in its linearized limit[4, 32]. These studies provide some support regarding the advantage of BSSN, while it is also shown
an example of an ill-posed solution in BSSN (as well in ADM) by Frittelli and Gomez [23].
As we discussed in [49], the stability of the BSSN formulation is due not only to the introductionsof new variables, but also to the replacement of terms in the evolution equations using the constraints.
Further, we will show several additional adjustments to the BSSN equations which are expected togive us more stable numerical simulations.
3

The box/column numbers in this subsection are numerated so as to ﬁt with [40].

Lecture Note at APCTP winter school 2003: H.Shinkai
2.5.5

26

Strategy 2: Hyperbolic reformulations

Deﬁnitions, properties, mathematical backgrounds The second eﬀort to re-formulate the Einstein equations is to make the evolution equations reveal a ﬁrst-order hyperbolic form explicitly. Thisis motivated by the expectation that the symmetric hyperbolic system has well-posed properties in
its Cauchy treatment in many systems and also that the boundary treatment can be improved if weknow the characteristic speed of the system.

Hyperbolic formulations
Box 2.5
We say that the system is a ﬁrst-order (quasi-linear) partial diﬀerential equation system, if acertain set of (complex-valued) variables uα (α = 1, · · · , n) forms∂t uα = Mlβ α (u) ∂l uβ + Nα (u),(2.82)

where M (the characteristic matrix) and N are functions of u but do not include any derivativesof u. Further we say the system is• a weakly hyperbolic system, if all the eigenvalues of the characteristic matrix are real.
• a strongly hyperbolic system (or a diagonalizable / symmetrizable hyperbolic system), if
the characteristic matrix is diagonalizable (has a complete set of eigenvectors) and has all
real eigenvalues.
• a symmetric hyperbolic system, if the characteristic matrix is a Hermitian matrix.

Writing the system in a hyperbolic form is a quite useful step in proving that the system is wellposed. The mathematical well-posedness of the system means (1◦ ) local existence (of at least one
solution u), (2◦ ) uniqueness (i.e., at most solutions), and (3◦ ) stability (or continuous dependence ofsolutions {u} on the Cauchy data) of the solutions. The resultant statement expresses the existenceof the energy inequality on its norm,||u(t)|| ≤ eατ ||u(t = 0)||,where 0 < τ < t,

α = const.

(2.83)

This indicates that the norm of u(t) is bounded by a certain function and the initial norm. Remark that
this mathematical boundness does not mean that the norm u(t) decreases along the time evolution.
The inclusion relation of the hyperbolicities is,symmetric hyperbolic ⊂ strongly hyperbolic ⊂ weakly hyperbolic.

(2.84)The Cauchy problem under weak hyperbolicity is not, in general, C ∞ well-posed. At the strongly
hyperbolic level, we can prove the ﬁniteness of the energy norm if the characteristic matrix is independent of u (cf [43]), that is one step deﬁnitely advanced over a weakly hyperbolic form. Similarly, thewell-posedness of the symmetric hyperbolic is guaranteed if the characteristic matrix is independentof u, while if it depends on u we have only limited proofs for the well-posedness.
From the point of numerical applications, to hyperbolize the evolution equations is quite attractive,
not only for its mathematically well-posed features. The expected additional advantages are thefollowing.
(a) It is well known that a certain ﬂux conservative hyperbolic system is taken as an essential formulation in the computational Newtonian hydrodynamics when we control shock wave formations
due to matter.

Lecture Note at APCTP winter school 2003: H.Shinkai

27(b) The characteristic speed (eigenvalues of the principal matrix) is supposed to be the propagation
speed of the information in that system. Therefore it is naturally imagined that these magnitudes
are equivalent to the physical information speed of the model to be simulated.
(c) The existence of the characteristic speed of the system is expected to give us an improved
treatment of the numerical boundary, and/or to give us a new well-deﬁned Cauchy problem
within a ﬁnite region (the so-called initial boundary value problem, IBVP).
These statements sound reasonable, but have not yet been generally conﬁrmed in actual numericalsimulations. But we are safe in saying that the formulations are not yet well developed to test theseissues.
Hyperbolic formulations of the Einstein equations Most physical systems can be expressed
as symmetric hyperbolic systems. In order to prove that the Einstein’s theory is a well-posed system,
to hyperbolize the Einstein equations is a long-standing research area in mathematical relativity.
The standard ADM system does not form a ﬁrst order hyperbolic system. This can be seen
immediately from the fact that the ADM evolution equation (2.11) has Ricci curvature in RHS. So
far, several ﬁrst order hyperbolic systems of the Einstein equation have been proposed. In constructing
hyperbolic systems, the essential procedures are (1◦ ) to introduce new variables, normally the spatiallyderivatived metric, (2◦ ) to adjust equations using constraints. Occasionally, (3◦ ) to restrict the gauge
conditions, and/or (4◦ ) to rescale some variables. Due to process (1◦ ), the number of fundamental
dynamical variables is always larger than that of ADM.
Due to the limitation of space, we can only list several hyperbolic systems of the Einstein equations.
• The Bona-Massó formulation [13, 14]• The Einstein-Ricci system [19, 1] / Einstein-Bianchi system [7]• The Einstein-Christoﬀel system [8]• The Ashtekar formulation [10]• The Frittelli-Reula formulation [24, 43]• The Conformal Field equations [21]• The Kidder-Scheel-Teukolsky (KST) formulation [26]Please refer [40] for each brief introductions.
Remarks When we discuss hyperbolic systems in the context of numerical stability, the following
questions should be considered:
Q From the point of the set of evolution equations, does hyperbolization actually contributeto numerical accuracy and stability? Under what conditions/situations will the advantagesof hyperbolic formulation be observed?
Unfortunately, we do not have conclusive answers to these questions, but many experiences are being
accumulated. Several earlier numerical comparisons reported the stability of hyperbolic formulations[14, 15, 33, 34]. But we have to remember that this statement went against the standard ADM
formulation, which has a constraint-violating mode for Schwarzschild spacetime as has been shown
recently[39].
These partial numerical successes encouraged the community to formulate various hyperbolic systems. Recently, Calabrese et al [17] reported there is a certain diﬀerences in the long-term convergencefeatures between weakly and strongly hyperbolic systems on the Minkowskii background space-time.
However, several numerical experiments also indicate that this direction is not a complete success.
Objections from numerical experiments• Above earlier numerical successes were also terminated with blow-ups.

Lecture Note at APCTP winter school 2003: H.Shinkai

28• If the gauge functions are evolved according to the hyperbolic equations, then their ﬁnite
propagation speeds may cause pathological shock formations in simulations [2, 3].
• There are no drastic diﬀerences in the evolution properties between hyperbolic systems(weakly, strongly and symmetric hyperbolicity) by systematic numerical studies by Hern[25] based on Frittelli-Reula formulation [24], and by the authors [38] based on Ashtekar’s
formulation [10, 46].
• Proposed symmetric hyperbolic systems were not always the best ones for numerical evolution. People are normally still required to reformulate them for suitable evolution. Sucheﬀorts are seen in the applications of the Einstein-Ricci system [34], the Einstein-Christoﬀel
system [12], and so on.
Of course, these statements only casted on a particular formulation, and therefore we have to becareful not to over-emphasize the results. In order to ﬁgure out the reasons for the above objections,
it is worth stating the following cautions:Remarks on hyperbolic formulations(a) Rigorous mathematical proofs of well-posedness of PDE are mostly for simple symmetric orstrongly hyperbolic systems. If the matrix components or coeﬃcients depend on dynamicalvariables (as in all any versions of hyperbolized Einstein equations), almost nothing wasproved in more general situations.
(b) The statement of “stability” in the discussion of well-posedness refers to the bounded
growth of the norm, and does not indicate a decay of the norm in time evolution.
(c) The discussion of hyperbolicity only uses the characteristic part of the evolution equations,and ignores the rest.
We think the origin of confusion in the community results from over-expectation on the above
issues. Mostly, point (c) is the biggest problem. The above numerical claims from Ashtekar andFrittelli-Reula formulations were mostly due to the contribution (or interposition) of non-principal
parts in evolution. Regarding this issue, the recent KST formulation ﬁnally opens the door. KST’s“kinematic” parameters enable us to reduce the non-principal part, so that numerical experiments are
hopefully expected to represent predicted evolution features from PDE theories. At this moment, theagreement between numerical behavior and theoretical prediction is not yet perfect but close [27].
If further studies reveal the direct correspondences between theories and numerical results, thenthe direction of hyperbolization will remain as the essential approach in numerical relativity, and therelated IBVP researches will become a main research subject in the future. Meanwhile, it will be
useful if we have an alternative procedure to predict stability including the eﬀects of the non-principal
parts of the equations. Our proposal of adjusted system in the next subsection may be one of them.
2.5.6

Strategy 3: Asymptotically constrained systems

The third strategy is to construct a robust system against the violation of constraints, such that theconstraint surface is an attractor. The idea was ﬁrst proposed as “λ-system” by Brodbeck et al [16],and then developed in more general situations as “adjusted system” by the authors [47].
The “λ-system” Brodbeck et al [16] proposed a system which has additional variables λ that obey
artiﬁcial dissipative equations. The variable λs are supposed to indicate the violation of constraints
and the target of the system is to get λ = 0 as its attractor.

Lecture Note at APCTP winter school 2003: H.Shinkai

29

The “λ-system” (Brodbeck-Frittelli-Hübner-Reula) [16]:Box 2.7
For a symmetric hyperbolic system, add additional variables λ and artiﬁcial force to reduce theviolation of constraints.
The procedure:
1.

Prepare a symmetric hyperbolic evolution system∂t u = M ∂i u + N

2.

Introduce λ as an indicator of violation of
constraint which obeys dissipative eqs. of motion∂t λ = αC − βλ(α = 0, β > 0)

3.

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

Modify evolution eqs so as to forma symmetric hyperbolic system

∂t




u
λ
u
λ









=

A
F

0
0

A
F

F̄
0





∂i


u
λ



∂i



u
λ



Since the total system is designed to have symmetric hyperbolicity, the evolution is supposed to be
unique. Brodbeck et al showed analytically that such a decay of λs can be seen for suﬃciently smallλ(> 0) with a choice of appropriate combinations of αs and βs.
Brodbeck et al presented a set of equations based on Frittelli-Reula’s symmetric hyperbolic formulation [24]. The version of Ashtekar’s variables was presented by the authors [37] for controlling theconstraints or reality conditions or both. The numerical tests of both the Maxwell-λ-system and the
Ashtekar-λ-system were performed [47], and conﬁrmed to work as expected. Although it is questionable whether the recovered solution is true evolution or not [41], we think the idea is quite attractive.
To enforce the decay of errors in its initial perturbative stage seems the key to the next improvements,
which are also developed in the next section on “adjusted systems”.
However, there is a high price to pay for constructing a λ-system. The λ-system can not be
introduced generally, because (i) the construction of λ-system requires the original evolution equations
to have a symmetric hyperbolic form, which is quite restrictive for the Einstein equations, (ii) the ﬁnal
system requires many additional variables and we also need to evaluate all the constraint equations at
every time step, which is a hard task in computation. Moreover, (iii) it is not clear that the λ-system
is robust enough for non-linear violation of constraints, or that λ-system can control constraints which
do not have any spatial diﬀerential terms.
The “adjusted system” Next, we propose an alternative system which also tries to control theviolation of constraint equations actively, which we named “adjusted system”. We think that thissystem is more practical and robust than the previous λ-system.

The Adjusted system (procedures):Box 2.91.

Prepare a set of evolution eqs.

∂t u = J∂i u + K

2.

Add constraints in RHS∂t u = J∂i u + K +κC 

Choose the coeﬀ. κ so as to make theeigenvalues of the homogenized adjusted∂t C eqs negative reals or pure imaginary.

∂t C = D∂i C + EC
∂t C = D∂i C + EC +F ∂i C + GC

3.







Lecture Note at APCTP winter school 2003: H.Shinkai

30

The process of adjusting equations is a common technique in other re-formulating eﬀorts as we reviewed. However, we try to employ the evaluation process of constraint ampliﬁcation factors as analternative guideline to hyperbolization of the system. We will explain these issues in the next section.
2.5.7

A uniﬁed treatment: Adjusted System

This section is devoted to present our idea of “asymptotically constrained system”. Original references
can be found in [47, 48, 39, 49].
Procedures : Constraint propagation equations and Proposals Suppose we have a set ofdynamical variables ua (xi , t), and their evolution equations,∂t ua = f (ua , ∂i ua , · · ·),(2.85)C α (ua , ∂i ua , · · ·) ≈ 0.

(2.86)and the (ﬁrst class) constraints,
Note that we do not require (2.85) forms a ﬁrst order hyperbolic form. We propose to investigate theevolution equation of C α (constraint propagation),∂t C α = g(C α , ∂i C α , · · ·),(2.87)

for predicting the violation behavior of constraints in time evolution. We do not mean to integrate(2.87) numerically together with the original evolution equations (2.85), but mean to evaluate them
analytically in advance in order to reformulate the equations (2.85).
There may be two major analyses of (2.87); (a) the hyperbolicity of (2.87) when (2.87) is a ﬁrst ordersystem, and (b) the eigenvalue analysis of the whole RHS in (2.87) after a suitable homogenization.
However, as we critically viewed the hyperbolization road in the previous section, we prefer to proceed
the road (b).
Ampliﬁcation Factors of Constraint Propagation equations:
We propose to homogenize (2.87) by a Fourier transformation, e.g.
∂t Ĉ α = ĝ(Ĉ α ) = M α β Ĉ β ,where C(x, t)ρ =

Ĉ(k, t)ρ exp(ik · x)d3 k,

Box 3.1(2.88)

then to analyze the set of eigenvalues, say Λs, of the coeﬃcient matrix, M α β , in (2.88).
We call Λs the constraint ampliﬁcation factors (CAFs) of (2.87).
The CAFs predict the evolution of constraint violations. We therefore can discuss the “distance” tothe constraint surface using the “norm” or “compactness” of the constraint violations (although we
do not have exact deﬁnitions of these “· · ·” words).
The next conjecture seems to be quite useful to predict the evolution feature of constraints:
Conjecture on Constraint Ampliﬁcation Factors (CAFs):

Box 3.2(A) If CAF has a negative real-part (the constraints are forced to be diminished), thenwe see more stable evolution than a system which has positive CAF.
(B) If CAF has a non-zero imaginary-part (the constraints are propagating away), thenwe see more stable evolution than a system which has zero CAF.

Lecture Note at APCTP winter school 2003: H.Shinkai

31

We found that the system becomes more stable when more Λs satisfy the above criteria. A generalfeature of the constraint propagation is reported in [50].
The above features of the constraint propagation, (2.87), will diﬀer when we modify the original
evolution equations. Suppose we add (adjust) the evolution equations using constraints∂t ua = f (ua , ∂i ua , · · ·) + F (C α , ∂i C α , · · ·),(2.89)

then (2.87) will also be modiﬁed as∂t C α = g(C α , ∂i C α , · · ·) + G(C α , ∂i C α , · · ·).

(2.90)

Therefore, the problem is how to adjust the evolution equations so that their constraint propagations
satisfy the above criteria as much as possible.
Applications For the Maxwell equation and the Ashtekar version of the Einstein equations, we
numerically found that this idea works to reduce the violation of constraints, and that the eﬀects aremuch better than by constructing its symmetric hyperbolic versions [38, 47].
Applications to ADM The idea was applied to the standard ADM formulation which is not
hyperbolic and several attractive adjustments were proposed [48, 39]. We made various predictions how
additional adjusted terms will change the constraint propagation. Systematic numerical comparisons
are also progressing, and we show two sample plots here.
Figure 4 (a) is a test numerical evolution of Detweiler-type adjustment [20] on the Minkowskii
background. We see the adjusted version gives convergence on to the constraint surface by arranging
the magnitude of the adjusting parameter, κ. Figure 4 (b) is obtained by a 3-dimensional numerical evolution of weak gravitational wave, the so-called Teukolsky wave [44]. The lines are of theoriginal/standard ADM evolution equations, Detweiler-type adjustment, and a part of Detweiler-typeadjustment. For a particular choice of κ, we observe again the L2 norm of constraint (violation of
constraints) is reduced than the standard ADM case, and can evolve longer than that.
Notion of Time Reversal Symmetry During the comparisons of adjustments, we found that itis necessary to create time asymmetric structure of evolution equations in order to force the evolution
on to the constraint surface. There are inﬁnite ways of adjusting equations, but we found that if we
follow the guideline Box 3.5, then such an adjustment will give us time asymmetric evolution.

Trick to obtain asymptotically constrained system:= Break the time reversal symmetry (TRS) of the evolution equation.

Box 3.5

1. Evaluate the parity of the evolution equation.
By reversing the time (∂t → −∂t ), there are variables which change their signatures (parity(−)) [e.g. Kij , ∂t γij , Mi , · · ·], while not (parity (+)) [e.g. gij , ∂t Kij , H, · · ·].
2. Add adjustments which have diﬀerent parity of that equation.
For example, for the parity (−) equation ∂t γij , add a parity (+) adjustment κH.

One of our criteria, the negative real CAFs, requires breaking the time-symmetric features of theoriginal evolution equations. Such CAFs are obtained by adjusting the terms which break the TRS ofthe evolution equations, and this is available even at the standard ADM system.

Lecture Note at APCTP winter school 2003: H.Shinkai
ADM and its adjusted versions-- Teukolsky wave evolution --Detweiler's adjustments
on Minkowskii spacetime
L= - 0.01
Log1 0 (L2 norm of constraints )10
L=-0.01
L=0.0
L=+0.01
L=+0.02
L=+0.03
L=+0.35

0.0

L= + 0.035-2.0

L=0 (standard ADM)

-4.0

-6.0

L= + 0.01-8.0

L= + 0.03

L2 norm of Hamiltonian constraint

2.0

32

0

10-1
Original ADM

10-2Standard ADM

10-3
Adjusted ADM (SimpDet)

10-4
Adjusted ADM (Det)

10-5
10-6

-10.0
0.0

0.5

1.0

time

1.5

2.0

0

50

100

150

200

250

300

time

Figure 3: Comparisons of numerical evolution between adjusted ADM systems. (a) Demonstration of
the Detweiler’s modiﬁed ADM system on Minkowskii background spacetime, 1-dimensional simulation.
The L2 norm of the constraints HADM and MADM is plotted in the function of time. Artiﬁcial
error was added at t = 0.25. L is the parameter in Detweiler’s adjustment. We see the evolution
is asymptotically constrained for small κ > 0. (b) L2 norm of the Hamiltonian constraint HADM
of evolution using ADM/adjusted ADM formulations for the case of Teukolsky wave, 3-dimensional
simulation.
Applications to BSSN This analysis was also applied to explain the advantages of the BSSN formulation, and again several alternative adjustments to BSSN equations were proposed [49]. Recently
Yo et al[45] reported their simulations of stationary rotating black hole, and mentioned that one of
our proposal was contributed to maintain their evolution of Kerr black hole (J/M up to 0.9M ) forlong time (t ∼ 6000M ). Their results also indicates that the evolved solution is closed to the exact
one, that is, the constrained surface.
2.5.8

OutlookWhat we have achieved• The constraint propagation features become diﬀerent by simply adding constraint terms to theoriginal evolution equations (we call this the adjustment of the evolution equations).
• There is a constraint-violating mode in the standard ADM evolution system when we apply it
to a single non-rotating black hole space-time, and its growth rate is larger near the black-holehorizon.
• Such a constraint-violating mode can be killed if we adjust the evolution equations with aparticular modiﬁcation using constraint terms (Box 2.7). An eﬀective guideline is to adjust
terms as they break the time-reversal symmetry of the equations (Box 3.5).
• Our expectations are borne out in simple numerical experiments using the Maxwell, Ashtekar,and ADM systems. However the modiﬁcations are not yet perfect to prevent non-linear growthof the constraint violation.
• We understand why the BSSN formulation works better than the ADM one in the limitedcase (perturbative analysis in the ﬂat background), and further we proposed modiﬁed evolution
equations along the lines of our previous procedure.

Lecture Note at APCTP winter school 2003: H.Shinkai

33

The common key to the problem is how to adjust the evolution equations with constraints. Any
adjusted systems are mathematically equivalent if the constraints are completely satisﬁed, but this is
not the case for numerical simulations. Replacing terms with constraints is one of the normal steps
when people hyperbolize equations. Our approach is to employ the evaluation process of constraint
ampliﬁcation factors for an alternative guideline to hyperbolization of the system.
Final remarks If we say the ﬁnal goal of this project is to ﬁnd a robust algorithm to obtain longterm accurate and stable time-evolution method, then the recipe should be a combination of (a)formulations of the evolution equations, (b) choice of gauge conditions, (c) treatment of boundary
conditions, and (d) numerical integration methods. We are in the stages of solving this mixed puzzle.
The ideal almighty algorithm may not exit, but we believe our accumulating experience will make theones we do have more robust and automatic.
I have written this review from the viewpoint that the general relativity is a constrained dynamicalsystem. This is not only a proper problem in general relativity, but also in many physical systems
such as electrodynamics, magnetohydrodynamics, molecular dynamics, mechanical dynamics, and soon. Therefore sharing the thoughts between diﬀerent ﬁeld will deﬁnitely accelerate the progress.

References[1] A. Abrahams, A. Anderson, Y. Choquet-Bruhat, and J.W. York, Jr., Phys. Rev. Lett. 75, 3377 (1995);
Class. Quant. Grav. 14, A9 (1997).
[2] M. Alcubierre, Phys. Rev. D 55, 5981 (1997).
[3] M. Alcubierre and J. Massó, Phys. Rev. D 57, R4511 (1998).
[4] M. Alcubierre, G. Allen, B. Brügmann, E. Seidel and W-M. Suen, Phys. Rev. D 62, 124011 (2000).
[5] M. Alcubierre, W. Benger, B. Brügmann, G. Lanfermann, L. Nerger, E. Seidel, and R. Takahashi, Phys.
Rev. Lett. 87 271103 (2001).
[6] M. Alcubierre, B. Brügmann, T. Dramlitsch, J.A. Font, P. Papadopoulos, E. Seidel, N. Stergioulas, and
R. Takahashi, Phys. Rev. D 62, 044034 (2000).
[7] A. Anderson, Y. Choquet-Bruhat, J.W. York, Jr., Topol. Methods in Nonlinear Analysis, 10, 353 (1997).
[8] A. Anderson and J. W. York, Jr, Phys. Rev. Lett. 82, 4384 (1999).
[9] R. Arnowitt, S. Deser and C.W. Misner, in Gravitation: An Introduction to Current Research, ed. by
L.Witten, (Wiley, New York, 1962).
[10] A. Ashtekar, Phys. Rev. Lett. 57, 2244 (1986); Phys. Rev. D36, 1587 (1987).
[11] T.W. Baumgarte and S.L. Shapiro, Phys. Rev. D 59, 024007 (1999).
[12] J.M. Bardeen, L.T. Buchman, Phys. Rev. D. 65, 064037 (2002).
[13] C. Bona, J. Massó, Phys. Rev. D 40, 1022 (1989); Phys. Rev. Lett. 68, 1097 (1992).
[14] C. Bona, J. Massó, E. Seidel and J. Stela, Phys. Rev. Lett. 75, 600 (1995); Phys. Rev. D 56, 3405 (1997).
[15] C. Bona, J. Massó, E. Seidel, and P. Walker, gr-qc/9804052.
[16] O. Brodbeck, S. Frittelli, P. Hübner, and O.A. Reula, J. Math. Phys. 40, 909 (1999).
[17] G. Calabrese, J. Pullin, O. Sarbach, and M. Tiglio, Phys. Rev. D 66, 041501 (2002).
[18] M.W. Choptuik, Phys. Rev. D 44, 3124 (1991).

Lecture Note at APCTP winter school 2003: H.Shinkai

34[19] Y. Choquet-Bruhat and J.W. York, Jr., C.R. Acad. Sc. Paris 321, Série I, 1089, (1995).
[20] S. Detweiler, Phys. Rev. D 35, 1095 (1987).
[21] H. Friedrich, Proc. Roy. Soc. A375, 169 (1981); ibid. A378, 401 (1981).
[22] S. Frittelli, Phys. Rev. D 55, 5992 (1997).
[23] S. Frittelli and R. Gomez, J. Math. Phys. 41, 5535 (2000).
[24] S. Frittelli and O.A. Reula, Phys. Rev. Lett. 76, 4667 (1996).
[25] S. D. Hern, PhD thesis, gr-qc/0004036.
[26] L.E. Kidder, M.A. Scheel, S.A. Teukolsky, Phys. Rev. D 64, 064017 (2001).
[27] L. Lindblom and M. Scheel, gr-qc/0206035[28] C.W. Misner, K.S. Thorne, J.A. Wheeler, Gravitation, (Freeman, N.Y., 1973).
[29] T. Nakamura, K. Oohara and Y. Kojima, Prog. Theor. Phys. Suppl. 90, 1 (1987).
[30] T. Nakamura and K. Oohara, in Frontiers in Numerical Relativity edited by C.R. Evans, L.S. Finn, andD.W. Hobill (Cambridge Univ. Press, Cambridge, England, 1989).
[31] Proceedings of Black Holes and Gravitational Waves ed. by T. Nakamura and H. Kodama, as the issueof Prog. Theor. Phys. Suppl. 136 (1999).
[32] O. Sarbach, G. Calabrese, J. Pullin, and M. Tiglio, Phys. Rev. D 66, 064002 (2002).
[33] M.A. Scheel, T.W. Baumgarte, G.B.Cook, S.L. Shapiro, S.A. Teukolsky, Phys. Rev. D 56, 6320 (1997).
[34] M.A. Scheel, T.W. Baumgarte, G.B. Cook, S.L. Shapiro, S.A. Teukolsky, Phys. Rev. D 58, 044020 (1998).
[35] M. Shibata and T. Nakamura, Phys. Rev. D52, 5428 (1995).
[36] M. Shibata and K. Uryu, Prog. Theor. Phys. 107, 265 (2002).
[37] H. Shinkai and G. Yoneda, Phys. Rev. D 60, 101502 (1999).
[38] H. Shinkai and G. Yoneda, Class. Quant. Grav. 17, 4799 (2000).
[39] H. Shinkai and G. Yoneda, Class. Quant. Grav. 19, 1027 (2002).
[40] H. Shinkai and G. Yoneda, in Progress in Astronomy and Astrophysics (Nova Science Publ) to be published. The manuscript is available as gr-qc/0209111.
[41] F. Siebel and P. Hübner, Phys. Rev. D 64, 024021 (2001).
[42] L. Smarr, J.W. York, Jr., Phys. Rev. D 17, 2529 (1978).
[43] J.M. Stewart, Class. Quant. Grav. 15, 2865 (1998).
[44] S.A. Teukolsky, Phys. Rev. D 26, 745 (1982).
[45] H-J. Yo, T.W. Baumgarte and S.L. Shapiro, Phys. Rev. D 66, 084026 (2002).
[46] G. Yoneda and H. Shinkai, Phys. Rev. Lett. 82, 263 (1999); Int. J. Mod. Phys. D. 9, 13 (2000).
[47] G. Yoneda and H. Shinkai, Class. Quant. Grav. 18, 441 (2001).
[48] G. Yoneda and H. Shinkai, Phys. Rev. D 63, 124019 (2001).
[49] G. Yoneda and H. Shinkai, Phys. Rev. D 66, 124003 (2002).
[50] G. Yoneda and H. Shinkai, to be published in Class. Quant. Grav. (gr-qc/0209106).
[51] J.W. York, Jr., in Sources of Gravitational Radiation, ed. by L.Smarr, (Cambridge, 1979).

Lecture Note at APCTP winter school 2003: H.Shinkai

3

35

Alternative Approaches to Numerical Relativity

In the pioneer days of numerical relativity (70s - 80s), people set the destination of a code so as it
simultaneously (1) avoids singularities, (2) handles black-holes, (3) maintains high accuracy, and (4)runs forever. This goal was also called “Holy Grail” of numerical relativity [1]. Various approaches
have been proposed and tested for these purposes. In this chapter, I introduce several alternative orcomplemental approaches to the standard 3+1 (ADM, or Cauchy) approach.

3.1

Full numerical, but diﬀerent foliations

3.1.1

Characteristic foliation

Cauchy versus Characteristic
Figure 4: Two major foliations for seeking evolution in general relativity.

Cauchy approachor ADM 3+1 formulation

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

time direction

Σ: Initial 3-dimensional Surface

ingoing
direction

outgoing
direction

S: Initial 2-dimensional Surface

pioneers

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

variables

easy to understand the concept oftime evolution

foliation

has Hamilton structure

initial data
evolutionneed to solve constraints
PDEs
need to avoid constraint violation
singularity
need to avoid by some method
disadvantages can not cover space-time globally

Characteristic (2+2) evolution
Bondi et al [2], Sachs [3], Penrose[4]Numerical works by Stewart et al [5]has geometrical meanings
1 complex function related to 2 GW
polarization modes
allows implementation of Penrose’s
space-time compactiﬁcation
no constraints
ODEs with consistent conditions
propagation eqs along the light rays
can truncate the griddiﬃculty in treating caustics
hard to treat matter

Table 1: Comparison of Cauchy and characteristic approaches. See reviews by Winicour [6].

Lecture Note at APCTP winter school 2003: H.Shinkai

36

A connection formula from ADM to Newman-Penrose The Newman-Penrose formulation[7, 8] has many advantages, especially for treating gravitational wave dynamics.
• Natural framework for calculations in radiative space-time.
• Variables have geometrical meanings.
• Practical advantages in treating Petrov type-D space-time.
• Closely related with spinor formalism.
Newman-Penrose’s variables are based on real-valued null vectors l, n and complex conjugate null
vectors m, m, which satisfy
l · n = la na = la na = 1,m · m = ma ma = ma ma = −1,

and else 0.

(3.1)

This set of null basis (la , na , ma , ma ) have relations with orthogonal tetrad basis (tâ , xâ , y â , z â )1la = oA oA = √ (tâ + z â ),21na = ιA ιA = √ (tâ − z â ),21
ma = oA ιA = √ (xâ − iy â )21
ma = ιA oA = √ (xâ + iy â ),
2

(3.2)
(3.3)

where I also put spinor expressions (oA , ιA ) of those. The indice rules arela = gab lb ,la = g ab lb ,gab = metric(3.4)

xâ = ηâb̂ xb̂ ,xâ = η âb̂ xb̂ ,ηâb̂ = (1, −1, −1, −1)(3.5)

gab = 2l(a nb) − 2m(a mb) ,g ab = 2l(a nb) − 2m(a mb)(3.6)

Metric gab will be recovered bygâb̂ = tâ tb̂ − xâ xb̂ − yâ yb̂ − zâ zb̂ .

(3.7)

The Weyl curvature Cabcd is deﬁned as1
Cabcd = Rabcd − ga[c Rd]b + gb[c Rd]a − Rga[c gd]b .
3

(3.8)

The 10 components of Weyl curvature are expressed by the following 5 complex scalars [9];
Ψ0 ≡ ψABCD oA oB oC oD = Cabcd la mb lc md ,na -directed transverse component, {4, 0}

Ψ1 ≡ ψABCD o o o ι = Cabcd l n l m ,a

A B C D

a b c

d

(3.9)

n -directed longitudinal component{2, 0} (3.10)

Ψ2 ≡ ψABCD o o ι ι = Cabcd l m m̄ n ,‘Coulomb’ component, {0, 0}

Ψ3 ≡ ψABCD o ι ι ι = Cabcd l n m̄ n ,l -directed longitudinal component, {−2, 0}(3.12)

Ψ4 ≡ ψABCD ι ι ι ι = Cabcd n m̄ n m̄ ,la -directed transverse component, {−4, 0} (3.13)

A B C DA B C DA B C D

a

a b
a

b

c d

c d

b c

d

(3.11)awhere {p, q} indicates spin- and boost-weighted type and prime-operation will be deﬁned later.
Gunnarsen-Shinkai-Maeda [11] derived a transformation formula of Weyl scalar Ψi from ADM
variables (γij , Kij ), motivated by an application to interpret numerically generated space-time. Here,we consider vacuum space-time. Let (M, ηab ) be real, 4-dimensional Lorentz vector space with volume
form εabcd ; εabcd εabcd = −4!. Let (ta , xa , y a , z a ) be orthonormal basis of (M, ηab ), and deﬁne
ta ta = +s (s = ±1),εabc = εabcd td ,(3.14)

Lecture Note at APCTP winter school 2003: H.Shinkai

37

where the tensor ﬁeld εabc = ε[abc] satisﬁes εabc εabc = 3! We formulate our equations in the signatures
both (+, −, −, −) and (−, +, +, +) by choosing s = 1 or −1, respectively4 , because the former notationis common in working with the spinors.
First, we deﬁne the Weyl curvature Cabcd by (3.8) and decompose those into its electric and amagnetic components,
Eab ≡ − Cambn tm tn , Bab ≡ − ∗ Cambn tn tm ,(3.15)
where ∗ Cabcd = 12 εabmn Cmncd is a dual of the Weyl tensor. These decomposed elements Eab and Bab
are also presented by the 3-metric γab and the extrinsic curvature Kab as [12]2(3.16)
Eab = Rab − Kam Kbm + KKab − Λγab ,3(3.17)
Bab = εamn Dm Knb .
This is why we emphasize that our inputs are ‘3+1’ elements. It follows from two constraint equations
that the ﬁelds Eab , Bab are both trace-free and symmetric. We can reconstruct the Weyl curvaturecompletely from Eab and Bab by
Cabcd = 4t[a Eb][c td] + 2εabm Bm[c td] + 2εcdm Bm[a tb] + εabm εcdn Emn .

(3.18)The next step is to choose a unit vector ﬁeld ẑ a on Σ, and to decompose Eab , Bab into components
along and perpendicular to ẑ a . We set
e = Eab ẑ a ẑ b ,ea = Ebc ẑ b (δac + sẑa ẑ c ),eab = Ecd (δac + sẑa ẑ c )(δb d + sẑb ẑ d ) + 12 esab ,b = Bab ẑ a ẑ b ,ba = Bbc ẑ b (δac + sẑa ẑ c ),(3.19)
1ccddbab = Bcd (δa + sẑa ẑ )(δb + sẑb ẑ ) + 2 bsab ,where sab = γab − ẑa ẑb . We note that Eab , Bab is again reconstructed from (3.19)
Eab = eẑa ẑb + 2e(a ẑb) + eab − (1/2)sab e.

(3.20)

Bab = bẑa ẑb + 2b(a ẑb) + bab − (1/2)sab b.

(3.21)Such decompositions will be useful to discuss the eﬀects of curvatures on the transversal plane to theẑ a direction.
We put a rotation operator on the plane spanned by x̂a and ŷa as,Jab ≡ εabcd ẑc td .

(3.22)

It is easy to check this mapping preserves sab , and is also easy to check Jac Jc b = −(δab + sẑa ẑ b ), whichshows us Jab has a complex structure, i.e., Jab lets us deﬁne complex multiples of vectors xa ∈ Pz ,according to the formula (m + in)xa = mxa + nJb a xb . In short, Jab expresses a rotation by 90 degrees
in the plane orthogonal to ẑ a .
By substituting (3.18) and (3.2, ??) into (3.9)-(3.13), we get Ψi using (3.19) and (3.22):
Ψ0 = −(eab + sJac bbc )ma mb ,√
Ψ1 = −(s/ 2) (ea + sJac bc )ma ,(3.23)

Ψ2 = −(1/2) (e + ib),√
Ψ3 = −(s/ 2)(ea − sJac bc )m̄a ,(3.25)(3.26)

Ψ4 = −(eab − sJa bbc )m̄ m̄ .

(3.27)

cab(3.24)

This relation has been applied to many groups’ numerical codes, and helps their simulation’s physicalunderstandings. Weyl scalars are also useful for evaluating Riemann (Kretchman) invariant as
Cabcd C abcd = Ψ4 Ψ0 − 4Ψ1 Ψ3 + 3Ψ22 .

(3.28)

Rabcd Rabcd = Cabcd C abcd + 2Rab Rab − (1/3)R2 .

(3.29)

Note that
4

That is, the metric is ηab = s(ta tb − xa xb − ya yb − za zb ).

Lecture Note at APCTP winter school 2003: H.Shinkai
3.1.2

38

Characteristic approach

As I described in Table 1, the characteristic approach is quite attractive unless the system does notmake a caustics in null hypersurface. Since the event horizon of black-hole is itself a characteristic
hypersurface, the characteristic technique is powerful tool as a stand-alone. It also allows us to express
compactiﬁed manifold, so that we can seek gravitational wave dynamics at time inﬁnity.
Although the applicability of characteristic foliation is limited, numerical codes are developed
extensively by by Pittsburgh group together with the Binary Black Hole Grand Challenge Alliance(1993-1998), http://www.npac.syr.edu/projects/bh/.
They adapted the Bondi-Sachs form of the metric for the null foliation. The coordinate areconstructed from a family of outgoing null hypersurfaces, emanating from a worldtube or a timelike
geodesic of topology S 2 × R, which is labelled with a parameter u. Each null ray on a speciﬁc
hypersurface is labelled with xA where (A = 2, 3), and let r be a surface area distance (i.e. surfaces at
r=constant have area 4πr2 ). Then, the resulting coordinate is xa = (u, r, xA ), and the Bondi-Sachs
form of the metric takes

ds2 = − e2βV− r2 hAB U A U B du2 − 2e2β dudr − 2r2 hAB U B dudxA + r2 hAB dxA dxB .
r(3.30)

This metric has six real ﬁeld variables, V, β, U A , and hAB . V can be understand in the analogueof the Newtonian potential, β represents the expansion of the light rays as they propagate radially.
hAB represents the conformal intrinsic geometry, which contains the two degrees of radiation freedom.
Note that on the r =const. timelike world tube, the intrinsic metric can be expressed similarly to the
Cauchy decomposition,(3)

ds2 = −e2βV 2du + r2 hAB (dxA − U A du)(dxB − U B du).
r

(3.31)For example, a Schwarzschild geometry in outgoing Eddington-Finkelstein coordinates is given by the
choice β = 0, V = r − 2m, U A = 0 and hAB be a unit sphere metric.
For a single black hole case the unlimited evolution was reported [13, 14, 15]. For the head-on
collision of black-holes, the detail analysis of the dynamics of the event horizon was reported [16].
Recently, the formulation is studied also in the direction of including matter dynamics, such as theimplementation of high resolution schemes [17] and the weak pressure ﬂuids [18].
3.1.3

Cauchy-characteristic matching approachThis idea is to combine Cauchy evolution (interior) with characteristic evolution (exterior), in orderto supply precise radiated waveform in the binary coalescence problem. Two evolution schemes arematched on the worldtube, and both sides of foliation supply the outer boundary values to the other.
Numerical codes were developed independently by Southampton group [19] and Pittsburgh group [20]in the middle 90s. A signiﬁcant advantage on the treatment of the outer boundary in Cauchy evolution
region was reported for the cases of pure gravitational wave problem.
3.1.4

Hyperbolical foliation, conformal ﬁeld equations

A series of works by Friedrich [21] attempted to construct a 3 + 1 formulation with hyperboloidal foliations (i.e. asymptotically null foliations), and with conformal compactiﬁcation. This is the ultimate
plan to remove the outer boundary problem in numerical simulation, and to provide a suitable foliation
for gravitational radiation problem. However, the current equations are rather quite complicated. In
its metric-based expression [22], the evolution variables are 57; γij , Kij , the connection coeﬃcientsγ a bc , projections (0,1) R̂a = nb γa c R̂bc and (1,1) R̂ab = γa c γb d R̂bd of 4-dimensional Ricci tensor R̂ab , theelectric and magnetic components of the rescaled Weyl tensor Cabc d , and the conformal factor Ω andLecture Note at APCTP winter school 2003: H.Shinkai

39

its related quantities Ω0 ≡ na ∇a Ω, ∇a Ω, ∇a ∇a Ω. By specifying suitable gauge functions (α, β a , R)
where R is the Ricci scalar, then the total system forms a symmetric hyperbolic system. Applications
to numerical relativity are in progress, but have not yet reached the stage of applying evolution in a
non-trivial metric. For more details, see reviews e.g. by Frauendiener [23] or by Husa [24].

References[1] S.L. Shapiro and Teukolsky, in Dynamical Space-times and Numerical Relativity, edited by J.
Centrella (Cambridge University Press, 1986).
[2] H. Bondi, M.J.G. van der Burg, and A.W.K. Metzner, Proc. R. Soc. London A 269, 21 (1962).
[3] R. Sachs, Proc. R. Soc. London A 270, 103 (1962).
[4] R. Penrose, Phys. Rev. Lett. 10, 66 (1963).
[5] J.M. Stewart, H.Friedrich, Proc. R. Soc. Lond. A 384, 427 (1982)
R.W. Corkill, J.M. Stewart, Proc. R. Soc. Lond. A 386, 373 (1983).
[6] J. Winicour, Livng Rev. Relativ. 2001-3 at http://www.livingreviews.org/
J. Winicour, Prog. Theor. Phys. Suppl. 136, 57 (1999).
J. Winicour, in Proceeding of CIMENICS 2000, The Vth International Congress on Numerical
Methods in Engineering and Applied Science (Puerto La Cruz, Venezuela, March 2000). available
as gr-qc/0003029.
[7] E. Newman and R. Penrose, J. Math. Phys. 3, 566 (1962); ibid. 4, 998 (1963) errata.
[8] R. Penrose and W. Rindler, Spinors and Space-time, vol.1 and 2, (Cambridge University Press,
1984, 1986)
S. Chandrasekhar, The Mathematical Theory of Black Holes, (Oxford Univ. Press, 1992).
J. Stewart, Advanced General Relativity, (Cambridge Univ. Press, 1990).
D. Kramer, H. Stephani, M. MacCallum, E. Herlt, Exact solutions of Einstein’s ﬁeld equations,(Cambridge University Press, 1980).
[9] R.K. Sachs, Proc. Roy. Soc. London, A264, 309 (1961); ibid A270, 103 (1962).
P. Szekeres, J. Math. Phys. 6, 1387 (1965).
[10] R. Gerosh, A. Held and R. Penrose, J. Math. Phys. 14, 874 (1973).
[11] L. Gunnarsen, H. Shinkai, and K. Maeda, Class. Quant. Grav. 12 (1995) 133[12] L. Smarr, Ann. N.Y. Acad. of Sci. 302, 569 (1977).
[13] N.T. Bishop, R. Gómez, L Lehner, and J. Winicour, Phys. Rev. D 54, 6153 (1996).
[14] N.T. Bishop, R. Gómez, L Lehner, M. Maharaj, and J. Winicour, Phys. Rev. D 56, 6298 (1997).
R. Gómez, L Lehner, M. Maharaj, and J. Winicour, Phys. Rev. D 57, 4778 (1998).
[15] The Binary Black Hole Grand Challenge Alliance, R. Gómez et al., Phys. Rev. Lett. 80, 3915
(1998).
[16] L Lehner, N.T. Bishop, R. Gómez, B. Szilagyi, and J. Winicour, Phys. Rev. D 60, 044005 (1999).
[17] P. Papadopoulos and J. Font, Phys. Rev. D 61, 024015 (2000).

Lecture Note at APCTP winter school 2003: H.Shinkai

40[18] N.T. Bishop, R. Gómez, L Lehner, M. Maharaj, and J. Winicour, Phys. Rev. D 60, 024005(1999).
[19] C. J. S. Clarke, R. A. d’Inverno, and J. A. Vickers, Phys. Rev. D 52, 6863 (1995).
M. Dubal, R. d’Inverno, and C. Clarke, Phys. Rev. D 52, 6868 (1995)
R. A. d’Inverno and J. A. Vickers, Phys. Rev. D 54, 4919 (1996).
R. A. d’Inverno and J. A. Vickers, Phys. Rev. D 56, 772 (1997).
M. R. Dubal, R. A. d’Inverno, and J. A. Vickers, Phys. Rev. D 58, 044019 (1998).
R.A. d’Inverno, M.R. Dubal, and E.A. Sarkies, Class. Quant. Grav. 17, 3157 (2000).
[20] R. Gómez, P. Laguna, P. Papadopoulos, and J. Winicour, Phys. Rev. D 54, 4719 (1996).
N.T. Bishop, R. Gómez, L. Lehner, and J. Winicour, Phys. Rev. D 54, 6153 (1996).
R. Gómez, R. Marsa, and J. Winicour, Phys. Rev. D 56, 6310 (1997).
N.T. Bishop, R. Gómez, P.R. Holvorcem, R.A. Matzner, P. Papadopoulos, and J. Winicour, J.
Comp. Phys. 136, 140 (1997).
[21] H. Friedrich, Proc. Roy. Soc. A375, 169 (1981); ibid. A378, 401 (1981).
[22] P. Hübner, Class. Quant. Grav. 16 2145 (1999), ibid. 16 2823 (1999), ibid. 18 1421 (2001), andibid. 18 1871 (2001).
[23] J. Frauendiener, Livng Rev. Relativ. 2000-4 at http://www.livingreviews.org/[24] S. Husa, in the Proceedings of the conference “The Conformal Structure of Spacetimes: Geometry,
Analysis, Numerics”, ed. by J. Frauendiener and H. Friedrich, by Springer Verlag, Lecture Notes
in Physics series), available as gr-qc/0204043;
in the Proceedings of the 2001 Spanish Relativity meeting, eds. L. Fernandez and L. Gonzalez,
to be published by Springer, Lecture Notes in Physics series), available as gr-qc/0204057.

3.2
3.2.1

Several approximations
Cauchy-perturbative matching: connection in spatial domain

This approximation intends to extract gravitational radiation information. and to provide stable outerboundary conditions for a Cauchy evolution numerical code. The procedure is to match the solutionof numerical simulation of non-linear Einstein equations to a set of one-dimensional linear equations
obtained through perturbation techniques over a curved background.
In early 90s, 1-dimensional version was implemented and tested 1-d [1]. The Binary Black HoleGrand Challenge Alliance and Rezzolla et al developed a 3-dimensional Cauchy-perturbative matching
technique [2, 3].
3.2.2

Close-limit approximation: connection in time domain

In binary black-hole problem, the ﬁnal outcome is a single black-hole which will ring down into
equilibrium. Using a perturbation theory, we know black-holes have quasi-normal modes that is
expected to be observed. The “close-limit approximation” is an extended idea to apply perturbation
theory just after one single common horizon around two black-holes formed [4]. On the validity of thisapproximation, i.e. from which regime and at what order of perturbation we need, was tested against
full numerical simulation using head-on collision of two black-holes [5]. Both radiated energy andwaveform agreed quite well even for large values of the momentum. (Numerical outputs are between
ﬁrst and second order perturbation outputs. )Recent development is one more step advanced. People try to evolve the system with perturbation
equations starting from a fully numerically evolved data. This would compensate the current limitation

Lecture Note at APCTP winter school 2003: H.Shinkai

41

of numerical relativity and derive more astrophysical conclusions. This idea was already tested forthe case of collapse of disks using Schwarzschild background [6], and also for the inspiralling binary
black-holes using Kerr background [7]. The latter project is named “Lazarus/Zorro”.
3.2.3

Quasi-spherical approximation

Hayward [8] proposed a new approximation scheme in a dual-null decomposition of space-time, with theaim of providing a computationally inexpensive estimate of the gravitational waveforms produced by ablack-hole or neutron-star collision, given a full numerical simulation up to (or close to) coalescence, oran analytical model thereof. The scheme truncates the Einstein equations by removing second-orderterms which would vanish in a spherically symmetric space-time.
Shinkai and Hayward [9] numerically implemented this scheme, testing it against angular momentum by applying it to Kerr black holes. As error measures, we take the conformal strain andspeciﬁc energy due to spurious gravitational radiation. The strain is found to be monotonic rather
than wavelike. The speciﬁc energy is found to be at least an order of magnitude smaller than the
1% level expected from typical black-hole collisions, for angular momentum up to at least 70% of the
maximum, for an initial surface as close as r = 3m .

References[1] A.M. Abrahams and C.R. Evans, Phys. Rev. D 37, 317 (1988), ibid. 42, 2585 (1990).
A.M. Abrahams et al., Phys. Rev. D 45, 3544 (1992).
A.M. Abrahams, R.H. Price, PRD 53, 1963 (1996), ibid. 53, 1963 (1996), ibid. 53, 1972 (1996).
[2] The Binary Black Hole Grand Challenge Alliance, A.M. Abrahams et al., Phys. Rev. Lett. 80,
1812 (1998).
[3] M. E. Rupright, A.M. Abrahams, and L. Rezzolla, Phys. Rev. D. 58, 044005 (1998).
L. Rezzolla, A.M. Abrahams, R.A. Matzner, M. E. Rupright, and S.L. Shapiro, Phys. Rev. D.
59, 064001 (1999).
[4] R. Price and J. Pullin, Phys. Rev. Lett. 72, 3297 (1994).
J. Pullin, Prog. Theor. Phys. Suppl. 136, 107 (1999).
[5] P. Anninos, R. Price, J. Pullin, E. Seidel, and W-M. Suen, Phys. Rev. D 52, 4462 (1995).
[6] A.M. Abrahams, S. Shapiro, and S. Teukolsky, Phys. Rev. D. 51, 4295 (1995).
[7] J. Baker, B. Brügmann, M. Campanelli, and C.O. Lousto, Class. Quant. Grav. 17, L149 (2000).
[8] S.A. Hayward, Phys. Rev. D 61, 101503 (2000).
[9] H. Shinkai and S.A. Hayward, Phys. Rev. D. 64, 044002 (2001).

Lecture Note at APCTP winter school 2003: H.Shinkai

442

Unsolved problems

I hope the statements below will inspire our future researches.

4.1

Gravitational Wave Physics and related problems

Please refer also Box 1.2 in §1 for numerical issues.
4.1.1

For extraction of physics• How binary behaves in the last stage of merger?
• What can we learn from waveform from the ﬁnal phase of binary merger?
• Can we determine equation of state of neutron star?
• Validity of alternative approaches?
• Validity of new approximations?
4.1.2

From numerical relativitic viewpoint• Physically reasonable initial data?
• Where to start the simulation? How connect from post-Newtonian evolution?
• How can we evolve the system stably?
• How can we achieve precise numerical simulations of coalescence of binary neutron stars and/or
black holes?
• How identify black-hole horizons?
• How to treat black hole singularity if it appears?
• How to extract gravitational wave?
• How can we manage the large-scale simulations?

4.2
4.2.1

Conjecture Hunting – reported and unreported issues –BH Uniqueness Theorems, No-hair Conjecture• Are colored BHs realistic?
• In higher dim.?
• Stable conﬁguration of Black String?
4.2.2

Cosmic Censorship Conjecture• Counter-examples?
• Strong version?

Lecture Note at APCTP winter school 2003: H.Shinkai
4.2.3

43

Gravitational Collapse and Hoop Conjecture• Deﬁnition of quasi-local mass?
• Validity?
• In higher dim.?
4.2.4

BH Thermodynamics• Why area, not volume?
• Under dynamical situation?
4.2.5

Dynamical Wormholes• Topology change in dynamical transition?
• New critical behavior for forming black-hole mass?
• Time-machine? (closed timelike curve? Chronological protection conjecture?)• Wormhole thermodynamics?
........... etc etc

4.3

Concluding Remark

There are many unsolved problems in general relativity. All realistic discussion requires numerical simulations. Our understanding for numerical procedures are accumulating, mature now. Computational
power is also suitable for actual researches.
Are you ready to go conjecture hunting?

Acknowledgment
I appreciate Prof. Sung-won Kim at Ewha Womans University, for inviting me at the Asia PaciﬁcCenter for Theoretical Physics Winter School on Gravitation and Cosmology, 2003. I am supportedby the special postdoctoral researchers program at RIKEN, and this work was supported partially by
the Grant-in-Aid for Scientiﬁc Research Fund of Japan Society of the Promotion of Science, (Kakenhi
Wakate B) No. 14740179.