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

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

5.3 Collisionless matter 171

matter source terms given in equations (5.226)–(5.229) below, cannot be determined from

F until the metric is known. The metric, in turn, is determined by solving Einstein’s ﬁeld

equations, which requires the stress-energy tensor. The task of constructing an equilibrium

conﬁguration therefore amounts to solving a set of coupled integro-differential equations

self-consistently. The equations are quite straightforward to solve in the case of spherical

symmetry, especially when the distribution of stellar velocities at each point in the cluster

is isotropic. In this case the equations can be made to look like the familiar OV equations

for ﬂuid stars in spherical equilibrium.

Exercise 5.42 Show that the velocity proﬁle of particles in a spherical cluster

in static equilibrium is isotropic as measured by a local observer whenever the

distribution function is independent of angular momentum, i.e., F = F(E, m).

The problem is conceptually similar, but quite a bit more complicated, in the case

of stationary axisymmetric equilibria. Nonspherical clusters with and without rotation

are included in this category. The simplest distribution functions that can give rise to

nonspherical, axisymmetric equilibria are of the form F = F(E, J

, m), where E is the

particle energy and J

is the particle angular momentum along the symmetry axis. By

suitable choice of the distribution in J

, one can construct models that are either prolate

or oblate. A depletion in the J

-distribution produces a prolate conﬁguration, while an

enhancement produces an oblate one.

We shall return to the construction of collisionless

equilibria and their utility as initial data in Chapter 8.2.

Now we wish to discuss how the spacetime for a relativistic cluster of collisionless

matter can be determined when the cluster is undergoing dynamical evolution, including

catastrophic collapse to a black hole. The basic goal is to solve the Liouville equation for

the matter simultaneously with the 3 +1 equations for the gravitational ﬁeld. There are

at least two different numerical strategies that have been employed successfully to track

the dynamical evolution of relativistic clusters. We describe both of them in the sections

immediately below. Applications of these techniques to problems involving collisionless

matter in strong gravitational ﬁelds are treated in Chapters 8 and 10.

Particle methods

Particle methods, sometimes referred to as “particle-mesh” or “particle-in-cell” methods,

are widely used in plasma physics and in Newtonian gravitation for the study of collisionless

systems.

The method has been extended to general relativity, yielding a robust mean-

ﬁeld, N -body evolution scheme in the context of the 3 + 1 equations that has been applied

to spherical and axisymmetric clusters.

In this approach, a statistical representation of

See Misner et al. (1973), Section 25.7.

See Shapiro and Teukolsky (1993a,b).

See, e.g., Hockney and Eastwood (1981)andSellwood (1987) for a discussion and review of applications to Newtonian

gravitation and plasma physics.

Shapiro and Teukolsky (1985a,b,c, 1986, 1992b) and references therein.

172 Chapter 5 Matter sources

the initial distribution function is constructed by specifying initial positions and velocities

for a large but ﬁnite number N of discrete particles. The motion of these particles is then

calculated by integrating simultaneously N geodesic equations in the mean gravitational

ﬁeld of the system. The source terms of the 3 + 1 equations are determined by smearing out

each particle over a small spatial volume, and adding up the contributions of all particles.

Working with collisionless matter via particle simulations has several advantages over

ﬂuid systems. The collisionless matter equations are ordinary differential equations (the

geodesic equations), while the hydrodynamic equations are more challenging partial dif-

ferential equations. Furthermore, collisionless matter is not subject to shocks or other

hydrodynamical discontinuities. These complications require special care or sophisticated

handling of the ﬂuid equations for adequate resolution, as we discussed in Section 5.2 of

this chapter.

Consider a collisionless gas sampled by a discrete number of particles. Group together

into distinct categories those particles at a point in spacetime that have the same mass and

4-velocity. Then the stress-energy tensor is



. (5.221)

Here m

is the rest-mass of category A particles, u

is their 4-velocity, and n

is their

comoving number density.

The quantity n

is proportional to a δ-function in the limit of

point particles, but it can be treated as a continuous, averaged quantity when the number

of particles is inﬁnite, the mass of each particle becomes inﬁnitesimal, and the particle

distribution is smooth.

The equation of motion of the particles is simply the geodesic equation

