Re-formulating the Einstein equations
for stable numerical simulations– Formulation Problem in 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 JapanGen Yoneda
yoneda@mse.waseda.ac.jp
Department of Mathematical Science Waseda University,
Ookubo, Shinjuku, Tokyo, 169-8555 Japan
December 1, 2002

Abstract
We review recent eﬀorts to re-formulate the Einstein equations for fully relativistic numericalsimulations. The so-called numerical relativity (computational simulations in general relativity) isa promising research ﬁeld matching with ongoing astrophysical observations such as gravitational
wave astronomy. 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 free evolution schemes. In this article, we ﬁrst review the eﬀorts of the community, categorizing them into the following three directions: (1) modiﬁcations of the standard
Arnowitt-Deser-Misner equations initiated by the Kyoto group, (2) rewriting of the evolution equations in hyperbolic form, and (3) construction of an “asymptotically constrained” system. We next
introduce our idea for explaining these evolution behaviors in a uniﬁed way using eigenvalue analysis of the constraint propagation equations. The modiﬁcations of (or adjustments to) the evolution
equations change the character of constraint propagation, and several particular adjustments using
constraints are expected to diminish the constraint-violating modes. We propose several new adjusted evolution equations, and include some numerical demonstrations. We conclude by discussing
some directions for future research.

This article is for a part of the book Progress in Astronomy and Astrophysics (Nova Science Publ.,
2003?). Also available as gr-qc/0209111.
1

Contents
1 Overview
1.1 Numerical Relativity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Formulation Problem in Numerical Relativity: Overview . . . . . . . . . . . . . .

3
3
4

2 The standard way and the three other roads
8
2.0 Strategy 0: The ADM formulation . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.1 Strategy 1: Modiﬁed ADM formulation by Nakamura et al (The BSSN formulation) 10
2.2 Strategy 2: Hyperbolic reformulations . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3 Strategy 3: Asymptotically constrained systems . . . . . . . . . . . . . . . . . . . 18
3 A uniﬁed treatment: Adjusted System
21
3.1 Procedures : Constraint propagation equations and Proposals . . . . . . . . . . . 21
3.2 Adjusted ADM formulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.3 Adjusted BSSN formulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4 Outlook
4.1 What we have achieved . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Next steps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Final remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36
36
37
37

A General expressions of ADM constraint propagation equations
A.1 The standard ADM equations and constraint propagations . . . . . . . . . . . . .
A.2 Constraint propagations for the adjusted ADM systems . . . . . . . . . . . . . .

39
39
40

B Numerical demonstrations using the Ashtekar formulation
41
B.1 The Ashtekar formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
B.2 Reformulate the Ashtekar evolution equations . . . . . . . . . . . . . . . . . . . . 42
B.3 Comparing numerical performance . . . . . . . . . . . . . . . . . . . . . . . . . . 44
References

47

2

1

Overview

1.1

Numerical Relativity

The theory of general relativity describes the nature of the strong gravitational ﬁeld. The Einsteinequation predicts quite unexpected phenomena such as gravitational collapse, gravitational waves,the expanding universe and so on, which are all attractive not only for researchers but also for the
public. The Einstein equation consists of 10 partial diﬀerential equations (elliptic and hyperbolic)for 10 metric components, and it is not easy to solve them for any particular situation. Over the
decades, people have tried to study the general-relativistic world by ﬁnding its exact solutions, bydeveloping approximation methods, or by simplifying the situations. Among these approaches, direct
numerical integration of the Einstein equations can be said to be the most robust way to study thestrong gravitational ﬁeld. This research ﬁeld is often called “numerical relativity”.

Numerical Relativity= Necessary for unveiling the nature of strong gravity. For example:

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, ...

Numerical relativity is now an essential ﬁeld in gravity research. The current mainstream innumerical relativity is to analyze the ﬁnal phase of compact binary objects (black holes and/or neutronstars) related to gravitational wave observations (see e.g. the conference proceedings [66]). Over the
past decades, many groups have developed their numerical simulations by trial and error. Simulations
require large-scale computational facilities, and long-time stable and accurate calculations. So far, wehave achieved certain successes in simulating the coalescence of binary neutron stars (see e.g. [75])and binary black holes (see e.g.[6]). However, people have still been faced with unreasonable numericalblow-ups at the end of simulations.
Diﬃculties in accurate/stable long-term evolution were supposed to be overcome by choosing
proper gauge conditions and boundary conditions. However, recent several numerical experiments
show that the (standard) Arnowitt-Deser-Misner (ADM) approach [12, 80, 98] is not the best formulation for numerics, and ﬁnding a better formulation has become one of the main research topics. A
majority of workers in the ﬁeld now believe in the existence of constraint-violating modes in most of
the formulations. Thus, the stability problem is now shedding light on the mathematical structure of
the Einstein equations.
The purpose of this article is to review the formulation problem in numerical relativity. Generally speaking, there are many open issues in numerical relativity, both theoretical (mathematical or
physical) and numerical. We list major topics in Box 1.2. More general and recent introductions tonumerical relativity are available, e.g. by d’Inverno (1996) [50], Seidel (1996/98/99) [73], Brügmann(2000) [25], Lehner (2001) [56], van Putten (2001) [88], and Baumgarte-Shapiro (2002) [16].

3Numerical Relativity – open issues

Box 1.2

0. How to select the foliation method of space-time
Cauchy (3 + 1), characteristic (2 + 2), or combined?
⇒ if the foliation is (3 + 1), then · · ·1. How to prepare the initial data
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?

⇐⇐⇐ this review

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?

1.2

Formulation Problem in Numerical Relativity: Overview

There are several diﬀerent approaches to simulating the Einstein equations. Among them the mostrobust way is to apply 3+1 (space + time) decomposition of space-time, as was ﬁrst formulated by
Arnowitt, Deser and Misner (ADM) [12] (we call this the “original ADM” system). 1
If we divide the space-time into 3+1 dimensions, the Einstein equations form a constrained system:
constraint equations and evolution equations. The system is quite similar to that of the Maxwell
equations (Box 1.3),
1
One alternative method of space-time foliation is the so-called characteristic approach (2 + 2 space-time decomposition). See reviews e.g. by d’Inverno (1996) [50], Winicour [89], Lehner (2001) [56]. Even in the 3+1 ADM approach, weconcentrate the standard ﬁnite diﬀerential scheme to express numerical expression of space-time. See e.g. Brewin [23]for a recent progress in a lattice method.

4

The Maxwell equations :The evolution equations: (∂t = ∂/∂t)Box 1.3∂t E = rot B − 4πj,and∂t B = −rot E(1.1)

div E = 4πρ,anddiv B = 0

(1.2)

Constraint equations:where people solve constraint equations on the initial data, and use evolution equations to follow the
dynamical behaviors.
In numerical relativity, this free-evolution approach is also the standard. This is because solving
the constraints (non-linear elliptic equations) is numerically expensive, and because free evolution
allows us to monitor the accuracy of numerical evolution. In black-hole treatments, recent “excision”
techniques do not require one to impose explicit boundary conditions on the horizon, which is alsoa reason to apply free evolution scheme. As we will show in the next section, the standard ADM
approach has two constraint equations; the Hamiltonian (or energy) and momentum constraints.
Up to a couple of years ago, the “standard ADM” decomposition [80, 98] of the Einstein equation
was taken as the standard formulation for numerical relativists. However, numerical simulations were
often interrupted by unexplained blow-ups (Figure.1). This was thought due to the lack of resolution,or inappropriate gauge choice, or the particular numerical scheme which was applied. However, after
the accumulation of much experience, people have noticed the importance of the formulation of theevolution equations, since there are apparent diﬀerences in numerical stability although the equations
are mathematically equivalent 2 .
At this moment, there are three major ways to obtain longer time evolutions. Of course, the ideas,
procedures, and problems are mingled with each other. The purpose of this article is to review all
three approaches and to introduce our idea to view them in a uniﬁed way. Table 1 is a list of references.
(1) The ﬁrst possibility is to use a modiﬁcation of the ADM system developed by the Kyoto group[63, 64] (often cited as Shibata and Nakamura [74]) and later re-introduced by Baumgarte andShapiro [15]. This is a combination of new variables, conformal decomposition, rescaling of
the conformal factor, and replacement of terms in the evolution equation using momentum
constraints (see §2.1).
(2) The second direction is to re-formulate the Einstein equations in a ﬁrst-order hyperbolic form.
This is motivated from 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 we know the characteristic speed of the system. In constructing hyperbolic systems,
the essential procedures are to adjust equations using constraints and to introduce new variables,
normally the spatially derivatived metric (see §2.2).
2

The word stability is used quite diﬀerent ways in the community.
• We mean by numerical stability a numerical simulation which continues without any blow-ups and in which dataremains on the constrained surface.
• Mathematical stability is deﬁned in terms of the well-posedness in the theory of partial diﬀerential equations, suchthat the norm of the variables is bounded by the initial data. See eq. (2.34) 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 to the original equations in its continuum limit) and stable (no error growing), then the simulation
represents the right (converging) solution. See [30] for the Einstein equations.

5errorBlow upt=0

Constrained Surface(satisfies
Einstein's constraints)

time

Figure 1: Origin of the problem for numerical relativists: Numerical evolutions depart from theconstraint surface.
(3) The third is to construct a system which is robust against the violation of constraints, such thatthe constraint surface is an attractor. The idea was ﬁrst proposed as a “λ-system” by Brodbeck
et al [24] in which they introduce artiﬁcial ﬂow to the constraint surface using a new variable
based on the symmetric hyperbolic system (see §2.3).
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 [94, 95, 78]. The method
is also applied to explain why the above approach (1) works, and also to propose alternative systems
based on the ADM [95, 78] and BSSN [96] equations. Section 3 is devoted to explain this idea withan analytical tool of the eigenvalue analysis of the constraint propagation.
We follow the notations of that of MTW[61], i.e. the signature of the space-time is (− + ++), andthe Riemann curvature is deﬁned asRµναβ ≡ ∂α Γµνβ − ∂β Γµνα + Γµσα Γσνβ − Γµσβ Γσνα
Rµν

≡

Rαµαν

(1.3)
(1.4)

We use µ, ν = 0, · · · , 3 and i, j = 1, · · · , 3 as space-time indices. The unit c = 1 is applied. The
discussion is mostly to the vacuum space-time, but the inclusion of matter is straightforward.

6

formulations(0) The standard ADM formulation
ADM
1962 Arnowitt-Deser-Misner [12, 80](1) The BSSN formulation
BSSN
1987 Nakamura et al [63, 64, 74]

1999 Baumgarte-Shapiro [15]1999 Alcubierre et al [8]1999 Frittelli-Reula [42]2002 Laguna-Shoemaker [55](2) The hyperbolic formulations
BM
1989 Bona-Massó [18, 19, 20]

KST1997 Bona et al [21]1999 Arbona et al [11]1995 Choquet-Bruhat and York [32]1995 Abrahams et al [1]1999 Anderson-York [10]1996 Frittelli-Reula [41]1996 Stewart [81]2001 Kidder-Scheel-Teukolsky [52]CFE

2002 Sarbach-Tiglio [69]1981 Friedrich[36]

CB-Y

FR

tetradAshtekar

1995 vanPutten-Eardley[86]1986 Ashtekar [13]1997 Iriondo et al [48]1999 Yoneda-Shinkai [92, 93](3) Asymptotically constrained formulations
λ-system to FR
1999 Brodbeck et al [24]to Ashtekar 1999 Shinkai-Yoneda [76]1987 Detweiler [33]adjusted to ADM
to ADM
2001 Shinkai-Yoneda [95, 78]
to BSSN
2002 Yoneda-Shinkai [96]numerical applications⇒ many⇒ 1987 Nakamura et al [63, 64]⇒ 1995 Shibata-Nakamura [74]
⇒ 2002 Shibata-Uryu [75] etc⇒ 1999 Baumgarte-Shapiro [15]⇒ 2000 Alcubierre et al [5, 7]⇒ 2001 Alcubierre et al [6] etc⇒ 2002 Laguna-Shoemaker [55]
⇒ 1995 Bona et al [20, 21, 22]⇒ 1997 Alcubierre, Massó [2, 4]⇒ 2002 Bardeen-Buchman [17]⇒ 1997 Scheel et al [70]
⇒ 1998 Scheel et al [71]⇒ 2002 Bardeen-Buchman [17]⇒ 2000 Hern [44]⇒ 2001 Kidder-Scheel-Teukolsky [52]
⇒ 2002 Calabrese et al [27]
⇒ 2002 Lindblom-Scheel [58]
⇒ 1998 Frauendiener [35]
⇒ 1999 Hübner [46]⇒ 1997 vanPutten [87]
⇒ 2000 Shinkai-Yoneda [77]
⇒ 2000 Shinkai-Yoneda [77, 94]⇒ 2001 Siebel-Hübner [79]
⇒ 2001 Yoneda-Shinkai [94]
⇒ 2001 Yoneda-Shinkai [95]
⇒ 2002 Mexico NR Workshop [59]
⇒ 2002 Mexico NR Workshop [59]
⇒ 2002 Yo-Baumgarte-Shapiro [90]Table 1: References to recent eﬀorts of reformulating the Einstein equations. We list mainly thosethat have been applied to actual numerical comparisons.

72

The standard way and the three other roads

2.0

Strategy 0: The ADM formulation

2.0.1

The original ADM formulation

The idea of space-time evolution was ﬁrst formulated by Arnowitt, Deser, and Misner (ADM) [12]. 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: such as
foliations of 3-dimensional hypersurface (Figure 2). This original ADM formulation was translated tonumerical relativists by Smarr [80] and York [98] in late 70s, with slightly diﬀerent notations. We referto the latter as the standard ADM formulation since this version is the starting point of the discussion.
The story begins by decomposing 4-dimensional space-time into 3 plus 1. The metric is expressedbyds2 = gµν dxµ dxν = −α2 dt2 + γij (dxi + β i dt)(dxj + β j dt),
(2.1)
where α and βj are deﬁned as α ≡ 1/ −g 00 and βj ≡ g0j , and called the lapse function and shift
vector, respectively. The projection operator or the intrinsic 3-metric gij is deﬁned as γµν = gµν +nµ nν ,where nµ = (−α, 0, 0, 0) [and nµ = g µν nν = (1/α, −β i /α)] is the unit normal vector of the spacelike
hypersurface, Σ (see Figure 2). By introducing the extrinsic curvature,Kij

1= − £n γij ,
2

(2.2)

and using the Gauss-Codacci relation, the Hamiltonian density of the Einstein equations can be writtenas
HGR = π ij γ̇ij − L,

where

L=

√

√
−gR = α γ[(3)R − K 2 + Kij K ij ],(2.3)

where π ij is the canonically conjugate momentum to γij ,π ij =∂L√= − γ(K ij − Kγ ij ),∂ γ̇ij(2.4)

omitting the boundary terms. The variation of HGR with respect to α and βi yields the constraints,
GRGR
and the dynamical equations are given by γ̇ij = δHand π̇ ij = − δHδhij .
δπ ijsurface normal line

shift vector, Ni
coordinate constant line

Nilapse function, Nt = constant hypersurface

Figure 2: 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.

8

2.0.2The standard ADM formulation

In the version of Smarr and York, Kij was used as a fundamental variable instead of the conjugate
momentum π ij (see also the footnote 3 ).

The Standard ADM formulation [80, 98]:
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∂t Kij= −2αKij + Di βj + Dj βi ,
= α

(2.5)

Rij + αKKij − 2αKik K j − Di Dj α

(3)k+(Di β k )Kkj + (Dj β k )Kki + β k Dk Kij(2.6)

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.7)
(2.8)

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

