1 Introduction

Numerical investigations of cosmological spacetimes can be categorized into two broad classes of calculations, distinguished by their computational (or even philosophical) goals: 1) geometrical and mathematical principles of cosmological models, and 2) physical and astrophysical cosmology. In the former, the emphasis is on the geometric framework in which astrophysical processes occur, namely the cosmological expansion, shear, and singularities of the many models allowed by the theory of general relativity. In the latter, the emphasis is on the cosmological and astrophysical processes in the real or observable Universe, and the quest to determine the model which best describes our Universe. The former is pure in the sense that it concerns the fundamental nonlinear behavior of the Einstein equations and the gravitational field. The latter is more complex, since it addresses the composition, organization and dynamics of the Universe from the small scales (fundamental particles and elements) to the large (galaxies and clusters of galaxies). However, the distinction is not always so clear, and geometric effects in the spacetime curvature can have significant consequences for the evolution and observation of matter distributions.

Any comprehensive model of cosmology must therefore include nonlinear interactions between different matter sources and spacetime curvature. A realistic model of the Universe must also cover large dynamical spatial and temporal scales, extreme temperature and density distributions, and highly dynamic atomic and molecular matter compositions. In addition, due to all the varied physical processes of cosmological significance, one must draw from many disciplines of physics to model curvature anisotropies, gravitational waves, electromagnetic fields, nucleosynthesis, particle physics, hydrodynamic fluids, etc. These phenomena are described in terms of coupled nonlinear partial differential equations and must be solved numerically for general inhomogeneous spacetimes. The situation appears extremely complex, even with current technological and computational advances. As a result, the codes and numerical methods that have been developed to date are designed to investigate very specific problems with either idealized symmetries or simplifying assumptions regarding the metric behavior, the matter distribution/composition or the interactions among the matter types and spacetime curvature.

It is the purpose of this article to review published numerical cosmological calculations addressing issues from the very early Universe to the present; from the purely geometrical dynamics of the initial singularity to the large scale structure of the Universe. There are three major sections: § 2 where a brief overview is presented of various defining events ocurring throughout the history of our Universe and in the context of the standard model; § 3 where brief summaries of early Universe and fully relativistic cosmological calculations are presented; and § 4 which focuses on structure formation in the post-recombination epoch and on testing cosmological models against observations.

2 Background

2.1 A Brief Chronology

With current observational constraints, the physical state of our Universe, as understood in the context of the standard, or Friedmann-Robertson-Lemaître-Walker (FLRW) model, can be crudely extrapolated back to the Planck epoch ∼ 10−43 seconds after the Big Bang, beyond which the classical theory of general relativity is invalid due to quantum corrections. At the earliest times, the Universe was a plasma of relativistic particles consisting of quarks, leptons, gauge bosons, and Higgs bosons represented by scalar fields with interaction and symmetry regulating potentials. It is believed that several spontaneous symmetry breaking (SSB) phase transitions transpired in the early Universe as it expanded and cooled, including: the grand unification transition (GUT) at ∼ 10−34 seconds after the Big Bang (Here, the strong nuclear force split off from the weak and electromagnetic forces. This also marks an era of inflationary expansion and the origin of matter-antimatter asymmetry through baryon, charge conjugation, and charge + parity violating interactions and nonequilibrium effects.); the electroweak (EW) SSB transition at ∼ 10−11 secs (when the weak nuclear force split from the electromagnetic force); and the chiral (QCD) symmetry breaking transition at ∼ 10−5 secs during which quarks condensed into hadrons. The most stable hadrons (baryons, or protons and neutrons comprised of three quarks) survived the subsequent period of baryon-antibaryon annihilations, which continued until the Universe cooled to the point at which new baryon-antibaryon pairs could no longer be produced. This resulted in a large number of photons and relatively few surviving baryons. A period of primordial nucleosynthesis followed from ∼ 10−2 to ∼ 102 secs during which light element abundances were synthesized to form 24% helium with trace amounts of deuterium, tritium, helium-3, and lithium.

By ∼ 1011 secs, the matter density became equal to the radiation density, identifying the start of the current matter-dominated era and the beginning of structure formation. Later, at ∼ 1013 secs (3 × 105 years), the free ions and electrons combine to form atoms, decoupling the matter from the radiation field as the Universe continued to expand and cool. This decoupling or post-recombination epoch marks the surface of last scattering and the boundary of the observable (via photons) Universe. Assuming a hierarchical (CDM-like) structure formation scenario, the subsequent development of our Universe is characterized by the growth of structures with increasing size. For example, the first stars are likely to form at t ∼ 108 years from molecular gas clouds when the Jeans mass of the background baryonic fluid is approximately 104 M, as indicated in Figure 1. This epoch of pop III star generation is followed by the formation of galaxies at t ∼ 109 years and then galaxy clusters.

Figure 1
figure 1

Schematic depicting the general sequence of events in the post-recombination Universe. The solid and dotted lines potentially track the Jeans mass of the average baryonic gas component from the recombination epoch at z ∼ 103 to the current time. A residual ionization fraction of nH+/nH ∼ 10−4 following recombination allows for Compton interactions with photons to z ∼ 200, during which the Jeans mass remains constant at 105 M. The Jeans mass then decreases as the Universe expands adiabatically until the first collapsed structures form sufficient amounts of hydrogen molecules to trigger a cooling instability and produce pop III stars at z ∼ 20. Star formation activity can then reheat the Universe and raise the mean Jeans mass to above 108 M. This reheating could affect the subsequent development of structures such as galaxies and the observed Lyman-alpha clouds.

2.2 Successes of the Standard Model

