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

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

4.5 Minimal distortion and variations 121

Here the factor of 3/4 is somewhat arbitrary, but leads to good numerical results, and the

parameter η is typically on the order of 1/(2M), where M measures the total mass of the

spacetime. The condition (4.83) is sometimes called the “nonshifting shift”, in contrast to

the “shifting shift” condition for which all time derivatives ∂

in equation (4.83) are replaced

with advective derivatives ∂

− β

∂

. This replacement has the effect of advecting various

unstable “gauge” modes off the computational grid. Various versions of these gauge condi-

tions have been adopted in simulations, and we refer the reader to the literature for compar-

ative studies and analysis.

The hyperbolic Gamma-driver condition (4.83), together with

the “1+log” slicing condition (4.51), have been extremely successful in dynamical “moving

puncture” black hole simulations, which we will discuss in detail in Chapter 13.1.3.

Box 4.1 Lapse and shift conditions: a sampler

Geodesic slicing,orGaussian normal coordinates,aredeﬁnedby

α = 1 β

= 0. (4.84)

Maximal slicing assumes K = 0 = ∂

K , which results in the lapse equation

α = α



+ 4π(ρ + S)



. (4.85)

Sometimes maximal slicing is imposed “approximately” via a K-driver,asinequation(4.39).

Harmonic coordinates satisfy

(4)



= 0, (4.86)

which results in conditions (4.44)and(4.45) for the lapse and shift. An important variation is

1+log slicing

(∂

− β

∂

)α =−2αK. (4.87)

For the minimal distortion gauge condition, the shift satisﬁes

(

β)

= 2A

α +

αγ

K + 16παS

. (4.88)

A related condition may be deﬁned in terms of the “conformal connection functions”



,which

leads to the hyperbolic Gamma-driver condition for the shift,

∂

,∂

= ∂



− ηB

, (4.89)

where η is a constant of order 1/(2M). Replacing the time derivatives ∂

in equation (4.89) with

advective derivatives ∂

− β

∂

is one useful variation.

We note in closing that in simulations of gravitational collapse to black holes, minimal

distortion and its close relatives often lead to poor coordinate resolution in the strong-ﬁeld

van Meter et al. (2006); see also Gundlach and Martin-Garcia (2006b) for a discussion of the mathematical properties

of the resulting evolution system.

122 Chapter 4 Choosing coordinates: the lapse and shift

central regions. The reason for this behavior can be understood in terms of a simpliﬁed

model equation

for condition (4.74),

= 16π S

, (4.90)

where we shall assume conformal ﬂatness ¯γ

= η

so that

is a ﬂat Laplace operator. If

matter is collapsing, the density is increasing, which, according to the continuity equation

(e.g., equation 5.12), leads to a negative divergence of the matter ﬂux S

.Takingthe

divergence of equation (4.90) we therefore ﬁnd

= 16π

< 0. (4.91)

Adopting boundary conditions to insure that β

vanishes assymptotically, we ﬁnd from

equation (4.91) that the divergence of the shift is positive,

> 0.

This means that the

shift moves coordinate observers away from the collapsing matter, thereby decreasing

the coordinate resolution of this matter. This response is unfortunate, since it is precisely

the strong ﬁeld that forms around the high-density, collapsing matter that leads to black hole

formation and hence this region is where we desire maximum resolution. To compensate

for this coordinate “blow-out” arising in minimal distortion (and related gauge choices), it

is possible to add an extra contribution to the shift vector that points inward towards the

collapsing region and increases the resolution there.

Shibata (1999b).

Viewed as Poisson’s equation

 =−4πρ

for the electromagnetic “potential”  ≡

, equation (4.91)hasa

positive “charge density” ρ

≡−4

on the right hand side, for which the potential must then be positive.

See Shibata (1999b); Duez et al. (2003).

Matter sources

In nonvacuum spacetimes, the stress-energy tensor T

contributes source terms in the

3 + 1 ﬁeld equations. The stress-energy tensor accounts for all sources of energy-

momentum in spacetime, excluding gravity. It thus arises from all forms of matter, electro-

magnetic ﬁelds, neutrinos, scalar ﬁelds, etc. in the Universe. For brevity, we shall sometimes

refer to these sources collectively as the “matter sources” and the terms that they contribute

in the 3 + 1 equations as the “matter source terms”. Matter source terms appear in the

Hamiltonian constraint equation (2.132), the momentum constraint equation (2.133), and

the evolution equation (2.135). The evolution equations for the “matter” sources are given

