There are many ways to write the Einstein equations. The most common is the ADM or 3 + 1 form  which decomposes spacetime into layers of three-dimensional space-like hypersurfaces, threaded by a time-like normal congruence , where we use Greek (Latin) indices to specify spacetime (spatial) quantities. The general spacetime metric is written as
It is worth noting that several alternative formulations of Einstein’s equations have been suggested, including hyperbolic systems  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 . Introducing with so that the determinant of is unity, and , evolution equations can be written in the conformal traceless system for , , , 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 , but can be less accurate than traditional ADM for short term evolutions .
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 , 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.
For cosmological simulations, one typically takes the shift vector to be zero, hence . 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  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 (), algebraic, and mean curvature slicing. The algebraic condition takes the form. 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 .
The mean curvature slicing equation is derived by taking the trace of the extrinsic curvature evolution equation (11),. It is also useful to consider replacing in Equation (17) with an exponentially driven form as suggested by Eppley , to reduce gauge drifting and numerical errors in maximal  and mean curvature  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 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 , Isenberg and Moncrief  proved that such foliations do exist and cover the entire spacetime.
A different approach to conventional (i.e., 3 + 1 ADM) techniques in numerical cosmology has been developed by Berger and Moncrief . For example, they consider Gowdy cosmologies on with the metric
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,. 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 ).
Symplectic integration methods are applicable to other spacetimes. For example, Berger et al.  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.
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  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  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.  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.
© Max Planck Society and the author(s)