Baumgarte T., Shapiro S. Numerical Relativity. Solving Einstein’s Equations on the Computer

Подождите немного. Документ загружается.

11.3 Generalized harmonic coordinates 381

which is symmetric hyperbolic. Interestingly, the gauge ﬁelds A

and  have disappeared

entirely from this equation.

If desired we can write the wave equation (11.16) as a system of coupled ﬁrst-order

equations by introducing a new set of ﬁrst-order variables. This is exactly how we found

the system (6.6) from the wave equation (6.3) in Section 6.1. We will discuss similar

approaches in general relativity in Section 11.4.

11.2.3 Auxiliary variables

A third approach to “ﬁxing” Maxwell’s equations is in some ways similar to the ﬁrst one.

We again use the auxiliary variable  deﬁned in equation (11.11). However, instead of

setting this variable to a predetermined function – given by the gauge source function –

we now treat this variable as a new, independent ﬁeld that we evolve. We can derive an

evolution equation for  from equation (11.4),

∂

 = ∂

= D

∂

=−D

− D

 =−D

 − 4πρ

. (11.17)

Note that we have used the constraint equation (11.6) in the last equality, similar to the way

we used the same constraint to arrive at the wave equation (11.16). Exercise 11.4 below

demonstrates how this step affects the properties of constraint violations.

In terms of , the evolution equation (11.5)forE

becomes

∂

=−D

+ D

 − 4π j

. (11.18)

Equations (11.4), (11.17) and (11.18) now constitute the evolution equations in this new

formulation, and equations (11.6) and (11.11) are the constraint equations. In this formu-

lation the mixed derivative term D

has been eliminated without using up any gauge

freedom, which is still imposed via the choice of . In Section 11.5 we will introduce an

analogous reformulation of the ADM equations.

11.3 Generalized harmonic coordinates

One way to avoid the problems associated with the standard ADM evolution equations

(11.7) and (11.8) is to abandon these equations completely. Instead of using a 3 + 1

decomposition of Einstein’s equations as a starting point, we could start with the original

4-dimensional version

= 8π T

(11.19)

To quote from Abrahams and York, Jr., (1997), “the dynamics of electromagnetism have been cleanly separated from

the gauge-dependent evolution of the vector and scalar potentials.”

382 Chapter 11 Recasting the evolution equations

in the form

(4)

= 8π (T

− (1/2)g

T ). Substituting the Ricci tensor as expressed by

equation (4.46) yields

∂

− g

c(a

∂

(4)



−

(4)



c (4)



(ab)c

− 2g

ed (4)



e(a

(4)



b)cd

− g

cd (4)



(4)



ecb

=−8π



−



, (11.20)

where we have reintroduced the deﬁnition (4.40),

(4)



≡ g

bc(4)



=−

|g|

1/2

∂

∂x



|g|

1/2



= g

∇

(11.21)

(see exercises 4.8 and 4.10 and the discussion that follows these exercises). The

(4)



are contractions of Christoffel symbols, and therefore do not transform like vectors under

coordinate transformations. We can now introduce a gauge by setting these quantities equal

to some given gauge source functions H

(4)



= H

(t, x

). (11.22)

This approach follows very closely our electromagnetic example of Section 11.2.1,and

equation (11.22) is the direct analog of (11.14). Just like the choice (11.14) led to the

elimination of the “mixed derivatives” terms in Maxwell’s equations (11.10), inserting

equation (11.22) into Einstein’s equations (11.20) yields a nonlinear wave equation for the

spacetime metric g

After some manipulations this equation can be brought into the

form

∂

+2∂

∂

b)d

+ 2H

(a,b)

− 2H

(4)



(4)



(4)



=−8π

(

− g

)

(11.23)

Here we have lowered the indices of H

with the spacetime metric g

, H

≡ g

. For the

special choice H

= 0 we recover the harmonic coordinates of Chapter 4.3. More generally

we refer to this approach as “generalized harmonic coordinates”.

This formalism was

adopted by Pretorius (2005a,b) in his simulations of binary black hole coalescence and

merger, which we will discuss in more detail in Chapter 13.

Identifying equations (11.21) and (11.22) imposes a new, 4-dimensional constraint

≡ H

− g

bc (4)



= 0. (11.24)