by ∇

= 0, which express the conservation of the total 4-momentum in spacetime.

Here T

is the total stress-energy tensor of the system. These conservation equations

must be solved simultaneously with the 3 +1 evolution equations for the gravitational

ﬁeld to determine the entire foliation of spacetime. Some of the quantities appearing in

the stress-energy tensor require auxiliary equations. These auxiliary equations include, for

example, the continuity equation and an equation of state in the case of hydrodynamic

matter, and Maxwell’s equations in the case of an electromagnetic ﬁeld, and so on.

AsshowninChapter2, the matter source terms ρ, S

and S

appearing in the 3 + 1

equations for the gravitational ﬁeld are projections of the stress-energy tensor into n

and

 and are given by

ρ = n

=−γ

= γ

(5.1)

The quantity ρ is the total mass-energy density as measured by a normal observer,

is the momentum density and S

is the stress. Finally, S is deﬁned as the trace

of S

S = γ

. (5.2)

In the following sections, we discuss some of the most important “matter” sources that

arise in astrophysical applications. These sources include hydrodynamic ﬂuids, magne-

tohydrodynamic plasmas threaded by magnetic ﬁelds, radiation gases (e.g., photon and

neutrino), collisionless matter, and scalar ﬁelds.

123

124 Chapter 5 Matter sources

5.1 Vacuum

At the risk of stating the obvious, vacuum spacetimes are characterized by the vacuum

stress-energy tensor

= 0. (5.3)

Spacetimes containing black holes and (or) gravitational waves, and nothing else, are char-

acterized by such a stress-energy tensor in Einstein’s ﬁeld equations. Vacuum spacetimes

are simpler to deal with numerically since they require no additional energy-momentum

conservation equations or auxiliary ﬁeld equations to solve. On the other hand, vacuum

black hole spacetimes, which are among the most important astrophysically, have all the

numerical complications in the gravitational ﬁeld sector that arise from the spacetime

singularities in the black hole interiors. Handling these complications will be a major issue

that we will address in subsequent chapters.

5.2 Hydrodynamics

Relativistic hydrodynamic matter is an important source of stress-energy in many astro-

physical applications. Loosely speaking, a hydrodynamic description of matter is appropri-

ate whenever the mean free path of a particle due to collisions with neighboring particles

is much shorter than the characteristic size or local scale length of the system. Many of the

global properties of relativistic stars, like neutron stars, are described by a hydrodynamic

ﬂuid. Also, the very early Universe is ﬁlled with hydrodynamic matter (e.g., baryons and

electrons) coupled to thermal radiation.

In this section we seek to cast the basic equations of relativistic hydrodynamics and mag-

netohydrodynamics into forms most suitable for building dynamical spacetimes numer-

ically, i.e., into forms best suited for numerical integration together with the Einstein

equations in 3 + 1 form for the gravitational ﬁeld. Although our discussion is reasonably

self-contained, readers of this section will ﬁnd it helpful to have been exposed previously to

a treatment of thermodynamics, hydrodynamics and electrodynamics in curved spacetime

at an elementary level.

5.2.1 Perfect gases

The stress-energy tensor of a perfect gas is given by

= ρ

+ Pg

. (5.4)

Here u

is the ﬂuid 4-velocity, ρ

is the rest-mass density, P is the pressure, and h is the

speciﬁc enthalpy

h = 1 +  + P/ρ

, (5.5)

See, e.g., Misner et al. (1973), Section 22.1–22.4.

5.2 Hydrodynamics 125

where  is the speciﬁc internal energy density. The total mass-energy density as measured

by an observer comoving with the ﬂuid is then given by ρ

∗

= ρ

(1 + ).

The equations of motion for the ﬂuid can be derived from the local conservation of

energy-momentum,

∇

= 0, (5.6)

and the conservation of rest-mass,

∇

(ρ

) = 0. (5.7)

Wilson scheme

The resulting equations can be cast in various forms, depending on the ﬂuid variables

chosen for numerical integration. One of the most straightforward schemes that readily

lends itself to numerical integration was introduced by Wilson.

His scheme was designed

to mimic the standard Eulerian equations of hydrodynamics in Newtonian theory, for which

the large body of expertise that has been acquired to solve the equations numerically can be

brought to tackle the relativistic system. To accomplish this goal, one deﬁnes a rest-mass

density variable

D ≡ ρ

W, (5.8)

an internal energy-density variable

E ≡ ρ

W, (5.9)