The isotropic and homogeneous FLRW cosmological model has been so successful in describing the observable Universe that it is commonly referred to as the “standard model”. Furthermore, and to its credit, the model is relatively simple; it allows for calculations and predictions to be made of the very early Universe, including primordial nucleosynthesis at 10−2 seconds after the Big Bang, and even particle interactions approaching the Planck scale at 10−43 seconds. At present, observational support for the standard model includes:

  • the expansion of the Universe as verified by the redshifts in galaxy spectra and quantified by measurements of the Hubble constant H0 = 100h km sec−1 Mpc−1 with Hubble parameter 0.5 < h < 1;

  • the deceleration parameter observed in distant galaxy spectra (although uncertainties about galactic evolution, intrinsic luminosities, and standard candles prevent even a crude estimate);

  • the large scale isotropy and homogeneity of the Universe based on temperature anisotropy measurements of the microwave background radiation and peculiar velocity fields of galaxies (although the light distribution from bright galaxies seems more tenuous);

  • the age of the Universe which yields roughly consistent estimates between the look-back time to the Big Bang in the FLRW model and observed data such as the oldest stars, radioactive elements, and cooling of white dwarf stars;

  • the cosmic microwave background radiation suggests that the Universe began from a hot Big Bang and the data is consistent with a black body at temperature 2.7 K;

  • the abundance of light elements such as 2H, 3He, 4He and 7Li, as predicted from the FLRW model, are consistent with observations and provides a bound on the baryon density and baryon-to-photon ratio;

  • the present mass density, as determined from measurements of luminous matter and galactic rotation curves, can be accounted for by the FLRW model with a single density parameter Ω to specify the metric topology;

  • the distribution of galaxies and larger scale structures can be reproduced by numerical simulations in the context of inhomogeneous perturbations of the FLRW models.

Because of these successes, most work in the field of physical cosmology (see § 4) has utilized the standard model as the background spacetime in which the large scale structure evolves, with the ambition to further constrain the cosmological parameters and structure formation scenarios through numerical simulations. The reader is referred to [34] for a more in-depth review of the standard model.

3 Relativistic Cosmology

This section is organized to track the chronological events in the history of the early or relativistic Universe; we focus mainly on four defining moments in our current understanding of the classical picture: 1) the Big Bang singularity and the dynamics of the very early Universe; 2) inflation and its generic nature; 3) QCD phase transitions; and 4) primordial nucleosynthesis and the freeze-out of the light elements. The late or post-recombination epoch is reserved to a separate section § 4.

3.1 Singularities

3.1.1 Mixmaster Dynamics

Belinsky, Lifshitz and Khalatnikov (BLK) [11, 12] and Misner [42] discovered that the Einstein equations in the vacuum homogeneous Bianchi type IX (or Mixmaster) cosmology exhibit complex behavior and are sensitive to initial conditions as the Big Bang singularity is approached. In particular, the solutions near the singularity are described qualitatively by a discrete map [10, 11] representing different sequences of Kasner spacetimes

$$d{s^2} = - d{t^2} + {t^{2{p_1}}}d{x^2} + {t^2}^{{p_2}}d{y^2} + {t^2}^{{p_3}}d{z^2},$$
(1)

with time changing exponents pi, but otherwise constrained by \({p_1} + {p_2} + {p_3} = p_1^2 + p_2^2 + p_3^2 = 1\). Because this discrete mapping of Kasner epochs is chaotic, the Mixmaster dynamics is presumed to be chaotic as well.

Mixmaster behavior can be studied in the context of Hamiltonian dynamics, with the Hamiltonian [43]

$$2{\mathcal H} = - p_\Omega ^2 + p_ + ^2 + p_ - ^2 + {e^{4\Omega}}(V - 1),$$
(2)

and a semi-bounded potential arising from the spatial scalar curvature (whose level curves are plotted in Figure 2)

$$V = 1 + {1 \over 3}{e^{- 8{\beta _ +}}} + {2 \over 3}{e^{4{\beta _ +}}}\left[ {\cosh (4\sqrt 3 {\beta _-}) - 1} \right] - {4 \over 3}{e^{- 2{\beta _ +}}}\cosh (2\sqrt 3 {\beta _ -}),$$
(3)

where eΩ and β± are the scale factor and anisotropies, and pΩ and p± are the corresponding conjugate variables.

Figure 2
figure 2

Contour plot of the Bianchi type IX potential V, where β± are the anisotropy canonical coordinates. Seven level surfaces are shown at equally spaced decades ranging from 10−1 to 105. For large isocontours (V > 1), the potential is open and exhibits a strong triangular symmetry with three narrow channels extending to spatial infinity. For V < 1, the potential closes and is approximately circular for β± ≪ 1.

A solution of this Hamiltonian system is an infinite sequence of Kasner epochs with parameters that change when the phase space trajectories bounce off the potential walls, which become exponentially steep as the system evolves towards the singularity. Several numerical calculations of the dynamical equations arising from (2) and (3) have indicated that the Liapunov exponents of the system vanish, in apparent contradiction with the discrete maps [17, 31]. However, it has since been shown that the usual definition of the Liapunov exponents is ambiguous in this case as it is not invariant under time reparametrizations, and that with a different time variable one obtains positive exponents [13, 27].

Although BLK conjectured that local Mixmaster oscillations might be the generic behavior for singularities in more general classes of spacetimes [12], this conjecture has yet to be established. However, Weaver et al. [56] have, through numerical investigations, established evidence that Mixmaster dynamics can occur in (at least a restricted class of) inhomogeneous spacetimes which generalize the Bianchi type VI0 with a magnetic field and two-torus symmetry.

3.1.2 AVTD Behavior

As noted in § 3.1.1, an interesting and (as yet) unresolved issue is whether or not the generic cosmological singularity is locally of a Mixmaster or BLK type, with a complex oscillatory behavior as the singularity is approached. Most of the other Bianchi models, including the Kasner solutions (1), are characterized by either open or no potentials in the Hamiltonian framework [48], and exhibit essentially monotonic or Asymptotically Velocity Term Dominated (AVTD) behavior, the opposite dynamics to the complex BLK oscillations.