Equation (11.23) can be integrated directly for the spacetime metric g

. To stabilize the

system, it is sometimes necessary to add linear combinations of the constraints (11.24)

to the evolution equations (11.23), i.e., it is necessary to add terms proportional to the

See Friedrich (1985); Garﬁnkle (2002). See also the related “Z4” formalism suggested by Bona et al. (2003).

See Pretorius (2005b).

This name is somewhat misleading, however, since it does not single out any particular family of coordinate systems.

Instead, any arbitrary coordinate system can be generated by (11.22) with a suitable choice of H

11.3 Generalized harmonic coordinates 383

quantities C

, which may not be identically zero due to numerical errors.

We will discuss

why adding constraints can stabilize evolution equations in Section 11.5.

The generalized harmonic approach to determing the metric in a dynamical spacetime

differs in several ways from the 3 +1 formalism, so that it is worthwhile to discuss

some aspects of this formalism in some more detail. Most importantly, equation (11.23)

is a second-order equation in time for the spacetime metric, while in the usual 3 +1

decomposition we integrate coupled equations for spatial metric and the extrinsic curvature

that are ﬁrst order in time. This difference effects both the initial data and the numerical

implementation.

Within the 3 + 1 decomposition, a set of initial data consists of values of the spatial

metric γ

and the extrinsic curvature K

that satisfy the constraint equations, e.g., (2.132)

and (2.133), at one instant of time. Given a choice for the lapse α and the shift β

, γ

and

can then be integrated forward in time with the evolution equations, e.g., (2.134)and

(2.135). Equation (11.23), on the other hand, requires the spacetime metric g

and its ﬁrst

time derivative at some instant of time t as initial data.

One way of constructing such initial data is the following. We could ﬁrst solve the

contraint equations (2.132) and (2.133), in, say, the conformal thin-sandwich formalism

described in Chapter 3.3. The freely speciﬁable variables then are the conformally related

metric ¯γ

and its time derivative, together with the trace of the extrinsic curvature K and

its time derivative. Solving the equations yields the conformal factor ψ, the lapse α and

the shift β

. This is all the information required to construct the spacetime metric g

(e.g., equation 2.131). To ﬁnd the time derivative of g

we can ﬁrst evaluate the evolution

equation (2.134), which yields ∂

at t = 0. We can ﬁnd the time derivatives of α and

from the condition

(4)



= H

, again evaluated at t = 0. For the special case H

= 0

this condition yields equations (4.44) and (4.45), which we can solve for the desired time

derivatives (see exercise 4.9). If H

= 0 these equations have additional source terms, but

the derivation is very similar. From the time derivatives of the spatial metric, the lapse

and the shift we can ﬁnally construct the time derivative of the spacetime metric, which

completes the initial data for equation (11.23).

The appearance of second-order time derivatives in equation (11.23) also poses some

numerical challenges. In ﬁnite difference applications these second derivatives can be

handled by choosing a three-level scheme and using a ﬁnite difference representation

similar to the one described by equation (6.22). The situation is complicated by the

presence of the mixed space-time derivatives, which couple the function values at the

new time level. To avoid an implicit scheme, the resulting system of equations can be

solved iteratively.

Alternatively, equation (11.23) can be recast into a system of coupled

ﬁrst-order equations.

To stabilize his simulations of binary black holes, Pretorius (2005a) added a combination of these constraints to the

evolution equations (11.23), as suggested by Gundlach et al. (2005).

Pretorius (2005b).

Lindblom et al. (2006).

384 Chapter 11 Recasting the evolution equations

Another important difference between this approach and 3 + 1 decompositions is the

way in which coordinates are imposed. In the 3 + 1 decomposition, we choose a coordinate

system with the help of the lapse α and the shift β

. These quantities are directly related

to the geometry of the spatial slices , and in Chapter 4 we have seen how geometric

considerations can guide the search for coordinate systems with particularly desirable

properties. In this “generalized harmonic coordinate” approach, on the other hand, the

coordinates are imposed by the gauge source functions H

. These quantities do not have a

direct geometrical meaning, and it is therefore much less clear how to construct a coordinate

system with desirable properties. Pretorius (2005a) chose H

to satisfy a wave equation

that includes the lapse α (computed from g

) as a source, and H

