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

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

10.4 Collapse of rotating toroidal clusters 371

event horizon as a toroid.

The initial matter distribution is based on a solution for a

relativistic toroidal cluster in dynamical equilibrium. A discussion of the construction of

initial data for rotating equilibrium systems, given a phase space distribution function, is

postponed to Chapter 14.1.3. The adopted distribution function for a toroidal cluster is

chosen to be of the form

f (E, J

) = g(E)h(J

), (10.18)

where E =−u

is the energy and J

= u

is the angular momentum of a particle, both

per unit mass. Toroidal clusters with net spin then can be generated most simply by taking

g(E) = K δ(E − E

max

), (10.19)

h( J

) = δ(J

− J

). (10.20)

Here the total angular momentum J = J

ADM

is related to J

by J = (M

/m)J

,where

is the total rest-mass of the conﬁguration and m is the rest-mass of one particle. The

quantity K is an arbitrary normalization constant and E

max

is chosen to be the maximum

energy of a particle in a spherical cluster of areal radius R

, E

max

= (1 − 2M/R

)

1/2

Here we will discuss the results for a toroidal cluster with an outer circumferential radius

of R

/M = 4.5. Collapse was induced by multiplying the angular velocity of each particle

by a factor of 0.5, resulting in a nonequilibrium cluster with a total angular momentum

of J/M

= 0.65. Re-solving the constraint equations after reducing the angular velocities

to get valid initial data, the collapse was followed numerically. The computational mesh

consisted of 200 radial zones and 16 angular zones (for one hemisphere) and 3000 particles.

The outer boundary of the mesh was placed at 50M, well outside the matter source.

The toroid initially collapses along the rotation axis to a thin hoop. Then, while undergo-

ing oscillations along the rotation axis, it collapses radially inwards. The ﬁnal conﬁguration

is a stationary Kerr black hole containing all of the matter. The matter distribution and

horizons are plotted at selected times in Figures 10.8 and 10.9. The event horizon is

determined by tracing null rays backwards in time from the end of the simulation at

t/M = 23.2, at which point the apparent horizon coincides with the event horizon.

The

topology of the event horizon is rather remarkable: it initially develops as a toroid, begin-

ning at t/M = 13.2. The event horizon is seen to form initially entirely in the vacuum

between the origin and the inner edge of the collapsing toroidal cluster. It then expands

to ﬁll up the the “doughnut hole”, becoming topologically spherical at t/M ≈ 13.5. At

this instant the outer edge of the event horizon has reached the inner edge of the matter

toroid. The spactime diagram for the collapse looks quite similar to the diagram plotted

in Figure 10.5. The line of crossover points at which light rays enter the horizon, and the

cusp formed by rays at the point at which the line of crossovers terminates, are all present

as in Figure 10.5.

Abrahams et al. (1994); Hughes et al. (1994).

Hughes et al. (1994); see the discussion in Chapter 7.2.

372 Chapter 10 Collapse of collisionless clusters in axisymmetry

Figure 10.8 Snapshots of the collapse of a rotating toroid of collisionless particles at selected times. The left

panels display meridional slices, the right panels show the equatorial plane. The solid line denotes the event

horizon, the dashed line the apparent horizon. The earliest times shown occur soon after the appearance of the

toroidal event horizon. By the ﬁnal time shown here, the event and apparent horizons coincide. The coordinate

labels (t, x , y, z) are in units of M.[FromHughes et al. (1994).]

Of course, some aspects of the collapse noted above are gauge-dependent; the collapse

could look quite different in a different time-slicing, for example. Nevertheless, when

the results of the numerical simulation were ﬁrst reported, it appeared that they might be

in conﬂict with black hole theorems regarding “topological censorship”. These theorems

permit a nonstationary black hole to have the topology of a 2-sphere or a torus,

but a

torus can persist only for a short time. According to topological censorship,

the hole in a

Gannon (1976).

Friedman et al. (1993); Jacobson and Venkataramani (1995).

10.4 Collapse of rotating toroidal clusters 373

Figure 10.9 Three-dimensional views of the collapse of the rotating toroid shown in Figure 10.8. The images