Considering inhomogeneous spacetimes, Isenberg and Moncrief [33] proved that the singularity in the polarized Gowdy model is AVTD type. This has also been conjectured to be the case for more general Gowdy models, and numerical simulations of one-dimensional plane symmetric Gowdy spacetimes support the notion. Furthermore, a symplectic numerical method has been applied to investigate the AVTD conjecture in even more general spacetimes, namely T3 × R spacetimes with U(1) symmetry [14]. Although there are concerns about the solutions due to difficulties in resolving the steep spatial gradients near the singularity, the numerical calculations find no evidence of BLK oscillations. Berger [14] attributes this to several possibilities: 1) the BLK conjecture is false; 2) the simulations have not been run long enough; 3) Mixmaster behavior is present but hidden in the variables; or 4) the initial data is insufficiently generic. In any case, further investigations are needed to confirm either the BLK or AVTD conjectures.

3.2 Inflation

The inflation paradigm is frequently invoked to explain the flatness (Ω ≈ 1) of the Universe, attributing it to an era of exponential expansion at about 10−34 seconds after the Big Bang. The mechanism of inflation is generally taken to be an effective cosmological constant from the dominant stress-energy of the inflaton scalar field that regulates GUT symmetry breaking, particle creation, and the reheating of the Universe. In order to study whether inflation can occur for arbitrary anisotropic and inhomogeneous data, many numerical simulations have been carried out with different symmetries, matter types and perturbations. A sample of such calculations are described in the following paragraphs.

3.2.1 Plane Symmetry

Kurki-Suonio et al. [35] extended the planar cosmology code of Centrella and Wilson [22, 23] (see § 3.5) to include a scalar field and simulate the onset of inflation in the early Universe with an inhomogeneous Higgs field and a perfect fluid with a radiation equation of state p = ρ/3. Their results suggest that whether inflation occurs or not can be sensitive to the shape of the potential ϕ. In particular, if the shape is flat enough and satisfies the slow-roll conditions (essentially upper bounds on ∂V/∂ϕ and 2V/∂ϕ2 [34] near the false vacuum ϕ ∼ 0), even large initial fluctuations of the Higgs field do not prevent inflation. They considered two different forms of the potential: a Coleman-Weinberg type

$$V(\phi) = \lambda {\phi ^4}\left[ {\ln \left({{{{\phi ^2}} \over {{\sigma ^2}}}} \right) - {1 \over 2}} \right] + {{\lambda {\sigma ^4}} \over 2}$$
(4)

which is very flat close to the false vacuum and does inflate; and a rounder “ϕ4” type

$$V(\phi) = \lambda {({\phi ^2} - {\sigma ^2})^2}$$
(5)

which, for their parameter combinations, does not.

3.2.2 Spherical Symmetry

Goldwirth and Piran [30] studied the onset of inflation with inhomogeneous initial conditions for closed, spherically symmetric spacetimes containing a massive scalar field and radiation field sources (described by a massless scalar field). In all the cases they considered, the radiation field damps quickly and only an inhomogeneous massive scalar field remains to inflate the Universe. They find that small inhomogeneities tend to reduce the amount of inflation and large initial inhomogeneities can even suppress the onset of inflation. Their calculations indicate that the scalar field must have “suitable” initial values over a domain of several horizon lengths in order for inflation to begin.

3.2.3 Bianchi V

Anninos et al. [4] investigated the simplest Bianchi model (type V) background that admits velocities or tilt in order to address the question of how the Universe can choose a uniform reference frame at the exit from inflation, since the deSitter metric does not have a preferred frame. They find that inflation does not isotropize the Universe in the short wavelength limit. However if inflation persists, the wave behavior eventually freezes in and all velocities go to zero at least as rapidly as tanh βR−1, where β is the relativistic tilt angle (a measure of velocity), and R is a typical scale associated with the radius of the Universe. Their results indicate that the velocities eventually go to zero as inflation carries all spatial variations outside the horizon, and that the answer to the posed question is that memory is retained and the Universe is never really deSitter.

3.2.4 Gravity Waves + Cosmological Constant

In addition to the inflaton field, one can consider other sources of inhomogeneity, such as gravitational waves. Although linear waves in deSitter space will decay exponentially and disappear, it is unclear what will happen if strong waves exist. Shinkai & Maeda [52] investigated the cosmic no-hair conjecture with gravitational waves and a cosmological constant in 1D plane symmetric vacuum spacetimes, setting up Gaussian pulse wave data with amplitudes \(0.02\Lambda \le \max (\sqrt I) \le 80\Lambda\) and widths 0.08lHl ≤ 2.5lH, where I is the invariant constructed from the 3-Riemann tensor and \({l_H} = \sqrt {3/\Lambda}\) is the horizon scale. They also considered colliding plane waves with amplitudes \(40\Lambda \le \max (\sqrt I) \le 125\Lambda\) and widths 0.08lHl ≤ 0.1lH. They find that for any large amplitude and/or small width with inhomogeneity in their parameter sets, the nonlinearity of gravity has little effect, and the spacetime always evolves into a deSitter spacetime.

3.2.5 Three Dimensional Inhomogeneous Spacetimes

The previous paragraphs discussed results from highly symmetric spacetimes, but the possibility of inflation remains to be established for more general inhomogeneous and nonperturbative data. To address this issue, Kurki-Suonio et al. [36] investigated fully three-dimensional inhomogeneous spacetimes with a chaotic inflationary potential V(ϕ) = Λϕ4/4. They considered basically two types of runs: small and large scale. In the small scale run, the grid length was initially set equal to the Hubble length so the inhomogeneities are well inside the horizon and the dynamical time scale is shorter than the expansion or Hubble time. As a result, the perturbations oscillate and damp while the field evolves and the spacetime inflates. In the large scale run, the inhomogeneities are outside the horizon and do not oscillate. They maintain their shape without damping and, because larger values of ϕ lead to faster expansion, the inhomogeneity in the expansion becomes steeper in time, since the regions of large ϕ and high inflation stay correlated. Both runs have sufficient inflation to solve the flatness problem.