and a momentum-density variable

≡ ρ

hWu

= (D + E + PW)u

. (5.10)

Here we have deﬁned W to be the Lorentz factor between normal and ﬂuid observers,

W ≡−n

= αu

. (5.11)

Exercise 5.1 (a) Determine the relation between the density ρ

∗

and the source term

ρ deﬁned in equation (5.1). Interpret your answer physically.

(b) Argue that ρ

∗

= ρ whenever a normal observer is a comoving observer.

∗

= ρ for a static metric, like the metric (1.73) for an Oppenheimer–

Volkoff equilibrium star.

(d) Argue that ρ

∗

= ρ = ρ

for the metric (1.88) describing Oppenheimer–Snyder

spherical dust collapse.

Exercise 5.2 Show that the spatial vector S

deﬁned by equation (5.10)istheﬂuid

contribution to the source term S

appearing in the 3 + 1 decomposition of Einstein’s

ﬁeld equations and deﬁned by equation (5.1).

Wilson (1972b, 1979); see also Hawley et al. (1984).

126 Chapter 5 Matter sources

In terms of these variables, the equation of continuity ∇

(ρ

) = 0 becomes

∂

(γ

1/2

D) + ∂

(γ

1/2

) = 0. (5.12)

Here v

= u

is the ﬂuid 3-velocity with respect to a coordinate observer, and γ is the

determinant of the spatial metric γ

. Contracting equation (5.6) with u

yields the energy

equation,

∂

(γ

1/2

E) +∂

(γ

1/2

) =−P



∂

(γ

1/2

W ) + ∂

(γ

1/2

W v

)



, (5.13)

while the spatial components of (5.6) yield the relativistic Euler equations,

∂

(γ

1/2

) + ∂

(γ

1/2

) =−αγ

1/2



∂

P +

2αS

∂



. (5.14)

In deriving equations (5.12), (5.13) and (5.14), we have used the identity (−g)

1/2

= αγ

1/2

(see exercise 2.25).

Exercise 5.3 Derive equations (5.12), (5.13)and(5.14).

Exercise 5.4 (a) Show that by contracting equation (5.6) with u

the energy equation

may be written in the alternative form

dρ

∗

dτ

(ρ

∗

+ P)

dρ

dτ

, (5.15)

where d/dτ ≡ u

∇

is the rate of change with proper time following a ﬂuid element.

(b) Show that by contracting equation (5.6) with the projection tensor P

= g

the Euler equation may be written in the form

(ρ

∗

+ P)u

∇

=∇

P − u

∇

P. (5.16)

motion is adiabatic, i.e., ds/dτ = 0, where s is the entropy per unit rest-mass.

Though entirely equivalent to equations (5.13)and(5.14), equations (5.15)and(5.16)

usually are not the most useful form of the equations of relativistic hydrodynamics

for numerical integration in an initial-value (evolution) problem.

For many purposes it is useful to employ a simple “-law” equation of state (EOS) of

the form

P = ( − 1)ρ

. (5.17)

Realistic applications involving relativistic objects are rarely described by EOSs obeying

this simple form. However, a -law EOS provides a computationally practical, albeit crude,

approximation that can be adapted to mimic the gross behavior of different states of matter

in many applications. For example, to model a stiff nuclear EOS in a neutron star, one

can adopt a moderately high value of  in a -law EOS, e.g.,  ≈ 2. By contrast, to

model a moderately soft, thermal radiation-dominated EOS governing a very massive or

5.2 Hydrodynamics 127

supermassive star, one can set  = 4/3. For isentropic ﬂow, a -law EOS is equivalent to

the equation of state of a polytrope,

P = Kρ



,= 1 + 1/n (5.18)

where n is the polytropic index and K is the gas constant. However, for nonisentropic ﬂow,

which is always the case when encountering a shock, K is no longer constant throughout

the ﬂuid, and equations (5.17) and equations (5.18) are no longer equivalent.

Exercise 5.5 Use the ﬁrst law of thermodynamics for isentropic ﬂow to show the

equivalence of equations (5.17) and equations (5.18).

When employing a -law EOS the source term on the right-hand side of the energy

equation (5.13) can be rearranged to yield

∂

(γ

1/2

∗

) + ∂

(γ

1/2

∗

) = 0, (5.19)

where we have introduced a new energy variable E

∗

deﬁned as

∗

≡ (ρ

)

1/

W. (5.20)

This simpliﬁcation has great computational advantages, since the time derivatives on the