are viewed from 45

◦

above the equatorial plane. The ﬁrst image shows the initial conﬁguration. The remaining

three images are at the same times as the last three images in Figure 10.8. The outermost shaded region is the

event horizon, the shaded region inside it is the apparent horizon. The scale of the last three images is larger by a

factor of 2 over the ﬁrst image. The clock on the lower right denotes the fraction of time that has elapsed since the

onset of the simulation. [From Hughes et al. (1994).]

374 Chapter 10 Collapse of collisionless clusters in axisymmetry

toroidal horizon must close up quickly, before a light ray can pass through. Moreover, the

black hole must be a 2-sphere if no null generators enter the horizon at later times.

Once the numerical simulations were studied in detail and the geometry of the transient

toroidal event horizon was analyzed by means of a spacetime diagram like Figure 10.5,the

results were shown to be completely consistent with the theorems concerning the topology

of black hole horizons.

At late times, when equilibrium has been reached, the topology

is spherical, in agreement with the theorem of Hawking.

At early times, the topology is

temporarily toroidal. However, the line of crossovers traced out by a point on the inner rim

of the torus is spacelike. That implies that the “hole” in the torus indeed closes up faster

than the speed of light, in compliance with topological censorship. Finally, at intermedidate

times, when the horizon is spanned by its full complement of null generators, the rotating

black hole has a spherical topology, in accord with the theorems of Browdy and Galloway.

Browdy and Galloway (1995).

Shapiro et al. (1995).

Hawking (1972).

Recasting the evolution

equations

At this point we might suspect that we are all set to carry out dynamical simulations involv-

ing relativistic gravitational ﬁelds in three spatial dimensions, i.e., dynamical simulations

in full 3 + 1 dimensional spacetimes. After all, we have derived the 3 + 1 evolution equa-

tions for the gravitational ﬁelds in Chapter 2, developed techniques for the construction

of initial ﬁeld data in Chapter 3, and discussed strategies for imposing suitable coordi-

nate conditions in Chapter 4. Should the spacetime in which we are interested contain

the most common matter sources, we can consult Chapter 5 for the matter source terms

and the matter equations of motion. Beyond assembling all of the relevant equations, we

have also sketched algorithms for solving them numerically in Chapter 6. If black holes

are present, we have derived methods to locate and measure their horizons in Chapter 7;

if gravitational waves are generated we have discussed how to extract radiation wave-

forms numerically in Chapter 9. Finally, we have evolved relativistic systems in spherical

symmetry (1 + 1 dimensions) in Chapter 8 and in axisymmetry (2 + 1 dimensions) in

Chapter 10. What else is left to do before going on to perform simulations in full 3 + 1

dimensions?

Suppose then we were to plunge ahead with the tools currently at our disposal to

build a dynamical spacetime in 3 + 1 dimensions. We could gain some conﬁdence and

computational experience by choosing to explore as our ﬁrst case a simple dynamical

system with a known analytic solution. A linear gravitational wave propagating in vacuum

would do very nicely. We could grab the initial data for such a wave from Chapter 9.1.2,

choose a reasonable set of coordinate conditions from the sample menu provided in

Chapter 4 and employ the standard 3 + 1 (ADM) evolution equations derived in Chapter 2

to evolve the wave.

What would we ﬁnd? Unfortunately, despite all of our methodical preparation and

best numerical implementation, most likely our code would crash after a rather short

time! The problem is illustrated in Figure 11.1, which shows the failure of our standard

evolution scheme to evolve such a linear wave in a stable fashion. What might we con-

clude? Simply the following: since even this simple, linear, vacuum spacetime encounters

numerical difﬁculties when evolved with the standard set of 3 + 1 equations, we can

expect that most other dynamical spacetime in 3 + 1 dimensions will encounter such

difﬁculties.

As it turns out, the standard evolution equations as derived in Chapter 2 are not yet in a

form that is suitable for stable numerical integration. Most unconstrained simulations that

375

376 Chapter 11 Recasting the evolution equations