3.3 Quark-Hadron Phase Transition

The standard picture of cosmology assumes that a phase transition (associated with chiral symmetry breaking after the electroweak transition) occurred at approximately 10−5 seconds after the Big Bang to convert a plasma of free quarks and gluons into hadrons. Although this transition can be of significant cosmological importance, it is not known with certainty whether it is of first order or higher, and what the astrophysical consequences might be on the subsequent state of the Universe. For example, the transition may give rise to significant baryon number inhomogeneities which can influence the outcome of primordial nucleosynthesis as evidenced in the distribution and averaged light element abundances. The QCD transition and baryon inhomogeneities may also play a significant and potentially observable role in the generation of primordial magnetic fields.

Rezolla et al. [46] considered a first order phase transition and the nucleation of hadronic bubbles in a supercooled quark-gluon plasma, solving the relativistic Lagrangian equations for disconnected and evaporating quark regions during the final stages of the phase transition. They numerically investigated a single isolated quark drop with an initial radius large enough so that surface effects can be neglected. The droplet evolves as a self-similar solution until it evaporates to a sufficiently small radius that surface effects break the similarity solution and increase the evaporation rate. Their simulations indicate that, in neglecting long-range energy and momentum transfer (by electromagnetically interacting particles) and assuming that the baryon number is transported with the hydrodynamical flux), the baryon number concentration is similar to what is predicted by chemical equilibrium calculations.

Kurki-Suonio and Laine [37] studied the growth of bubbles and the decay of droplets using a spherically symmetric code that accounts for a phenomenological model of the microscopic entropy generated at the phase transition surface. Incorporating the small scale effects of the finite wall width and surface tension, but neglecting entropy and baryon flow through the droplet wall, they demonstrate the dynamics of nucleated bubble growth and quark droplet decay. They also find that evaporating droplets do not leave behind a global rarefaction wave to dissipate any previously generated baryon number inhomogeneity.

3.4 Nucleosynthesis

Observations of the light elements produced during Big Bang nucleosynthesis following the quark/hadron phase transition (roughly 10−2 − 102 seconds after the Big Bang) are in good agreement with the standard model of our Universe (see § 2.2). However, it is interesting to investigate other more general models to assert the role of shear and curvature on the nucleosynthesis process.

Rothman and Matzner [47] considered primordial nucleosynthesis in anisotropic cosmologies, solving the strong reaction equations leading to 4He. They find that the concentration of 4He increases with increasing shear; this is due to time scale effects and the competition between dissipation and enhanced reaction rates from photon heating and neutrino blue shifts. Their results have been used to place a limit on anisotropy at the epoch of nucleosynthesis. Kurki-Suonio and Matzner [38] extended this work to include 30 strong 2-particle reactions involving nuclei with mass numbers A ≤ 7, and to demonstrate the effects of anisotropy on the cosmologically significant isotopes 2H, 3He, 4He and 7Li as a function of the baryon to photon ratio. They conclude that the effect of anisotropy on 2H and 3He is not significant, and the abundances of 4He and 7Li increase with anisotropy in accord with [47].

Furthermore, it is possible that neutron diffusion, the process whereby neutrons diffuse out from regions of very high baryon density just before nucleosynthesis, can affect the neutron to proton ratio in such a way as to enhance deuterium and reduce 4He compared to a homogeneous model. However, plane symmetric, general relativistic simulations with neutron diffusion [39] show that the neutrons diffuse back into the high density regions once nucleosynthesis begins there — thereby wiping out the effect. As a result, although inhomogeneities influence the element abundances, they do so at a much smaller degree then previously speculated. The numerical simulations also demonstrate that, because of the back diffusion, a cosmological model with a critical baryon density cannot be made consistent with helium and deuterium observations, even with substantial baryon inhomogeneities and the anticipated neutron diffusion effect.

3.5 Plane Symmetric Gravitational Waves

Gravitational waves are an inevitable product of the general Einstein equations, and one can expect a wide spectrum of wave signals propagating throughout our Universe due to shear anisotropies, primordial metric and matter fluctuations, collapsing matter structures, ringing black holes, and colliding neutron stars, for example. The discussion here is restricted to the pure vacuum field dynamics and the fundamental nonlinear behavior of gravitational waves in numerically generated cosmological spacetimes.

Centrella and Matzner [20, 21] studied a class of plane symmetric cosmologies representing gravitational inhomogeneities in the form of shocks or discontinuities separating two vacuum expanding Kasner cosmologies (1). By a suitable choice of parameters, the constraint equations can be satisfied at the initial time with an Euclidean 3-surface and an algebraic matching of parameters across the different Kasner regions that gives rise to a discontinuous extrinsic curvature tensor. They performed both numerical calculations and analytical estimates using a Green’s function analysis to establish and verify (despite the numerical difficulties in evolving discontinuous data) certain aspects of the solutions, including gravitational wave interactions, the formation of tails, and the singularity behavior of colliding waves in expanding vacuum cosmologies.

Shortly thereafter, Centrella and Wilson [22, 23] developed a more general plane symmetric code for cosmology, adding also hydrodynamic sources. In order to simplify the resulting differential equations, they adopted a diagonal 3-metric of the form

$${\gamma _{ij}} = {A^2}(t, \, z)\,\,\left({\begin{array}{*{20}c} 1 & 0 & 0 \\ 0 & {{h^2}(t, \, z)} & 0 \\ 0 & 0 & 1 \\ \end{array}} \right),$$
(6)

which is maintained in time with a proper choice of shift vector. The metric (6) allows an overall conformal factor A to simplify the initial value problem, and a dynamical transverse wave component in the variable h. The hydrodynamic equations are solved using artificial viscosity methods for shock capturing and Barton’s method for monotonic transport [57]. The evolutions are fully constrained (solving both the momentum and Hamiltonian constraints at each time step) and use the mean curvature slicing condition. This work was subsequently extended by Anninos et al. [1, 3], implementing more robust numerical methods and an improved parametric treatment of the initial value problem.