= 0. A priori there is

not much reason to expect this choice to lead to a well-behaved coordinate system, but in

binary black hole merger simulations it evidently does.

While the wave-like structure of equation (11.23) is very appealing, the uncertainty

over how to impose suitable coordinate conditions has served to discourage many

researchers from using this formalism. Instead, a number of other formulations of

Einstein’s equations have been constructed with similarly desirable mathematical proper-

ties as equation (11.23) but which, by contrast, impose the gauge conditions via a lapse and

shift.

11.4 First-order symmetric hyperbolic formulations

As we discussed in Section 11.1, it could prove quite desirable to cast the evolution

equations into a ﬁrst-order hyperbolic form.

The ﬁrst such formulations were based on

harmonic coordinates

and used as a starting point the formalism of Section 11.3 with

= 0. The ﬁrst formulations that departed from the assumption of harmonic coordinates

were based on a spin-frame formalism.

Other formulations introduced partial derivatives

of the metric and other quantities as new independent variables.

In analogy to the electro-

magnetic example in Section 11.2.2, it is also possible to take a time derivative of equation

(11.8) to derive what is sometimes referred to as the “Einstein–Ricci” system of equa-

tions.

Another system of equations, sometimes called the “Einstein–Bianchi” system,

can be derived from the Bianchi identities.

The “Einstein–Christoffel” system is con-

structed by introducing additional “connection” variables.

We also mention a so-called

“λ”-system that embeds Einstein’s equations into a larger symmetric hyperbolic system

with the constraint surface of Einstein’s equations acting as an attractor of the evolution.

See Reula (1998) for a more extensive survey than provided in this section.

Choquet-Bruhat (1952, 1962); Fischer and Marsden (1972).

Friedrich (1981, 1985)

Bona and Mass

o (1992); Frittelli and Reula (1994).

Choquet-Bruhat and Ruggeri (1983); Abrahams et al. (1995); Abrahams and York, Jr., (1997).

Friedrich (1996); Anderson et al. (1997).

Anderson and York, Jr. (1999).

Brodbeck et al. (1999).

11.4 First-order symmetric hyperbolic formulations 385

It is beyond the scope of this volume to review all of these formulations. Some have

features that are not very desirable numerically, in that they restrict the gauge freedom,

introduce extra derivatives of the matter variables, or introduce a large number of auxil-

iary variables. Only few of these formulations have been implemented numerically, and

many of these implementations adopted simplifying symmetry conditions (e.g., spherical

symmetry). Some of these implementations display advantages over the standard 3 + 1

formalism,

but others reveal additional problems. For example, a particular equation in

the “Einstein–Ricci” system turns out to produce an exponentially growing mode, which

can be removed in spherical symmetry, but not in more general 3 + 1 simulations.

Only

a few of these systems have been implemented in 3 + 1 dimensions, including the one of

Bona and Mass

o (1992), and versions of the “Einstein–Christoffel” system of Anderson

and York, Jr. (1999).

Since this “Einstein–Christoffel” formulation is particularly ele-

gant, we provide a brief summary of this system as an example of a ﬁrst-order hyperbolic

formulation.

Starting with the standard 3 + 1 or ADM formalism of Chapter 2 we deﬁne the new

variables

kij

≡ 

(ij)k

+ γ



[lj]m

+ γ



[li]m

. (11.25)

These functions are now promoted to independent functions. It can then be shown that the

evolution equations (11.7) and (11.8) can be rewritten as the coupled system

=−2α K

+ αγ

∂

kij

= α M

kij

+ α∂

= α N

kij

(11.26)

Here we have introduced the abbreviation

≡ ∂

− L

, (11.27)

and the source terms M

and N

ijk

given by

= γ

− 2K

) + γ

(4 f

kmi

[ln] j

+4 f

km[n

l]ij

− f

ikm

jln

+ 8 f

(ij)k

[ln]m

+ 4 f

km(i

j)ln

−8 f

kli

mnj

+ 20 f

kl(i

j)mn

− 13 f

ikl

jmn

)

−∂

∂

ln ¯α − (∂

ln ¯α)(∂

ln ¯α) + 2γ

( f

kmn

∂

ln ¯α

− f

kml

∂

ln ¯α) + γ