implement these equations, at least in three spatial dimensions, turn out to be unstable.

The

failure of these equations can be understood in terms of their mathematical properties, as we

will brieﬂy discuss in Section 11.1. Following this summary, we will discuss reformulations

of the evolution equations that avoid these problems and have proven more robust and

successful numerically. But ﬁrst, in order to model the shortcomings of the standard

evolution equations of Chapter 2 in a simple way and to obtain some guidance on how

to ﬁx them, we shall return to Maxwell’s equations in Section 11.2. Fortiﬁed by the

insight provided by this example, we will then present several different reformulations

of the 3 + 1 evolution equations that have proven very successful in many numerical

simulations.

11.1 Notions of hyperbolicity

Recall from Chapter 6 that we can cast a hyperbolic partial differential equation – for

example the simple wave equation (6.3) – in the ﬁrst-order form

∂

u + A · ∂

u = S, (11.1)

where u is a solution vector, S = S(u) is a source vector, and where we have called the

matrix A the velocity matrix (see equation 6.7). Here we have assumed that the solution

depends on t and x only; in more than one spatial dimension we may generalize equation

(11.1) to read

∂

u + A

· ∂

u = S. (11.2)

If the solution vector u has n components, then each matrix A

is an n × n matrix. For the

purposes of this section we may ignore the source term and set S = 0.

We call a problem well-posed if we can deﬁne some norm ... so that the norm of the

solution vector satisﬁes

u(t, x

)≤ke

αt

u(0, x

) (11.3)

for all times t ≥ 0. Here k and α are two constants that are independent of the initial data

u(0, x

). Stated differently, solutions of a well-posed problem cannot increase more rapidly

than exponentially. Clearly, this is a very desirable property, but it is not guaranteed for all

hyperbolic systems.

As it turns out, there exist different types of hyperbolicity, and only for some of these

types are the equations well-posed. To analyze these hyperbolicity properties, we consider

an arbitrary unit vector n

and construct the matrix P = A

, which is sometimes referred

The high degree of spatial symmetry in spherical symmetry and even axisymmetry allows for a number of different

strategies to obtain stable evolutions, including fully or partially constrained evolution, the use of special coordinate

systems, or the introduction of new variables.

Much of this section follows the arguments of Alcubierre (2008a); see also Kreiss and Lorenz (1989)andAlcubierre

(2008b).

11.1 Notions of hyperbolicity 377

to as the principal symbol or the characteristic matrix of the system. Based on the properties

of P we then distinguish different notions of hyperbolicity. In particular, we call the system

symmetric hyperbolic if P can by symmetrized in a way that is independent of n

strongly hyperbolic if, for all unit vectors n

, P has real eigenvalues and a complete set

of eigenvectors,

and weakly hyperbolic if P has real eigenvalues but not a complete set of eigenvectors.

Simple wave equations, for example, are symmetric hyperbolic. Symmetric hyperbolic

systems are automatically strongly hyperbolic. The key result for our purposes states

that strongly hyperbolic equations are well-posed, while systems that are only weakly

hyperbolic are not.

We are not quite ready yet to analyze the hyperbolicity of the ADM evolution equations

(2.134) and (2.135), since our analysis above applies to ﬁrst-order equations, while the

ADM equations are second order in space. One approach to deal with this is to introduce a

so-called ﬁrst-order reduction of the equations, achieved by writing all second derivatives

in terms of ﬁrst derivatives of a new set of auxiliary functions that contain ﬁrst derivatives

of the original variables.

An analysis of a ﬁrst-order reduction of the ADM evolution

equations shows that these equations are in fact only weakly hyperbolic.

As a consequence,

the evolution problem is not well-posed, and we have no reason to expect the solutions –

or numerical implementations – to be well-behaved.

Clearly, then, we should try to recast these equations in a form that is strongly hyperbolic.

How can we do that without modifying Einstein’s equations? The answer lies in the

constraint equations. We ﬁrst observe that the constraints vanish, at least analytically.

are therefore entitled to add multiples of the constraint equations to the evolution equations.

Furthermore, the constraints contain up to second derivatives of the gravitational ﬁelds,