right-hand side of (5.13) are difﬁcult to handle in strongly relativistic ﬂuid ﬂow.

Exercise 5.6 Use the ﬁrst law of thermodynamics to show that the left-hand side

of equation (5.19) actually describes the evolution of s, the entropy per unit mass,

following a ﬂuid element:

∂

(γ

1/2

∗

) + ∂

(γ

1/2

∗

) =

αγ

1/2



∗



(1−)

dτ

, (5.21)

where d/dτ = u

∇

is the rate of change with proper time following the ﬂuid.

Hence equation (5.19) guarantees that the ﬂow is adiabatic, which is always true for

a perfect gas in the absence of shocks.

For given values of u

, W can be found from the normalization relation u

=−1,

W = αu



1 + γ



1/2

, (5.22)

and v

from

= αγ

/W − β

. (5.23)

Equations (5.12), (5.14) and (5.19), expressing the conservation of rest-mass, momen-

tum and energy, form a convenient set that can integrated simultaneously to determine the

evolution of a perfect gas in the absence of shocks. However, this particular set of equa-

tions must be suitably modiﬁed to handle the appearance of shock discontinuities. Shocks

occur whenever ﬂuid elements collide supersonically, which can occur during stellar col-

lapse, accretion ﬂows, stellar collisions, etc. Mathematically, shocks arise whenever ﬂuid

128 Chapter 5 Matter sources

characteristic curves cross, leading to sharp discontinuities in the ﬂuid variables. For realis-

tic ﬂuids containing viscosity, the shock transition region has a ﬁnite thickness amounting

to several particle collision mean-free-paths, and the ﬂuid variables vary continuously

across this ﬁnite interval. Two strategies are commonly adopted to treat shocks when the

details of the transition region are of no physical consequence and the perfect gas equa-

tions are used throughout, whereby the thickness of the transition region is allowed to be

arbitrarily thin. The more traditional approach is to add an artiﬁcial viscosity term to the

equations.

This term mimics the effect of physical viscosity, except that it is employed

only in the vicinity of a shock and serves to spread the shock transition region over a few

spatial grid spacings in ﬁnite-difference codes.

For this purpose, the artiﬁcial viscosity

term P

vis

is nonzero only where the ﬂuid is compressed and is added to the pressure on the

right-hand sides of both the energy equation (5.13) and the Euler equation (5.14). Such a

term has the approximate form

vis



vis

(δv)

for δv < 0,

0otherwise,

(5.24)

where δv = ∂

x, x is the local spatial grid size and C

vis

is a dimensionless constant of

order unity. Adding such a term allows the ﬂuid to satisfy the Rankine–Hugoniot “jump”,

or junction, conditions across the shock, which we shall derive below. These conditions

ensure the continuity of rest-mass, momentum and energy ﬂux across any surface in the

ﬂuid, including a shock front. For a ﬂuid element traversing a shock, the junction conditions

serve to convert some of the bulk kinetic energy into internal energy and to increase the

entropy. Artiﬁcial viscosity schemes have the virtue of being quite robust and very easy

to implement. For shocks occurring in Newtonian ﬂuids with modest Mach numbers,

artiﬁcial viscosity enables the ﬂuid to satisfy the Rankine–Hugoniot jump conditions to

reasonable accuracy. Artiﬁcial viscosity has also been used successfully in many relativistic

applications,

but it can lead to less satisfactory results for ultrarelativistic ﬂows or high

Mach numbers.

Exercise 5.7 Show that the addition of P

vis

to the pressure in the stress-energy

tensor modiﬁes equation (5.19) according to

∂

(γ

1/2

∗

) + ∂

(γ

1/2

∗

) =−



∗



(1−)

vis



∂

(W γ

1/2

), (5.25)

where v

= u

. Comparing with equation (5.21), we see that the role of P

vis

to generate the entropy jump required across a shock discontinuity.

von Neumann and Richtmeyer (1950).

Finite-difference techniques for integrating partial differential equations are discussed in Chapter 6.2.

May and White (1966); Wilson (1972b); Shapiro and Teukolsky (1980); Hawley et al. (1984); Shibata (1999a); Duez

et al. (2003).

Winkler and Norman (1986).

5.2 Hydrodynamics 129

High-resolution shock-capturing (HRSC) schemes

An alternative approach to handling shocks in ﬁnite-difference algorithms involves recast-

ing the equations in a “ﬂux-conservative form” and adopting a so-called “high-resolution

shock-capturing scheme” or HRSC scheme.