The formulation has 12 free ﬁrst-order dynamical variables (γij , Kij ), with 4 freedom of gauge choice(α, βi ) and with 4 constraint equations, (2.7) and (2.8). The rest freedom expresses 2 modes ofgravitational waves.
The constraint propagation equations, which are the time evolution equations of the Hamiltonian
constraint (2.7) and the momentum constraints (2.8), can be written as

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.9)∂t Mi = −(1/2)α(∂i H) − (∂i α)H + β (∂j Mi )j

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

(2.10)

Further expressions of these constraint propagations are in Appendix A.
From these equations, we know that if the constraints are satisﬁed on the initial slice Σ, then theconstraints are satisﬁed throughout evolution. The normal numerical scheme is to solve the elliptic
3

We remark that there is one replacement in (2.6) using (2.7) in the process of conversion from the original ADM to
the standard ADM equations. This is the key issue in the later discussion, and we shall be back this point in §3.2.

9

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.
2.0.3

Remarks

The ADM formulation was the standard formulation for numerical relativity up to the middle 90s.
Numerous successful simulations were obtained for the problems of gravitational collapse, criticalbehavior, cosmology, and so on. However, stability problems have arisen for the simulations suchas the gravitational radiation from compact binary coalescence, because the models require quite along-term time evolution.
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. From
the next subsection we shall look back of them by summarizing “three strategies”. We then unify
these three roads as “adjusted systems”, and as its by-product we show in §3.2 that the standard
ADM equations has a constraint violating mode in its constraint propagation equations even for asingle black-hole (Schwarzschild) spacetime [78].

2.1
2.1.1

Strategy 1: Modiﬁed ADM formulation by Nakamura et al (The BSSN formulation)
Introduction

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. [63, 64, 74]. The usefulness of this reformulation was re-introduced by Baumgarte and Shapiro [15], then was conﬁrmed by other groups to show
a long-term stable numerical evolution [5, 7]. The procedure is to apply conformally decompositionof the ADM variables and to implement their dynamical equations with several replacements.
2.1.2

Basic variables and equations

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

γij ,

(2.12)

K = γ Kij ,(2.13)

γ̃ij= e(2.11)

ij

Ãij
Γ̃

i

−4ϕ

= e
=

(Kij − (1/3)γij K),(2.14)

Γ̃ijk γ̃ jk .

(2.15)The new variable Γ̃i was introduced in order to calculate Ricci curvature more accurately. Γ̃i also
contributes to making system re-produce wave equations in its linear limit.
In BSSN formulation, Ricci curvature is not calculated as
ADM= ∂k Γkij − ∂i Γkkj + Γlij Γklk − Γlkj Γkli ,Rij(2.16)

but
BSSN
Rij

ϕ
= Rij
+ R̃ij ,(2.17)ϕRij= −2D̃i D̃j ϕ − 2γ̃ij D̃ D̃k ϕ + 4(D̃i ϕ)(D̃j ϕ) − 4γ̃ij (D̃ ϕ)(D̃k ϕ),(2.18)

R̃ij= −(1/2)γ̃ lk ∂l ∂k γ̃ij + γ̃k(i ∂j) Γ̃k + Γ̃k Γ̃(ij)k + 2γ̃ lm Γ̃kl(i Γ̃j)km + γ̃ lm Γ̃kim Γ̃klj ,(2.19)

k

k

10

BSSN
where D̃i is covariant derivative associated with γ̃ij . These are approximately equivalent, but Rijdoes have wave operator apparently in the ﬂat background 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,

(2.20)

during evolution. This is a kind of deﬁnition, but can also be treated as a constraint.

The BSSN formulation [63, 64, 74, 15]:
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.21)= −2αÃij + γ̃ik (∂j β ) + γ̃jk (∂i β ) − (2/3)γ̃ij (∂k β ) + β (∂k γ̃ij ), (2.22)
k

k

k

k

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

−4ϕ

(2.23)−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.24)−∂j (β k (∂k γ̃ ij ) − γ̃ kj (∂k β i ) − γ̃ ki (∂k β j ) + (2/3)γ̃ ij (∂k β k )).

(2.25)

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.26)(2.27)(2.28)(2.29)(2.30)(2.26) and (2.27) are the Hamiltonian and momentum constraints (the “kinematic” constraints), while the latter three are “algebraic” constraints due to the requirements of
BSSN formulation.

Hereafter we will write HBSSN and MBSSN simply as H and M respectively.
Taking careful account of these constraints, (2.26) and (2.27) can be expressed directly as
H = e−4ϕ R̃ − 8e−4ϕ D̃j D̃j ϕ − 8e−4ϕ (D̃j ϕ)(D̃j ϕ) + (2/3)K 2 − Ãij Ãij − (2/3)AK,(2.31)Mi = 6Ã i (D̃j ϕ) − 2A(D̃i ϕ) − (2/3)(D̃i K) + γ̃ (D̃j Ãki ).

(2.32)jkj

In summary, the fundamental dynamical variables in BSSN are (ϕ, γ̃ij , K,Ãij ,Γ̃i ), 17 in all. The
gauge quantities are (α, β i ) of which there are 4, and the constraints are (H, Mi , G i , A, S), i.e. 9components. As a result, 4 (2 by 2) components are left which correspond to two gravitationalpolarization modes.

112.1.3

Remarks

Why BSSN is better than the standard ADM? Together with numerical comparisons with the standard
ADM case[7, 57], this question has been studied by many groups using diﬀerent approaches. Using
numerical test evolution, Alcubierre et al [5] found that the essential improvement is in the process of
replacing terms by the momentum constraints. They also pointed out that the eigenvalues of BSSN
evolution equations have fewer “zero eigenvalues” than those of ADM, and they conjectured that the
instability might be caused by these “zero eigenvalues”. Miller [60] applied von Neumann’s stability
analysis to the plane wave propagation, and reported that BSSN has a wider range of parameters,
which produces stable evolution. Analogical conformal decomposition of the Maxwell equations arealso reported [53]. An eﬀort was made to understand the advantage of BSSN from the point ofhyperbolization of the equations in its linearized limit [5, 68]. 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 [40]. (Inspired by the BSSN’s conformal decomposition,
several related hyperbolic formulations have also been proposed [11, 42, 8].)
As we shall discuss in §3.3, 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.
Recently, Laguna and Shoemaker [55] modiﬁed the BSSN slightly, and found a great improvementin simulating a Schwarzschild black-hole.

The Laguna-Shoemaker version of BSSN [55]:Modiﬁcations to BSSN are

Box 2.4• to introduce conformal scalings also to K, Aij and N as
(LS)

K̂ = e6 n ϕ K,(LS)

Âi j = e6 n ϕ Ai j ,N = e−6 n ϕ α.

where γ̂ij = e−4ϕ γij and n is a parameter (n = 0 recovers the BSSN variables)• to use a mixed indices form (LS) Âi j rather than Aij .
• to use a densitized lapse N rather than α.

There is no explicit explanation why these small changes work better than before, but we expect that
our method can also be applied to ﬁnding the reason.

2.2
2.2.1

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. This is motivated by the expectation that the symmetric
hyperbolic system has well-posed properties in its Cauchy treatment in many systems and also thatthe boundary treatment can be improved if we know the characteristic speed of the system. As acomprehensive review of the hyperbolic formulation, we refer to those by Choquet-Bruhat and York(1980) [31], Geroch (1996) [43], Reula (1998) [67], and Friedrich-Rendall (2000) [37],We use the following deﬁnition:

12

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.33)

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.

We treat Mlβ α as a n × n matrix when the l-index is ﬁxed. The following properties of these matrices
apply to every basis of l-index. We say λl is an eigenvalue of Mlβ α when the characteristic equation,det(Mlβ α − λl δ β α ) = 0, is satisﬁed. The eigenvectors, pα , are given by solving Mlα β plβ = λl plα .
The strong hyperbolicity is identiﬁed, e.g. by judging whether the multiplicity of its eigenvalue, nλ ,satisﬁes rank(Mlβ α − λl δ β α ) = n − nλ for all λ. In order to deﬁne the symmetric hyperbolic system,
we need to declare an inner product u|u to determine whether Mlβ α is Hermitian. In other words,
we are required to deﬁne the way of lowering the index α of uα . We say Mlβ α is Hermitian withrespect to this index rule, when Ml βα = M̄l αβ for every l, where the overhead bar denotes complex
conjugate. To avoid this requirement of the deﬁnition of the inner product, people sometimes use
the word symmetrizable, when the characteristic matrix becomes Hermitian by a certain symmetrizer(positive deﬁnite matrix). In our classiﬁcation, this is only equivalent to the strongly hyperbolic
system. 4
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.34)

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.35)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 [81]), that is one step deﬁnitely advanced over a weakly hyperbolic form. Similarly, the4Several groups use a slightly diﬀerent deﬁnition for a symmetric hyperbolic system: deﬁning when the principal
symbol of the system iMlβ α kl is anti-Hermitian for arbitrary vector kl . The two deﬁnitions are equivalent when thevector kl is real-valued, but diﬀerent eigenvalues. In our deﬁnition, all eigenvalues are real-valued, while in the otherthey are all pure imaginary.

13

well-posedness of the symmetric hyperbolic is guaranteed if the characteristic matrix is independent ofu, while if it depends on u we have only limited proofs for the well-posedness. From the mathematicalpoint of view, proving well-posedness with less strict conditions is an old but still active research
problem. Therefore we have to recall that even if we construct a symmetric hyperbolic system ingeneral relativity, that equation does not necessarily guarantee numerical stability.
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 (e.g. [45]).
(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. For example, IBVP studies are preliminary yet, and most works are based on a particularsymmetric hyperbolic and in a limited space-time symmetry [26, 29, 38, 49, 81, 82, 83, 84]. We will
come back to this issue in §2.2.3 or §4, but meanwhile let us view the hyperbolic formulations fromthe comparisons of pure evolution equations.
We note that rich mathematical theories on partial diﬀerencial equations are obtained in a ﬁrstorder form, while there is a study on linearized ADM equations in a second-order form [54].
2.2.2Hyperbolic formulations of the Einstein equations

As was discussed by Geroch [43], 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 theEinstein equations is a long-standing research area in mathematical relativity. As we mentioned in
the introduction, numerical relativity shed light on this mathematical problem.
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.6) 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.
In the following discussion, we brieﬂy review several hyperbolic systems of the Einstein equations. We only mention the systems which applied numerical comparisons. See Table. 1 for a morecomprehensive list.
• The Bona-Massó formulation [18, 19, 20, 21]This was the ﬁrst active eﬀort to apply hyperbolized equations to numerical relativity. They
introduced auxiliary variables to reduce the system ﬁrst order in space: Ak = ∂k ln α, Bki =(1/2)∂k β i , and Dkij = (1/2)∂k gij , and construct a ﬁrst order ﬂux conservative system. In theirlatest formulation [21], the set of dynamical variables is (α, γij , Kij , Ai , Drij , Vi ≡ 2D[ir] r ), 37
14

functions in total, and the lapse function is restricted to a certain functional condition. The
system is a symmetrizable hyperbolic. They observed improved numerical stability compared tothat of the standard ADM system in spherical symmetric spacetime evolution [20]. An advantageof having characteristic speed is also applied to improve the treatment of the outer boundary
condition [22]. However, the appearance of shock formation was also reported unless the lapse
function is determined by solving the elliptic or parabolic equations [2, 4].
• The Einstein-Ricci system [32, 1] / Einstein-Bianchi system [9]Series of works by Choquet-Bruhat, York, and their colleagues developed hyperbolic systems inslightly diﬀerent ways. They intended to construct wave-type equations (2φ = S) for (γij , Kij ),that require the introduction of new variables (45 or 66 variables in total) and the use of Bianchi
identities. The resultant system is a third-order system in time for γij since they use the equation∂t Rij , but has only physical characteristic speed (zero or light speed). Scheel et al [70, 71] developed a numerial code with this formulation, and reported that they could evolve Schwarzschild
black hole (with 1-dimensional code) quite long time (∼ 104 M ), but futher modiﬁcations toequations were necesarrily.
• The Einstein-Christoﬀel system [10]Anderson and York [10] constructed a symmetrizable hyperbolic system which only has physicalcharacteristic speed. The variables are fundamentally (γij , Kij , Γkij ), and they also derived a setof equations with (γij , Kij , Gkij ≡ ∂k γij ), 30 functions. Their formulation diﬀers from the BonaMassó formulation in the way the momentum constraint is used to make the system hyperbolic.
Using a model of plane wave propagation, Bardeen and Buchman [17] numerically compared thisformulation with ADM and Bona-Massó with slight changes in variables. Their experiments areof plane-symmetric wave propagations, and they studies the boundary treatment in detail. We
will mention their work later.
• The Ashtekar formulation [13]Ashtekar’s reformulation of space-time was introduced to provide a new non-perturbative approach to quantum gravity. The new basic variables are the densitized inverse triad, Ẽai , andthe SO(3,C) self-dual connection, Aai , where the indices i, j, · · · indicate the 3-spacetime, anda, b, · · · are for SO(3) space. The formulation requires additional constraints and reality conditions in order to describe the classical Lorentzian space-time evolution, but the evolution
system itself forms a weakly hyperbolic system. By keeping the number of variables the same[(Ẽai , Aai )=total 18 (minimum ever)], we can construct both strongly and symmetric hyperbolic
systems by adjusting equations with constraints and/or restricting gauge conditions [92, 93].
The authors made numerical comparisons between the hyperbolicity of the systems with planesymmetric gravitational wave propagation using a periodic boundary condition [77, 94]. The
outcome is that there are no drastic diﬀerences in numerical stability between the three levelsof the hyperbolic equations. We will describe the details in Appendix B.
• The Frittelli-Reula formulation [41, 81]
This is a one-parameter family of a symmetric hyperbolic system. The procedure is to deﬁne
new variables, to adjust evolution equations with constraints (with one parameter), and todensitize the lapse function (with one parameter). The variables are (γij , M ijk ≡ 12 (γ ij ,k +aγ ij γrs γ rs,k ), P ij ≡ K ij + bγ ij K) where a, b are two other parameters, and make a total of 30.
They deﬁne a symmetric hyperbolic system by non-diverging energy norm, and restricted freeparameter spaces. A numerical comparison was also made for Gowdy spacetime [44], with quite
similar conclusions as in the Ashtekar version.
• The Conformal Field equations [36]A series of works by Friedrich [36] attempted to construct a 3 + 1 formulation with hyperboloidal
15

foliations (i.e. asymptotically null foliations), and with comformal compactiﬁcation. This is theultimate plan to remove the outer boundary problem in numerical simulation, and to provide asuitable foliation for gravitational radiation problem. However, the current equations are rather
quite complicated. In its metric-based expression [46], 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 4dimensional Ricci tensor R̂ab , the electric and magnetic components of the rescaled Weyl tensor
Cabc d , and the comformal factor Ω and 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 systemforms a symmetric hyperbolic system. Applications to numerical relativity are in progress, buthave not yet reached the stage of applying evolution in a non-trivial metric. For more details,
see reviews e.g. by Frauendiener [34] or by Husa [47].
• The Kidder-Scheel-Teukolsky (KST) formulation [52]This set of equations is a generalized formulation of the previous ones. It has 12 free parameters,
and includes both Anderson-York[10] and Frittelli-Reula[41] formulations as a subset. Therefore
it might be useful to discuss further hyperbolization methods starting with this KST formulation.
We brieﬂy summarize it in Box 2.6.

The Kidder-Scheel-Teukolsky (KST) formulation [52]:Box 2.6• Starting from a set of (gij , Kij , dkij ≡ ∂k gij ), the fundamental dynamical variables aredeﬁned (gij , Pij , Mkij ), as
Pij
Mkij