as do the evolution equations. Adding the constraints to the evolution equations therefore

affects the appearance of the highest order derivatives, which in turn affects the principal

symbol P and hence the hyperbolicity. We will see that many reformulations of the evolution

equations also involve the introduction of new variables that absorb some of the ﬁrst

derivatives of the gravitational ﬁelds.

Before proceeding we add a word of caution. The notion of well-posedness rules out

modes that grow faster than exponentially with time.

From a numerical perspective,

Note that different authors use slightly different conventions.

See, e.g., Kreiss and Lorenz (1989). Note that this analysis focuses on the principal part of the equations only; the

source term may also lead to faster-than-exponential growth.

For an alternative approach see, e.g., Gundlach and Martin-Garcia (2006a), who introduce a notion of hyperbolicity

that applies to systems that are second order in space.

See, e.g., Kidder et al. (2001).

Here we mean that if we move the source terms appearing in equations (3.1)and(3.2) to the left-hand sides, the

combination of terms on the left-hand sides must then add up to zero at all times; see, e.g., equations (11.48)and

(11.49).

At least in the absence of source terms; see footnote above.

378 Chapter 11 Recasting the evolution equations

however, exponentially growing modes are still very bad, and can easily terminate a

simulation after only a short time. In fact, many of the symmetric or strongly hyperbolic

systems that have been introduced over the years lead to such exponentially growing modes

in numerical implementations. Unless these modes can be controlled,

these systems are not

very useful for numerical simulations. From a numerical perspective, then, well-posedness

is a necessary but not a sufﬁcient condition.

In the following we will discuss several different approaches that lead to different

reformulations of the evolution equations that are both strongly hyperbolic and that have

been used successfully in numerical simulations. In the remainder of this chapter we

will pursue a more heuristic approach than the one sketched in this section and will

marshal mathematical intuition and exploit analogies to motivate the different approaches.

In fact, some of the systems discussed below were developed in exactly this way and were

discovered to have desirable numerical properties empirically before they were analyzed

mathematically and revealed to be strongly hyperbolic. To illustrate some of these heuristic

arguments for a simple and familiar example we will begin by discussing Maxwell’s

equations in Section 11.2.

11.2 Recasting Maxwell’s equations

In Chapter 2.2 we have seen that we can bring Maxwell’s evolution equations in a

Minkowski spacetime into the form

∂

=−E

− D

 (11.4)

∂

=−D

+ D

− 4π j

(11.5)

(see equations 2.11 and 2.12). Recall that A

is the three-vector potential, the magnetic

ﬁeld satisﬁes B

= 

ijk

and therefore is automatically divergence-free, φ is a gauge

potential, and the electric ﬁeld E

has to satisfy the constraint equation (2.5),

= 4πρ

. (11.6)

In the above equations ρ

is the electric charge density and j

the current density.

In Chapter 2.7 we have discussed some of the similarities of Maxwell’s equations in the

above form with the ADM evolution equations (2.134) and (2.135), namely

∂

=−2α K

+ D

(11.7)

and

∂

= α(R

− 2K

+ KK

) − D

α − 8πα(S

−

(S − ρ))

+β

∂

+ K

∂

+ K

∂

(11.8)

See, e.g., Scheel et al. (1998); Kidder et al. (2001) for examples.

11.2 Recasting Maxwell’s equations 379

If we identify the vector potential A

with the spatial metric γ

and the electric ﬁeld E

with the extrinsic curvature K

, we see that the right-hand sides of both equations (11.4)

and (11.7) contain a ﬁeld variable and a spatial derivative of a gauge variable, while the

right-hand sides of both equations (11.5) and (11.8) involve matter sources as well as

second spatial derivatives of the second ﬁeld variable. In equation (11.8) these second

derivatives are hidden in the Ricci tensor R

, which we can write, for example, as



∂

+ ∂

∂

− ∂

∂

− ∂

∂



+ γ





mkj

− 



mkl



. (11.9)

We can now exploit these similarities by focusing on the simpler Maxwell system of

equations to identify some of the computational shortcomings of these forms of the