In applications of these codes, Centrella [19] investigated nonlinear gravity waves in Minkowski space and compared the full numerical solutions against a first order perturbation solution to benchmark certain numerical issues such as numerical damping and dispersion. A second order perturbation analysis was also used to model the transition into the nonlinear regime. Anninos et al. [2] considered small and large perturbations in the two degenerate Kasner models: p1 = p3 = 0 or 2/3, and p2 = 1 or −1/3 respectively, where pi are parameters in the Kasner metric (1). Carrying out a second order perturbation expansion and computing the Newman-Penrose (NP) scalars, Riemann invariants and Bel-Robinson vector, they demonstrated, for their particular class of spacetimes, that the nonlinear behavior is in the Coulomb (or background) part represented by the leading order term in the NP scalar and not in the gravitational wave component. For standing-wave perturbations, the dominant second order effects in their variables are an enhanced monotonic increase in the background expansion rate, and the generation of oscillatory behavior in the background spacetime with frequencies equal to the harmonics of the first order standing-wave solution.

3.6 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 [29] 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 [58] where the conformal metric, trace of the extrinsic curvature, and momentum variables are all freely specifiable.

Although additional work is needed to apply (and develop) Regge Calculus techniques to more general spacetimes, the early results of Gentle and Miller are promising. In particular, their evolutions have reproduced the continuum Kasner solution, achieved second order convergence, and sustained numerical stability.

4 Physical Cosmology

The phrase “physical cosmology” is generally associated with the large (galaxy and cluster) scale structure of the post-recombination epoch where gravitational effects are modeled approximately by Newtonian physics on an uniformly expanding, matter dominated FLRW background. A discussion of the large scale structure is included in this review since any viable model of our Universe which allows a regime where strongly general relativistic effects are important must match onto the weakly relativistic (or Newtonian) regime. Also, since certain aspects of this regime are directly observable, one can hope to constrain or rule out various cosmological models and/or parameters, including the density (Ω0), Hubble (H0 = 100h km sec−1 Mpc−1), and cosmological (Λ) constants.

Due to the vast body of literature on numerical simulations of the post-recombination epoch, it is possible to mention only a small fraction of all the published papers. Hence, the following summary is limited to cover just a few aspects of computational physical cosmology, and in particular those that can potentially be used to discriminate between cosmological model parameters, even within the realm of the standard model.

4.1 Microwave Background

The Cosmic Microwave Background Radiation (CMBR), which is a direct relic of the early Universe, currently provides the deepest probe of cosmological structures and imposes severe constraints on the various proposed matter evolution scenarios and cosmological parameters. Although the CMBR is a unique and deep probe of both the thermal history of the early Universe and the primordial perturbations in the matter distribution, the associated anisotropies are not exclusively primordial in nature. Important modifications to the CMBR spectrum can arise from large scale coherent structures, even well after the photons decouple from the matter at redshift z ∼ 103, due to the gravitational redshifting of the photons through the Sachs-Wolfe effect arising from potential gradients [49, 5]

$${{\Delta T} \over T} = {\phi _e} - {\phi _r} - \int\nolimits_r^e {{{2\overrightarrow l \cdot \nabla \phi} \over a}dt,}$$
(7)

where the integral is evaluated from the emission (e) to reception (r) points along the spatial photon paths \(\vec l,\phi\) is the gravitational potential, ΔT/T defines the temperature fluctuations, and a(t) is the cosmological scale factor in the standard FLRW metric. Also, if the intergalactic medium (IGM) reionizes sometime after the decoupling, say from an early generation of stars, the increased rate of Thomson scattering off the free electrons will erase sub-horizon scale temperature anisotropies, while creating secondary Doppler shift anisotropies. To make meaningful comparisons between numerical models and observed data, all of these effects (and others, see for example § 4.1.3) must be incorporated self-consistently into the numerical models and to high accuracy in order to resolve the weak signals.

4.1.1 Ray-Tracing Methodology

Many computational analyses based on linear perturbation theory have been carried out to estimate the temperature anisotropies in the sky (for example see [40] and the references cited in [32]). Although such linearized approaches yield reasonable results, they are not well-suited to discussing the expected imaging of the developing nonlinear structures in the microwave background. An alternative ray-tracing approach has been developed by Anninos et al. [5] to introduce and propagate individual photons through the evolving nonlinear matter structures. They solve the geodesic equations of motion and subject the photons to Thomson scattering in a probabilistic way and at a rate determined by the local density of free electrons in the model. Since the temperature fluctuations remain small, the equations of motion for the photons are treated as in the linearized limit, and the anisotropies are computed according to

$${{\Delta T} \over T} = {{\delta z} \over {1 + z}},$$
(8)

where

$$1 + z = {{{{({k^\mu}{u_\mu})}_e}} \over {{{({k^\mu}{u_\mu})}_r}}},$$
(9)

and the photon wave vector kμ and matter rest frame four-velocity uμ are evaluated at the emission (e) and reception (r) points. Applying their procedure to a Hot Dark Matter (HDM) model of structure formation, Anninos et al. [5] find the parameters for this model are severely constrained by COBE data such that Ω0h2 ≈ 1, where Ω0 and h are the density and Hubble parameters.

4.1.2 Effects of Reionization

In models where the IGM does not reionize, the probability of scattering after the photon-matter decoupling epoch is low, and the Sachs-Wolfe effect dominates the anisotropies at angular scales larger than a few degrees. However if reionization occurs, the scattering probability increases substantially and the matter structures, which develop large bulk motions relative to the comoving background, induce Doppler shifts on the scattered CMBR photons and leave an imprint of the surface of last scattering. The induced fluctuations on subhorizon scales in reionization scenarios can be a significant fraction of the primordial anisotropies, as observed by Tuluie et al. [55]. They considered two possible scenarios of reionization: A model that suffers early and gradual (EG) reionization of the IGM as caused by the photoionizing UV radiation emitted by decaying neutrinos, and the late and sudden (LS) scenario as might be applicable to the case of an early generation of star formation activity at high redshifts. Considering the HDM model with Ω0 = 1 and h = 0.55, which produces CMBR anisotropies above current COBE limits when no reionization is included (see § 4.1.1), they find that the EG scenario effectively reduces the anisotropies to the levels observed by COBE and generates smaller Doppler shift anisotropies than the LS model, as demonstrated in Figure 3.