∇

= 0, (5.222)

where we drop the label A when referring to any one of the particles.

Exercise 5.43 Show that equation (5.222) follows from the conservation of stress-

energy, ∇

= 0, and the conservation of particles, ∇

(

)

= 0.

Exercise 5.44 Show that the equations describing the motion of a particle in phase

space can be written in terms of the 3 + 1metricas

=−αu

∂

α + u

∂

−

∂

, (5.223)

= γ

− β

, (5.224)

αu



1 + γ



1/2

. (5.225)

Exercise 5.45 Derive the Newtonian limit of the equations of motion in exer-

cise 5.44.

See, e.g., Misner et al. (1973) Section 5.4.

5.3 Collisionless matter 173

Hint: In the Newtonian limit, γ

→ η

, α → 1 + ,where is the gravitational

potential, and β

→ 0.

Thus, given the metric at any time t, equations (5.223)–(5.225) may be integrated for the

new positions and 4-velocities of the particles at t + t. Once the new particle positions

and velocities are determined, the new source terms needed in the ﬁeld equations can be

computed as a discrete sum over the particles using equations (5.207)–(5.210):

ρ =



, (5.226)



(5.227)



(5.228)

S = ρ −



. (5.229)

We regard each particle as being its own category A, so that the sum over A is a sum over

particles. Then



W γ

1/2

x

. (5.230)

Here 

is the invariant 3-volume as measured in a frame comoving with particle A (see

equation 5.201). In writing the last equality in equation (5.230) we divided up coordinate

space into bins of coordinate volume x

x

and we treated the particle as smeared

out over the entire bin in which it resides. Note that more sophisticated binning and

averaging algorithms can be constructed which assign fractions of particles to neighboring

bins in order to calculate smoother and statistically more reliable source terms from the

discrete particle sample.

Once the source terms have been assembled, the 3 + 1 equations can be integrated

to give the mean gravitational ﬁeld quantities at the new time. The whole cycle can be

repeated to evolve the particles to the next time step in the mean ﬁeld, and so on.

Phase space methods

Sampling the phase space distribution function f by a discrete set of particles does not

determine f , but only moments of f , like the source terms of the ﬁeld equations.

Some-

times it is important to acquire the full knowledge of f to understand the underlying

dynamical state of a collisionless system. Many complicated phenomena of great astro-

physical interest, such as “violent relaxation”, can only be understood by a study of the

detailed evolution of f . Violent relaxation, which can lead to the virialization of a cluster

The discussion in this section is patterned after Rasio et al. (1989b).

174 Chapter 5 Matter sources

that is initially far from dynamical equilibrium, is the result of collisionless damping of

perturbations by “phase mixing”.

Studying violent relaxation requires that the phase

space density be tracked in detail, including all distortions due to phase mixing.

A phase space method does not rely on any statistical representation of the system by

particles. Instead, a phase space method explicitly constructs the smooth distribution func-

tion of matter in phase space. The source terms of the ﬁeld equations are obtained by direct

numerical quadratures of f over the momentum (or velocity) space, as in equation (5.212).

This has the great advantage of eliminating random statistical ﬂuctuations in the data

while at the same time providing the full distribution function of the system. However,

very few phase space methods have so far been developed successfully, even in the much

simpler framework of Newtonian gravity. Part of the problem is the extreme complexity

of working in phase space instead of physical space. The large number of dimensions

in phase space (already three in spherical symmetry where physical space has only one)

would already discourage many attempts. In addition, distribution functions often have

irregular structures that can be hard to represent accurately on a numerical grid in phase

space. Such irregular structures can arise from discontinuities in the initial data, but even

in the case of very smooth initial data, phase mixing will usually produce increasingly

intricate “ﬁne-grained” structures.

One phase space method developed for general relativistic systems seems to work well,

as least for the spherical systems to which it has been applied.

The method exploits

Liouville’s theorem to determine the evolution of the phase space distribution function

directly. The method consists of three basic steps to propagate the distribution function

f from time t

to t

> t

: (1) compute the source terms for the ﬁeld equations (e.g.,

equation 5.212) by integrating f at t

over momentum space; (2) evolve the 3 + 1ﬁeld

equations to determine the metric at t