evolution equations.

We ﬁrst note that equations (11.4) and (11.5) almost can be combined to yield a wave

equation, which would make the system symmetric hyperbolic. To see this, take a time

derivative of equation (11.4) and insert equation (11.5) to form a single equation for the

vector potential A

−∂

+ D

− D

= D

∂

 − 4π j

. (11.10)

On the left-hand side, the second time derivative combines with the Laplace operator

to form a wave operator (d’Alembertian). Equations (11.4) and (11.5) would then

constitute a wave equation for the components A

if it weren’t for the mixed derivative

term D

In general relativity the situation is very similar. The Ricci tensor R

on the right hand

side of equation (11.8) contains three mixed derivative terms in addition to the term with a

Laplace-like operator acting on γ

, i.e., γ

∂

. Without these mixed derivative terms

the standard ADM equations could be written as a set of wave equations for the components

of the spatial metric, which would make them symmetric hyperbolic. As we discussed in

Section 11.1 we would like to evolve a system that is well-posed, but the presence of these

additional terms may spoil this property.

These considerations suggest that it would be desirable to eliminate the mixed derivative

terms. In electrodynamics, three different approaches can be taken to eliminate the D

term: one can make a special gauge choice; one can bring Maxwell’s equations into a ﬁrst-

order symmetric hyperbolic form; or, one can introduce an auxiliary variable.Inthe

remainder of this section we will discuss brieﬂy each one of these three strategies.

11.2.1 Generalized Coulomb gauge

The most straightforward approach to eliminating the undesirable terms is to choose a

gauge so that the term D

disappears. Deﬁne the quantity

 ≡ D

. (11.11)

380 Chapter 11 Recasting the evolution equations

To eliminate the offending term, we may then set

 = 0. (11.12)

This choice is the familiar Coulomb gauge condition.

Exercise 11.1 Show that in the Coulomb gauge the gauge potential  satisﬁes

Poisson’s elliptic equation

 =−4πρ

. (11.13)

In general relativity, an analogous approach can be taken by choosing harmonic coor-

dinates, which bring the equations into the form of a wave equation.

(The reader may

wish to review Chapter 4.3,and,especially,exercise4.10.) We will return to this strategy

in Section 11.3 below.

This approach has disadvantages, too. As exercise 11.1 demonstrates, the choice (11.12)

forces us to satisfy a constraint that may be inconvenient to solve. In general relativity, the

analogous harmonic gauge condition,

(4)



= 0, similarly imposes a coordinate choice –

the lapse and shift now have to satisfy equations (4.44) and (4.45). These coordinates may

not be well-suited to the problem at hand, and may lead to coordinate singularities. One

way to regain coordinate freedom is not to set  = 0asinequation(11.12), but instead to

set  equal to some yet-to-be-chosen gauge source function H (t, x

 = H (t, x

). (11.14)

Retracing the steps in exercise 11.1 we could ﬁnd out how a particular choice of H affects

the gauge potential . Similarly, in general relativity we can relax the harmonic gauge

condition with the help of a 4-dimensional gauge source function,

(4)



= H

(t, x

), as we

will see in Section 11.3. This approach is commonly referred to as “generalized harmonic

coordinates”, and, in analogy, we could refer to the condition (11.14) imposed on  as the

“generalized Coulomb gauge”.

11.2.2 First-order hyperbolic formulations

An alternative, gauge-covariant approach to bringing Maxwell’s equations into a symmetric

hyperbolic form is to take a time derivative of equation (11.5) instead of equation (11.4),

which then yields

∂

= D

(−E

− D

) − D

(−E

− D

) − 4π∂

. (11.15)

Using the constraint (11.6) we can eliminate the ﬁrst term and ﬁnd a wave equation for E

−∂

+ D

= 4π (∂

+ D

), (11.16)

This property was ﬁrst realized by De Donder (1921)andLanczos (1922). Many of the early hyperbolic formulations

of Einstein’s equations were based on this gauge choice, e.g., Choquet-Bruhat (1952, 1962); Fischer and Marsden

(1972).