Figure 3
figure 3

The top two images represent temperature fluctuations (i.e., ΔT/T) due to the Sachs-Wolfe effect and Doppler shifts in a standard critically closed HDM model with no reionization and baryon fractions 0.02 (plate 1, 4° × 4°, rms=2.8 × 10−5) and 0.2 (plate 2, 8° × 8°, rms=3.4 × 10−5). The bottom two plates image fluctuations in an “early and gradual” reionization scenario of decaying neutrinos with baryon fraction 0.02 (plate 3, 4° × 4°, rms=1.3 × 10−5; and plate 4, 8° × 8°, rms = 1.4 × 10−5).

The LS scenario of reionization is not able to reduce the anisotropy levels below the COBE limits, and can even give rise to greater Doppler shifts than expected at decoupling.

4.1.3 Secondary Anisotropies

Additional sources of CMBR anisotropy can arise from the interactions of photons with dynamically evolving matter structures and nonstatic gravitational potentials. Tuluie et al. [54] considered the impact of nonlinear matter condensations on the CMBR in Ω0 ≤ 1 Cold Dark Matter (CDM) models, focusing on the relative importance of secondary temperature anisotropies due to three different effects: 1) time-dependent variations in the gravitational potential of nonlinear structures as a result of collapse or expansion; 2) proper motion of nonlinear structures such as clusters and superclusters across the sky; and 3) the decaying gravitational potential effect from the evolution of perturbations in open models. They applied the ray-tracing procedure of [5] to explore the relative importance of these secondary anisotropies as a function of the density parameter Ω0 and the scale of matter distributions. They find that the secondary temperature anisotropies are dominated by the decaying potential effect at large scales, but that all three sources of anisotropy can produce signatures of order ΔT/T ∼ 10−6 as shown in Figure 4.

Figure 4
figure 4

The top two images represent the proper motion and Rees-Sciama effects in the CMBR for a critically closed CDM model (upper left), together with the corresponding column density of voids and clusters over the same region (upper right). The bottom two images show the secondary anisotropies dominated here by the decaying potential effect in an open cosmological model (bottom left), together with the corresponding gravitational potential over the same region (bottom right). The rms fluctuations in both cases are on the order of ±5 × 10−7, though the open model carries a somewhat larger signature.

In addition to the effects discussed in the previous paragraphs, many other sources of secondary anisotropies (such as gravitational lensing, the Vishniac effect accounting for matter velocities and flows into local potential wells, and the Sunyaev-Zel’dovich distortions from the Compton scattering of CMB photons in the hot cluster medium) can also be significant. See reference [32] for a more complete list and thorough discussion of the different sources of CMBR anisotropies.

4.2 Gravitational Lensing

Observations of gravitational lenses [50] provide measures of the abundance and strength of nonlinear potential fluctuations along the lines of sight to distant objects. Since these calculations are sensitive to the gravitational potential, they may be more reliable than galaxy and velocity field measurements as they are not subject to the same ambiguities associated with biasing effects. Also, since different cosmological models predict different mass distributions, especially at the higher redshifts, lensing calculations can potentially be used to confirm or discard competing cosmological models.

As an alternative to solving the computationally demanding lens equations, Cen et al. [18] developed a more efficient scheme to identify regions with surface densities capable of generating multiple images accurately for splittings larger than three arcseconds. They applied this technique to a standard CDM model with Ω0 = 1, and found that this model predicts more large angle splittings (> 8″) than are known to exist in the observed Universe. Their results suggest that the standard CDM model should be excluded as a viable model of our Universe. A similar analysis for a flat low density CDM model with a cosmological constant \(({\Omega _0} = 0.3,\,\,\Lambda/3H_0^2 = 0{.}7)\) suggests a lower and perhaps acceptable number of lensing events. However, an uncertainty in their studies is the nature of the lenses at and below the resolution of the numerical grid. They model the lensing structures as simplified Singular Isothermal Spheres (SIS) with no distinctive cores.

Large angle splittings are generally attributed to larger structures such as clusters of galaxies, and there are indications that clusters have small but finite core radii rcore ∼ 20 − 30h−1 kpc. Core radii of this size can have an important effect on the probability of multiple imaging. Flores and Primack [28] considered the effects of assuming two different kinds of splitting sources: isothermal spheres with small but finite core radii \(\rho \propto {({r^2} + r_{core}^2)^{- 1}}\), and spheres with a Hernquist density profile ρr−1(r + a)−3, where rcore ∼ 20 − 30h−1 kpc and a ∼ 300h−1 kpc. They find that the computed frequency of large-angle splittings, when using the nonsingular profiles, can potentially decrease by more than an order of magnitude relative to the SIS case and can bring the standard CDM model into better agreement with observations. They stress that lensing events are sensitive to both the cosmological model (essentially the number density of lenses) and to the inner lens structure, making it difficult to probe the models until the structure of the lenses, both observationally and numerically, is better known.

4.3 First Star Formation

In CDM cosmogonies, the fluctuation spectrum at small wavelengths has a logarithmic dependence at mass scales smaller than 108 solar masses, which indicates that all small scale fluctuations in this model collapse nearly simultaneously in time. This leads to very complex dynamics during the formation of these first structures. Furthermore, the cooling in these fluctuations is dominated by the rotational/vibrational modes of hydrogen molecules that were able to form using the free electrons left over from recombination and the ones produced by strong shock waves as catalysts. The first structures to collapse may be capable of producing population III stars and have a substantial influence on the subsequent thermal evolution of the intergalactic medium, as suggested by Figure 1, due to the radiation emitted by the first generation stars as well as supernova driven winds. To know the subsequent fate of the Universe and which structures will survive or be destroyed by the UV background, it is first necessary to know when the first stars formed.