In such a scheme, one divides up the spatial

domain into a discrete set of contiguous cells. At the center of each cell is a grid point where

one keeps track of the ﬂuid variables as they evolve in time. One treats all ﬂuid variables

as constant in each grid cell. The discontinuous ﬂuid variables at the grid interfaces serve

as initial conditions for a local Riemann shock tube problem, which can be treated either

exactly or approximately. The general Riemann problem involves the evolution of a gas in

which, initially, the ﬂuid variables are everywhere constant on either side of an interface,

but discontinuous across the interface. The solution to this idealized problem is known

and provides the basis for constructing an HRSC scheme. Allowing for discontinuities,

including shocks, lies at the core of such schemes, and does not require any additional

artiﬁcial viscosity. Constructing Riemann solvers for HRSC schemes requires knowledge

of the local characteristic structure of the equations to be solved. This has motived the

development of several ﬂux-conservative hydrodynamics schemes for which this charac-

teristic structure can be determined.

A key feature of the hydrodynamical equations in

ﬂux-conservative form is that they do not contain any derivatives of the ﬂuid variables in

the source terms on the right-hand sides, in contrast to equation (5.14), which contains a

pressure gradient.

The ﬂux-conservative equations of relativistic hydrodynamics take on the general form

∂

U + ∂

= S (5.26)

where

U is the state vector of conserved variables built out of the so-called primitive ﬂuid

variables

P = (ρ

, P), the F

are the ﬂux vectors (one for each spatial dimension,

i), and where the source vector

S does not contain any derivatives of the primitive ﬂuid

variables. A particularly useful choice which can be written in this way is

U =







˜τ













1/2

Wρ

1/2

αT

1/2

−







, (5.27)

where

D and

differ from D and S

deﬁned in equations (5.8) and (5.10) by a factor of

1/2

. The ﬂux vectors F

are then given by







αγ

1/2

−







, (5.28)

See the text by Toro (1999), or the review by Mart

ıandM

uller (1999), and references therein.

See, e.g., Courant and Friedrichs (1948).

See, e.g., review by Font (2000).

Font et al. (2000, 2002).

130 Chapter 5 Matter sources

and the source vector S by

S =







αγ

1/2

ab, j

αγ

1/2

∂

α −

(4)



α)







. (5.29)

The ﬁrst row of this equation is just the continuity equation (5.12) expressing rest-mass

conservation. The second row arises from the equation of energy-momentum conservation

∇

= 0, or

∂

(

√

−gT

) + ∂

(

√

−gT

) =

√

−gT

∂

. (5.30)

The second row is simply the a = j ,or j-momentum, component of equation (5.30).

The third row arises from projecting the equations of motion along the normal vector:

∇

= 0, and subtracting off the continuity equation (5.12) to achieve the desired

form. The subtraction removes the conserved rest-mass contribution from the energy

density so that the internal energy density, which equals the (typically small) difference

between the total and rest-mass energy densities, can be integrated directly without being

corrupted by numerical cancelation errors.

Exercise 5.8 On ﬁrst glance it appears that the source term in the energy equation

(i.e., the third row in

S), which we will call s

˜τ

, contains explicit time derivatives of

the lapse and shift, as well as the spatial metric.

(a) Show that by expanding out the source term s

˜τ

and using the deﬁnition of the

Christoffel symbol, the time derivatives of the lapse and shift miraculously cancel

out.

(b) Now let us recast the source term to get a convenient expression free of time

derivatives. First use n

∇

= 0 to show that

˜τ

=−αγ

1/2

∇

= αγ

1/2

+ a

), (5.31)

where K

is the extrinsic curvature and a

= n

∇

= D

(ln α) (see exercise 2.13)

is the 4-acceleration of the normal observer.

α = 0thatD

α = β

∂

α.

(d) Finally, show that K

= β

and K

= β

. Conclude that the source

term can be written as

˜τ

= αγ

1/2

[(T

+ 2T

+ T

− (T

+ T

)∂

α], (5.32)

which contains no explicit time derivatives.

The basic strategy is to integrate the coupled set of equations (5.26) for the conservative

variables (

, ˜τ) from time t to t + t, and then combine the conservative variables at

the new time to solve algebraically for the new set of primitive variables (ρ

, P). The

procedure is repeated for each subsequent time step. For a self-consistent determination of

the background gravitational ﬁeld, Einstein’s equations in 3 +1 form must be integrated

in time simultaneously to determine the spacetime metric.