The Einstein equations

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

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,

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,

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 [147 , 31 ] 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 [31 ], but can be less accurate than traditional ADM for short term evolutions [6 ].

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 [31 , 6 , 13 ] (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 [60 , 61 ]. 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

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 [5 , 86 ]. 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 [7

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

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 (17

) for the quantity $α∕∂tK$ [61 ]. It is also useful to consider replacing $∂ K t$ in Equation (17

) with an exponentially driven form as suggested by Eppley [71 ], to reduce gauge drifting and numerical errors in maximal [25 ] and mean curvature [7

] 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 [40 ]. For example, they consider Gowdy cosmologies on $T 3 × R$ with the metric

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 [40 , 121 ] founded on recognizing that the second order equations can be obtained from the variation of a Hamiltonian decomposed into kinetic and potential subhamiltonians,

The symplectic method approximates the evolution operator by

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 [163 ] 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.