≡ Kij + ẑgij K,(2.36)ˆ j) )],≡ (1/2)[k̂dkij + êd(ij)k + gij (âdk + b̂bk ) + gk(i (ĉdj) + db

(2.37)ˆ ê, k̂, ẑ) are “kinematical” parameters.
where dk = g ab dkab , bk = g ab dabk , and (â, b̂, ĉ, d,• The 3-hypersurface Σ is foliated with gauge functions, (α, β i ), the lapse and shift vector.
However, the densitized lapse, Q = log(αg −σ ), (with a parameter σ) is actually used.
• The evolution equations are adjusted with constraints [in a version of (gij , Kij , dkij )]∂ˆ0 gij∂ˆ0 Kij= −2αKij ,(2.38)= (· · ·) + γαgij H + ζαg ab Ca(ij)b ,(2.39)∂ˆ0 dkij= (· · ·) + ηαgk(i Mj) + χαgij Mk ,(2.40)

where ∂ˆ0 = ∂t − £β and (γ, ζ, η, χ) are parameters. The terms (· · ·) are original terms
derived from the ADM equations and are available as (2.14) and (2.24) in [52].
• Constraints are (H, Mi , Cklij ), where Cklij ≡ ∂[k dl]ij .

In short, the number of dynamical variables and constraints are 30 and 22 and there are 12 freeparameters. Although there is no speciﬁc method to specify 12 free parameters KST showed numericalexamples of the Schwarzschild black hole (mass M ) evolution, which runs quite a longer time (t ∼6000M ) [52, 72].

16

We think the essential advantage of the KST system is the introduction of “kinematical” parameters. These 6 parameters• do not change the eigenvalues of evolution eqs.,• do not aﬀect the principal part of the constraint evolution eqs.,• do aﬀect the eigenvectors of the evolution system, and• do aﬀect the nonlinear terms of evolution eqs/constraint evolution eqs.
As Calabrese et al [27] pointed out, KST equations at linearized level on the ﬂat spacetime have nonon-principal terms, and these “kinematical” parameters ﬁnally enable us to discuss the features ofhyperbolicity in numerical experiments. Several variations of the KST formulation are presented in[69]. Recently, Lindblom and Scheel [58] tried to explain the relation between the numerical error
growing rate and the predicted error growing rate from the characteristic matrix of KST evolution
equations. Their trial did not succeed in matching the two completely, but at least began to reveal asimilar order before non-linear blow-up begins.
2.2.3

Remarks

When we discuss hyperbolic systems in the context of numerical stability, the following questions
should be considered:
Questions to hyperbolic formulations on its applications to numerics
Q(A) 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?
Q(B) If the answer to Q(A) is aﬃrmative, which level of hyperbolicity is practically useful fornumerical applications? Strongly hyperbolic, symmetric hyperbolic, or other?
Q(C) If the answer to Q(A) is negative, then can we ﬁnd other practical advantages to hyperbolization? Treatment of boundary conditions, or other?
Unfortunately, we do not have conclusive answers to these questions, but many experiences arebeing accumulated. Several earlier numerical comparisons reported the stability of hyperbolic formulations [20, 21, 22, 70, 71]. 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 beenshown recently (see §3.2).
These partial numerical successes encouraged the community to formulate various hyperbolic systems. Recently, Calabrese et al [28] 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.
• If the gauge functions are evolved according to the hyperbolic equations, then their ﬁnite
propagation speeds may cause pathological shock formations in simulations [2, 4].
• There are no drastic diﬀerences in the evolution properties between hyperbolic systems(weakly, strongly and symmetric hyperbolicity) by systematic numerical studies by Hern[44] based on Frittelli-Reula formulation [41], and by the authors [77] based on Ashtekar’s
formulation [13, 92, 93].
• Proposed symmetric hyperbolic systems were not always the best ones for numerical evolution. People are normally still required to reformulate them for suitable evolution. Such17

eﬀorts are seen in the applications of the Einstein-Ricci system [71], the Einstein-Christoﬀel
system [17], the conformal ﬁeld equations [47], and so on.
• Bardeen and Buchmann [17] conﬁrmed the usefulness of hyperbolicity when they treatednumerical boundary conditions in plane-symmetric wave propagation problem, but theyalso mentioned that the same techniques can not be applied in 2 or 3-dimensional cases.
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, as was already pointed out in several places [94, 56, 78].
The above numerical claims from Ashtekar and Frittelli-Reula formulations were mostly due to thecontribution (or interposition) of non-principal parts in evolution.
Regarding this issue, the recent KST formulation ﬁnally opens the door. As we saw, 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,
the agreement between numerical behavior and theoretical prediction is not yet perfect but close [58].
If further studies reveal the direct correspondences between theories and numerical results, then the
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 theeﬀects of the non-principal parts of the equations, which are neglected in the discussion of hyperbolicity.
Our proposal of adjusted system in the next subsection may be one of them.

2.3

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 (Figure 3). The idea was ﬁrst proposed as “λ-system” by Brodbeck
et al [24], and then developed in more general situations as “adjusted system” by the authors [94].
2.3.1

The “λ-system”

Brodbeck et al [24] proposed a system which has additional variables λ that obey artiﬁcial dissipativeequations. The variable λs are supposed to indicate the violation of constraints and the target of thesystem is to get λ = 0 as its attractor. Their proposal can be summarized as in Box 2.7.

18

The “λ-system” (Brodbeck-Frittelli-Hübner-Reula) [24]: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
λ



The main ideas are to introduce additional variables, λs, to impose dissipative dynamical equations for
λs, and to construct a symmetric hyperbolic system for both original variables and λs. 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 choiceof appropriate combinations of αs and βs.
Brodbeck et al presented a set of equations based on Frittelli-Reula’s symmetric hyperbolic formulation [41]. The version of Ashtekar’s variables was presented by the authors [76] for controllingthe constraints or reality conditions or both (see §B.2.2). The numerical tests of both the Maxwellλ-system and the Ashtekar-λ-system were performed [94], and conﬁrmed to work as expected (see§B.3.3). Although it is questionable whether the recovered solution is true evolution or not [79], 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

Blow up

error

t=0

?
Stabilize?

Constrained Surface(satisfies
Einstein's constraints)

time

Figure 3: Schematic picture of “asymptotically constrained system”.
19do not have any spatial diﬀerential terms.
2.3.2

The “adjusted system”

Next, we propose an alternative system which also tries to control the violation of constraint equations
actively, which we named “adjusted system”. We think that this system is more practical and robust
than the previous λ-system.

The Adjusted system (essentials) [94]:Box 2.8

Purpose:Control the violation of constraints by reformulating the system so as tohave a constrained surface as attractor.

Procedure:

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

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

Available even if the base system is not symmetric hyperbolic.

Advantages:

Keeps the number of the variables the same as in the original system.

We will describe the details in the next section, but summarize the procedures in advance:

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.
(See Box 3.2 and 3.3)

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

3.







The details are in §3.
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.
For the Maxwell equation and the Ashtekar version of the Einstein equations, we numericallyfound that this idea works to reduce the violation of constraints, and that the eﬀects are much better
than by constructing its symmetric hyperbolic versions [94] (see §B.3.4). The idea was applied to thestandard ADM formulation which is not hyperbolic and several attractive adjustments were proposed[95, 78] (see §3.2). This analysis was also applied to explain the advantages of the BSSN formulation,
and again several alternative adjustments to BSSN equations were proposed [96] (see §3.3). We will
explain these issues in the next section.
20

3

A uniﬁed treatment: Adjusted System

This section is devoted to present our idea of “asymptotically constrained system”, which was brieﬂy
introduced in Box 2.8 and 2.9 in the previous section. We begin with an overview of the adjustingprocedure and the idea of background structure in §3.1. Next, we show the applications both to ADM(§3.2) and BSSN (§3.3) formulations. Original references can be found in [94, 95, 78, 96].

3.1

Procedures : Constraint propagation equations and Proposals

Suppose we have a set of dynamical variables ua (xi , t), and their evolution equations,∂t ua = f (ua , ∂i ua , · · ·),(3.1)

C α (ua , ∂i ua , · · ·) ≈ 0.

(3.2)

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

for predicting the violation behavior of constraints in time evolution. We do not mean to integrate(3.3) numerically together with the original evolution equations (3.1), but mean to evaluate them
analytically in advance in order to reformulate the equations (3.1).
There may be two major analyses of (3.3); (a) the hyperbolicity of (3.3) when (3.3) is a ﬁrst ordersystem, and (b) the eigenvalue analysis of the whole RHS in (3.3) after a suitable homogenization.
(a) Hyperbolicity analysis of (3.3)
If (3.3) forms a ﬁrst-order system, the standard hyperbolic PDE analysis is applicable. As weviewed in §2.2, the analysis is mainly to identify the level of hyperbolicity and to calculate thecharacteristic speed of the system, from eigenvalues of the principal matrix.
For example, the evolution equations of the ADM formulations, (3.1) [(2.5) and (2.6)], do notform a ﬁrst-order system, while their constraint propagation equations, (3.3) [(2.9) and (2.10)],
do form. Therefore, we can apply the classiﬁcation on the hyperbolicity (weakly, strongly orsymmetric) to the constraint propagation equations. However, if one adjusts the ADM equations
with constraints, then this ﬁrst-order characters will not be guaranteed.
As we mentioned in §2.2.3, another big problem in the hyperbolic analysis is that it only discusses
the principal part of the system. If there is a method to characterize the non-principal part,
then that will help to clarify our understanding of evolution behavior.
(b) Eigenvalue analysis of the whole RHS in (3.3) after a suitable homogenizing procedure.
This analysis may compensate for the above problems.
Ampliﬁcation Factors of Constraint Propagation equations:
We propose to homogenize (3.3) 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(3.4)

then to analyze the set of eigenvalues, say Λs, of the coeﬃcient matrix, M α β , in (3.4). We
call Λs the constraint ampliﬁcation factors (CAFs) of (3.3).

21

The CAFs predict the evolution of constraint violations. We therefore can discuss the “distance” to the 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.

We found that the system becomes more stable when more Λs satisfy the above criteria. (Theﬁrst observations were in the Maxwell and Ashtekar formulations [77, 94]). Actually, supporting
mathematical proofs are available when we classify the fate of constraint propagations as follows.

Classiﬁcation of Constraint propagations:Box 3.3
If we assume that avoiding the divergence of constraint norm is related to the numericalstability, the next classiﬁcation would be useful:(C1) Asymptotically constrained : All the constraints decay and converge to zero.
This case can be obtained if and only if all the real parts of CAFs are negative.
(C2) Asymptotically bounded : All the constraints are bounded at a certain value. (thisincludes the above asymptotically constrained case.)
This case can be obtained if and only if (a) all the real parts of CAFs are not positive
and the constraint propagation matrix M α β is diagonalizable, or (b) all the real partsof CAFs are not positive and the real part of the degenerated CAFs is not zero.
(C3) Diverge : At least one constraint will diverge.
The details are shown in [97].

This classiﬁcation roughly indicates the heuristic statements (A) in Box 3.2. A practical procedure for this classiﬁcation is drawn in Figure 4.
We remark that this eigenvalue analysis requires the ﬁxing of a particular background spacetime,since the CAFs depend on the dynamical variables, ua .
The above features of the constraint propagation, (3.3), 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 α , · · ·),(3.5)

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

(3.6)

Therefore, the problem is how to adjust the evolution equations so that their constraint propagations
satisfy the above criteria as much as possible.
22Q1: Is there a CAF which real part is positive?
Diverge

NO / YES
Q2: Are all the real parts of CAFs negative?

AsymptoticallyConstrained

NO / YES

Q3: Is the constraint propagation matrix diagonalizable?
AsymptoticallyBoundedNO / YES

Q4: Is a real part of the degenerated CAFs is zero?
YES / NO

AsymptoticallyBounded

Q5: Is the associated Jordan matrix diagonal?
AsymptoticallyBoundedNO / YES

Diverge

Figure 4: A ﬂowchart to classify the constraint propagations (Box 3.3).

3.2

Adjusted ADM formulations

We show an application of the above idea, to the standard ADM system, Box 2.1.
3.2.1

Adjustments to ADM equations and its eﬀects on constraint propagations

Generally, we can write the adjustment terms to (2.5) and (2.6) using (2.7) and (2.8) by the following
combinations (using up to the ﬁrst derivatives of constraints for simplicity):

The adjusted ADM formulation [78]:Modify the evolution eqs of (γij , Kij ) by constraints H and Mi ,adjustment terms of ∂t γij :Box 3.4+Pij H + Qk ij Mk + pk ij (∇k H) + q kl ij (∇k Ml ),adjustment terms of ∂t Kij : +Rij H + S ij Mk + r ij (∇k H) + s
k

k

kl

(3.7)

ij (∇k Ml ),(3.8)

where P, Q, R, S and p, q, r, s are multipliers. That is,∂t γij= (2.5) + Pij H + Qk ij Mk + pk ij (∇k H) + q kl ij (∇k Ml ),(3.9)∂t Kij= (2.6) + Rij H + S ij Mk + r ij (∇k H) + s(3.10)kkkl

ij (∇k Ml ).

According to this adjustment, the constraint propagation equations are also modiﬁed as∂t H = (2.9) + H1mn (3.7) + H2imn ∂i (3.7) + H3ijmn ∂i ∂j (3.7) + H4mn (3.8),∂t Mi = (2.10) + M1i

mn(3.7) + M2i

jmn∂j (3.7) + M3imn(3.8) + M4i

jmn(3.11)∂j (3.8). (3.12)

with appropriate changes in indices. (See Appendix A. H1 , · · · , M1 , · · · are deﬁned there.)
We show two examples of adjustments here. We will list several others later in Table 2.
23

1. The standard ADM vs. original ADM
The ﬁrst comparison is to show the diﬀerences between the standard ADM [98] and the original
ADM system [12] (see §2.0). In the notation of (3.7) and (3.8), the adjustment,
Rij = κF αγij ,(and set the other multipliers zero)(3.13)

will distinguish two, where κF is a constant. Here κF = 0 corresponds to the standard ADM(no adjustment), and κF = −1/4 to the original ADM (without any adjustment to the canonicalformulation by ADM). As one can check by (3.11) and (3.12) adding Rij term keeps the constraint
propagation in a ﬁrst-order form. Frittelli [39] (see also [95]) pointed out that the hyperbolicity
of constraint propagation equations is better in the standard ADM system.
2. Detweiler type
Detweiler [33] found that with a particular combination, the evolution of the energy norm of theconstraints, H2 + M2 , can be negative deﬁnite when we apply the maximal slicing condition,
K = 0. His adjustment can be written in our notation in (3.7) and (3.8), asPij= −κL α3 γij ,(3.14)Rij= κL α3 (Kij − (1/3)Kγij ),(3.15)

k

S ij

=

kl

=

s

ij

k
κL α [3(∂(i α)δj)− (∂l α)γij γ kl ],k l
κL α3 [δ(iδj) − (1/3)γij γ kl ],2(3.16)(3.17)

everything else is zero, where κL is a constant. Detweiler’s adjustment, (3.14)-(3.17), does not
put constraint propagation equation to a ﬁrst order form, so we cannot discuss hyperbolicity
or the characteristic speed of the constraints. We conﬁrmed numerically, using perturbation on
Minkowskii and Schwarzschild spacetime, that Detweiler’s system provides better accuracy than
the standard ADM, but only for small positive κL .
3.2.2

Procedure to evaluate constraint ampliﬁcation factors in spherically symmetricspacetime

Before we compare between particular adjustment examples, we describe our procedure to evaluate
CAFs in spherically symmetric spacetime. According to our motivation, the actual procedure to
analyze the adjustments is to substitute the perturbed metric to the (adjusted) evolution equations ﬁrst
and to evaluate the according perturbative errors in the (adjusted) constraint propagation equations.
However, for the simplicity, we apply the perturbation to the pair of constraints directly and analyze
the eﬀects of adjustments in its propagation equations. The latter, we think, presents the feature ofconstraint propagation more clearly for our purposes.
The discussion becomes clear if we expand the constraint Cµ := (H, Mi )T using vector harmonics,
Alm alm + B lm blm + C lm clm + Dlm dlm ,Cµ =(3.18)

l,m

where we choose the basis asalm (θ, ϕ) = (Ylm , 0, 0, 0)T ,(3.19)

T(3.20)blm (θ, ϕ) = (0, Ylm , 0, 0) ,rclm (θ, ϕ) = (0, 0, ∂θ Ylm , ∂ϕ Ylm )T ,l(l + 1)r
1dlm (θ, ϕ) = (0, 0, −∂ϕ Ylm , sin θ ∂θ Ylm )T ,sin θl(l + 1)24

(3.21)(3.22)

(b)(a)

original ADM (κ = - 1/4)
F

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

0.5

0

-0.5

-1
0

5

10

r

15

0.5

0

-0.5
0

20

5

10

r

sch

15

20

sch

Figure 5: Ampliﬁcation factors (CAFs, eigenvalues of homogenized constraint propagation equations)
are shown for the standard Schwarzschild coordinate, with (a) no adjustments, i.e., standard ADM,
and (b) original ADM (κF = −1/4). [see eq. (3.13)]. The solid lines and the dotted lines withcircles are real parts and imaginary parts, respectively. They are four lines each, but actually the two
eigenvalues are zero for all cases. Plotting range is 2 < r ≤ 20 using Schwarzschild radial coordinate.
c
We set k = 1, l = 2, and m = 2 throughout the article. (Reprinted from [78], APS
2002)
and the coeﬃcients Alm , · · · , Dlm are functions of (t, r). Here Ylm is the spherical harmonic function,
Ylm (θ, ϕ) = (−1)(m+|m|)/2(2l + 1) (l − |m|)! m
P (cos θ)eimϕ .
4π (l + |m|)! l(3.23)The basis (3.19)-(3.22) are normalized so that they satisfy
Cµ , Cν  =

 2π

 π

dϕ
0

0

Cµ∗ Cρ η µρ sin θdθ,(3.24)

where η µρ is Minkowskii metric and the asterisk denotes the complex conjugate. Therefore
Alm = alm(ν) , Cν ,∂t Alm = alm(ν) , ∂t Cν ,

etc.

(3.25)In order to analyze the radial dependences, we also express these evolution equations using theFourier expansion on the radial coordinate,
Alm =

ikr
Âlm
(k) (t) e

etc.

(3.26)klm
T
So that we can obtain the RHS of the evolution equations for (Âlm(k) (t), · · · , D̂(k) (t)) in a homogeneous
form.

3.2.3

Constraint ampliﬁcation factors in Schwarzschild spacetime

We present our CAF analysis in Schwarzschild black hole spacetime, which metric is
ds2 = −(1 −

2M 2
dr2
)dt +
+ r2 dΩ2 ,r
1 − 2M/r(the standard expression)(3.27)

where M is the mass of a black hole. For numerical relativists, evolving a single black hole is theessential test problem, though it is a trivial at ﬁrst sight. The standard expression, (3.27), has a
25

(a)

(b)

Detweiler type, κ = + 1/2

Detweiler type, κ = - 1/2
L1

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

L

0.5

0

-0.5

-1
0

5

10

r

15

1

0.5

0

-0.5

20

-1
0

5

10

r

sch

15

20

sch

Figure 6: Ampliﬁcation factors of the standard Schwarzschild coordinate, with Detweiler type adjustments, (3.14)-(3.17). Multipliers used in the plot are (a) κL = +1/2, and (b) κL = −1/2. Plottingcdetails are the same as Fig.5. (Reprinted from [78], APS
2002)
coordinate singularity at r = 2M , so that we need to move another coordinate for actual numericaltime integrations.
The ingoing Eddington-Finkelstein (iEF) coordinate has become popular in numerical relativity, inorder to excise black hole singularity, since iEF penetrates the horizon without an irregular coordinate.
The expression is,
ds2 = −(1 −

2M 2
4M
2M
)dtiEF +dtiEF dr + (1 +)dr2 + r2 dΩ2rrr(the iEF expression)(3.28)

which is given by tiEF = t + 2M log(r − 2M ) and the radial coordinate is common to (3.27).
We show CAFs in the following cases. More examples are available in [78].
1. in the standard Schwarzschild metric expression (3.27)• For the adjustment (3.13), CAFs are obtained as√√
Λi = (0, 0, a, − a),(3.29)
2
2
2
4M k r (r − M ) + 2M (2r − M ) + l(l + 1)r(r − 2M ) + ikr(2r − 3M r − 2M 2 )a = −k 2 +r4
for the choice of κ1 = 0 (the standard ADM), while they are√√Λi = (0, 0, b, − b),b=

M (2r − M ) + irkM (2M − r)
r4

(3.30)

for the choice of κ1 = −1/4 (the original ADM).
These are plotted in Fig.5. The solid lines and dotted lines with circles are real parts andimaginary parts of CAFs, respectively. They are four lines each, but as we showed in (3.29),
two of them are zero. The 
plotting range is 2 < r ≤ 20 in Schwarzschild radial coordinate.
The CAFs at r = 2 are ± 3/8 and 0. The existence of this positive real CAF near thehorizon is an important result. We show only the cases with l = 2 and k = 1, because we
judged that the plots of l = 0 and other ks are qualitatively the same.
The adjustment (3.13) with κF = −1/4 returns the system back to the original ADM.
CAFs are (3.30) and we plot them in Fig.5(b). We can see that the imaginary parts areapparently diﬀerent from those of the standard ADM [Fig.5(a)]. This is the same feature as
26

(a)

(b)

iEF coordinate, no adjustments (standard ADM)

iEF coordinate, Detweiler type κ =+0.5
LReal / Imaginary parts of Eigenvalues (AF)Real / Imaginary parts of Eigenvalues

2

1.5

1

0.5

0

-0.5

-1
0

5

10

r

15

2

1.5

1

0.5

0

-0.5

-1
0

20

5

10

r

sch

15

20

sch

Figure 7: CAFs in the iEF coordinate (3.28) on t = 0 slice for the standard ADM formulation (i.e. no
adjustments) [See Fig.5(a) for the standard Schwarzschild coordinate.] and for Detweiler adjustmentswith κL = +1/2 [See Fig.6(a) for the standard Schwarzschild coordinate.]. The solid four lines andthe dotted four lines with circles are real parts and imaginary parts, respectively. (Reprinted fromc[78], APS
2002)
in the case of the ﬂat background [95]. According to our conjecture, the non-zero imaginary
values are better than zeros, so we expect that the standard ADM has a better evolution
property than the original ADM system. Negative κF always makes the asymptotical realvalues ﬁnite.
• The Detweiler-type adjustment (3.14)-(3.17) makes the feature completely diﬀerent.
Fig.6(a) and (b) are the cases of κL = ±1/2. A great improvement can be seen in thepositive κL case where all real parts become negative in large r. Moreover all imaginaryparts are apart from zero. These are the desired features according to our conjecture.
Therefore we expect the Detweiler adjustment has good stability properties except near theblack hole region. The CAF near the horizon has a positive real component. This is notcontradictory with the Detweiler’s original idea. His idea came from suppressing the total
L2 norm of constraints on the spatial slice, while our plot indicates the existence of a local
violation mode. The change of signature of κL can be understood just by changing thesignature of CAFs, and this fact can also be seen to the other plot. In [78] we reported
that a partial adjustment (apply only (3.14) or (3.17)) is also eﬀective.
2. in the iEF coordinates (3.28)• For the adjustment (3.13), we plotted CAFs in Figure 7. We Figure 7(a) is qualitatively
diﬀerent from Fig.5(a). This is because the iEF expression is asymmetric to time, i.e. has
non-zero extrinsic curvature. We notice that while some CAFs in iEF remain positive inlarge r region, that their nature changes due to the adjustments.
• For the Detweiler-type adjustment (3.14)-(3.17), CAFs are as in Figure 7(b). Interestingly, all plots indicate that all real parts of CAFs are negative, and imaginary parts are
non-zero (again except near the black hole region). By arranging the multiplier parameter,
there is a chance to get all negative real CAFs outside the black hole horizon. For example,
all the real-part goes negative outside the black hole horizon if κL > 3.1, while large κL
may introduce another instability problem [94].
Such kinds of test can be done with other combinations. In Table 2, we listed our results for more
examples. We deﬁned the adjustment terms so that their positive multiplier parameter, κ > 0, makes
27

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

adjustment
–
Pij
Pij
Pij
Pij
Pij
Qk ij
Qk ij
Qk ij
Qk ij
Rij
Rij
Rij
S k ijS k ij
pk ij
pk ij
pk ij
q kl ij
q kl ijr k ijr k ijr k ijskl ijskl ijskl ijno adjustments−κL α3 γij−κL αγij
Prr = −κ
−κγij
−κγrr
καβ k γij
Qr rr = κQr ij = κγij
Qr rr = κγrr
κF αγij
Rrr = −κµ α
Rrr = −κγrr
S k ij = (3.16)
καγ lk (∂l γij )pr ij = −καγij
pr rr = καpr rr = καγrr
q rr ij = καγij
q rr rr = −καγrr
r r ij = καγij
r r rr = −καr r rr = −καγrr
skl ij = (3.17)
srr ij = −καγij
srr rr = −καγrr

1st?
yes
no
no
no
no
no
no
no
no
no
yes
yes
yes
yes
yes
no
no
no
no
no
no
no
no
no
no
no

TRS
–
no
no
no
no
no
no
yes
yes
yes
yes
no
no
no
no
no
no
no
no
yes
yes
yes
yes
no
no
no

Sch coord.
real.
imag.
–
–
makes 2 Neg.
not app.
makes 2 Neg.
not app.
slightly enl.Neg.
not app.
makes 2 Neg.
not app.
red. Pos./enl.Neg.
not app.
N/A
N/A
red. abs vals.
not app.
red. abs vals.
not app.
red. abs vals.
not app.
κF = −1/4 min. abs vals.
not app.
not app.
enl. vals.
not app.
not app.
not app.
makes 2 Neg.
not app.
red. Pos.
red. vals.
red. Pos.
red. vals.
makes 2 Neg.
enl.vals.
κ = 1/2 min.vals.
red.vals.
red. abs vals.
not app.
not app.
not app.
red. abs vals.
enl. vals.
red. abs vals.
enl. vals.
makes 4 Neg.
not app.
makes 2 Neg.
red. vals.
makes 2 Neg.
red. vals.

iEF coord.
real.
imag.
–
–
makes 2 Neg.
not app.
makes 2 Neg.
not app.
slightly enl.Neg.
not app.
makes 2 Neg.
not app.
red.Pos./enl.Neg.
not app.
κ ∼ 1.35 min.vals.
not app.
red. abs vals.
not app.
enl.Neg.
enl. vals.
red. abs vals.
not app.
κF = −1/4 min. vals.
red.Pos./enl.Neg.
enl. vals.
red.Pos./enl.Neg.
enl. vals.
not app.
not app.
makes 2 Neg.
not app.
red. Pos.
enl. vals.
red.Pos/enl.Neg.
enl. vals.
red.Pos.vals.
red.vals.
not app.
enl.vals.
not app.
not app.
not app.
enl. vals.
red. abs vals.
enl. vals.
red. abs vals.
enl. vals.
makes 4 Neg.
not app.
makes 2 Neg.
red. vals.
makes 2 Neg.
red. vals.

Table 2: List of ADM adjustments we tested in the Schwarzschild spacetime. The column of adjustments are nonzero multipliers in terms of (3.7) and (3.8). The column ‘TRS’ indicates whether eachadjusting term satisﬁes the time reversal symmetry or not on the standard Schwarzschild coordinate.
(‘No’ means a candidate that makes asymmetric CAFs.) The column ‘1st?’ indicates whether eachadjusting term breaks the ﬁrst-order feature of the standard constraint propagation equation, (2.9)
and (2.10). (‘Yes’ keeps the system ﬁrst-order, ‘No’ may break hyperbolicity of constraint propagation
depending a choice of κ.) The eﬀects to CAFs (when κ > 0) are commented for both coordinates andboth real/imaginary parts, respectively. The ‘N/A’ means that there is no eﬀect due to the coordinate
properties; ‘not app.’ (not apparent) means the adjustment does not change the CAFs eﬀectively
according to our conjecture; ‘enl./red./min.’ means enlarge/reduce/minimize, and ‘Pos./Neg.’ means
positive/negative, respectively. These judgements are made at the r ∼ O(10M ) region on their t = 0
slice. See more detail in [78].
the system better in stability according to our conjecture. (Here better means in accordance with ourconjecture in Box 3.2 and 3.3). The table includes the above results and is intended to extract thecontributions of each term in (3.7) and (3.8). The eﬀects of adjustments (of each κ > 0 case) to CAFs
are commented upon for each coordinate system and for real/imaginary parts of CAFs, respectively.
These judgements are made at the r ∼ O(10M ) region on their t = 0 slice. Among them, No. R-2
in Table 2 explains why a particular adjustment by PennState group [51] gives better stability than
before.
3.2.4

Remarks

Numerical demonstrations The above analyses are only predictions, and supporting numericaldemonstrations are necessary for the next steps. Systematic numerical comparisons are progressing,
and we show two sample plots here.
Figure 8 (a) is a test numerical evolution of Detweiler-type adjustment on the Minkowskii background. We see the adjusted version gives convergence on to the constraint surface by arranging the28

ADM and its adjusted versions-- Teukolsky wave evolution --Detweiler's adjustments
on Minkowskii spacetime
2.0

0.0

L= + 0.035

L2 norm of Hamiltonian constraint

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-2.0

L=0 (standard ADM)

-4.0

-6.0

L= + 0.01-8.0

L= + 0.03

010-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

1.5

2.0

0

50

100

150

200

250

300

time

time

Figure 8: 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 used in (3.14)-(3.16). We see the evolution iscasymptotically constrained for small κ > 0. (Reprinted from [95], APS2001) (b) L2 norm of the
ADM
Hamiltonian constraint Hof evolution using ADM/adjusted ADM formulations for the case of
Teukolsky wave, 3-dimensional simulation. Cactus-based original ADM module (CactusGR) was used.
(Shinkai-Yoneda, in preparation).
magnitude of the adjusting parameter, κ.
Figure 8 (b) is an example from the project of “Comparison of formulations of the Einstein equations for numerical relativity” [59]. (We will describe the project in §4). The plot was obtained by a
3-dimensional numerical evolution of weak gravitational wave, the so-called Teukolsky wave [85]. Thelines are of the original/standard ADM evolution equations, Detweiler-type adjustment, and a part of
Detweiler-type adjustment (actually used only (3.14)). For a particular choice of κ, we observe again
the L2 norm of constraint (violation of constraints) is reduced than the standard ADM case, and canevolve 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.

29

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. The TRS features (onthe standard Schwarzschild metric, which is time-symmetric background) are also denoted in Table 2.
The criteria of Box 3.5 will be applied also to the BSSN equation (see §3.3.3).
Diﬀerences with Detweiler’s requirement We comment on the diﬀerences between Detweiler’s
criteria [33] and ours. Detweiler calculated the L2 norm of the constraints, Cρ , over the 3-hypersurface
and imposed the negative deﬁniteness of its evolution,
Detweiler’s criteria ⇔ ∂tCρ C ρ dV < 0, ∀ non zero Cρ .

(3.31)

where Cρ C ρ =: Gρσ Cρ Cσ , and Gρσ = diag[1, γij ] for the pair of Cρ = (H, Mi ).
Assuming the constraint propagation to be ∂t Ĉρ = Aρ σ Ĉσ in the Fourier components, the time
derivative of the L2 norm can be written as¯∂t (Ĉρ Ĉ ρ ) = (Aρσ + Āσρ + ∂t Ḡρσ )Ĉρ Ĉ σ .

(3.32)

Together with the fact that the L2 norm is preserved by Fourier transform, we can say, for the caseof static background metric,
Detweiler’s criteria ⇔ eigenvalues of (A + A† ) are all negative ∀k.

(3.33)

Our criteria ⇔ eigenvalues of A are all negative ∀k.

(3.34)

Therefore for the case of static background, Detweiler’s criterion is stronger than ours. For example,
the matrix
−1 awhere a is constant,(3.35)
A=
0 −1√
for the evolution system (Ĉ1 , Ĉ2 ) satisﬁes our criterion but not Detweiler’s when |a| ≥ 2. Thismatrix however gives asymptotical decay for (Ĉ1 , Ĉ2 ). Therefore we may say that Detweiler requires
the monotonic decay of the constraints, while we assume only asymptotical decay.
We remark that Detweiler’s truncations on higher order terms in C-norm corresponds to our
perturbational analysis; both are based on the idea that the deviations from constraint surface (theerrors expressed non-zero constraint value) are initially small.
Ranges of eﬀective κ We do not discuss the ranges of the eﬀective multiplier parameter, κ, since
the range depends on the characteristic speeds of the models and numerical integration schemes as we
observed in [94]. At this moment, we can only estimate the maximum value for adjusting parameterκs from von Neumann’s stability analysis, while we do not have a background theoretical explanation
to predict the optimized value of κ for numerical evolution. We will comment this point later againin §4.

3.3

Adjusted BSSN formulations

Next we apply the idea of adjusted system to the BSSN system, Box 2.3.
3.3.1

Constraint propagation analysis of the BSSN equations

The BSSN has 5 constraint equations, (H, Mi , G i , A, S), see (2.26)-(2.30). We begin identifying where
the BSSN evolution equations were adjusted already in its standard notation, (2.21)-(2.25).
30

Taking careful account of these constraints, (2.26) and (2.27) can be expressed directly as
H = e−4ϕ R̃ − 8e−4ϕ D̃j D̃j ϕ − 8e−4ϕ (D̃j ϕ)(D̃j ϕ) + (2/3)K 2 − Ãij Ãij − (2/3)AK,(3.36)Mi = 6Ã i (D̃j ϕ) − 2A(D̃i ϕ) − (2/3)(D̃i K) + γ̃ (D̃j Ãki ).

(3.37)jkj

By a straightforward calculation, we get:∂tB ϕ = ∂tA ϕ + (1/6)αA − (1/12)γ̃ −1 (∂j S)β j ,∂tB γ̃ij∂tB K∂tB Ãij===(3.38)∂tA γ̃ij − (2/3)αγ̃ij A + (1/3)γ̃ −1 (∂k S)β k γ̃ij ,∂tA K − (2/3)αKA − αH + αe−4ϕ (D̃j G j ),∂tA Ãij + ((1/3)αγ̃ij K − (2/3)αÃij )A+((1/2)αe−4ϕ (∂k γ̃ij ) − (1/6)αe−4ϕ γ̃ij γ̃ −1 (∂k S))G k+αe−4ϕ γ̃k(i (∂j) G k ) − (1/3)αe−4ϕ γ̃ij (∂k G k ),(3.39)(3.40)(3.41)∂tB Γ̃i = ∂tA Γ̃i + ( − (2/3)(∂j α)γ̃ ji − (2/3)α(∂j γ̃ ji ) − (1/3)αγ̃ ji γ̃ −1 (∂j S) + 4αγ̃ ij (∂j ϕ))A−(2/3)αγ̃ ji (∂j A) + 2αγ̃ ij Mj − (1/2)(∂k β i )γ̃ kj γ̃ −1 (∂j S) + (1/6)(∂j β k )γ̃ ij γ̃ −1 (∂k S)+(1/3)(∂k β k )γ̃ ij γ̃ −1 (∂j S) + (5/6)β k γ̃ −2 γ̃ ij (∂k S)(∂j S) + (1/2)β k γ̃ −1 (∂k γ̃ ij )(∂j S)+(1/3)β k γ̃ −1 (∂j γ̃ ji )(∂k S).

(3.42)

where ∂tA denotes the part of no replacements, i.e. the terms only use the standard ADM evolution
equations in its time derivatives.
From (3.38)-(3.42), we understand that all BSSN evolution equations are adjusted using constraints.
However, from the viewpoint of time reversal symmetry (Box 3.5), all the above adjustments in (3.38)(3.42) unfortunately keep the time reversal symmetry. Therefore we can not expect direct decays of
constraint violation in the present form.
The set of the constraint propagation equations, ∂t (H, Mi , G i , A, S)T , turns to be not a ﬁrst-order
hyperbolic and includes many non-linear terms, (see the Appendix in [96]). In order to understandthe fundamental structure, we show an analysis on the ﬂat spacetime background.
For the ﬂat background metric gµν = ηµν , the ﬁrst order perturbation equations of (3.38)-(3.42)
can be written as∂t(1)ϕ = −(1/6)(1)K + (1/6)(κϕ − 1)(1)A,∂t(1)γ̃ij∂t(1)K∂t(1)Ãij===(3.43)−2 Ãij − (2/3)(κγ̃ − 1)δij(1)A,(3.44)
(1)
(1) j
(1)
−(∂j ∂j α) + (κK1 − 1)∂j G − (κK2 − 1) H,(3.45)(1) BSSN T F
(1)
TF
(1) k
(1) k
(Rij
) − (D̃i D̃j α) + (κA1 − 1)δk(i (∂j) G ) − (1/3)(κA2 − 1)δij (∂k G (3.46)),(1)∂t(1)Γ̃i = −(4/3)(∂i(1)K) − (2/3)(κΓ̃1 − 1)(∂i(1)A) + 2(κΓ̃2 − 1)(1)Mi ,(3.47)

where we introduced parameters κs, all κ = 0 reproduce no adjustment case from the standard ADM
equations, and all κ = 1 correspond to the BSSN equations. We express them asκadj:= (κϕ , κγ̃ , κK1 , κK2 , κA1 , κA2 , κΓ̃1 , κΓ̃2 ).

(3.48)

Constraint propagation equations at the ﬁrst order in the ﬂat spacetime, then, become:∂t(1)H =κγ̃ − (2/3)κΓ̃1 − (4/3)κϕ + 2 ∂j ∂j(1)A + 2(κΓ̃2 − 1)(∂j(1)Mj ),∂t(1)Mi(−(2/3)κK1 + (1/2)κA1 − (1/3)κA2 + (1/2)) ∂i ∂j(1)G j+(1/2)κA1 ∂j ∂j(1)G i + ((2/3)κK2 − (1/2)) ∂i(1)H,2κΓ̃2(1)Mi + (−(2/3)κΓ̃1 − (1/3)κγ̃ )(∂i(1)A),−2κγ̃(1)A,(κA1 − κA2 )(∂j(1)G j ).

∂t(1)G i∂t(1)S∂t(1)A

=
=
=
=

31

(3.49)(3.50)(3.51)(3.52)(3.53)

3.3.2The origin of the advantages of the BSSN equations

We next discuss CAFs of (3.49)-(3.53). Hereafter we let k 2 = kx2 + ky2 + kz2 for Fourier wave numbers.
1. The no-adjustment case, κadj =(all zeros). This is the starting point of the discussion. In thiscase,CAFs = (0 (×7), ± −k 2 ),i.e., (0 (×7), ±pure imaginary (1 pair)). In the standard ADM formulation, which uses (γij , Kij ),CAFs are (0, 0, ±Pure Imaginary) [95]. Therefore if we do not apply adjustments in BSSN
equations the constraint propagation structure is quite similar to that of the standard ADM.
2. For the BSSN equations, κadj =(all 1s),

CAFs = (0 (×3), ± −k 2 (3 pairs)),i.e., (0 (×3), ±Pure Imaginary (3 pairs)). The number of pure imaginary CAFs is increased overthat of No.1, and we conclude this is the advantage of adjustments used in BSSN equations.
3. No S-adjustment case. All the numerical experiments so far apply the scaling condition S forthe conformal factor ϕ. The S-originated terms appear many places in BSSN equations (2.21)(2.25), so that we guess non-zero S is a kind of source of the constraint violation. However, since
all S-originated terms do not appear in the ﬂat spacetime background analysis, [no adjusted
terms in (3.43)-(3.47)], our case can not say any eﬀects due to S-constraint.
4. No A-adjustment case. The trace (or traceout) condition for the variables is also considerednecessary (e.g. [3]). This can be checked with κadj = (κ, κ, 1, 1, 1, 1, κ, 1), and we getCAFs = (0 (×3), ± −k 2 (3 pairs)),
independent of κ. Therefore the eﬀect of A-adjustment is not apparent from this analysis.
5. No G i -adjustment case. The introduction of Γi is the key in the BSSN system. If we do not
apply adjustments by G i , (κadj = (1, 1, 0, 1, 0, 0, 1, 1)) then we getCAFs = (0(×7), ± −k 2 ),which is the same with No.1. That is, adjustments due to G i terms are eﬀective to make aprogress from ADM.
6. No Mi -adjustment case. This can be checked with κadj = (1, 1, 1, 1, 1, 1, 1, κ), and we getCAFs = (0, ± −κk 2 (2 pairs),± −k 2 (−1 + 4κ + |1 − 4κ|)/6,

± −k 2 (−1 + 4κ − |1 − 4κ|)/6).



If κ = 0, then (0(×7), ± k 2 /3), which is (0(×7), ±real value). Interestingly, these real values
indicate the existence of the error growing mode together with the decaying mode. Alcubierre
et al. [5] found that the adjustment due to the momentum constraint is crucial for obtaining
stability. We think that they picked up this error growing mode. Fortunately at the BSSN limit(κ = 1), this error growing mode disappears and turns into a propagation mode.
7. No H-adjustment case. The set κadj = (1, 1, 1, κ, 1, 1, 1, 1) gives
CAFs = (0 (×3), ± −k 2 (3 pairs)),
independently to κ. Therefore the eﬀect of H-adjustment is not apparent from this analysis.
32

No.

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

S (1)

diag?

CAFs
in Minkowskii background

0.
1.
2.

standard ADM
BSSN no adjustment
the BSSNuse
use
use+adjuse
use
use+adjuse
use+adjuse
use+adjuse
use+adj

yes
yes
no

(0, 0, , )(0, 0, 0, 0, 0, 0, 0, , )(0, 0, 0, , , , , , )

3.
4.
5.
6.
7.

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

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

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

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

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

no
no
no
no
no

no diﬀerence in ﬂat background(0, 0, 0, , , , , , )(0, 0, 0, 0, 0, 0, 0, , )(0, 0, 0, 0, 0, 0, 0, , )(0, 0, 0, , , , , , )

Table 3: Contributions of adjustments terms and eﬀects of introductions of new constraints in the
BSSN system. The center column indicates whether each constraints are taken as a component ofconstraints in each constraint propagation analysis (‘use’), and whether each adjustments are on(‘adj’). The column ‘diag?’ indicates diagonalizability of the constraint propagation matrix. The
right column shows ampliﬁcation factors, where  and  means pure imaginary and real eigenvalue,
respectively.
These tests are on the eﬀects of adjustments that are already in the BSSN equations. We will consider
whether much better adjustments are possible in the next section.
We list the above results in Table 3. The most characteristic points of the above are No. 5 andNo.6 that denote the contribution of the momentum constraint adjustment and the importance ofthe new variable Γ̃i . It is quite interesting that the unadjusted BSSN equations (case 2) does nothave apparent advantages from the ADM system. As we showed in the case 5 and 6, if we missed aparticular adjustment, then the expected stability behaviour occationally gets worse than the startingADM system. Therefore we conclude that the better stability of the BSSN formulation is obtained bytheir adjustments in the equations, and the combination of the adjustments is in a good balance.
3.3.3

Proposals of the modiﬁed BSSN equations

We next consider the possibility whether we can obtain a system which has much better properties;
whether more pure imaginary CAFs or negative real CAFs.
Heuristic examples(A) A system which has 8 pure imaginary CAFs:
One direction is to seek a set of equations which make fewer zero CAFs than the standard BSSN
case. Using the same set of adjustments in (3.43)-(3.47), CAFs are written asCAFs =

0, ± −k 2 κA1 κΓ̃2 (2 pairs), ±complicated expression, ±complicated expression .

The terms in the ﬁrst line certainly give four pure imaginary CAFs (two positive and negative
real pairs) if κA1 κΓ̃2 > 0 (< 0). Keeping this in mind, by choosing κadj = (1, 1, 1, 1, 1, κ, 1, 1), we
ﬁnd
CAFs =



0, ± −k 2 (2 pairs), ± −k 2 (2 + κ + |κ − 4|)/6, ± −k 2 (2 + κ − |κ − 4|)/6, .

√Therefore the adjustment κadj = (1, 1, 1, 1, 1, 4, 1, 1) gives CAFs = 0, ± −k 2 (4 pairs) , which
is one step advanced from BSSN according our guidelines.
We note that such a system can be obtained in many ways, e.g. κadj = (0, 0, 1, 0, 2, 1, 0, 1/2) also
gives four pairs of pure imaginary CAFs.
33

(B) A system which has negative real CAF:
One criterion to obtain a decaying constraint mode (i.e. an asymptotically constrained system)
is to adjust an evolution equation as it breaks time reversal symmetry (Box 3.5). For example,
we consider an additional adjustment to the BSSN equation as∂t γ̃ij = ∂tB γ̃ij + κSD αγ̃ij H,(3.54)

which is a similar adjustment of the simpliﬁed Detweiler-type [33]. The constaint ampliﬁcation
factors becomeCAFs = (0 (×2), ± −k 2 (3 pairs), (3/2)k 2 κSD ),in which the last one becomes negative real if κSD < 0.
(C) Combination of above (A) and (B)Naturally we next consider both adjustments:∂t γ̃ij∂t Ãij= ∂tB γ̃ij + κSD αγ̃ij H=(3.55)∂tB Ãij − κ8 αe−4ϕ γ̃ij ∂k G k(3.56)

where the second one produces the 8 pure imaginary CAFs. We then obtain
CAFs =

0, ± −k 2 (3 pairs), (3/4)k 2 κSD ±

k 2 (−κ8 + (9/16)k 2 κSD )which reproduces case (A) when κSD = 0, κ8 = 1, and case (B) when κ8 = 0. These CAFs canbecome (0, pure imaginary (3 pairs), complex numbers with a negative real part (1 pair)), with
an appropriate combination of κ8 and κSD .
Possible adjustments In order to break time reversal symmetry of the evolution equations (Box3.5), the possible simple adjustments are (1) to add H, S or G i terms to the equations of ∂t φ, ∂t γ̃ij ,or ∂t Γ̃i , and/or (2) to add Mi or A terms to ∂t K or ∂t Ãij . We write them generally, including theabove proposal (B), as∂t φ = ∂tB φ + κφH αH + κφG αD̃k G k ,∂t γ̃ij=(3.58)=+ κγ̃S1 αγ̃ij S + κγ̃S2 αD̃i D̃j S,
B∂t K + κKM αγ̃ jk (D̃j Mk ),∂tB Ãij + κAM1 αγ̃ij (D̃k Mk ) + κAM2 α(D̃(i Mj) )(3.60)=+ κAA1 αγ̃ij A + κAA2 αD̃i D̃j A,
B i∂t Γ̃ + κΓ̃H αD̃i H + κΓ̃G1 αG i + κΓ̃G2 αD̃j D̃j G i + κΓ̃G3 αD̃i D̃j G j ,∂t K =∂t Ãij∂t Γ̃

i(3.57)∂tB γ̃ij + κSD αγ̃ij H + κγ̃G1 αγ̃ij D̃k G k + κγ̃G2 αγ̃k(i D̃j) G k(3.59)(3.61)

where κs are possible multipliers (all κ = 0 reduce the system the standard BSSN system).
We show the eﬀects of each terms in Table 4. The CAFs in the table are on the ﬂat spacebackground. We see several terms make negative CAFs, which might improve the stability than theprevious system. For the readers convenience, we list up several best candidates here.
(D) A system which has 7 negative CAFs
Simply adding D̃(i Mj) term to ∂t Ãij equation, say∂t Ãij = ∂tBSSN Ãij + κAM2 α(D̃(i Mj) )with κAM2 > 0, CAFs on the ﬂat background are 7 negative real CAFs.

34(3.62)

∂t φ∂t φ∂t γ̃ij∂t γ̃ij∂t γ̃ij∂t γ̃ij∂t γ̃ij∂t K∂t Ãij∂t Ãij∂t Ãij∂t Ãij∂t Γ̃i∂t Γ̃i∂t Γ̃i∂t Γ̃iadjustment
κφH αHκφG αD̃k G k
κSD αγ̃ij H
κγ̃G1 αγ̃ij D̃k G k
κγ̃G2 αγ̃k(i D̃j) G k
κγ̃S1 αγ̃ij S
κγ̃S2 αD̃i D̃j S
κKM αγ̃ jk (D̃j Mk )
κAM1 αγ̃ij (D̃k Mk )
κAM2 α(D̃(i Mj) )κAA1 αγ̃ij A
κAA2 αD̃i D̃j A
κΓ̃H αD̃i H
κΓ̃G1 αG iκΓ̃G2 αD̃j D̃j G iκΓ̃G3 αD̃i D̃j G j

CAFs√(0, 0, ±√−k2 (∗3), 8κφH k2 )(0, 0, ±√−k2 (∗2), long expressions)(0, 0, ±√−k2 (∗3), (3/2)κSD k2 )(0, 0, ± −k2 (∗2), long expressions)(0,0, long√ expressions)(0, 0, ±√−k2 (∗3), 3κγ̃S1 )2
2(0, 0, ± −k
√ (∗3), −κγ̃S2 k )(0, 0, 0,√± −k2 (∗2), long expressions)(0, 0, ± −k2 (∗3), −κAM1 k2 )(0,0, long√ expressions)(0, 0, ±√−k2 (∗3), 3κAA1 )(0, 0, ±√−k2 (∗3), −κAA2 k2 )(0, 0, ± −k2 (∗3), −κAA2 k2 )(0,0, long expressions)(0,0, long expressions)(0,0, long expressions)

diag?
no
yes
yes
yes
yes
no
no
no
yes
yes
yes
yes
no
yes
yes
yes

eﬀect of the adjustment
κφH < 0 makes 1 Neg.
κφG < 0 makes 2 Neg. 1 Pos.
κSD < 0 makes 1 Neg.
Case (B)κγ̃G1 > 0 makes 1 Neg.
κγ̃G2 < 0 makes 6 Neg. 1 Pos.
Case (E1)
κγ̃S1 < 0 makes 1 Neg.
κγ̃S2 > 0 makes 1 Neg.
κKM < 0 makes 2 Neg.
κAM1 > 0 makes 1 Neg.
κAM2 > 0 makes 7 Neg
Case (D)κAA1 < 0 makes 1 Neg.
κAA2 > 0 makes 1 Neg.
κΓ̃H > 0 makes 1 Neg.
κΓ̃G1 < 0 makes 6 Neg. 1 Pos.
Case (E2)κΓ̃G2 > 0 makes 2 Neg. 1 Pos.
κΓ̃G3 > 0 makes 2 Neg. 1 Pos.

Table 4: Possible adjustements which make a real-part constraint ampliﬁcation factors negative. The
column of adjustments are nonzero multipliers in terms of (3.57)-(3.61), which all violate time reversal
symmetry of the equation. The column ‘diag?’ indicates diagonalizability of the constraint propagation
matrix. Neg./Pos. means negative/positive respectively.
(E) A system which has 6 negative and 1 positive CAFs
The below two adjustments will make 6 negative real CAFs, while they also produce one positive
real CAF (a constraint violating mode). The eﬀectiveness is not clear at this moment, but wethink they are worth to be tested in numerical experiments.

3.3.4

(E1)

∂t γ̃ij = ∂tBSSN γ̃ij + κγ̃G2 αγ̃k(i D̃j) G k ,withκγ̃G2 < 0.

(3.63)(E2)∂t Γ̃i = ∂tBSSN Γ̃i + κΓ̃G2 αD̃j D̃j G i ,withκΓ̃G2 < 0.

(3.64)

Remarks

We have studied step-by-step where the replacements in the equations aﬀect and/or newly added
constraints work, by checking whether the error of constraints (if it exists) will decay or propagate
away. Alcubierre et al [5] pointed out the importance of the replacement (adjustment) of terms in
the evolution equation due to the momentum constraint, and our analysis clearly explain why they
concluded this is the key. Not only this adjustment, we found, but also other adjustments and otherintroductions of new constraints also contribute to making the evolution system more stable. We
found that if we missed a particular adjustment, then the expected stability behaviour occationallygets worse than the ADM system. We further propose other adjustments of the set of equations which
may have better features for numerical treatments.
The discussion was only in the ﬂat background spacetime, and further analysis is in progress.
However, we rather believe that the general fundamental aspects of constraint propagation analysis
are already revealed here. This is because, for the ADM and its adjusted formulation cases, we foundthat the better formulations in the ﬂat background are also better in the Schwarzschild spacetime, while
there are diﬀerences on the eﬀective adjusting multipliers or the eﬀective coordinate ranges [78, 95].
Actually, recently Yo, Baumgarte and Shapiro [90] reported their simulations of stationary rotating
black hole, and mentioned that the above proposal (B) was contributed to maintain their evolution of
Kerr black hole (J/M up to 0.9M ) for long time (t ∼ 6000M ). Their results also indicates that theevolved solution is closed to the exact one, that is, the constrained surface.

35

4

Outlook

4.1

What we have achieved

Let us summarize our story ﬁrst. The beginnings of the study are:
• General relativistic numerical simulations are quite important in astrophysical studies, but we do
not have a deﬁnite recipe to integrate the Einstein equations in a long-term stable and accurate
manner.
• The most standard approach is to decompose the space-time into 3-space and time (3 + 1
decomposition), to solve the constraints for obtaining the initial data, and to evolve the spacetime applying a free evolution scheme.
• Over a period of decades, the community has been observing the violation of constraints in theirsimulations, and their consensus is that we have to ﬁnd a better formulation of the Einstein
evolution equations.
We reviewed recent eﬀorts on this problem by categorizing them into(0) The standard ADM formulation (§2.0, Box 2.1),(1) The modiﬁed ADM (so-called BSSN) formulation (§2.1, Box 2.3),(2) Hyperbolic formulations (§2.2, Box 2.5), and(3) Asymptotically constrained formulations (§2.3).
Among them, the approach (2) is perhaps best justiﬁed on mathematical grounds. However, as we
critically reviewed in §2.2.3, the practical advantages may not be available unless we kill the lower-orderterms in its hyperbolized equation, as in KST’s formulation (Box 2.6).
We therefore proceeded in the direction (3). Our approach, which we term adjusted system, is toconstruct a system that has its constraint surface as an attractor. Our uniﬁed view is to understandthe evolution system by evaluating its constraint propagation. Especially we proposed to analyze theconstraint ampliﬁcation factors (Box 3.1) which are the eigenvalues of the homogenized constraint
propagation equations. We analyzed the system based on our conjecture (Box 3.2/3.3) whether theconstraint ampliﬁcation factors suggest the constraint to decay/propagate or not. We concluded that
• 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.
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 is36

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.

4.2

Next steps

Here are some directions for future researches in this formulation problem of the Einstein equations.
Generalize the adjusted systems: We have tried to unify all the eﬀorts of the community applying
the idea of “adjusted systems”. Our newly modiﬁed equations are working as desired up to the current
numerical tests, but we see also that the eﬀect is not yet perfect to remove non-linear error growththat terminates numerical simulations. We suspect that this is due to the current eigenvalue analysis
based on a perturbative analysis with ﬁxing adjusting multiplier, κ. One remedy is to generalize thedetermination process of κ, say, dynamically and automatically under a suitable principle. We are
now working on a method to control the violation of each constraint independently, or an additional
supporting mathematical criteria to realize more robust stabilizations.
More on hyperbolic formulations: We have already pointed out several directions on hyperbolic
eﬀorts in §2.2.3. We here only point out the links to the initial-boundary value problem (IBVP). In
order to avoid unphysical incoming information from the boundary of computation, a proper treatmentof the boundary is quite important problem in numerical simulations. If we can treat boundary
condition as a part of mathematical framework, that would be a great step to promote researches.
The IBVP approach to the Einstein equations is quite new and has only been studied in a restricted
cases. The current proposals are based on a particular symmetric hyperbolic formulation or under acertain assumption to the symmetry of space-time. More general treatments on IBVP are expected.
Alternative new ideas?: We have to be open to develop an alternative approach to this formulation
problem. If our goal is to obtain a stable system on its constrained surface, then there may be at least
three fundamental approaches: (1) construct an improved evolution system (as we discussed most),(2) develop a “maintenance method” for the system, or (3) redeﬁne or classify the stability of thesystem. An idea of automatic adjustment of multipliers in adjusted system is a sort of “maintenance”.
We guess there might be a key in existing theories such as control theories, optimization methods(convex functional theories), mathematical programming methods, or others.
Numerical comparisons of formulations: Fortunately, an eﬀort toward systematic numericalcomparisons of diﬀerent formulations of the Einstein equations has been organized and started recently.
The workshop “Comparisons of Formulations of Einstein’s equations for Numerical Relativity” washeld at Mexico City in May 2002 5 , and more than 20 people attended this 2-week workshop. This
project intends to provide a format for comparisons in a common numerical framework. That is,comparing diﬀerent formulations using the same initial data, the same resolution, the same integration
scheme, the same boundary treatment, and the same output. Comparisons are now in progress for
vacuum and regular space-time, and we hope to extend them to black-hole space-time next. The ﬁrst
reports are in preparation [59].

4.3

Final remarks

If we say the ﬁnal goal of this project is to ﬁnd a robust algorithm to obtain long-term accurate
and stable time-evolution method, then the recipe should be a combination of (a) formulations of the5

Follow-up information is available at http://www.nuclecu.unam.mx/˜gravit/main rn.html.

37

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 the ones we do havemore robust and automatic.
We have written this review from the viewpoint that the general relativity is a constrained dynamical system. 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.
When we discuss asymptotically constrained manifolds, we implicitly assume that the dynamics
could be expressed on a wider manifold. Recently we found that such a proposal is quite similar
with techniques in e.g. molecular dynamics. For example, people assume an extended environment
when simulating molecular dynamics under constant pressure (with a potential piston) or a constanttemperature (with a potential thermostat) (see e.g. [65]). We have also noticed that a dynamical
adjusting method of Lagrange multipliers has been developed in multi-body mechanical dynamics (seee.g. [62]). We are now trying to apply these ideas back into numerical relativity.
In such a way, communication and interaction between diﬀerent ﬁelds is encouraged. Cooperation
between numerical and mathematical scientists is necessary. By interchanging ideas, we hope we willreach our goal in next few years, and obtain interesting physical results and predictions. It is our
personal view that exciting revolutions in numerical relativity are coming soon.

38A

General expressions of ADM constraint propagation equations

For the reader’s convenience, we express here the constraint propagation equations generally, considering the adjustments to the evolution equations.

A.1

The standard ADM equations and constraint propagations

We start by analyzing the standard ADM system, that is, with evolution equations (2.5) and (2.6)
and constraint equations (2.7) and (2.8).
The constraint propagation equations, which are the time evolution equations of the Hamiltonian
constraint (2.7) and the momentum constraints (2.8).
Expression using H and Mi The constraint propagation equations can be written as [these are
the same with (2.9) and (2.10) ]∂t H = β j (∂j H) + 2αKH − 2αγ ij (∂i Mj )+α(∂l γmk )(2γ ml γ kj − γ mk γ lj )Mj − 4γ ij (∂j α)Mi ,(A.1)∂t Mi = −(1/2)α(∂i H) − (∂i α)H + β (∂j Mi )j

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

(A.2)

This is a suitable form to discuss hyperbolicity of the system. The simplest derivation of (A.1) and(A.2) is by using the Bianchi identity, which can be seen in Frittelli [39].
A shorter expression is available, e.g.
√∂t H = β l ∂l H + 2αKH − 2αγ −1/2 ∂l ( γMl ) − 4(∂l α)Ml
= β l ∇l H + 2αKH − 2α(∇l Ml ) − 4(∇l α)Ml ,(A.3)∂t Mi = −(1/2)α(∂i H) − (∂i α)H + β ∇l Mi + αKMi + (∇i βl )Mll= −(1/2)α(∇i H) − (∇i α)H + β l ∇l Mi + αKMi + (∇i βl )Ml ,(A.4)

or by using Lie derivatives along αnµ ,√£αn µ H = 2αKH − 2αγ −1/2 ∂l ( γMl ) − 4(∂l α)Ml ,£αn µ Mi(A.5)

= −(1/2)α(∂i H) − (∂i α)H + αKMi .

(A.6)

Expression using γij and Kij In order to check the eﬀects of the adjustments in (2.5) and (2.6)
to constraint propagation, it is useful to re-express (A.1) and (A.2) using γij and Kij . By a straightforward calculation, we obtain an expression as∂t H = H1mn (∂t γmn ) + H2imn ∂i (∂t γmn ) + H3ijmn ∂i ∂j (∂t γmn ) + H4mn (∂t Kmn ),(A.7)∂t Mi = M1i∂j (∂t Kmn ),(A.8)+γ ij γ np (∂i γ mk )(∂j γkp ) − γ mp γ ni (∂i γ kj )(∂j γkp ) − 2KK mn + 2K n j K mj ,

(A.9)

mn

(∂t γmn ) + M2i

jmn∂j (∂t γmn ) + M3i

mn(∂t Kmn ) + M4i

jmn

where
H1mn := −2R(3)mn − Γpkj Γkpi γ mi γ nj + Γm Γn
H2imn
H3ijmnH4mn
M1i mn:= −2γΓ − (3/2)γ (∂j γ

mi n:= −γ γij mn:= 2(Kγnj

mn:= γ (∂i K

ij

in mj

+γ γ
−K

m

mn

),

j) − γ

mj

mn

)+γ

mj

in

(∂j γ ) + γ

mn i

Γ,

(A.10),(A.11)(A.12)
n(∂j K i ) + (1/2)(∂j γ
39

mn

j

n

)K i + Γ K

mi,(A.13)

M2i jmn := −γ mj K n i + (1/2)γ mn K j i + (1/2)K mn δij ,
mn

:=

jmn

:=

M3i
M4i

(A.14)−δin Γm − (1/2)(∂i γ mn ),γ mj δin − γ mn δij ,(A.15)(A.16)

ij
where we expressed Γm = Γmij γ .

A.2

Constraint propagations for the adjusted ADM systems

Generally, we here write the adjustment terms to (2.5) and (2.6) using (2.7) and (2.8) by the following
combinations, [these are the same with (3.7) and (3.8)]
adjustment term of∂t γij :adjustment term of ∂t Kij :+Pij H + Qk ij Mk + pk ij (∇k H) + q kl ij (∇k Ml ),(A.17)+Rij H + S ij Mk + r ij (∇k H) + s(A.18)kkkl

ij (∇k Ml ),where P, Q, R, S and p, q, r, s are multipliers (please do not confuse Rij with three Ricci curvature(3)
that we write as Rij ). We adjust them only using up to the ﬁrst derivatives in order to make thediscussion simple.
By substituting the above adjustments into (A.7) and (A.8), we can write the adjusted constraint
propagation equations as∂t H = (original terms)+H1mn [Pmn H + Qkmn Mk + pk mn (∇k H) + q kl mn (∇k Ml )]+H2imn ∂i [Pmn H + Qkmn Mk + pk mn (∇k H) + q kl mn (∇k Ml )]+H3ijmn ∂i ∂j [Pmn H + Qkmn Mk + pk mn (∇k H) + q kl mn (∇k Ml )]k
Mk + rk mn (∇k H) + skl mn (∇k Ml )],+H4mn [Rmn H + Smn(A.19)∂t Mi = (original terms)+M1i mn [Pmn H + Qkmn Mk + pk mn (∇k H) + q kl mn (∇k Ml )]+M2i jmn ∂j [Pmn H + Qkmn Mk + pk mn (∇k H) + q kl mn (∇k Ml )]k
Mk + rk mn (∇k H) + skl mn (∇k Ml )]+M3i mn [Rmn H + Smnk
Mk + rk mn (∇k H) + skl mn (∇k Ml )].
+M4i jmn ∂j [Rmn H + Smn(A.20)

Here the “original terms” can be understood either as (A.1) and (A.2), or as (A.7) and (A.8). Therefore, for example, we can see that adjustments to ∂t γij do not always keep the constraint propagation
equations in the ﬁrst order form, due to their contribution in the third adjusted term in (A.19).
We note that these expressions of constraint propagation equations are equivalent when we include
the cosmological constant and/or matter terms.

40

B

Numerical demonstrations using the Ashtekar formulation

This appendix is devoted to introduce our numerical comparisons between three levels of hyperbolicity
using Ashtekar variables. Details are available in [77, 94].

B.1

The Ashtekar formulation

The key feature of Ashtekar’s formulation of general relativity [13] is the introduction of a self-dual
connection as one of the basic dynamical variables. Let us write the metric gµν using the tetrad EµI
as gµν = EµI EνJ ηIJ 6 . Deﬁne its inverse, EIµ , by EIµ := EνJ g µν ηIJ and we impose Ea0 = 0 as the gauge
condition. We deﬁne SO(3,C) connections ±Aaµ := ωµ0a ∓ (i/2)a bc ωµbc , where ωµIJ is a spin connection
1-form (Ricci connection), ωµIJ := E Iν ∇µ EνJ . Ashtekar’s plan is to use only the self-dual part of theconnection +Aaµ and to use its spatial part +Aai as a dynamical variable. Hereafter, we simply denote+Aa as Aa .
µµ
The lapse function, N , and shift vector, N i , both of which we treat as real-valued functions7 , areexpressed as E0µ = (1/N, −N i /N ). This allows us to think of E0µ as a normal vector ﬁeld to Σ spannedby the condition t = x0 =const., which plays the same role as that of the ADM formulation. Ashtekar
treated the set (Ẽai , Aai ) as basic dynamical variables, where Ẽai is an inverse of the densitized triad
deﬁned by Ẽai := eEai , where e := det Eia is a density 8 . This pair forms the canonical set. In the caseof pure gravitational spacetime, the Hilbert action takes the formS =d4 x[(∂t Aai )Ẽai + (i/2)NẼai Ẽbj Fijc ab c − N i Fija Ẽaj + Aa0 Di Ẽai ],∼(B.21)

a := 2∂ Aa − ia Ab Ac is the curvature 2-form, D Ẽ j := ∂ Ẽ j − i c Ab Ẽ j .
:= e−1 N , Fµν
where Ni ai abc µ ν
ab
[µ ν]
i c
∼

The action (B.21) gives us the following evolution equations and constraints:

The Ashtekar formulation [13]:The dynamical variables are (Ẽai , Aai ).
The evolution equations for a set of (Ẽai , Aai ) areBox B.1

Ẽcj Ẽbi ) + 2Dj (N [j Ẽai] ) + iAb0 ab c Ẽci ,∂t Ẽai = −iDj (cb a N∼(B.22)∂t Aai = −iab c N
Ẽbj Fijc + N j Fjia + Di Aa0 ,∼(B.23)

where Dj Xaji := ∂j Xaji − iab c Abj Xcji , and Fija := 2∂[i Aaj] − ia bc Abi Acj .
Constraint equations: (Hamiltonian, momentum and Gauss constraints)
ASH:= (i/2)ab c Ẽai Ẽbj Fijc ≈ 0,
CH

(B.24)

ASH:= −Fija Ẽaj ≈ 0,CMi

(B.25)

ASH
CGa(B.26):=

Di Ẽai ≈ 0,

We use I, J = (0), · · · , (3) and a, b = (1), · · · , (3) are SO(1, 3), SO(3) indices respectively. We raise and lowerµ, ν, · · · by g µν and gµν (the Lorentzian metric); I, J, · · · by η IJ = diag(−1, 1, 1, 1) and ηIJ ; i, j, · · · by γ ij and γij (thethree-metric); a, b, · · · by δ ab and δab . We also use volume forms abc : abc abc = 3!.
7
N and N i are the same with α and β i in the previous text. We follow this conventional notation in this appendix.
8
For later convenience, e2 = det Ẽai = (det Eia )2 = (1/6)abc ∼ ijk Ẽai Ẽbj Ẽck , where ijk := abc Eia Ejb Ekc and ∼ ijk :=6

e−1 ijk . When (i, j, k) = (1, 2, 3), we have ijk = e, ∼ ijk = 1, ijk = e−1 , and ˜ijk = 1.

41

We have to consider the reality conditions when we use this formalism to describe the classical
Lorentzian spacetime. Fortunately, the metric will remain on its real-valued constraint surface during
time evolution automatically if we prepare initial data which satisﬁes the reality condition. Morepractically, we further require that triad is real-valued. But again this reality condition appears asa gauge restriction on Aa0 [91], which can be imposed at every time step. In our actual simulation,
we prepare our initial data using the standard ADM approach, so that we have no diﬃculties inmaintaining these reality conditions.

B.2

Reformulate the Ashtekar evolution equations

B.2.1

Strongly and Symmetric Hyperbolic systems

The authors’ recent studies showed the following:(a) The original set of dynamical equations (B.22) and (B.23) [the original equations] already forms
a weakly hyperbolic system [93]. So that we regard the mathematical structure of the original
equations as one step advanced from the standard ADM.
(b) Further, we can construct higher levels of hyperbolic systems by restricting the gauge condition
ASH , C ASH and C ASH , to the original equations.
and/or by adding constraint terms, CH
Mi
Ga– by requiring additional gauge conditions or adding constraints to the dynamical equations,
we can obtain a strongly hyperbolic system [93],– by requiring additional gauge conditions and adding constraints to the dynamical equations,
we can obtain a symmetric hyperbolic system [92, 93].
(c) Based on the above symmetric hyperbolic system, we can construct an Ashtekar version of the λsystem [24] which is robust against perturbative errors for both constraints and reality conditions[76].
In order to obtain a symmetric hyperbolic system 9 , we add constraint terms to the right-hand-sideof (B.22) and (B.23). The adjusted dynamical equations,
ASH b,Ẽcj Ẽbi ) + 2Dj (N [j Ẽai] ) + iAb0 ab c Ẽci + κ1 P i ab CG∂t Ẽai = −iDj (cb a N
∼

where

 c Ẽci ,P i ab ≡ N i δab + iN
∼ ab

ASH
ASH
∂t Aai = −iab c N+ κ3 Ri ja CM
Ẽbj Fijc + N j Fjia + Di Aa0 + κ2 Qai CH
j ,
∼

where

(B.27)

Ẽia ,
Qai ≡ e−2 N∼(B.28)

Ri ja ≡ ie−2 N
ac b Ẽib Ẽcj∼

form a symmetric hyperbolicity if we further require κ1 = κ2 = κ3 = 1 and the gauge conditions,
Aa0 = Aai N i ,∂i N = 0.

(B.29)We remark that the adjusted coeﬃcients, P i ab , Qai , Ri ja , for constructing the symmetric hyperbolic
ASH , C ASH for ∂ Ẽ i ,system are uniquely determined, and there are no other additional terms (say, no CH
t a
M
ASH
a
no CG for ∂t Ai ) [93]. The gauge conditions, (B.29), are consequences of the consistency with (triad)
reality conditions.
We can also construct a strongly (or diagonalizable) hyperbolic system by restricting to a gauge
N l = 0, ±N γ ll (where γ ll is the three-metric and we do not sum indices here) for the original
9

Iriondo et al [48] presented a symmetric hyperbolic expression in a diﬀerent form. The diﬀerences between ours andtheirs are discussed in [92, 93].

42

systemvariables

Eqs of motion

I
II
III
adj

Ashtekar (weakly hyp.)Ashtekar (strongly hyp.)Ashtekar (symmetric hyp.)Ashtekar (adjusted)(B.22), (B.23) (original)(B.27), (B.28) (with κ = 1)
(B.27), (B.28) (with κ = 1)
(B.27), (B.28) (with κ = 1)λAshtekar-λ-system(Ẽai , Aai )(Ẽai , Aai )(Ẽai , Aai )(Ẽai , Aai )(Ẽai , Aai ,remarkcontrols CH , CM i , CGa(B.37)

λ, λi , λa )“original” eqs.
(B.30) required(B.29) required

Table 5: List of systems that we compare in this appendix.

equations (B.22), (B.23). Or we can also construct from the adjusted equations, (B.27) and (B.28),
together with the gauge condition(B.30)
Aa0 = Aai N i .
As for the strongly hyperbolic system, we hereafter take the latter expression.
In Table 5, we summarized the equations to be used for our comparisons.
B.2.2

Ashtekar-λ-systemIn §2.3.1, we introduced an idea to construct a robust evolution system against a perturbative error,
named “λ-system” [24]. Based on the above symmetric hyperbolic equations, we constructed the
Ashtekar version of the “λ-system” [76]. Here we present our system which evolves the spacetime tothe constraint surface, CH ≈ CM i ≈ CGa ≈ 0 as the attractor. In [76], we also presented a system
which controls the perturbative violation of the reality condition.
We introduce new variables (λ, λi , λa ), as they obey the dissipative evolution equations∂t λ = α1 CH − β1 λ,∂t λi = α2 C˜M i − β2 λi ,(B.32)∂t λa = α3 CGa − β3 λa ,(B.33)(B.31)

where αi = 0 (allowed to be complex numbers) and βi > 0 (real numbers) are constants.
(DL)
If we take uα= (Ẽai , Aai , λ, λi , λa ) as a set of dynamical variables, then the principal part of(B.31)-(B.33) can be written as∂t λ ∼= −iα1 bcd Ẽcj Ẽdl (∂l Abj ),∼ α2 [−eδ l Ẽ j + eδ j Ẽ l ](∂l Ab ),∂t λi =

(B.35)∂t λa ∼= α3 ∂l Ẽal .

(B.36)

i

b

i

b

j

(B.34)(DL)

The characteristic matrix of the system uαdoes not form a Hermitian matrix. However, if we
modify the right-hand-side of the evolution equation of (Ẽai , Aai ), then the set becomes a symmetric hyperbolic system. This is done by adding ᾱ3 γ il (∂l λa ) to the equation of ∂t Ẽai , and by addingiᾱ1 a c d Ẽic Ẽdl (∂l λ) + ᾱ2 (−eγ lm Ẽia + eδim Ẽ la )(∂l λm ) to the equation of ∂t Aai . The ﬁnal principal part,then, is written as  Al bi
Ẽai
a m
0
 Aai   
0∂t  λ  ∼
=

 

λi
0
λa
α3 δa b δ l m
00

Dla ib miᾱ1 a c d Ẽic Ẽdl
cd m l−iα1 b Ẽc Ẽd00
α2 e(δim Ẽbl − δil Ẽbm )0043

  m
ᾱ3 γ il δa b
Ẽb
  Abm ᾱ2 e(δim Ẽ la − γ lm Ẽia )0  ∂l  λ 
0
0
.
 


λm
0
0
λb
0
0
(B.37)
0Clearly, the solution (Ẽai , Aai , λ, λi , λa ) = (Ẽai , Aai , 0, 0, 0) represents the original solution of the
Ashtekar system. If the λs decay to zero after the evolution, then the solution also describes the(DL)
original solution of the Ashtekar system in that stage. Since the dynamical system of uα , (B.37),
constitutes a symmetric hyperbolic form, the solutions to the λ-system are unique. Therefore, thedynamical system, (B.37), is useful for stabilizing numerical simulations from the point that it recovers
the constraint surface automatically.
B.2.3

Adjusted Ashtekar system

We also try to compare a set of evolution system, which we proposed as “adjusted-system”. The
fundamental equations that we will demonstrate are the same with (B.27) and (B.28), but here thereal-valued constant multipliers κs are not necessary equals to unity. We set κ ≡ κ1 = κ2 = κ3 forsimplicity. Apparently the set of (B.27) and (B.28) becomes the original weakly hyperbolic system if
κ = 0, becomes the symmetric hyperbolic system if κ = 1 and N = const.. The set remains strongly
hyperbolic systems for other choices of κ except κ = 1/2 which only forms weakly hyperbolic system.

B.3

Comparing numerical performance

B.3.1

Model and Numerical method

The model we present here is gravitational wave propagation in a planar spacetime under periodic
boundary condition. We performed a full numerical simulation using Ashtekar’s variables. We prepare
two +-mode strong pulse waves initially by solving the ADM Hamiltonian constraint equation, using
York-O’Murchadha’s conformal approach. Then we transform the initial Cauchy data (3-metric andextrinsic curvature) into the connection variables, (Ẽai , Aai ), and evolve them using the dynamical
equations. For the presentation in this article, we apply the geodesic slicing condition (ADM lapse
N = 1 or densitized lapse N= 1, with zero shift and zero triad lapse). We have used both the∼Brailovskaya integration scheme, which is a second order predictor-corrector method, and the socalled iterative Crank-Nicholson integration scheme for numerical time evolution. The details of thenumerical method are described in [77].
More speciﬁcally, we set our initial guess 3-metric as


1
0
2
2
γ̂ij =  sym. 1 + K(e−(x−L) + e−(x+L) )
sym.
sym.

0
0


,

(B.38)

2
2
1 − K(e−(x−L) + e−(x+L) )in the periodically bounded region x = [−5, +5]. Here K and L are constants and we set K = 0.3 andL = 2.5 for the plots.
In order to show the expected “stabilization behavior” clearly, we artiﬁcially add an error in the
middle of the time evolution. We added an artiﬁcial inconsistent rescaling once at time t = 6 for theA2y component as A2y → A2y (1 + error).
B.3.2

Diﬀerences between three levels of hyperbolicity

We have performed comparisons of stability and/or accuracy between weakly and strongly hyperbolic
systems, and between weakly and symmetric hyperbolic systems[77]. (We can not compare strongly
and symmetric hyperbolic systems directly, because these two requires diﬀerent gauge conditions.)
We omit ﬁgures in this report, but one can see a part of results in Fig.9 and Fig.10. We may
conclude that higher level hyperbolic system gives us slightly accurate evolution. However, if we
evaluate the magnitude of L2 norms, then we also conclude that there is no measurable diﬀerences
between strongly and symmetric hyperbolicities. This last fact will be supported more aﬃrmatively
in the next experiments.
44

(a)

0.035

0.030

0.025

L2 norm of CH(ASH)

L2 norm of CH(ASH)

0.030

weakly hyp.
strongly hyp.
λ-system (α=1, β=1)λ-system (α=1, β=5)λ-system (α=1, β=10)

(b)

0.035

no errorerror +10%
error +20%

0.020
0.015
0.010
0.005

0.025
0.020
0.015
0.010
0.005

artificial error added at time=6
0.000
0.0

5.0

10.0

0.000

15.0

0.0

5.0

10.0

15.0

time

time

Figure 9: Demonstration of the Ashtekar-λ-system for the cases of plane wave propagation under the periodic
boundary. We plot the L2 norm of the Hamiltonian constraint equation, CH . Fig. (a) shows how the system
goes bad depending on the amplitude of artiﬁcial error at t = 6. All the lines are of the evolution of Ashtekar’s
original equation (no λ-system). Fig. (b) shows the eﬀect of λ-system. All the lines are 20% error amplitude,
but shows the diﬀerence of evolution equations. The solid line is for Ashtekar’s original equation (the same as inFig.(a)), the dotted line is for the strongly hyperbolic Ashtekar’s equation. Other lines are of λ-systems, whichcproduces better performance than that of the strongly hyperbolic system. (Reprinted from [94], IOP
2001)

B.3.3

Demonstrating “λ-system”

Next, we show a result of the “λ-system” [94]. Fig.9 (a) shows how the violation of the Hamiltonian
constraint equation, CH , become worse depending on the term error. The oscillation of the L2 norm
CH in the ﬁgure due to the pulse waves collide periodically in the numerical region. We, then, ﬁx theerror term as a 20% spike, and try to evolve the same data in diﬀerent equations of motion, i.e., theoriginal Ashtekar’s equation [solid line in Fig.9 (b)], strongly hyperbolic version of Ashtekar’s equation(dotted line) and the above λ-system equation (other lines) with diﬀerent βs but the same α. As we
expected, all the λ-system cases result in reducing the Hamiltonian constraint errors.
(b)(a)κ=−1.0κ=0.0 (weakly hyp.)κ=1.0 (symmetric hyp.)κ=+2.0

0.050

0.035

κ=−1.0κ= 0.0 (weakly hyp.)κ=1.0 (symmetric hyp.)κ=+2.0

0.030

0.040

L2 norm of CH(ASH)0.020

κ=−1.0at time=6, we added
an artificial error

0.015

κ=0.0

L2 norm of CM(ASH)κ=−1.0
0.025

κ=0.0

0.030

0.020

κ=1.0

0.010

0.010

κ=1.0
0.005

κ=2.0
κ=2.0

0.000

0.000
0.0

4.0

8.0

12.0

0.0

16.0

4.0

8.0

12.0

16.0

time

time

Figure 10: Demonstration of the adjusted-Ashtekar systems. The experiments as in Fig.9(b). Fig. (a) and (b)are L2 norm of the Hamiltonian constraint equation, CH , and momentum constraint equation, CM x , respectively.
The solid line is the case of κ = 0, that is the case of “no adjusted” original Ashtekar equation (weakly hyperbolic
system). The dotted line is for κ = 1, equivalent to the symmetric hyperbolic system. We see other line (κ = 2.0)cshows better performance than the symmetric hyperbolic case. (Reprinted from [94], IOP
2001)

45

B.3.4

Demonstrating “adjusted system”

We here examine how the adjusted multipliers contribute to the system’s stability [94]. In Fig.10, weshow the results of this experiment. We plot the violation of the constraint equations both CH andCM x . An artiﬁcial error term was added in the same way as above. The solid line is the case of κ = 0,
that is the case of “no adjusted” original Ashtekar equation (weakly hyperbolic system). The dottedline is for κ = 1, equivalent to the symmetric hyperbolic system. We see other line (κ = 2.0) shows
better performance than the symmetric hyperbolic case.

Acknowledgment
HS thanks M. Alcubierre, A. Ashtekar, T. Baumgarte, C. Bona, G. Calabrese, S. Detweiler, N. Dorband, S. Husa, P. Laguna, C. Lechner, L.Lehner, L. Lindblom, V. Moncrief, T. Nakamura, D. Pollney,
O. Sarbach, M. Scheel, E. Seidel, M. Shibata, D. Shoemaker, and J. Winicour for scientiﬁc discussion
in the past year. He also thanks Albert Einstein Institute Cactus terms for technical helps at the
Mexico numerical relativity workshop. He also thanks T. Ebisuzaki for continuous encouragement.
HS is supported by 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 Promotionof Science, No. 14740179. We thank J. Overduin for careful reading a part of the manuscript. GY
thanks H. Shinkai for writing all the manuscript.

46

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 B. Brügmann, Phys. Rev. D 63, 104006 (2001).
[4] M. Alcubierre and J. Massó, Phys. Rev. D 57, R4511 (1998).
[5] M. Alcubierre, G. Allen, B. Brügmann, E. Seidel and W-M. Suen, Phys. Rev. D 62, 124011(2000).
[6] M. Alcubierre, W. Benger, B. Brügmann, G. Lanfermann, L. Nerger, E. Seidel, and R. Takahashi,
Phys. Rev. Lett. 87 271103 (2001).
[7] 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).
[8] M. Alcubierre, B. Brügmann, M. Miller and W-M. Suen, Phys. Rev. D 60, 064017 (1999).
[9] A. Anderson, Y. Choquet-Bruhat and J.W. York, Jr., Topol. Methods in Nonlinear Analysis,
10, 353 (1997), available as gr-qc/9710041.
[10] A. Anderson and J. W. York, Jr, Phys. Rev. Lett. 82, 4384 (1999).
[11] A. Arbona, C. Bona, J. Massó, and J. Stela, Phys. Rev. D 60, 104014 (1999).
[12] R. Arnowitt, S. Deser and C.W. Misner, “The Dynamics of General Relativity”, in Gravitation:
An Introduction to Current Research, ed. by L.Witten, (Wiley, New York, 1962).
[13] A. Ashtekar, Phys. Rev. Lett. 57, 2244 (1986); Phys. Rev. D36, 1587 (1987).
[14] A. Ashtekar, Lectures on Non-Perturbative Canonical Gravity (World Scientiﬁc, Singapore,1991).
[15] T.W. Baumgarte and S.L. Shapiro, Phys. Rev. D 59, 024007 (1999).
[16] T.W. Baumgarte and S.L. Shapiro, gr-qc/0211028.
[17] J.M. Bardeen, L.T. Buchman, Phys. Rev. D. 65, 064037 (2002).
[18] C. Bona, J. Massó, Phys. Rev. D 40, 1022 (1989).
[19] C. Bona, J. Massó, Phys. Rev. Lett. 68, 1097 (1992).
[20] C. Bona, J. Massó, E. Seidel and J. Stela, Phys. Rev. Lett. 75, 600 (1995).
[21] C. Bona, J. Massó, E. Seidel and J. Stela, Phys. Rev. D 56, 3405 (1997).
[22] C. Bona, J. Massó, E. Seidel, and P. Walker, gr-qc/9804052.
[23] L. Brewin, Class. Quant. Grav. 19, 429 (2002).
[24] O. Brodbeck, S. Frittelli, P. Hübner, and O.A. Reula, J. Math. Phys. 40, 909 (1999).
[25] B. Brügmann, Annalen Phys. 9, 227 (2000), available as gr-qc/9912009.
[26] G. Calabrese, L. Lehner, and M. Tiglio, Phys. Rev. D 65, 104031 (2002).
[27] G. Calabrese, J. Pullin, O. Sarbach, and M. Tiglio, Phys. Rev. D 66 064011 (2002).
[28] G. Calabrese, J. Pullin, O. Sarbach, and M. Tiglio, Phys. Rev. D 66, 041501 (2002).
47

[29] G. Calabrese, J. Pullin, O. Sarbach, M. Tiglio, and O. Reula, gr-qc/0209017.
[30] M.W. Choptuik, Phys. Rev. D 44, 3124 (1991).
[31] Y. Choquet-Bruhat and J.W. York. Jr., in General Relativity and Gravitation , vol. 1, ed. by A.
Held, (Plenum, New York, 1980).
[32] Y. Choquet-Bruhat and J.W. York, Jr., C.R. Acad. Sc. Paris 321, Série I, 1089, (1995), available
as gr-qc/9506071.
[33] S. Detweiler, Phys. Rev. D 35, 1095 (1987).
[34] J. Frauendiener, Livng Rev. Relativ. 2000-4 at http://www.livingreviews.org/[35] J. Frauendiener, Phys. Rev. D 58 064002 (1998); ibid. 58 064003 (1998); Class. Quant. Grav.
17 373 (2000).
[36] H. Friedrich, Proc. Roy. Soc. A375, 169 (1981); ibid. A378, 401 (1981).
[37] H. Friedrich and A. Rendall, in Einstein’s ﬁeld equations and their physical interpretation, ed.
by B.G.Schmidt (Springer, Berlin 2000), available as gr-qc/0002074.
[38] H. Friedrich and G. Nagy, Comm. Math. Phys. 201, 619 (1999).
[39] S. Frittelli, Phys. Rev. D 55, 5992 (1997).
[40] S. Frittelli and R. Gomez, J. Math. Phys. 41, 5535 (2000).
[41] S. Frittelli and O.A. Reula, Phys. Rev. Lett. 76, 4667 (1996).
[42] S. Frittelli and O.A. Reula, J. Math. Phys. 40, 5143 (1999).
[43] R. Geroch, in General Relativity, ed. by G.S. Hall and J.R. Pulham (Institute of Physics Publishing, 1996), available as gr-qc/9602055.
[44] S. D. Hern, PhD thesis, gr-qc/0004036.
[45] C. Hirsch, Numerical Computation of Internal and External Flows Vol I, (John Wiley & Sons,
1988).
[46] P. Hübner, Class. Quant. Grav. 16 2145 (1999), ibid. 16 2823 (1999), ibid. 18 1421 (2001), andibid. 18 1871 (1999).
[47] 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.
[48] M.S. Iriondo, E.O. Leguizamón and O.A. Reula, Phys. Rev. Lett. 79, 4732 (1997), Adv. Theor.
Math. Phys. 2, 1075 (1998).
[49] M.S. Iriondo and O.A. Reula, Phys.Rev. D 65, 044024 (2002).
[50] R. d’Inverno, in General Relativity, ed. by G.S. Hall and J.R. Pulham (Institute of Physics
Publishing, 1996)[51] B. Kelly, P. Laguna, K. Lockitch, J. Pullin, E. Schnetter, D. Shoemaker, and M. Tiglio, Phys.
Rev. D 64, 084013 (2001).
[52] L.E. Kidder, M.A. Scheel, S.A. Teukolsky, Phys. Rev. D 64, 064017 (2001).
[53] A.M. Knapp, E.J. Walker, T.W. Baumgarte Phys. Rev. D 65, 064031 (2002).
48

[54] H-O. Kreiss, O. E. Ortiz, gr-qc/0106085.
[55] P. Laguna and D. Shoemaker, Class. Quant. Grav. 19 (2002) 3679.
[56] L. Lehner, Class. Quant. Grav. 18, R25 (2001); Proceedings of GR16, available as gr-qc/0202055.
[57] L. Lehner, M. Huq, D. Garrison, Phys. Rev. D 62, 084016 (2000).
[58] L. Lindblom and M. Scheel, gr-qc/0206035[59] Mexico Numerical Relativity Workshop 2002 participants, in preparation.
[60] M. Miller, gr-qc/0008017.
[61] C.W. Misner, K.S. Thorne, J.A. Wheeler, Gravitation, (Freeman, N.Y., 1973).
[62] T. Nagata, in Proceedings of the ﬁrst Asian Conference on Multibody Dynamics 2002
(ACMD02), held at Fukushima, Japan.
[63] T. Nakamura, K. Oohara and Y. Kojima, Prog. Theor. Phys. Suppl. 90, 1 (1987).
[64] T. Nakamura and K. Oohara, in Frontiers in Numerical Relativity edited by C.R. Evans, L.S.
Finn, and D.W. Hobill (Cambridge Univ. Press, Cambridge, England, 1989).
[65] S. Nosé, Prog. Theor. Phys. 103, 1 (1991).
[66] Proceedings of Black Holes and Gravitational Waves ed. by T. Nakamura and H. Kodama, as
the issue of Prog. Theor. Phys. Suppl. 136 (1999).
[67] O.A. Reula, Livng Rev. Relativ. 1998-3 at http://www.livingreviews.org/[68] O. Sarbach, G. Calabrese, J. Pullin, and M. Tiglio, Phys. Rev. D 66 064002.
[69] O. Sarbach and M. Tiglio, gr-qc/0205086.
[70] M.A. Scheel, T.W. Baumgarte, G.B.Cook, S.L. Shapiro, and S.A. Teukolsky, Phys. Rev. D 56,
6320 (1997).
[71] M.A. Scheel, T.W. Baumgarte, G.B. Cook, S.L. Shapiro and S.A. Teukolsky, Phys. Rev. D 58,
044020 (1998).
[72] M.A. Scheel, L.E. Kidder, L. Lindblom, H.P. Pfeiﬀer and S.A. Teukolsky, gr-qc/0209115.
[73] 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.
[74] M. Shibata and T. Nakamura, Phys. Rev. D52, 5428 (1995).
[75] M. Shibata and K. Uryu, Prog. Theor. Phys. 107, 265 (2002).
[76] H. Shinkai and G. Yoneda, Phys. Rev. D 60, 101502 (1999).
[77] H. Shinkai and G. Yoneda, Class. Quant. Grav. 17, 4799 (2000).
[78] H. Shinkai and G. Yoneda, Class. Quant. Grav. 19, 1027 (2002).
[79] F. Siebel and P. Hübner, Phys. Rev. D 64, 024021 (2001).
[80] L. Smarr, J.W. York, Jr., Phys. Rev. D 17, 2529 (1978).

49

[81] J.M. Stewart, Class. Quant. Grav. 15, 2865 (1998).
[82] B. Szilágyi, R. Gomez, N.T. Bishop and J. Winicour, Phys. Rev. D 62, 104006 (2000).
[83] B. Szilágyi, B. Schmidt and J. Winicour, Phys. Rev. D 65, 064015 (2002).
[84] B. Szilágyi and J. Winicour, gr-qc/0205044.
[85] S.A. Teukolsky, Phys. Rev. D 26, 745 (1982).
[86] M.H.P.M. van Putten and D.M. Eardley, Phys. Rev. D 53, 3056 (1996).
[87] M.H.P.M. van Putten, Phys. Rev. D 55, 4705 (1997).
[88] 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.
[89] J. Winicour, Livng Rev. Relativ. 2001-3 at http://www.livingreviews.org/[90] H-J. Yo, T.W. Baumgarte and S.L. Shapiro, Phys. Rev. D 66 084026 (2002).
[91] G. Yoneda and H. Shinkai, Class. Quantum Grav. 13, 783 (1996).
[92] G. Yoneda and H. Shinkai, Phys. Rev. Lett. 82, 263 (1999).
[93] G. Yoneda and H. Shinkai, Int. J. Mod. Phys. D. 9, 13 (2000).
[94] G. Yoneda and H. Shinkai, Class. Quant. Grav. 18, 441 (2001).
[95] G. Yoneda and H. Shinkai, Phys. Rev. D 63, 124019 (2001).
[96] G. Yoneda and H. Shinkai, gr-qc/0204002 (to be published in Phys. Rev. D).
[97] G. Yoneda and H. Shinkai, gr-qc/0209106.
[98] J.W. York, Jr., “Kinematics and Dynamics of General Relativity”, in Sources of Gravitational
Radiation, ed. by L.Smarr, (Cambridge, 1979).

50