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 as 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 time 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 various published numerical cosmological calculations; from the very early Universe to the present; from the purely geometrical dynamics of the vacuum field to the quark-hadron phase transition and the large scale structure of the Universe. There are two major sections: § 2 where brief summaries of early Universe and fully relativistic cosmological calculations are presented; and § 3 which focuses on structure formation in the post-combination epoch and on testing cosmological models against observations.

2 Relativistic Cosmology

The following sections are organized to track the chronological events in the history of our Universe (see figure 1), focusing on five 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) the quark-hadron phase transition; 4) Big Bang nucleosynthesis and the freeze-out of the light elements; and 5) the decoupling or post-combination epoch and large scale structure formation. This particular section summarizes some results on the early Universe (the first four subjects) and a few other topics in relativistic cosmology. The fifth item on the late or post-combination epoch is reserved to a separate section § 3.

Figure 1
figure 1

A chronological timeline of events in the Universe. The main cosmological epochs discussed in this review are the Big Bang singularity at t = 0, inflation at t ∼ 10−34 secs, the quark-hadron phase transition at t ∼ 10−5 secs, Big Bang nucleosynthesis at t ∼ 10−2–102 secs, and the post-combination Universe where matter decouples from radiation at t ∼ 106 years and structures such as galaxies and clusters begin to form at t ∼ 109 years. (Image courtesy of Barry Sanders, LCA/NCSA.)

Full size image

2.1 Singularities

2.1.1 Mixmaster Dynamics

Belinsky, Lifshitz and Khalatnikov (BLK) [12, 13] and Misner [48] 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 [11, 12] representing different sequences of Kasner spacetimes and, because this map is chaotic, the Mixmaster dynamics is presumed to be chaotic as well.

Mixmaster behavior can be studied in the context of Hamiltonian dynamics with a semi-bounded potential (figure 2) arising from the spatial scalar curvature [49]. 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 Liapunov exponents indicated that the Mixmaster system is nonchaotic and in apparent contradiction with the discrete maps [19, 35]. However, it has since been suggested that the usual definition of the Liapunov exponents is not appropriate in this case as it is not invariant under time reparametrizations, and that with a more natural time variable one obtains positive exponents [14, 30], confirming the chaotic nature of the anisotropy bounces.

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.

Full size image

BLK conjectured that local Mixmaster oscillations might be the generic behavior for singularities in more general classes of spacetimes [13]. However, this conjecture has yet to be established.

2.1.2 AVTD behavior

As noted in § 2.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, are characterized by either open or no potentials in the Hamiltonian framework [55], 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 [38] 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 [16, 50] has been applied to investigate the AVTD conjecture in even more general spacetimes, namely T3 × R spacetimes with U (1) symmetry [15]. 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 [15] 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.

2.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 inflation 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.

2.2.1 Plane Symmetry

Kurki-Suonio et al. [41] extended the planar cosmology code of Centrella and Wilson [24, 25] (see § 2.5) to include a scalar field and simulate the onset of inflation in the early Universe starting with an inhomogeneous Higgs field. 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, even large initial fluctuations of the Higgs field do not prevent inflation. They considered two different forms of the potential: a flat Coleman-Weinberg type which is very flat close to the false vacuum and does inflate; and a rounder “ϕ4” type which, for their parameter combinations, does not.

2.2.2 Spherical Symmetry

Goldwirth and Piran [33] 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.

2.2.3 Bianchi V

Anninos et al. [5] 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 fast 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.

2.2.4 Λ + Gravity Waves

In addition to the inflation 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 [58] 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Λ ≤ max(√I) ≤ 80Λ and widths 0.08lH ≤ l ≤ 2.5lH, where I is the invariant constructed from the 3-Riemann tensor and lH = √3/Λ is the horizon scale. They also considered colliding plane waves with amplitudes 40Λ ≤ max(√I) ≤ 125Λ and widths 0.08lH ≤ l ≤ 0.1lH. They find that for any large amplitude and/or small width inhomogeneity in their parameter sets, the nonlinearity of gravity has little effect and the spacetime always evolves into a deSitter spacetime.

2.2.5 Three-D 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 this end, Kurki-Suonio et al. [42] 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 they 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.

2.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 of such a transition might be on the formation of baryon number inhomogeneities, primordial nucleosynthesis, and the generation of primordial magnetic fields.

Rezolla et al. [53] 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 [43] 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.

2.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 generally consistent with the standard model of our Universe (see §3.1). 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 [54] 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 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 [44] 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 [54].

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

2.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 through 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 [22, 23] studied a class of plane symmetric cosmologies representing gravitational inhomogeneities in the form of shocks or discontinuities separating two vacuum expanding Kasner cosmologies. 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 [24, 25] 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 γij = A(t, z)2[1, h(t, z)2, 1], which is maintained in time with a proper choice of shift vector. The hydrodynamic equations are solved using artificial viscosity methods for shock capturing and Barton’s method for monotonic transport [62]. 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. [2, 4], implementing more robust numerical methods and an improved parametric treatment of the initial value problem.

In applications of these codes, Centrella [21] 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. [3] considered small and large perturbations in the axisymmetric Kasner models. 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 NP scalar Ψ2, 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.