Ostriker and Gnedin [45] have carried out high resolution numerical simulations of the reheating and reionization of the Universe due to star formation bursts triggered by molecular hydrogen cooling. Accounting for the chemistry of the primeval hydrogen/helium plasma, self-shielding of the gas, radiative cooling, and a phenomenological model of star formation, they find that two distinct star populations form: the first generation pop III from H2 cooling prior to reheating at redshift z ≥ 14; and the second generation pop II at z < 10 when the virial temperature of the gas clumps reaches 104 K and hydrogen line cooling becomes efficient. Star formation slows in the intermittent epoch due to the depletion of H2 by photo-destruction and reheating. In addition, the objects which formed pop III stars also initiate pop II sequences when their virial temperatures reach 104 K through continued mass accretion.

4.4 Lyman-alpha Forest

The Lyman-alpha forest represents the optically thin (at the Lyman edge) component of Quasar Absorption Systems (QAS), a collection of absorption features in quasar spectra extending back to high redshifts. QAS are effective probes of the matter distribution and the physical state of the Universe at early epochs when structures such as galaxies are still forming and evolving. Although stringent observational constraints have been placed on competing cosmological models at large scales by the COBE satellite and over the smaller scales of our local Universe by observations of galaxies and clusters, there remains sufficient flexibility in the cosmological parameters that no single model has been established conclusively. The relative lack of constraining observational data at the intermediate to high redshifts (0 < z < 5), where differences between competing cosmological models are more pronounced, suggests that QAS can potentially yield valuable and discriminating observational data.

Several combined N-body and hydrodynamic numerical simulations of the Lyman forest have been performed recently [26, 41, 60], and all have been able to fit the observations remarkably well, including the column density and Doppler width distributions, the size of absorbers [24], and the line number evolution. Despite the fact that the cosmological models and parameters are different in each case, the simulations give similar results, provided that the proper ionization bias is used (bion ≡ (ΩBh2)2/Γ, where ΩB is the baryonic density parameter, h is the Hubble parameter and Γ is the photoionization rate at the hydrogen Lyman edge). A theoretical paradigm has thus emerged from these calculations in which Lyman-alpha absorption lines originate from the relatively small scale structure in pregalactic or intergalactic gas through the bottom-up hierarchical formation picture in CDM-like universes. The absorption features originate in structures exhibiting a variety of morphologies commonly found in numerical simulations (see Figure 5), including fluctuations in underdense regions, spheroidal minihalos, and filaments extending over scales of a few megaparsecs.

Figure 5
figure 5

Distribution of the gas density at redshift z = 3 from a numerical hydrodynamics simulation of the Lyman-alpha forest with a CDM spectrum normalized to second year COBE observations, a Hubble parameter of h = 0.5, a comoving box size of 9.6 Mpc, and baryonic density of Ωb = 0.06 composed of 76% hydrogen and 24% helium. The region shown is 2.4 Mpc (proper) on a side. The isosurfaces represent baryons at ten times the mean density and are color coded to the gas temperature (dark blue = 3 × 104 K, light blue = 3 × 105 K). The higher density contours trace out isolated spherical structures typically found at the intersections of the filaments. A single random slice through the cube is also shown, with the baryonic overdensity represented by a rainbow-like color map changing from black (minimum) to red (maximum). The He+ mass fraction is shown with a wire mesh in this same slice. To emphasize fine structure in the minivoids, the mass fraction in the overdense regions has been rescaled by the gas overdensity wherever it exceeds unity.

4.5 Galaxy Clusters

Clusters of galaxies are the largest gravitationally bound systems known to be in quasi-equilibrium. This allows for reliable estimates to be made of their mass as well as their dynamical and thermal attributes. The richest clusters, arising from 3σ density fluctuations, can be as massive as 1014 − 1015 solar masses, and the environment in these structures is composed of shock heated gas with temperatures of order 107 − 108 degrees Kelvin which emits thermal bremsstrahlung and line radiation at X-ray energies. Also, because of their spatial size ∼ 1h−1 Mpc and separations of order 50h−1 Mpc, they provide a measure of nonlinearity on scales close to the perturbation normalization scale 8h−1 Mpc. Observations of the substructure, distribution, luminosity, and evolution of galaxy clusters are therefore likely to provide signatures of the underlying cosmology of our Universe, and can be used as cosmological probes in the easily observable redshift range 0 ≤ z ≤ 1.

4.5.1 Internal Structure

Thomas et al. [53] have investigated the internal structure of galaxy clusters formed in high resolution N-body simulations of four different cosmological models, including standard, open, and flat but low density universes. They find that the structure of relaxed clusters is similar in the critical and low density universes, although the critical density models contain relatively more disordered clusters due to the freeze-out of fluctuations in open universes at late times. The profiles of relaxed clusters are very similar in the different simulations, since most clusters are in a quasi-equilibrium state inside the virial radius and follow the universal density profile of Navarro et al. [44]. There does not appear to be a strong cosmological dependence in the profiles, as suggested by previous studies of clusters formed from pure power law initial density fluctuations [25]. However, because more young and dynamically evolving clusters are found in critical density universes, Thomas et al. suggest that it may be possible to discriminate among the density parameters by looking for multiple cores in the substructure of the dynamic cluster population. They note that a statistical population of 20 clusters can distinguish between open and critically closed universes.

4.5.2 Number Density Evolution

