6.1 The Einstein equations

6.1.1 ADM formalism

There are many ways to write the Einstein equations. The most common is the ADM or 3 + 1 form [23] which decomposes spacetime into layers of three-dimensional space-like hypersurfaces, threaded by a time-like normal congruence μ i n = (1,− β )∕α, where we use Greek (Latin) indices to specify spacetime (spatial) quantities. The general spacetime metric is written as

ds2 = (− α2 + βiβi) dt2 + 2βi dxidt + γij dxidxj, (8 )
where α and βi are the lapse function and shift vector respectively, and γij is the spatial 3-metric. The lapse defines the proper time between consecutive layers of spatial hypersurfaces, the shift propagates the coordinate system from 3-surface to 3-surface, and the induced 3-metric is related to the 4-metric via γμν = gμν + nμnν. The Einstein equations are written as four constraint equations,
(3)R + K2 − KijKij = 16 πG ρH, (9 ) ( ij ij ) j ∇i K − γ K = 8πGs , (10 )
twelve evolution equations,
∂tγij = ℒ βγij − 2αKij, [ ( )] (11 ) ∂tKij = ℒ βKij − ∇i ∇j α + α (3)Rij − 2KikKk + KKij − 8πG sij − 1sγij + 1ρHγij , j 2 2
and four kinematical or coordinate conditions for the lapse function and shift vector that can be specified arbitrarily (see Section 6.1.2). Here Kij is the extrinsic curvature describing how the 3-metric evolves along a time-like normal vector. It is formally defined as the Lie derivative (ℒn) of the spatial metric along the vector nμ, Kij ≡ − ℒn γij∕2. Also,
ℒ γ = ∇ β + ∇ β, β ij i j j i (12 ) ℒ βKij = βk ∇kKij + Kik∇j βk + Kkj ∇iβk,
where ∇i is the spatial covariant derivative with respect to γij, (3) Rij is the spatial Ricci tensor, K is the trace of the extrinsic curvature Kij, and G is the gravitational constant. The matter source terms ρH, sj, sij and s = sii as seen by the observers at rest in the time slices are obtained from the appropriate projections
ρ = nμn νT , (13 ) H μ μν si = − γi nνT μν, (14 ) s = γμγ νT (15 ) ij i j μν
for the energy density, momentum density and spatial stresses, respectively. Here c = 1, and Greek (Latin) indices refer to 4(3)-dimensional quantities.

It is worth noting that several alternative formulations of Einstein’s equations have been suggested, including hyperbolic systems [136] which have nice mathematical properties, and conformal traceless systems [14731Jump To The Next Citation Point] which make use of a conformal decomposition of the 3-metric and trace-free part of the extrinsic curvature Aij = Kij − γij K ∕3. Introducing &tidle;γij = e−4ψγij with e4ψ = γ1∕3 so that the determinant of &tidle;γij is unity, and &tidle;Aij = e−4ψAij, evolution equations can be written in the conformal traceless system for ψ, &tidle;γij, K, &tidle; Aij and the conformal connection functions, though not all of these variables are independent. However, it is not yet entirely clear which of these methods is best suited for generic problems. For example, hyperbolic forms are easier to characterize mathematically than ADM and may potentially be more stable, but can suffer from greater inaccuracies by introducing additional equations or higher order derivatives. Conformal treatments are considered to be generally more stable [31Jump To The Next Citation Point], but can be less accurate than traditional ADM for short term evolutions [6Jump To The Next Citation Point].

Many numerical methods have been used to solve the Einstein equations, including variants of the leapfrog scheme, the method of McCormack, the two-step Lax–Wendroff method, and the iterative Crank–Nicholson scheme, among others. For a discussion and comparison of the different methods, the reader is referred to [43], where a systematic study was carried out on spherically symmetric black hole spacetimes using traditional ADM, and to [31613] (and references therein) which discuss the stability and accuracy of hyperbolic and conformal treatments.

6.1.2 Kinematic conditions

For cosmological simulations, one typically takes the shift vector to be zero, hence ℒ βγij = ℒ βKij = 0. However, the shift can be used advantageously in deriving conditions to maintain the 3-metric in a particular form, and to simplify the resulting differential equations [6061Jump To The Next Citation Point]. See also [146] describing an approximate minimum distortion gauge condition used to help stabilize simulations of general relativistic binary clusters and neutron stars.

Several options can be implemented for the lapse function, including geodesic (α = 1), algebraic, and mean curvature slicing. The algebraic condition takes the form

μ α = F1(x )F2 (γ ), (16 )
where F (xμ) 1 is an arbitrary function of coordinates xμ, and F (γ) 2 is a dynamic function of the determinant of the 3-metric. This choice is computationally cheap, simple to implement, and certain choices of F2 (i.e., 1 + lnγ) can mimic maximal slicing in its singularity avoidance properties [8]. On the other hand, numerical solutions derived from harmonically-sliced foliations can exhibit pathological gauge behavior in the form of coordinate “shocks” or singularities which will affect the accuracy, convergence and stability of solutions [586]. Also, evolutions in which the lapse function is fixed by some analytically prescribed method (either geodesic or near-geodesic) can be unstable, especially for sub-horizon scale perturbations [7Jump To The Next Citation Point].