; (3) compute f at time t

, assuming that it satisﬁes

the Liouville equation (5.216). The key idea is to use Liouville’s theorem for step (3),

which can be written as

f (t

, x

, p

) = f (t

, x

, p

), (5.231)

where (x

, p

) represents the position in phase space at time t

of a test particle that will

actually reach the position (x

, p

) at time t

. By interpolating the metric in the interval

, t

) one can integrate the equations of motion, equations (5.223)–(5.225), to construct

the actual trajectory of such a test particle. Since f is known at all times t ≤ t

one can

actually determine f (t

, x

, p

)forall(x

, p

The original scheme along these lines was developed in Newtonian theory.

However,

it was soon discovered that in many cases the scheme developed numerical instabilities,

rendering the results inaccurate.

The reason was that in the original scheme f was

Lynden-Bell (1967).

Rasio et al. (1989b).

Fujiwara (1981, 1983).

Inagaki et al. (1984); Nishida (1986).

5.4 Scalar ﬁelds 175

constructed on a grid of points in phase space. In employing equation (5.231) to determine

f at t

by placing a particle at (x

, p

) and integrating backward in time from t

to t

,the

old position of the particle (x

, p

) is not a grid point at t

in general. As a result, one must

perform a (multidimensional) interpolation on the grid to get f (t

, x

, p

) and thereby the

new grid point value f (t

, x

, p

). The error in the interpolation on the grid at t

propagates

to the new grid at t

, where it becomes an error on the values of f at the grid points. This

leads to an ampliﬁcation of the error, which is reﬂected in the failure to conserve rest mass

and other conserved quantities, and may even lead to the violation of the positivity of the

distribution function. The problem is most severe when the distribution function exhibits

discontinuities, or even a mild degree of phase mixing.

One cure adopted in the relativistic scheme that eliminates interpolation error is to

extend the value of t

in equation (5.231)tot = 0, where the distribution function is given

by the initial data and is therefore known everywhere to arbitrary accuracy. When f is

required at some point in phase space (e.g., in performing the quadratures for the source

terms), this point is simply tracked along a dynamical path all the way back to t = 0, where

f can be accurately evaluated from the initial data. For problems involving discontinuous

distribution functions or large degrees of phase mixing, more points can be added in phase

space to do the quadratures to the required accuracy. There is never any need to introduce

any grid at all in phase space, since intermediate values of f are never needed and therefore

need not be stored. The only disadvantage, of course, is that the computational time per

time step increases dramatically with time, since longer and longer trajectories have to be

constructed to evaluate f at a given point.

5.4 Scalar ﬁelds

A classical real scalar ﬁeld provides one of the simplest sources of stress-energy in

Einstein’s equations; a classical complex scalar ﬁeld furnishes another example, only

slightly more complicated. Upon quantization, the momentum eigenstates of a scalar ﬁeld

are observable as particles. Scalar ﬁelds give rise to particles of spin 0, while vector ﬁelds

(like the electromagnetic ﬁeld) give rise to particles of spin 1 (like the photon, in the

case of electromagnetism), and tensor ﬁelds of rank two or higher give rise to higher-spin

particles. A complex scalar ﬁeld has two degrees of freedom instead of just one, and it can

be interpreted as a particle and an antiparticle. Real ﬁelds are their own antiparticles. A

neutral π-meson is an example of a real scalar ﬁeld, while the charged π

-andπ

−

-mesons

are described by complex scalar ﬁelds.

Interest in scalar ﬁelds has been stimulated in recent years by theoretical developments in

particle physics and cosmology. Both the standard model of elementary particles as well as

their superstring extensions involve scalar ﬁelds, although the existence of a fundamental

elementary scalar particle has yet to be conﬁrmed by an accelerator experiment. For

example, the Higgs boson, which serves to generate masses for the W

and Z

gauge

vector bosons in the electroweak theory of Weinberg, Salam and Glashow, is a particle

176 Chapter 5 Matter sources

of ﬁnite mass described by a neutral spin 0 scalar ﬁeld. Its discovery is a top priority of

experimental particle physics. There are many attempts at unifying the standard model

with gravitation at the quantum level, like string theory. Typically, these theories give rise

to 4-dimensional “effective” models in which the usual spacetime metric g

of gravity is

accompanied by one or more scalar ﬁelds. The massive scalar dilaton and the (pseudo-)

scalar axion are two examples.

In cosmology, the inﬂationary phase in the early Universe can be driven by the vacuum

energy density provided by the potential of a time-dependent, but slowly varying scalar

ﬁeld ϕ(t) called an “inﬂaton”. It is also possbile that the existence of “dark energy” in

the Universe, which manifests itself as a nonzero cosmological constant in the standard

Big Bang model, may be represented by the vacuum energy density associated with the

potential of yet another scalar ﬁeld at a much lower characteristic energy (“quintessence”).

Candidates for “dark matter” in the Universe include bosonic, as well as fermionic

particles. It is an issue of wide speculation whether or not “boson stars” could actually arise

during the gravitational condensation of bosonic dark matter in the early Universe. Boson

stars are self-gravitating, stationary equilibrium conﬁgurations constructed from complex

scalar ﬁelds in asymptotically ﬂat spacetimes.

They are macroscopic quantum states

that are supported against gravitational collapse by the Heisenberg uncertainty principle;

they can be modeled by classical scalar ﬁelds. Like neutron stars, boson stars can have

highly relativistic gravitational ﬁelds and yet are nonsingular and have no event horizons.

However, if the scalar ﬁelds have self-interactions, then, unlike neutron stars, boson stars

can be very massive. It is therefore not surprising that boson stars are sometimes invoked

as alternatives to black holes to model massive, compact stars. By contrast with boson

stars, “soliton stars”, which are constructed from real scalar ﬁelds, are not stationary but

periodic, both in the spacetime geometry and the matter ﬁeld. Nonsingular, self-gravitating

stationary solutions do not exist for real, massive scalar ﬁelds.

Apart from their possible physical signiﬁcance, scalar ﬁelds serve as very useful tools

for probing strong gravitational ﬁeld phenomena and for learning how to do numeri-

cal relativity. The dynamical equation governing a scalar ﬁeld is the simple, classical

Klein–Gordon equation (see discussion below). In contrast to the partial differential equa-

tions describing hydrodynamic or magnetohydrodynamic matter, which can exhibit shock

waves and other discontinuities that require special handling (see Section 5.2), the Klein–

Gordon equation does not tend to develop discontinuities from smooth initial data and

is thus straightforward to integrate. As a result, integrating the scalar wave equation is

often chosen as the ﬁrst test of a new numerical evolution scheme, a choice that can

prove instructive even when the equation is integrated in ﬂat spacetime. Another conse-

quence of the simple, well-behaved nature of the scalar wave equation is that it is fairly

easy to implement advanced numerical techniques such as “adaptive mesh reﬁnement”

See, e.g., Kaup (1968); Rufﬁni and Bonazzola (1969); Colpi et al. (1986); Lee and Pang (1992); Seidel and Suen

(1990, 1991); Yu an et al. (2004); Schunck and Mielke (2003).

5.4 Scalar ﬁelds 177

(AMR)

for high-resolution integrations of the equation, particularly in 1 + 1 dimen-

sional simulation (e.g., spherical symmetry). The detailed study of gravitational collapse

and black hole formation using AMR with a scalar wave source has lead to the discovery

of black hole “critical pheonomena”.

Here the behavior observed in dynamical simula-

tions of conﬁgurations at the onset of black formation along a one-parameter family of

initial data assumes many of the features of a phase transition in a statistical mechanical

system.

Finally, scalar ﬁelds generate radiation ﬁelds – waves which travel at the speed of

light locally, have amplitudes that fall-off with radius like r

−1

at large distances from a

central source, and carry away mass-energy – even in spherical symmetry! By contrast,

gravitational waves, which are an important yield of numerical simulations in general

relativity, cannot arise in spherical symmetry, as we know from Birkhoff’s theorem. So

working with scalar ﬁelds in spherical symmetry allows one to gain experience tackling

some of the computational challenges of wave generation and propagation in dynamical

spacetimes without having to integrate in more than one spatial dimension.

The stress-energy tensor of a real scalar ﬁeld, ϕ(x

), is given by

=∇

ϕ∇

ϕ −

∇

ϕ − g

V (ϕ), (5.232)

where V (ϕ) is the potential. The potential may be decomposed according to

V (ϕ) =

+ V

int

, (5.233)

where m is the mass of the ﬁeld (i.e., the mass of the momentum eigenstates, or “particles”,

when the ﬁeld is quantized), and V

int

is an interaction potential.

When m = 0theﬁeld

is massless; when V

int

= 0 the ﬁeld is noninteracting, apart from gravitation (“minimal

coupling”). The equations of motion (i.e., the scalar ﬁeld equations) follow from ∇

0, which gives the Klein–Gordon equation with a potential term,

∇

ϕ − m

ϕ −

int

dϕ

= 0. (5.234)

Exercise 5.46 Calculate the total energy density for a real scalar ﬁeld measured in

the rest frame of an observer in Minkowski spacetime and give a physical interpre-

tation for each of the contributions in your expression.

Consider the case of a massless ﬁeld in the absence of an interaction potential (i.e., V =

0), for which equation (5.234) reduces to the massless Klein–Gordon equation in curved

spacetime,

∇

ϕ = 0. (5.235)

Berger and Oliger (1984); see Chapter 6.2.5.

Choptuik (1993). See Chapter 8 for a discussion.

Here we set

h = 1 in addition to c = 1 = G, so that terms like (mc/

become m

178 Chapter 5 Matter sources

In the special case of a stationary black hole background, the only nonsingular, time-

independent solution to equation (5.235)isϕ = 0. If the initial value of ϕ varies in space

or time, it will evolve to this solution eventually. This result was demonstrated by Price

and is consistent with black-hole uniqueness (“no-hair”) theorems.

Exercise 5.47 Consider the propagation of a massless, noninteracting scalar ﬁeld

in a static spherical spacetime,

=−e

2(r)

+ e

2(r)

(dθ

+ sin

θdφ

). (5.236)

Decompose the ﬁeld into spherical harmonics

ϕ(t, r,θ,φ) =

(θ,φ). (5.237)

Show that the scalar ﬁeld equation reduces to

−∂

∂

ψ + ∂

∗

∂

∗

ψ = V

(r)ψ, (5.238)

where r

∗

is a “generalized tortoise coordinate” satisfying

∗

= e

−

, (5.239)

and where V

(r) is an effective potential given by

(r) = e

2



l(l + 1)

−2

∂

(

 − 

)



. (5.240)

Note: In the case of ﬂat spacetime, the solutions to the above radial wave equation

are known analytically for all l.

For a static vacuum Schwarzschild spacetime,

the numerical solutions exhibit three separate phases: an initial burst, a quasinor-

mal ringing phase, and a power-law tail phase. During the late-time tail phase, all

multipoles of the scalar wave decay at ﬁnite r

∗

according to ψ ∼ t

−(2l+3)

as t →∞

(“Price’s theorem” for scalar waves). The same radiative behavior is also observed

for scalar wave propagation in the static spacetime of a spherical equilibrium (OV)

star.

In exercise (5.48) we construct a conserved energy integral that provides a useful check

on numerical integrations of the scalar wave equation in static spacetimes, as derived in

exercise (5.47).

Exercise 5.48 In a static spacetime ξ

= ∂/∂t is a Killing vector and J

= T

is a conserved current, i.e., ∇

= 0.

(a) Prove the conservation law

(

δ



= 0, (5.241)

where δ is the closed 3-surface enclosing a 4-volume . The surface element d



is given by d



= (1/3!)

abcd

,where

abcd

is the Levi-Civita tensor.

Price (1972a,b).

See, e.g., Wald (1984).

Burke (1971).

See, e.g., Kokkotas and Schmidt (1999); Pavlidou et al. (2000), and references therein.

5.4 Scalar ﬁelds 179

(Note: here δ is a 3-surface that lies inside the worldtube of the matter depicted in

Figure 3.5.)

(b) Adopt spherical polar coordinates and choose  to be the volume conﬁned by

the spherical radii r

and r

and the times t

and t

. Show that the conservation law

(5.241) can be written as

E(t

) − E(t

) = J (r

) − J (r

), (5.242)

where

E(t) ≡



r=r



θ=0



2π

φ=0

√

−gdrdθ dφ (5.243)

is the energy of the matter ﬁeld contained between r

and r

at time t,and

J (r) ≡



t=t



θ=0



2π

φ=0

√

−gdtdθdφ (5.244)

is the radial ﬂux across r, integrated over time between t

and t

. Discuss the meaning

of equation (5.242).

An interesting mathematical connection exists between the massless Klein–Gordon

equation (5.235) and the evolution equation for a perfect, irrotational ﬂuid. As shown in

exercise (5.49), the hydrodynamic evolution of a relativistic ﬂuid obeying the special EOS

P = ρ

∗

can actually be determined by integrating the wave equation for a real, massless,

noninteracing scalar ﬁeld!

Exercise 5.49 The relativistic vorticity tensor is deﬁned as

= P

[

∇

(hu

) −∇

(hu

)

]

, (5.245)

where u

is the 4-velocity, h = (ρ

∗

+ P)/ρ

is the speciﬁc enthalpy, ρ

∗

= ρ

(1 + )

is the total mass-energy density, ρ

is the rest-mass density, and P

= g

+ u

is the projection tensor.

(a) Show that for a perfect ﬂuid, equation (5.62) may be used to recast equa-

tion (5.245)as

[

∇

(hu

) −∇

(hu

)

]

. (5.246)

(b) Argue that if the vorticity is zero, the quantity hu

can be expressed as the

gradient of a potential,

=∇

ϕ. (5.247)

What kind of ﬂow does this correspond to physically?

∇

[

(ρ

/ h)∇

]

= 0. (5.248)

(d) The equation of state relates ρ

to h,andh is found from the normalization

condition h = (−∇

ϕ∇

ϕ)

1/2

. Relate ρ

to h for a polytropic gas where P is given

by equation (5.18).

180 Chapter 5 Matter sources

(e) Consider the extreme equation of state P = ρ

∗

. Find  and show that for this

equation of state, equation (5.248) simplies to yield

∇

(

∇

)

= 0. (5.249)

What is the sound speed in this gas?

To integrate the Klein–Gordon equation (5.235) numerically it is useful to cast it into

ﬁrst-order form. For this purpose, introduce the following new variables:

 ≡

−1



∂

ϕ − β

∂



, (5.250)

≡ ∂

ϕ. (5.251)

For a general 3 + 1 metric, the Klein–Gordon equation (5.235) becomes

∂

ϕ = β

∂

ϕ − α, (5.252)

∂

 = β

∂

 − αg

∂

+ αg



− g

∂

α + αK , (5.253)

∂

= β

∂

+ ψ

∂

− α∂

 − ∂

α. (5.254)

It is noteworthy that the deﬁnition of ψ

, equation (5.251),turnsintoasetofthree

constraints,

≡ ∂

ϕ − ψ

= 0, (5.255)

that must be satisﬁed at all times. When solved exactly, the system of equations (5.252)–

(5.254) automatically preserves the constraints C

= 0, provided they are satisﬁed initially.

However, in a numerical simulation, truncation errors and boundary errors can cause C

wander away from zero. Hence it is useful to monitor the evolution of C

to help calibrate

the accuracy of the simulation. In the same way, monitoring how well the 3 + 1 constraints

are maintained serves to calibrate the accuracy of the gravitational ﬁeld integrations.

The Klein–Gordon system of ﬁrst-order, symmetric hyperbolic equations thus provides a

simple arena for exploring the growth of constraint violations in numerical simulations of

hyperbolic systems and devising methods for controlling them.

Evaluating the scalar ﬁeld equation (5.234) for a Robertson–Walker (RW) metric (see

exercise 4.1), assuming the ﬁeld is everywhere homogeneous (∇

ϕ = 0) gives a simple

second-order ordinary differential equation for the evolution of ϕ(t),

¨ϕ + 3H ˙ϕ + m

ϕ +

int

dϕ

= 0. (5.256)

In equation (5.256) the dot denotes time differentiation, H (t) =

a/a is Hubble’s constant,

and a(t) is the expansion parameter appearing in the RW metric. This equation is the

usual starting point of most discussions of inﬂation in the early Universe, in which case

Scheel et al. (2004).