The evolution of the number density of rich clusters of galaxies can be used to compute Ω0 and σ8 (the power spectrum normalization on scales of 8h−1 Mpc) when numerical simulation results are combined with the constraint \({\sigma _8}\Omega _0^{0.5} \approx 0.5\), derived from observed present-day abundances of rich clusters. Bahcall et al. [9] computed the evolution of the cluster mass function in five different cosmological model simulations and found that the number of high mass (Comalike) clusters in flat, low σ8 models (i.e., the standard CDM model with σ8 ≈ 0.5) decreases dramatically by a factor of approximately 103 from z = 0 to z ≈ 0.5. For low Ω0, high σ8 models, the data results in a much slower decrease in the number density of clusters over the same redshift interval. Comparing these results to observations of rich clusters in the real Universe, which indicate only a slight evolution of cluster abundances to redshifts z ≈ 0.5 − 1, they conclude that critically closed standard CDM and Mixed Dark Matter (MDM) models are not consistent with the observed data. The models which best fit the data are the open models with low bias (Ω0 = 0.3 ± 0.1 and σ8 = 0.85 ± 0.5), and flat low density models with a cosmological constant (Ω0 = 0.34 ± 0.13 and Ω0 + Λ = 1).

4.5.3 X-Ray Luminosity Function

The evolution of the X-ray luminosity function, and the size and temperature distribution of rich clusters of galaxies are all potentially important discriminants of cosmological models. Bryan et al. [16] investigated these properties in a high resolution numerical simulation of a standard CDM model normalized to COBE. Although the results are highly sensitive to grid resolution (see [6] for a discussion of the effects from resolution constraints on the properties of rich clusters), their primary conclusion, that the standard CDM model predicts too many bright X-ray emitting clusters and too much integrated X-ray intensity, is robust, since an increase in resolution will only exaggerate these problems.

4.6 Cosmological Sheets

Cosmological sheets, or pancakes, form as overdense regions collapse preferentially along one axis. Originally studied by Zel’dovich [59] in the context of neutrino-dominated cosmologies, sheets are ubiquitous features in nonlinear structure formation simulations of CDM-like models with baryonic fluid, and manifest on a spectrum of length scales and formation epochs. Gas collapses gravitationally into flattened sheet structures, forming two plane parallel shock fronts that propagate in opposite directions, heating the infalling gas. The heated gas between the shocks then cools radiatively and condenses into galactic structures. Sheets are characterized by essentially five distinct components: the preshock inflow, the postshock heated gas, the strongly cooling/recombination front separating the hot gas from the cold, the cooled postshocked gas, and the unshocked adiabatically compressed gas at the center. Several numerical calculations [15, 51, 8] have been performed of these systems which include the baryonic fluid with hydrodynamical shock heating, ionization, recombination, dark matter, thermal conductivity, and radiative cooling (Compton, bremstrahlung, and atomic line cooling), in both one and two spatial dimensions to ascert the significance of each physical process and to compute the fragmentation scale.

In addition, it is well known that gas which cools to 1 eV through hydrogen line cooling will likely cool faster than it can recombine. This nonequilibrium cooling increases the number of electrons and ions (compared to the equilibrium case) which, in turn, increases the concentrations of H and H+, the intermediaries that produce hydrogen molecules H2. If large concentrations of molecules form, excitations of the vibrational/rotational modes of the molecules can efficiently cool the gas to well below 1 eV, the minimum temperature expected from atomic hydrogen line cooling. Because the gas cools isobarically, the reduction in temperature results in an even greater reduction in the Jeans mass, and the bound objects which form from the fragmentation of H2 cooled cosmological sheets may be associated with massive stars or star clusters. Anninos and Norman [7] have carried out 1D and 2D high resolution numerical calculations to investigate the role of hydrogen molecules in the cooling instability and fragmentation of cosmological sheets, considering the collapse of perturbation wavelengths from 1 Mpc to 10 Mpc. They find that for the more energetic (long wavelength) cases, the mass fraction of hydrogen molecules reaches nH2/nH ∼ 3 × 10−3, which cools the gas to 4 × 10−3 eV and results in a fragmentation scale of 9 × 103 M. This represents reductions of 50 and 103 in temperature and Jeans mass respectively when compared, as in Figure 6, to the equivalent case in which hydrogen molecules were neglected.

Figure 6
figure 6

Two different model simulations of cosmological sheets are presented: (a) a six species model including only atomic line cooling, and (b) a nine species model including also hydrogen molecules. The evolution sequences in the images show the baryonic overdensity and gas temperature at three redshifts following the initial collapse at z = 5. In each figure, the vertical axis is 32 kpc long (parallel to the plane of collapse) and the horizontal axis extends to 4 Mpc on a logarithmic scale to emphasize the central structures. Differences in the two cases are observed in the cold pancake layer and the cooling flows between the shock front and the cold central layer. When the central layer fragments, the thickness of the cold gas layer in the six (nine) species case grows to 3 (0.3) kpc and the surface density evolves with a dominant transverse mode corresponding to a scale of approximately 8(1) kpc. Assuming a symmetric distribution of matter along the second transverse direction, the fragment masses are approximately 107 (105) solar masses.

5 Conclusion

This review is intended to provide a flavor of the variety of numerical cosmological calculations that have been performed to investigate different phenomena throughout the history of our Universe. The topics discussed range from the strong field dynamical behavior of spacetime geometry at early times near the Big Bang singularity and the epoch of inflation, to the late time evolution of large scale matter fluctuations and the formation of clusters of galaxies. Although many problems have been addressed, the complexity of the varied physical processes in the extreme time, spatial, and dynamical scales of our Universe guarantees that many more interesting and unresolved issues remain. Indeed, even the background cosmological model which best describes our Universe is essentially unknown. The topology of our Universe, the generic singularity behavior, the fundamental nonlinear field and gravitational wave interactions, the existence and role of a cosmological constant, the correct structure formation scenario, and the origin and spectrum of primordial fluctuations, for example, are uncertain. However the field of numerical cosmology has matured, in the development of computational techniques, the modeling of microphysics, and in taking advantage of current computing technologies, to the point that it is now possible to perform reliable comparisons of numerical models with observed data.