2.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 [32] 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 [63] 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.

3 Physical Cosmology

The phrase “physical cosmology” is generally associated with the large (galaxy and cluster) scale structure of the post-combination epoch where gravitational effects are modeled approximately by Newtonian physics on an uniformly expanding, matter dominated Friedmann-Robertson-Lemaître-Walker (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.

Due to the vast body of literature on numerical simulations of the post-combination 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.

3.1 The Standard Model

The isotropic and homogeneous FLRW cosmological model, although certainly incomplete, 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 so that 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,

  • 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.75 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 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 [40] for a more in-depth review of the standard model.

3.2 Newtonian Order Perturbations

The large (galaxy and cluster) scale inhomogeneities are generally described by Newtonian order, scalar type perturbations of the FLRW spacetimes, assuming spatial scales much smaller than the horizon radius, peculiar velocities small compared to the speed of light, and a gravitational potential that is both much smaller than unity (in geometric units) and slowly varying in time. A comprehensive review of the theory of cosmological perturbations can be found in [51].

Many numerical techniques have been developed to solve the hydrodynamic and collisionless particle equations in this limit. For the hydrodynamic equations, the methods range from Eulerian finite difference techniques [8, 18, 56] to unstructured particle methods [29, 34]. For the collisionless particle or dark matter equations, the canonical choices are treecodes or PM and P3M methods [36], although many variants have been developed to optimize computational performance and accuracy, including grid and particle refinement methods. The reader is referred to [39] for a comparison of several different numerical methods applied to a common problem of structure formation. A method for solving non-equilibrium, multi-species chemical reaction flows together with the hydrodynamic equations in a background FLRW model is described in [1, 9].

3.3 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 due to the gravitational redshifting of the photons through the Sachs-Wolfe effect, even well after the photons decouple from the matter at redshift z ∼ 103. Also, if the intergalactic medium reionizes sometime after the decoupling, say from an early generation of stars, the increased rate of Thomson scattering off the free electrons will erase subhorizon 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 §3.3.3) must be incorporated self-consistently into the numerical models and to high accuracy in order to resolve the weak signals.

3.3.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 [46] and the references cited in [37]). 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. [6] 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 ΔT/T = δz/(1 + z), where 1 + z = (kuuμ)e/(kuuμ)r, and the photon wave vector ku 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. [6] find the parameters for this model are severely constrained by the COBE data such that Ω0h2 ≈ 1, where Ω0 and h are the density and Hubble parameters.

3.3.2 Effects of Reionization

In models where the InterGalactic Medium (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. [61]. 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 the decaying HDM, 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, 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. 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.

3.3.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. [60] 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 separate 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 [6] 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 and are therefore important contributors to the composite images (see figure 3 for a visual example of secondary anisotropy effects).

Figure 3
figure 3

Temperature map from a 240h−1 Mpc simulation of the large scale structure in a critically closed universe using 105 photons in a 80 × 80 window. The signal is displayed at a redshift z = 0. 425 with Δz = 0.025, and shows the secondary anisotropy from the intrinsic Rees-Sciama effect and the proper motion of a cluster of galaxies moving across the sky. The minimum (maximum) values in the image are ΔT/T = −4.3 × 10−7 (4.5 × 10−7). The proper motion effect leaves a clear signature in the center of the image, forming a dipolar pattern with the clusters at the center.

Full size image

In addition to the effects discussed in the previous paragraphs, many other sources of secondary anisotropies (such as gravitational lensing, the Vishniac effect and the Sunyaev-Zel’dovich effect) can also be significant. See reference [37] for a more complete list and thorough discussion of the different sources of CMBR anisotropies.

3.4 Gravitational Lensing

Observations of gravitational lenses [57] 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. [20] 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, where h is the Hubble parameter. Core radii of this size can have an important effect on the probability of multiple imaging. Flores and Primack [31] 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.

3.5 Lyα Forest

The Lyα forest represents the optically thin (at the Lyman edge) component of Quasar Absorption Systems (QAS), a collection of absorption features in QSO 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 Lyα forest have been performed recently [28, 47, 64], and all have been able to fit the observations remarkably well, including the column density and Doppler width distributions, the size of absorbers [26], 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 Lyα 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, including fluctuations in under-dense regions, spheroidal minihalos, and filaments extending over scales of a few megaparsecs (figure 4). However, it is not yet clear what effect different cosmological models have on these systems.

Figure 4
figure 4

Distribution of the gas density at redshift z = 3 from a numerical hydrodynamics simulation of the Lyα forest. The simulation adopted a CDM spectrum of primordial density fluctuations, normalized to the 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 cosmic density (characteristic of typical filamentary structures) 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. Notice that there is fine structure everywhere. 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.

Full size image

3.6 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.

3.6.1 Internal Structure

Thomas et al. [59] 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. [52]. 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 [27]. 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.

3.6.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. [10] computed the evolution of the cluster mass function in five different cosmological model simulations and find that the number of high mass (Comalike) clusters in flat, low σ8 models (ie. 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).

3.6.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. [17] 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 [7] 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 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.