The mean curvature slicing equation is derived by taking the trace of the extrinsic curvature evolution equation (11View Equation),

[ ] ∇i∇i α = α KijKij + 4πG (ρH + s) + βi∇iK − ∂tK, (17 )
and assuming K = K (t), which can either be specified arbitrarily or determined by imposing a boundary condition on the lapse function after solving Equation (17View Equation) for the quantity α∕∂tK [61]. It is also useful to consider replacing ∂ K t in Equation (17View Equation) with an exponentially driven form as suggested by Eppley [71], to reduce gauge drifting and numerical errors in maximal [25] and mean curvature [7Jump To The Next Citation Point] sliced spacetimes. The mean curvature slicing condition is the most natural one for cosmology as it foliates homogeneous cosmological spacetimes with surfaces of homogeneity. Also, since K represents the convergence of coordinate curves from one slice to the next and if it is constant, then localized caustics (crossing of coordinate curves) and true curvature singularities can be avoided. Whether general inhomogeneous spacetimes can be foliated with constant mean curvature surfaces remains unknown. However, for Gowdy spacetimes with two Killing fields and topology T 3 × R, Isenberg and Moncrief [97] proved that such foliations do exist and cover the entire spacetime.

6.1.3 Symplectic formalism

A different approach to conventional (i.e., 3 + 1 ADM) techniques in numerical cosmology has been developed by Berger and Moncrief [40Jump To The Next Citation Point]. For example, they consider Gowdy cosmologies on T 3 × R with the metric

2 −λ∕2 τ∕2( −2τ 2 2) −τ [ P 2 P ( P 2 −P) 2] ds = e e − e d τ + d 𝜃 + e e dσ + 2e Q dσ dδ + e Q + e dδ , (18 )
where λ, P and Q are functions of 𝜃 and τ, and the coordinates are bounded by 0 ≤ (𝜃,σ, δ) ≤ 2π. The singularity corresponds to the limit τ → ∞. For small amplitudes, P and Q may be identified with + and × polarized gravitational wave components and λ with the background cosmology through which they propagate. An advantage of this formalism is that the initial value problem becomes trivial since P, Q and their first derivatives may be specified arbitrarily (although it is not quite so trivial in more general spacetimes).

Although the resulting Einstein equations can be solved in the usual spacetime discretization fashion, an interesting alternative method of solution is the symplectic operator splitting formulation [40121] founded on recognizing that the second order equations can be obtained from the variation of a Hamiltonian decomposed into kinetic and potential subhamiltonians,

∮ ∮ 1- 2π ( 2 −2P 2) 1- 2π −2τ ( 2 2P 2) H = H1 + H2 = 2 0 d𝜃 πP + e πQ + 2 0 d𝜃 e P ,𝜃 + e Q ,𝜃 . (19 )
The symplectic method approximates the evolution operator by
HΔτ H2 Δτ∕2 H1Δτ H2Δ τ∕2 3 e = e e e + 𝒪 (Δτ ), (20 )
although higher order representations are possible. If the two Hamiltonian components H1 and H2 are each integrable, their solutions can be substituted directly into the numerical evolution to provide potentially more accurate solutions with fewer time steps [36]. This approach is well-suited for studies of singularities if the asymptotic behavior is determined primarily by the kinetic subhamiltonian, a behavior referred to as Asymptotically Velocity Term Dominated (see Section 3.1.2 and [37]).

Symplectic integration methods are applicable to other spacetimes. For example, Berger et al. [39] developed a variation of this approach to explicitly take advantage of exact solutions for scattering between Kasner epochs in Mixmaster models. Their algorithm evolves Mixmaster spacetimes more accurately with larger time steps than previous methods.

6.1.4 Regge calculus model

A unique approach to numerical cosmology (and numerical relativity in general) is the method of Regge Calculus in which spacetime is represented as a complex of 4-dimensional, geometrically flat simplices. The principles of Einstein’s theory are applied directly to the simplicial geometry to form the curvature, action, and field equations, in contrast to the finite difference approach where the continuum field equations are differenced on a discrete mesh.

A 3-dimensional code implementing Regge Calculus techniques was developed recently by Gentle and Miller [81] and applied to the Kasner cosmological model. They also describe a procedure to solve the constraint equations for time asymmetric initial data on two spacelike hypersurfaces constructed from tetrahedra, since full 4-dimensional regions or lattices are used. The new method is analogous to York’s procedure (see [163Jump To The Next Citation Point] and Section 6.3) where the conformal metric, trace of the extrinsic curvature, and momentum variables are all freely specifiable. These early results are promising in that they have reproduced the continuum Kasner solution, achieved second order convergence, and sustained numerical stability. Also, Barnett et al. [29] discuss an implicit evolution scheme that allows local (vertex) calculations for efficient parallelism. However, the Regge Calculus approach remains to be developed and applied to more general spacetimes with complex topologies, extended degrees of freedom, and general source terms.

  Go to previous page Go up Go to next page