(2 f

(ij)k

− f

kij

)∂

ln ¯α

+4 f

kl(i

∂

ln ¯α − 3( f

ikl

∂

ln ¯α + f

jkl

∂

ln ¯α)



−8π S

+ 4πγ

(11.28)

e.g., Bona and Mass

o (1992).

Scheel et al. (1997, 1998).

The latter have been implemented numerically by Kidder et al. (2001, see also Kidder et al. (2000)) using spectral

methods.

Here we adopt the notation of Kidder et al. (2000).

386 Chapter 11 Recasting the evolution equations

and

kij

= γ



k(i

j)mn

− 4 f

mn(i

j)k

+ K

(2 f

mnk

− 3 f

kmn

)



+2γ



(γ

k(i

j)qn

− 2 f

qn(i

j)k

)

+γ

k(i

j)m

(8 f

npq

− 6 f

pqn

) + K

(4 f

pq(i

j)k

− 5γ

k(i

j) pq

)



− K

∂

ln ¯α + 2γ

m(i

j)k

∂

ln ¯α − K

k(i

∂

ln ¯α)

+16πγ

k(i

(11.29)

We have also used the “densitized” lapse function

¯α = γ

−1/2

α (11.30)

(see equation 3.100) and, in addition to the matter source terms deﬁned in Chapter 2,the

4-dimensional trace of the stress energy tensor,

T = g

. (11.31)

The ﬁrst-order, symmetric hyperbolic (“FOSH”) system (11.26) is equivalent to the

original set of evolution equations (11.7) and (11.8). Since the f

kij

are evolved as indepen-

dent functions, the deﬁning relations (11.25) can be considered as a new set of constraint

equations in addition to the usual ones given by equations (2.132) and (2.133). Note also

that the source terms M

and N

ijk

on the right-hand sides do not contain any derivatives of

the fundamental variables (other than of the arbitrary lapse function ¯α). Equations (11.26)

can be combined to yield a wave equation for the components of the spatial metric γ

which the right-hand sides appear as source terms.

Evolutions of a single black hole using a spectral implementation of this system are still

unstable, but the lifetime of these simulation can be extended to late times by a general-

ization of the equations.

This generalization involves a redeﬁnition of the independent

variables and the addition of new constraints. The above equations can be embedded into

a 12-parameter family of strongly hyperbolic formulations. The stability properties of the

system depend sensitively on the choice of the free parameters, which can be understood

analytically in terms of energy arguments.

Given this need for ﬁne-tuning, and given the

successes of both the generalized coordinate formalism of Section 11.3 and the BSSN for-

mulation discussed below, this “Einstein–Christoffel” system has lost some of its appeal.

11.5 The BSSN formulation

In our illustration of Section 11.2 employing electromagnetism we have seen that we can

also eliminate the “mixed second derivatives” in Maxwell’s equations (11.10) by promot-

ing the contraction  deﬁned in equation (11.11) to be a new, independent function. In

Kidder et al. (2001).

See, e.g., Lindblom and Scheel (2002). Also see exercise 11.4 below and the related discussion, which illustrate how

the addition of constraints affects the properties of evolution systems.

11.5 The BSSN formulation 387

Sections 4.3 and 11.3 we have seen similarly that we can absorb the mixed second deriva-

tives in the Ricci tensor with the help of the connection functions (11.21). The BSSN

formalism

adopts a similar strategy to simplify the three-dimensional, spatial Ricci ten-

sor.

In addition, the conformal factor and the trace of the extrinsic curvature are evolved

separately in the BSSN formalism, which follows the philosophy of separating transverse

from longitudinal, or, equivalently, radiative from nonradiative, degrees of freedom. The

later idea was also the basis of the conformal decompositions discussed in Chapter 3 for

the construction of initial data.

To derive this formulation, we begin by writing the conformal factor ψ as ψ = e

that we have

¯γ

= e

−4φ

. (11.32)

We then require that the determinant of the conformally related metric ¯γ

be equal to that

of the ﬂat metric η

in whatever coordinate system we are using, i.e.,

φ =





. (11.33)

In the following we will adopt a Cartesian coordinate system, so that ¯γ = η = 1.

As in equation (3.31) we split off from the extrinsic curvature its trace and conformally

rescale the remaining traceless piece A

. However, we choose a conformal rescaling that

is different from equation (3.35) and, instead, rescale A

the same way as we did for the

metric itself,

= e

−4φ

. (11.34)

We will use tildes as opposed to the bars used in Chapter 3 to distinguish between these

different rescalings. Indices of

will be raised and lowered with the conformal metric

¯γ

, so that

= e

4φ

Evolution equations for φ and K can now be found by taking the trace of the evolution

equations (2.136) and (2.137), which yields

∂

φ =−

αK + β

∂

φ +

∂

(11.35)

and

∂

K =−γ

α + α(

) + 4πα(ρ + S) + β

∂

K. (11.36)

Subtracting these equations from the evolution equations (11.7) and (11.8)leavesthe

traceless parts of the evolution equations for ¯γ

and

according to

∂

¯γ

=−2α

+ β

∂

¯γ

+ ¯γ

∂

+ ¯γ

∂

−

¯γ

∂

, (11.37)

Shibata and Nakamura (1995); Baumgarte and Shapiro (1998).

See also Nakamura et al. (1987).

388 Chapter 11 Recasting the evolution equations

and

∂

= e

−4φ



−(D

α)

+ α(R

− 8π S

)



+ α(K

− 2

)

+β

∂

−

∂

(11.38)

In the last equation, the superscript TFdenotes the trace-free part of a tensor, e.g., R

− γ

R/3. Note that in equations (11.35) through (11.38) the shift terms arise from

Lie derivatives

of the respective evolution variable appearing on the left-hand side. The

divergence of the shift, ∂

, appears in the Lie derivative because the choice ¯γ = 1 makes

φ a tensor density of weight 1/6, and ¯γ

and

tensor densities of weight −2/3(see

Section A.3).

Exercise 11.2 Derive equations (11.35) through (11.38).

According to equation (3.10) we can split the Ricci tensor into two terms

+ R

, (11.39)

where only R

depends on the conformal function φ. We can identify the form of R

inserting φ = ln ψ into equation (3.10). We could compute the conformally related Ricci

tensor

by inserting ¯γ

into equation (2.143), but that would again introduce the mixed

second derivatives that we are trying to avoid. Analogously to the way we introduced a new

variable  in equation (11.11) to eliminate the mixed derivatives in Maxwell’s evolution

equations, we can now deﬁne “conformal connection functions”



≡ ¯γ



=−∂

¯γ

(11.40)

to accomplish the same task in the above evolution equations for the gravitational ﬁeld.

Here the



are the connection coefﬁcients associated with ¯γ

, and the last equality holds

in Cartesian coordinates when ¯γ = 1. In terms of these conformal connection functions

we can now write the Ricci tensor as

=−

¯γ

∂

¯γ

+ ¯γ

k(i

∂



(ij)k

+ ¯γ





l(i



j)km



klj



. (11.41)

The only explicit second-derivative operator acting on ¯γ

in this expression involves

a Laplacian, ¯γ

∂

– all other second derivatives are absorbed in ﬁrst derivatives of



. This derivative structure is similar to the one we encountered in the 4-dimensional

Ricci tensor in Sections 4.3 and 11.3, where the only explicit second derivatives on

the metric appearing in that quantity managed to form a convenient wave operator

(d’Alembertian).

In the generalized harmonic formalism of Section 11.3 we set the 4-dimensional analogs

(4)



of the conformal connection functions



equal to some speciﬁed gauge source func-

tion. Here we do not follow this approach but instead promote the



to new independent

functions, in complete analogy to our treatment of Maxwell’s equations in Section 11.2.3.

Adopting this approach requires us to derive separate evolution equations for the



.By

11.5 The BSSN formulation 389

analogy with the derivation of equation (11.17) we interchange a partial time and space

derivative in the deﬁnition (11.40) to obtain

∂



=−∂



2α

− 2¯γ

m( j

∂

¯γ

∂

+ β

∂

¯γ



. (11.42)

We can now eliminate the divergence of the extrinsic curvature with the help of the

momentum constraint (2.133), which then yields the desired evolution equation,

∂



=−2

∂

α + 2α





−

¯γ

∂

K − 8π ¯γ

+ 6

∂



+β

∂



−



∂



∂

¯γ

∂

+ ¯γ

∂

(11.43)

Equations (11.35) through (11.38), together with (11.43), form a new system of evolution

equations that is equivalent to equations (11.7) and (11.8). Since the



are evolved as

independent functions, the deﬁning relation (11.40) serves as a new constraint equation, in

addition to equations (2.132) and (2.133). We summarize the resulting BSSN formalism

in Box 10.1.

Exercise 11.3 Show that the shift terms in equation (11.43) arise from the Lie

derivative of



along β



= β

∂



−



∂



∂

¯γ

∂

+ ¯γ

∂

. (11.44)

Hint: First show that the conformal connection coefﬁcients transform according to





= J

−W

∂x



∂x



+ J

−W

¯γ

∂x



∂x



∂x



∂x

∂

∂x



∂x



−

−W

¯γ

∂x



∂x

∂

(ln J

(11.45)

where J is the Jacobian of the transformation and W is the weight of ¯γ

.Theﬁrst

term is the usual transformation term for a tensor (except for the Jacobian factor), the

second term arises because the connection coefﬁcients do not transform like tensors,

and the third term appears because ¯γ

is a tensor density. Then use deﬁnition (A.7)

to ﬁnd equation (11.44)forW =−2/3(cf. equations A.36 and A.37).

While the different formulations are equivalent analytically, the difference in perfor-

mance of numerical implementations is striking. In Figure 11.1 we return to the example

of linear gravitational wave propagation cited at the beginning of this chapter. In this exam-

ple,

a small amplitude, time-symmetric, even-parity l = 2, m = 0 gravitational wave of

the kind discussed in Chapter 9.1.2 is evolved with harmonic slicing (see Chapter 4.3), zero

shift, and a simple outgoing wave boundary condition. Both the standard 3 + 1 and BSSN

systems give very similar results early on, but the standard system crashes very soon,

while the BSSN system remains stable. Similar improvements have been found for many

other applications that employ a BSSN scheme, including the propagation of nonlinear

Baumgarte and Shapiro (1998); see also a similar example in Shibata and Nakamura (1995).

390 Chapter 11 Recasting the evolution equations

Box 11.1 The BSSN equations

In the BSSN formulation of the 3 + 1 equations the spatial metric γ

is decomposed into a

conformally related metric ¯γ

with determinant ¯γ = 1 (assuming Cartesian coordinates) and a

conformal factor e

= e

4φ

¯γ

. (11.46)

We also decompose the extrinsic curvature into its trace and traceless parts and conformally

transform the traceless part as we do the metric,

= e

4φ

K. (11.47)

In terms of these variables the Hamiltonian constraint (2.132) becomes

0 =

H = ¯γ

−

R +

5φ

−

5φ

+ 2πe

5φ

ρ, (11.48)

while the momentum constraint (2.133) becomes

0 =

6φ

) −

6φ

K − 8π e

6φ

. (11.49)

The evolution equation (2.136)forγ

splits into two equations,

∂

φ =−

αK + β

∂

φ +

∂

, (11.50)

∂

¯γ

=−2α

+ β

∂

¯γ

+ ¯γ

∂

+ ¯γ

∂

−

¯γ

∂

, (11.51)

while the evolution equation (2.135)forK

splits into the two equations

∂

K =−γ

α + α(

) + 4πα(ρ + S) +β

∂

K, (11.52)

∂

= e

−4φ



−(D

α)

+ α(R

− 8π S

)



+ α(K

− 2

)

+β

∂

−

∂

(11.53)

In the last equation the superscript TF denotes the trace-free part of a tensor, e.g., R

− γ

R/3. We also split the Ricci tensor into R

+ R

,whereR

can be found by

inserting φ = ln ψ into equation (3.10). We express

in terms of the conformal connection

functions



≡ ¯γ



=−∂

¯γ

, which yields

=−

¯γ

∂

¯γ

+ ¯γ

k(i

∂



(ij)k

+ ¯γ





l(i



j)km



klj



. (11.54)

The



are now treated as independent functions that satisfy their own evolution equations,

∂



=−2

∂

α + 2α





−

¯γ

∂

K − 8π ¯γ

+ 6

∂



+β

∂



−



∂



∂

¯γ

∂

+ ¯γ

∂

(11.55)