Computational Cosmology: From the Early Universe to the Large Scale Structure

Abstract

In order to account for the observable Universe, any comprehensive theory or model of cosmology must draw from many disciplines of physics, including gauge theories of strong and weak interactions, the hydrodynamics and microphysics of baryonic matter, electromagnetic fields, and spacetime curvature, for example. Although it is difficult to incorporate all these physical elements into a single complete model of our Universe, advances in computing methods and technologies have contributed significantly towards our understanding of cosmological models, the Universe, and astrophysical processes within them. A sample of numerical calculations (and numerical methods applied to specific issues in cosmology are reviewed in this article: from the Big Bang singularity dynamics to the fundamental interactions of gravitational waves; from the quark-hadron phase transition to the large scale structure of the Universe. The emphasis, although not exclusively, is on those calculations designed to test different models of cosmology against the observed Universe.

1 Introduction

Numerical investigations of cosmological spacetimes can be categorized into two broad classes of calculations, distinguished by their computational goals: (1) geometrical and mathematical principles of cosmological models, and (ii) physical and astrophysical cosmology. In the former, the emphasis is on the geometric framework in which astrophysical processes occur, for example cosmological expansion, topological singularities, geometrodynamics in general, and classification characteristics or invariants 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 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 problems 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: Section 2 where a brief overview is presented of various defining events occurring throughout the history of our Universe and in the context of the standard model, Section 3 where reviews of early Universe and relativistic cosmological calculations are presented, and Section 4 which focuses on structure formation in the post-recombination epoch and on testing cosmological models against observations. Following the summary paragraphs in Section 5, an appendix in Section 6 presents the basic Einstein equations, kinematic considerations, matter source equations with curvature, and the equations of perturbative physical cosmology on background isotropic models. References to numerical methods are also supplied and reviewed for each case.

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 Friedman-Lemaitre-Robertson-Walker (FLRW) model, can be crudely extrapolated back to ∼ 10⁻⁴³ seconds after the Big Bang, before which the classical description of general relativity is expected to give way to a quantum theory of gravity. As the time-line in Figure 1 shows, the Universe was a plasma of relativistic particles at the earliest times consisting of quarks, leptons, gauge bosons, and Higgs bosons represented by scalar fields with interaction and symmetry regulating potentials.

Figure 1

A historical time-line showing the major evolutionary stages of our Universe according to the standard model, from the earliest moments of the Planck era to the present. The horizontal axis represents logarithmic time in seconds (or equivalently energy in electron-Volts or temperature in Kelvin, and the solid red line roughly models the radius of the Universe, showing the different rates of expansion at different times: exponential during inflation, shallow power law during the radiation dominated era, and a somewhat steeper power law during the current matter dominated phase.

It is believed that several spontaneous symmetry breaking (SSB) phase transitions occured in the early Universe as it expanded and cooled, including the grand unification transition (GUT) at ∼ 10–34 s after the Big Bang in which 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 s when the weak nuclear force split from the electromagnetic force; and the chiral or quantum chromodynamic (QCD) symmetry breaking transition at ∼ 10⁻⁵ s 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. Topological defects, defined as stable configurations of matter in the symmetric (high temperature) phase, may persist after any of the phase transitions described above to influence the subsequent evolution of matter structures. The nature of the defects is determined by the phase transition and the symmetry properties of the matter, and some examples include domain walls, cosmic strings, monopoles, and textures.

A period of primordial nucleosynthesis followed from ∼ 10⁻² to ∼ 10² during which light element abundances were synthesized to form 24% helium with trace amounts of deuterium, tritium, helium-3, and lithium. Observations of these relative abundances represent the earliest confirmation of the standard model. It is also during this stage that neutrinos (produced from proton-proton and proton-photon interactions, and from the collapse or quantum evaporation/annihilation of topological defects) stopped interacting with other matter, such as neutrons, protons, and photons. Neutrinos that existed at this time separated from these other forms of matter and traveled freely through the Universe at very high velocities, near the speed of light.

By ∼ 10¹¹ s, the matter density became equal to the radiation density as the Universe continued to expand, identifying the start of the current matter-dominated era and the beginning of structure formation. Later, at ∼ 10¹³ s (3 x 10⁵ yr), the free ions and electrons combined to form atoms, effectively decoupling the matter from the radiation field as the Universe cooled. This decoupling or post-recombination epoch marks the surface of last scattering and the boundary of the observable (via photons) Universe, and plays an important role in the history of the Cosmic Microwave Background Radiation (CMBR). Assuming a hierarchical Cold Dark Matter (CDM) 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 have formed at t ∼ 10⁸y from molecular gas clouds when the Jeans mass of the background baryonic fluid was approximately 10⁴ M_⊙, as indicated in Figure 2. This epoch of pop III star generation is followed by the formation of galaxies at t ∼ 10⁹ yr and subsequently galaxy clusters. Though somewhat controversial, estimates of the current age of our Universe range from 10 to 20 Gy, with a present-day linear structure scale radius of about 8 h⁻¹ Mpc, where h is the Hubble parameter (compared to 2–3 Mpc typical for the virial radius of rich galaxy clusters).

Figure 2

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 ∼ 10³ to the current time. A residual ionization fraction of n_H+/n_H ∼ 10⁻⁴ following recombination allows for Compton interactions with photons to z ∼ 200, during which the Jeans mass remains constant at 10⁵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 10⁸M_⊙. This reheating could affect the subsequent development of structures such as galaxies and the observed Ly_α clouds.

2.2 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 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 s. 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 H₀ = 100h km s⁻¹ Mpc⁻¹, where 0.5 ≤ h ≤ 1 is the Hubble constant;

• the deceleration parameter observed in distant galaxy spectra (although uncertainties about galactic evolution, intrinsic luminosities, and standard candles prevent an accurate 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 is somewhat contradictory);

• 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 mostly isotropic model and a black body at temperature 2.7 K;

• CMBR precision measurements suggest best fit cosmological parameters in accord with the critical density standard model;

• the abundance of light elements such as ²H, ³He, ⁴ He, and ⁷Li, as predicted from the FLRW model, is consistent with observations, provides a bound on the baryon density and baryon-to-photon ratio, and is the earliest confirmation of the standard model;

• 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;

• the detection of dark energy from observations of supernovae is generally consistent with accepted FLRW model parameters and cold dark matter + cosmological constant numerical structure formation models.

Because of these remarkable agreements between observation and theory, most work in the field of physical cosmology (see Section 4) has utilized the standard model as the background spacetime in which the large scale structure evolves, with the ambition to further constrain parameters and structure formation scenarios through numerical simulations. The most widely accepted form of the model is described by a set of dimensionless density parameters which sum to

$${\Omega _b} + {\Omega _d} + {\Omega _\gamma } + {\Omega _\Lambda } = {\Omega _0},$$

((1))

where the different components measure the present mean baryon density Ω_b, the dark matter density Ω_d, the radiation energy Ω_γ, and the dark energy Ω_Λ. The relative contributions of each source and their sum Ω₀ (which determines the topological curvature of the model) remains one of the most important issues in modern computational and observational cosmology. The reader is referred to [104] for a more in-depth review of the standard model, and to [128, 154] for a summary of observed cosmological parameter constraints and best fit “concordance” models. Peebles and Ratra [133] provide a comprehensive literature survey and an excellent review of the standard model, cosmological tests, and the evidence for dark energy and the cosmological constant.

However, some important unanswered questions about the standard model concern the nature of the special conditions that produced an essentially geometrically flat Universe that is also homogeneous and isotropic to a high degree over large scales. In an affort to address these questions, it should be noted that many other cosmological models can be constructed with a late time behavior similar enough to the standard model that it is difficult to exclude them with absolute certainty. Consider, for example, the collection of homogeneous but arbitrarily anisotropic vacuum spacetimes known as the Bianchi models [141, 69]. There are nine unique models in this family of cosmologies, ranging from simple Bianchi I models representing the Kasner class of spacetimes (the flat FLRW model, sometimes referred to as Type I-homogeneous, belongs to this group), to the more complex and chaotic Bianchi IX or Mixmaster model (which also includes the closed FLRW model, Type IX-homogeneous). Several of these models will be discussed in the first section on relativistic cosmology (Section 3) dealing pre-dominately with the early Universe, where the models differ the most.

Anisotropic solutions, as well as more general (and in some cases exact) inhomogeneous cosmological models with initial singularities, can isotropize through anisotropic damping mechanisms and adiabatic cooling by the expansion of the Universe to resemble the standard FLRW model at late times. Furthermore, if matter is included in these spacetimes, the effective energy of anisotropy, which generally dominates matter energy at early times, tends to become less important over time as the Universe expands. The geometry in these matter-filled anisotropic spacetimes thus evolves towards an isotropic state. Quantum mechanical effects have also been proposed as a possible anisotropy damping mechanism that takes place in the early Universe to convert vacuum geometric energy to radiation energy and create particles. All of this suggests that the early time behavior and effects of local and global geometry are highly uncertain, despite the fact that the standard FLRW model is generally accepted as accurate enough for the late time description of our Universe.

Further detailed information on homogeneous (including Blanch!) universes, as well as more general classes of inhomogeneous cosmological models can be found in [105, 158, 70].

3 Relativistic Cosmology

This section is organized to track the chronological events in the history of the early or relativistic Universe, focusing mainly on four defining moments: (1) the Big Bang singularity and the dynamics of the very early Universe, (ii) inflation and its generic nature, (iii) QCD phase transitions, and (iv) 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) [32, 33] and Misner [119] 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 [30, 32] 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_2}}}d{z^2},$$

((2))

with time changing exponents p_i, but otherwise constrained by p₁ + p₂ + p₃ = p ²₁ + p ²₂ + p ²₃ - 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 a Hamiltonian [120]

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

((3))

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

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

((4))

where e^α and β_± are the scale factor and anisotropies, and p_α and p_± are the corresponding conjugate variables. 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.

Figure 3

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⁻¹ to 10⁵. 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.

Some of the earliest numerical simulations of this dynamical system were performed by Matzner, Shepley, and Warren [116], and Moser, Matzner and Ryan [123] who followed phase space trajectories and provided examples of solutions for various initial conditions and special cases. Several, more recent, numerical calculations of the equations arising from Equations (3) and (4) have indicated that the Liapunov exponents of the system vanish, in apparent contradiction with the discrete maps [53, 89], and putting into question the characterization of Mixmaster dynamics as chaotic. 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 [35, 73]. Also, coordinate independent methods using fractal basin boundaries to map basins of attraction in the space of initial conditions indicates Mixmaster spacetimes to be chaotic [64].

Although BLK conjectured that local Mixmaster oscillations might be the generic behavior for singularities in more general classes of spacetimes [33], it is only recently that this conjecture has begun to be supported by numerical evidence (see Section 3.1.2 and [37]).

3.1.2 AVTD vs. BLK oscillatory behavior

As noted in Section 3.1.1, an interesting and important issue in classical cosmology is whether or not the generic Big Bang singularity is locally of a Mixmaster or BLK type, with complex oscillatory behavior as the singularity is approached. Many of the Bianchi models, including the Kasner solutions (2), are characterized by either open or no potentials in the Hamiltonian framework [141], and exhibit essentially monotonic or Asymptotically Velocity Term Dominated (AVTD) behavior.

Considering inhomogeneous spacetimes, Isenberg and Moncrief [98] proved that the singularity in the polarized Gowdy model is AVTD type, as are more general polarized T² symmetric cosmologies [38]. Early numerical studies using symplectic methods confirmed AVTD behavior and found no evidence of BLK oscillations, even in T³ × R spacetimes with U(1) symmetry [36] (although there were concerns about the solutions due to difficulties in resolving steep spatial gradients near the singularity [36], which were addressed later by Hern and Stewart [87] for the Gowdy T³ models).

However, Weaver et al. [160] established the first evidence through numerical simulations that Mixmaster dynamics can occur in a class of inhomogeneous spacetimes which generalize the Bianchi type VI₀ model with a magnetic field and two-torus symmetry. Berger and Moncrief [41, 42] also demonstrated that U(1) symmetric vacuum cosmologies exhibit local Mixmaster dynamics consistent with the BLK conjecture, despite numerical difficulties in resolving steep gradients (which they managed by enforcing the Hamiltonian constraint and spatially averaging the solutions). Another more recent example supporting the BLK conjecture is provided by Garfinkle [79], who finds local oscillating behavior approaching the singularity in closed vacuum (but otherwise generic) spacetimes with no assumed symmetry in the initial data.

3.2 Inflation

The inflation paradigm is frequently invoked to explain the flatness (Ω₀ ≈ 1 in the context of the FLRW model) and nearly isotropic nature of the Universe at large scales, attributing them to an era of exponential expansion at about 10⁻³⁴ s after the Big Bang. This expansion acts as a strong dampening mechanism to random curvature or density fluctuations, and may be a generic attractor in the space of initial conditions. An essential component needed to trigger this inflationary phase is a scalar or inflaton field φ representing spin zero particles. The vacuum energy of the field acts as an effective cosmological constant that regulates GUT symmetry breaking, particle creation, and the reheating of the Universe through an interaction potential V (φ) derived from the form of symmetry breaking that occurs as the temperature of the Universe decreases.

Early analytic studies focused on simplified models, treating the interaction potential as flat near its local maximum where the field does not evolve significantly and where the formal analogy to an effective cosmological constant approximation is more precise. However, to properly account for the complexity of inflaton fields, the full dynamical equations (see Section 6.2.2) must be considered together with consistent curvature, matter and other scalar field couplings. Also, many different theories of inflation and vacuum potentials have been proposed (see, for example, a recent review by Lyth and Motto [113] and an introductory article by Liddle [111]), and it is not likely that simplified models can elucidate the full nonlinear complexity of scalar fields (see Section 3.3) nor the generic nature of inflation.

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 is described in the following paragraphs.

3.2.1 Plane symmetry

Kurki-Suonio et al. [106] extended the planar cosmological code of Centrella and Wilson [59, 60] (see Section 3.6.1) 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, where p is the pressure and ρ is the energy density. 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 ∂²V/∂φ² [104] 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 with interaction strength λ and distance between true and false vacua σ

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

((5))

which is very flat near the false vacuum and does inflate; and a rounder “φ⁴” type

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

((6))

which, for their parameter combinations, does not.

3.2.2 Spherical symmetry

Goldwirth and Piran [83] 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 (local conditions for which an equivalent homogeneous Universe will inflate) over a domain of several horizon lengths in order to trigger inflation.

3.2.3 Bianchi V

Anninos et al. [14] 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 de Sitter 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⁻¹, 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 de Sitter.

3.2.4 Gravitational waves + cosmological constant

In addition to the inflaton field, one can consider other sources of inhomogeneity, such as gravitational waves. Although linear waves in de Sitter space will decay exponentially and disappear, it is unclear what will happen if strong waves exist. Shinkai and Maeda [148] 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 \leqslant max (\sqrt I ) \leqslant 80\Lambda$ and widths 0.08l_H ≤ l ≤ 2.5l_H, 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 \leqslant \max (\sqrt I ) \leqslant 125\Lambda$ and widths 0.08l_H ≤ l ≤ 0.1l_H. They find that for any large amplitude or small width inhomogeneity in their parameter sets, the nonlinearity of gravity has little effect and the spacetime always evolves towards de Sitter.

3.2.5 3D 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. In an effort to address this issue, Kurki-Suonio et al. [107] investigated fully three-dimensional inhomogeneous spacetimes with a chaotic inflationary potential V(φ) = λφ⁴/4. They considered basically two types of runs: small and large scale. For the small scale runs, the grid dimensions were 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. For the large scale runs, 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 produce enough inflation to solve the flatness problem.

3.3 Chaotic scalar field dynamics

Many studies of cosmological models generally reduce complex physical systems to simplified or even analytically integrable systems. In sufficiently simple models the dynamical behavior (or fate) of the Universe can be predicted from a given set of initial conditions. However, the Universe is composed of many different nonlinear interacting fields including the inflaton field which can exhibit complex chaotic behavior. For example, Cornish and Levin [63] consider a homogenous isotropic closed FLRW model with various conformal and minimally coupled scalar fields (see Section 6.2.2). They find that even these relatively simple models exhibit chaotic transients in their early preinflationary evolution. This behavior in exiting the Planck era is characterized by fractal basins of attraction, with attractor states being to (1) inflate forever, (ii) inflate over a short period of time then collapse, or (iii) collapse without inflating. Monerat et al. [122] investigated the dynamics of the pre-inflationary phase of the Universe and its exit to inflation in a closed FLRW model with radiation and a minimally coupled scalar field. They observe complex behavior associated with saddle-type critical points in phase space that give rise to deSitter attractors with multiple chaotic exits to inflation that depend on the structure of the scalar field potential. These results suggest that distinctions between exits to inflation may be manifested in the process of reheating and as a selected spectrum of inhomogeneous perturbations influenced by resonance mechanisms in curvature oscillations. This could possibly lead to fractal patterns in the density spectrum, gravitational waves, cosmic microwave background radiation (CMBR) field, or galaxy distribution that depend on specific details including the number of fields, the nature of the fields, and their interaction potentials.

Chaotic behavior can also be found in more general applications of scalar field dynamics. Anninos et al. [20] investigated the nonlinear behavior of colliding kink-antikink solitons or domain walls described by a single real scalar field with self-interaction potential λ(φ² -1 )². Domain walls can form as topological defects during the spontaneous symmetry breaking period associated with phase transitions, and can seed density fluctuations in the large scale structure. For collisional time scales much smaller than the cosmological expansion, they find that whether a kink-antikink collision results in a bound state or a two-soliton solution depends on a fractal structure observed in the impact velocity parameter space. The fractal structure arises from a resonance condition associated with energy exchanges between translational modes and internal shape-mode oscillations. The impact parameter space is a complex self-similar fractal composed of sequences of different n-bounce (the number of bounces or oscillations in the transient semi-coherent state) reflection windows separated by regions of oscillating bion states (see Figure 4). They compute the Lyapunov exponents for parameters in which a bound state forms to demonstrate the chaotic nature of the bion oscillations.

Figure 4

Fractal structure of the transition between reflected and captured states for colliding kinkantikink solitons in the parameter space of impact velocity for a λ(φ²-1)² scalar field potential. The top image (a) shows the 2-bounce windows in dark with the rightmost region (v/c > 0.25) representing the single-bounce regime above which no captured state exists, and the leftmost white region (v/c < 0.19) representing the captured state below which no reflection windows exist. Between these two marker velocities, there are 2-bounce reflection states of decreasing widths separated by regions of bion formation. Zooming in on the domain outlined by the dashed box, a self-similar structure is apparent in the middle image (b), where now the dark regions represent 3-bounce windows of decreasing widths. Zooming in once again on the boundaries of these 3-bounce windows, a similar structure is found as shown in the bottom image (c) but with 4-bounce reflection windows. This pattern of self-similarity characterized by n-bounce windows is observed at all scales investigated numerically.

3.4 Quark-hadron phase transition

The standard picture of cosmology assumes that a phase transition (associated with chiral symmetry breaking following the electroweak transition) occurred at approximately 10⁻⁵ s 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 play a potentially observable role in the generation of primordial magnetic fields. The QCD transition may also give rise to important baryon number inhomogeneities which can affect the distribution of light element abundances from primordial Big Bang nucleosynthesis. The distribution of baryons may be influenced hydrodynamically by the competing effects of phase mixing and phase separation, which arise from any inherent instability of the interface surfaces separating regions of different phase. Unstable modes, if they exist, will distort phase boundaries and induce mixing and diffusive homogenization through hydrodynamic turbulence [102, 112, 95, 4, 137].

In an effort to support and expand theoretical studies, a number of one-dimensional numerical simulations have been carried out to explore the behavior of growing hadron bubbles and decaying quark droplets in simplified and isolated geometries. For example, Rezolla et al. [138] considered a first order phase transition and the nucleation of hadronic bubbles in a supercooled quarkgluon plasma, solving the relativistic Lagrangian equations for disconnected and evaporating quark regions during the final stages of the phase transition. They investigated numerically 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 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 [108] studied the growth of bubbles and the decay of droplets using a one-dimensional 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 finite wall width and surface tension, but neglecting entropy and baryon flow through the droplet wall, they simulate the process by which nucleating bubbles grow and evolve to a similarity solution. They also compute the evaporation of quark droplets as they deviate from similarity solutions at late times due to surface tension and wall effects.

Ignatius et al. [96] carried out parameter studies of bubble growth for both the QCD and electroweak transitions in planar symmetry, demonstrating that hadron bubbles reach a stationary similarity state after a short time when bubbles grow at constant velocity. They investigated the stationary state using numerical and analytic methods, accounting also for preheating caused by shock fronts.

Fragile and Anninos [76] performed two-dimensional simulations of first order QCD transitions to explore the nature of interface boundaries beyond linear stability analysis, and determine if they are stable when the full nonlinearities of the relativistic scalar field and hydrodynamic system of equations are accounted for. They used results from linear perturbation theory to define initial fluctuations on either side of the phase fronts and evolved the data numerically in time for both deflagration and detonation configurations. No evidence of mixing instabilities or hydrodynamic turbulence was found in any of the cases they considered, despite the fact that they investigated the parameter space predicted to be potentially unstable according to linear analysis. They also investigated whether phase mixing can occur through a turbulence-type mechanism triggered by shock proximity or disruption of phase fronts. They considered three basic cases (see image sequences in Figures 5, 6, and 7 below): interactions between planar and spherical deflagration bubbles, collisions between planar and spherical detonation bubbles, and a third case simulating the interaction between both deflagration and detonation systems initially at two different thermal states. Their results are consistent with the standard picture of cosmological phase transitions in which hadron bubbles expand as spherical condensation fronts, undergoing regular (non-turbulent) coalescence, and eventually leading to collapsing spherical quark droplets in a medium of hadrons. This is generally true even in the detonation cases which are complicated by greater entropy heating from shock interactions contributing to the irregular destruction of hadrons and the creation of quark nuggets.

Figure 5

Image sequence of the scalar field from a 2D calculation showing the interaction of two deflagration systems (one planar wall propagating from the right side, and one spherical bubble nucleating from the center. The physical size of the grid is set to 1000 x 1000 fm and resolved by 512 x 512 zones. The run time of the simulation is about two sound crossing times, where the sound speed is $c/\sqrt 3$, so the shock fronts leading the condensing phase fronts travel across the grid twice. The hot quark (cold hadron) phases have smaller (larger) scalar field values and are represented by black (color) in the colormap.

Figure 6

Image sequence of the scalar field from a 2D calculation showing the interaction of two detonation systems (one planar wall propagating from the right side, and one spherical bubble nucleating from the center. The physical size of the grid is set to 1000 x 1000 fm and resolved by 1024 x 1024 zones. The run time of the simulation is about two sound crossing times.

Figure 7

Image sequence of the scalar field from a 2D calculation showing the interaction of shock and rarefaction waves with a deflagration wall (initiated at the left side and a detonation wall (starting from the right. A shock and rarefaction wave travel to the right and left, respectively, from the temperature discontinuity located initially at the grid center (the right half of the grid is at a higher temperature. The physical size of the domain is set to 1806.1 x 451.53 fm and resolved by 2048 x 512 zones. The run time of the simulation is about two sound crossing times.

However, Fragile and Anninos also note a deflagration ‘instability’ or acceleration mechanism evident in their third case for which they assume an initial thermal discontinuity in space separating different regions of nucleating hadron bubbles. The passage of a rarefaction wave (generated at the thermal discontinuity) through a slowly propagating deflagration can significantly accelerate the condensation process, suggesting that the dominant modes of condensation in an early Universe which super-cools at different rates within causally connected domains may be through supersonic detonations or fast moving (nearly sonic) deflagrations. A similar speculation was made by Kamionkowski and Freese [102] who suggested that deflagrations become unstable to perturbations and are converted to detonations by turbulent surface distortion effects. However, in the simulations, deflagrations are accelerated not from turbulent mixing and surface distortion, but from enhanced super-cooling by rarefaction waves. In multi-dimensions, the acceleration mechanism can be exaggerated further by upwind phase mergers due to transverse flow, surface distortion, and mode dissipation effects, a combination that may result in supersonic front propagation speeds, even if the nucleation process began as a slowly propagating deflagration.

3.5 Nucleosynthesis

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

Rothman and Matzner [140] considered primordial nucleosynthesis in anisotropic cosmologies, solving the strong reaction equations leading to ⁴He. They find that the concentration of ⁴He 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 [109] 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 ²H, ³He, ⁴He and ⁷Li as a function of the baryon to photon ratio. They conclude that the effect of anisotropy on ²H and ³He is not significant, and the abundances of ⁴He and ⁷Li increase with anisotropy in accord with [140].

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 ⁴He compared to a homogeneous model. However, plane symmetric, general relativistic simulations with neutron diffusion [110] 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.6 Cosmological gravitational waves

Gravitational waves are an inevitable product of the Einstein equations, and one can expect a wide spectrum of wave signals propagating throughout our Universe due to anisotropic and inhomogeneous 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.

3.6.1 Planar symmetry

Centrella and Matzner [57, 58] studied a class of plane symmetric cosmologies representing gravitational inhomogeneities in the form of shocks or discontinuities separating two vacuum expanding Kasner cosmologies (2). By a suitable choice of parameters, the constraint equations can be satisfied at the initial time with a 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 [59, 60] developed a polarized plane symmetric code for cosmology, adding also hydrodynamic sources with artificial viscosity methods for shock capturing and Barton’s method for monotonic transport [162]. 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. [9, 11, 7], implementing more robust numerical methods, an improved parametric treatment of the initial value problem, and generic unpolarized metrics.

In applications of these codes, Centrella [61] investigated nonlinear gravitational 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 used to model the transition into the nonlinear regime. Anninos et al. [10] considered small and large perturbations in the two degenerate Kasner models: p₁ = p₃ = 0 or 2/3, and p₂ = 1 or −1/3, respectively, where p2 are parameters in the Kasner metric (2). 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.

Expanding these investigations of the Coulomb nonlinearity, Anninos and McKinney [16] used a gauge invariant perturbation formalism to construct constrained initial data for general relativistic cosmological sheets formed from the gravitational collapse of an ideal gas in a critically closed FLRW “background” model. They compared results to the Newtonian Zel’dovich [165] solution over a broad range of field strengths and flows, and showed that the enhanced growth rates of nonlinear modes (in both the gas density and Riemann curvature invariants) accelerate the collapse process significantly compared to Newtonian and perturbation theory. They also computed the back-reaction of these structures to the mean cosmological expansion rate and found only a small effect, even for cases with long wavelengths and large amplitudes. These structures were determined to produce time-dependent gravitational potential signatures in the CMBR (essentially fully relativistic Rees-Sciama effects) comparable to, but still dominated by, the large scale Sachs-Wolfe anisotropies. This confirmed, and is consistent with, the assumptions built into Newtonian calculations of this effect.

3.6.2 Multi-dimensional vacuum cosmologies

Two additional examples of general relativistic codes developed for the purpose of investigating dynamical behaviors in non-flat, vacuum, cosmological topologies are attributed to Holcomb [91] and Ove [129]. Holcomb considered vacuum axisymmetric models to study the structure of General Relativity and the properties of gravitational waves in non-asymptotically flat spacetimes. The code was based on the ADM 3 + 1 formalism and used Kasner matching conditions at the outer edges of the mesh, mean curvature slicing, and a shift vector to enforce the isothermal gauge in order to simplify the metric and to put it in a form that resembles quasi-isotropic coordinates. However, a numerical instability was observed in cases where the mesh domain exceeded the horizon size. This was attributed to the particular gauge chosen, which does not appear well-suited to the Kasner metric as it results in super-luminal coordinate velocities beyond the horizon scale.

Ove developed an independent code based on the ADM formalism to study cosmic censorship issues, including the nature of singular behavior allowed by the Einstein equations, the role of symmetry in the creation of singularities, the stability of Cauchy horizons, and whether black holes or a ring singularity can be formed by the collision of strong gravitational waves. Ove adopted periodic boundary conditions with 3-torus topology and a single Killing field, and therefore generalizes to two dimensions the planar codes discussed in the previous section. This code also used a variant of constant mean curvature slicing, was fully constrained at each time cycle, and the shift vector was chosen to put the metric into a (time-dependent) conformally flat form at each spatial hypersurface.

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 (Ω₀), Hubble (H₀ = 100 h km s⁻¹ Mpc⁻¹), and cosmological (A) constants.

Due to the vast body of literature on numerical simulations dealing with the post-recombination epoch, only a very small fraction of published work can be reviewed in this paper. 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.

For a general overview of theoretical and observational issues associated with structure formation, the reader is referred to [132, 131], and to [45] for a broad review of numerical simulations (and methods) of structure formation.

4.1 Cosmic microwave background

The Cosmic Microwave Background Radiation (CMBR) is a direct relic of the early Universe, and currently provides the deepest probe of evolving cosmological structures. Although the CMBR is primarily a uniform black body spectrum throughout all space, fluctuations or anisotropies in the spectrum can be observed at very small levels to correlate with the matter density distribution. Comparisons between observations and simulations generally support the mostly isotropic, standard Big Bang model, and can be used to constrain the various proposed matter evolution scenarios and cosmological parameters. For example, sky survey satellite observations [34, 149] suggest a flat A-dominated Universe with scale-invariant Gaussian fluctuations that is consistent with numerical simulations of large sale structure formation (e.g., galaxy clusters, Lya forest).

As shown in the timeline of Figure 8, CMBR signatures can be generally classified into two main components: primary and secondary anisotropies, separated by a Surface of Last Scattering (SoLS). Both of these components include contributions from two distinctive phases: a surface marking the threshold of decoupling of ions and electrons from hydrogen atoms in primary signals, and a surface of reionization marking the start of multiphase secondary contributions through nonlinear structure evolution, star formation, and radiative feedback from the small scales to the large.

Figure 8

Historical time-line of the cosmic microwave background radiation showing the start of photon/nuclei combination, the surface of last scattering (SoLS), and the epoch of reionization due to early star formation. The times are represented in years (to the right and redshift (to the left. Primary anisotropies are collectively attributed to the early effects at the last scattering surface and the large scale Sachs-Wolfe effect. Secondary anisotropies arise from path integration effects, reionization smearing, and higher order interactions with the evolving nonlinear structures at relatively low redshifts.

4.1.1 Primordial black body effects

The black body spectrum of the isotropic background is essentially due to thermal equilibrium prior to the decoupling of ions and electrons, and few photon-matter interactions after that. At sufficiently high temperatures, prior to the decoupling epoch, matter was completely ionized into free protons, neutrons, and electrons. The CMB photons easily scatter off electrons, and frequent scattering produces a blackbody spectrum of photons through three main processes that occur faster than the Universe expands:

• Compton scattering in which photons transfer their momentum and energy to electrons if they have significant energy in the electron’s rest frame. This is approximated by Thomson scattering if the photon’s energy is much less than the rest mass. Inverse Compton scattering is also possible in which sufficiently energetic (relativistic) electrons can blueshift photons.

• Double Compton scattering can both produce and absorb photons, and thus is able to thermalize photons to a Planck spectrum (unlike Compton scattering which conserves photon number, and, although it preserves a Planck spectrum, relaxes to a Bose-Einstein distribution).

• Bremsstrahlung emission of electromagnetic radiation due to the acceleration of electrons in the vicinity of ions. This also occurs in reverse (free-free absorption) since charged particles can absorb photons. In contrast to Coulomb scattering, which maintains thermal equilibrium among baryons without affecting photons, Bremsstrahlung tends to relax photons to a Planck distribution.

Although the CMBR is a unique and deep probe of both the thermal history of the early Universe and primordial perturbations in the matter distribution, the associated anisotropies are not exclusively primordial in nature. Important modifications to the CMBR spectrum, from both primary and secondary components, can arise from large scale coherent structures, even well after the photons decouple from the matter at redshift z ∼ 10³, due to gravitational redshifting, lensing, and scattering effects.

4.1.2 Primary anisotropies

The most important contributions to primary anisotropies between the start of decoupling and the surface of last scattering include the following effects:

• Sachs-Wolfe (SW) effect: Gravitational redshift of photons between potentials at the SoLS and the present. It is the dominant effect at large angular scales comparable to the horizon size at decoupling (θ ∼ 2° Ω¹/²).

• Doppler effect: Dipolar patterns are imprinted in the energy distribution from the peculiar velocities of the matter structures acting as the last scatterers of the photons.

• Acoustic peaks: Anisotropies at intermediate angular scales (0.1° < 0 < 2°) are atttributed to small scale processes occurring until decoupling, including acoustic oscillations of the baryon-photon fluid from primordial density inhomogeneities. This gives rise to acoustic peaks in the thermal spectrum representing the sound horizon scale at decoupling.

• SoLS damping: The electron capture rate is only slightly faster than photodissociation at the start of decoupling, causing the SoLS to have a finite thickness (Δz ∼ 100). Scattering over this interval damps fluctuations on scales smaller than the SoLS depth (θ < 10′ Ω^1/2).

• Silk damping: Photons can diffuse out of overdense regions, dragging baryons with them over small angular scales. This tends to suppress both density and radiation fluctuations.

All of these mechanisms perturb the black body background radiation since thermalization processes are not efficient at redshifts smaller than ∼ 10⁷.

4.1.3 Secondary anisotropies

Secondary anisotropies consist of two principal effects, gravitational and scattering. Some of the more important gravitational contributions to the CMB include:

• Early ISW effect: Photon contributions to the energy density of the Universe may be nonnegligible compared to ordinary matter (dark or baryonic) at the last scattering. The decreasing contribution of photons in time results in a decay of the potential, producing the early Integrated Sachs-Wolfe (ISW) effect.

• Late ISW effect: In open cosmological models or models with a cosmological constant, the gravitational potential decays at late times due to a greater rate of expansion compared to flat spacetimes, producing the late ISW effect on large angular scales.

• Rees-Sciama effect: Evolving nonlinear strucutures (e.g., galaxies and clusters) generate time-varying potentials which can seed asymmetric energy shifts in photons crossing potential wells from the SoLS to the present.

• Lensing: In contrast to ISW effects which change the energy but not directions of the photons, gravitational lensing deflects the paths without changing the energy. This effectively smears out the imaging of the SoLS.

• Proper motion: Compact objects such as galaxy clusters can imprint a dipolar pattern in the CMB as they move across the sky.

• Gravitational waves: Perturbations in the spacetime fabric affect photon paths, energies, and polarizations, predominantly at scales larger than the horizon at decoupling.

Secondary scattering effects are associated with reionization and their significance depends on when and over what scales it takes place. Early reionization leads to large optical depths and greater damping due to secondary scattering. Over large scales, reionization has little effect since these scales are not in causal contact. At small scales, primordial anisotropies can be wiped out entirely and replaced by secondary ones. Some of the more important secondary scattering effects include:

• Thomson scattering: Photons are scattered by free electrons at sufficiently large optical depths achieved when the Universe undergoes a global reionization at late times. This damps out fluctuations since energies are averaged over different directions in space.

• Vishniac effect: In a reionized Universe, high order coupling between the bulk flow of electrons and their density fluctuations generates new anisotropies at small angles.

• Thermal Sunyaev-Zel’dovich effect: Inverse Compton scattering of the CMB by hot electrons in the intracluster gas of a cluster of galaxies distorts the black body spectrum of the CMB. Low frequency photons will be shifted to high frequencies.

• Kinetic Sunyaev-Zel’dovich effect: The peculiar velocities of clusters produces anisotropies via a Doppler effect to shift the temperature without distorting the spectral form. Its effect is proportional to the product of velocity and optical depth.

• Polarization: Scattering of anisotropic radiation affects polarization due to the angular dependence of scattering. Polarization in turn affects anisotropies through a similar dependency and tends to damp anisotropies.

To make meaningful comparisons between numerical models and observed data, all of these (low and high order) effects from both the primary and secondary contributions (see for example Section 4.1.4 and [94, 101]) must be incorporated self-consistently into any numerical model, and to high accuracy in order to resolve and distinguish amongst the various weak signals. The following sections describe some work focused on incorporating many of these effects into a variety of large-scale numerical cosmological models.

4.1.4 Computing CMBR anisotropies with ray-tracing methods

Many efforts based on linear perturbation theory have been carried out to estimate temperature anisotropies in our Universe (for example see [114] and references cited in [131, 94]). 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. Also, because photons are intrinsically coupled to the baryon and dark matter thermal and gravitational states at all spatial scales, a fully self-consistent treatment is needed to accurately resolve the more subtle features of the CMBR. This can be achieved with a ray-tracing approach based on Monte-Carlo methods to track individual photons and their interactions through the evolving matter distributions. A fairly complete simulation involves solving the geodesic equations of motion for the collisionless dark matter which dominate potential interactions, the hydrodynamic equations for baryonic matter with high Mach number shock capturing capability, the transport equations for photon trajectories, a reionization model to reheat the Universe at late times, the chemical kinetics equations for the ion and electron concentrations of the dominant hydrogen and helium gases, and the photon-matter interaction terms describing scattering, redshifting, depletion, lensing, and Doppler effects.

Such an approach has been developed by Anninos et al. [15], and applied to a Hot Dark Matter (HDM) model of structure formation. In order to match both the observed galaxy-galaxy correlation function and LOBE measurements of the CMBR, they find, for that model and neglecting reionization, the cosmological parameters are severely constrained to Ω₂h² ≈ 1, where Ω₀ and h are the density and Hubble parameters respectively.

In models where the IGM does not reionize, the probability of scattering after the photonmatter 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. [157] also using ray-tracing methods. 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 Ω₀ = 1 and h = 0.55, which produces CMBR anisotropies above current LOBE limits when no reionization is included (see Section 4.1.4), they find that the EG scenario effectively reduces the anisotropies to the levels observed by LOBE and generates smaller Doppler shift anisotropies than the LS model, as demonstrated in Figure 9. The LS scenario of reionization is not able to reduce the anisotropy levels below the LOBE limits, and can even give rise to greater Doppler shifts than expected at decoupling.

Figure 9

Temperature fluctuations (ΔT/T) in the CMBR due to the primary Sachs-Wolfe (SW) effect and secondary integrated SW, Doppler, and Thomson scattering effects in a critically closed model. The top two plates are results with no reionization and baryon fractions 0.02 (plate 1, 4° × 4°, ΔT/T |_rms = 2.8 × 10⁻⁵), and 0.2 (plate 2, 8° x 8°, ΔT/T |_rms = 3.4 x 10⁻⁵). The bottom two plates are results from an “early and gradual” reionization scenario of decaying neutrinos with baryon fraction 0.02 (plate 3, 4° × 4°, ΔT/T|rms = 1.3 × 10⁻⁵; and plate 4, 8° × 8°, AT/T |_rms = 1.4 × 10⁻⁵). If reionization occurs, the scattering probability increases and anisotropies are damped with each scattering event. At the same time, matter structures develop large bulk motions relative to the comoving background and induce Doppler shifts on the CMB. The imprint of this effect from last scattering can be a significant fraction of primary 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. [156] considered the impact of nonlinear matter condensations on the CMBR in Ω₀ ≤ 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 (the Rees-Sciama effect), (ii) proper motion of nonlinear structures such as clusters and superclusters across the sky, and (iii) the decaying gravitational potential effect from the evolution of perturbations in open models. They applied the ray-tracing procedure of [15] to explore the relative importance of these secondary anisotropies as a function of the density parameter S2o and the scale of matter distributions. They find that 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|_rms ∼ 10⁻⁶ as shown in Figure 10.

Figure 10

Secondary anisotropies from the proper motion of galaxy clusters across the sky and Rees-Sciama effects are presented in the upper-left image over 8° × 8° in a critically closed Cold Dark Matter model. The corresponding column density of matter over the same region (z = 0.43, Δz = 0.025) is displayed in the upper-right, clearly showing the dipolar nature of the proper motion effect. Anisotropies arising from decaying potentials in an open Ω = 0.3 model over a scale of 8° × 8° are shown in the bottom left image, along with the gravitational potential over the same region (z = 0.33, Δz = 0.03) in the bottom right, demonstrating a clear anti-correlation. Maximum temperature fluctuations in each simulation are ΔT/T = (5 × 10⁻⁷, 1.0 × 10⁻⁶) respectvely. Secondary anisotropies are dominated by decaying potentials at large scales, but all three sources (decaying potential, proper motion, and R-S) produce signatures of order 10⁻⁶.

In addition to the effects discussed in this section, many other sources of secondary anisotropies (as mentioned in Section 4.1, including gravitational lensing, the Vishniac effect accounting for matter velocities and flows into local potential wells, and the Sunyaev-Zel’dovich (SZ) (Section 4.5.4) distortions from the Compton scattering of CMB photons by electrons in the hot cluster medium) can also be fairly significant. See [94, 152, 28, 80, 93] for more thorough discussions of the different sources of CMBR anisotropies.

4.2 Gravitational lensing

Observations of gravitational lenses [143] 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, because 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. [55] developed an 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 Ω₀ = 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 (Ω₀ = 0.3, Λ/3H ²₀ = 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 more r_core ∼ 20–30h⁻¹ kpc. Core radii of this size can have an important effect on the probability of multiple imaging. Flores and Primack [74] considered the effects of assuming two different kinds of splitting sources: isothermal spheres with small but finite core radii and radial density profiles ρ ∝ (r² + r ²_core )⁻¹, and spheres with a Hernquist density profile ρ ∝ r⁻¹ (r + a)⁻³, where r_core ∼ 20–30 h⁻¹ kpc and a ∼ 300 h⁻¹ 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 those produced by strong shock waves as catalysts. The first structures to collapse may be capable of producing pop III stars and have a substantial influence on the subsequent thermal evolution of the intergalactic medium, as suggested by Figure 2, 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 and how the first stars formed.

Ostriker and Gnedin [127] 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 H₂ 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 10⁴ K and hydrogen line cooling becomes efficient. Star formation slows in the intermittent epoch due to the depletion of H₂ by photo-destruction and reheating. In addition, the objects which formed pop III stars also initiate pop II sequences when their virial temperatures reach 10⁴ K through continued mass accretion.

In resolving the details of a single star forming region in a CDM Universe, Abel et al. [2, 3] implemented a non-equilibrium radiative cooling and chemistry model [1, 21] together with the hydrodynamics and dark matter equations, evolving nine separate atomic and molecular species (H, H⁺, He, He⁺, He⁺⁺, H^-, H ^-₂ , H₂, and e^-, according to the reactive network described in Section 6.4.1) on nested and adaptively refined numerical grids. They follow the collapse and fragmentation of primordial clouds over many decades in mass and spatial dynamical range, finding a core of mass ∼ 200 M_⊙ forms from a halo of about ∼ 10⁵ M_⊙ (where a significant number fraction of hydrogen molecules are created) after less than one percent of the halo gas cools by molecular line emission. Bromm et al. [48] use a different Smoothed Particle Hydrodynamics (SPH) technique and a six species model (H, H⁺, H^-, H ⁺₂ , H₂, and e^-) to investigate the initial mass function of the first generation pop III stars. They evolve an isolated 3σ peak of mass 2 × 10⁶M_⊙ which collapses at redshift z ∼ 30 and forms clumps of mass 10² – 10³ M_⊙ which then grow by accretion and merging, suggesting that the very first stars were massive and in agreement with [3].

The implications of an early era of massive star populations on the thermal and chemical state of the intergalactic medium was investigated by Yoshida et al. [164]. They considered the effects of feedback and radiation transfer in early structure formation simulations to show that a significant fraction of the IGM can be ionized and polluted by metals from the first stars to form and become supernovae by z ∼ 15, thus affecting subsequent stellar populations. They also argue that observed elemental abundances in the intracluster medium are not affected by metals originating from the first stars.

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

Many complex multi-component numerical simulations have been performed of the Lyman forest, which include the effects of dark matter (N-body), baryons (hydrodynamics), chemical composition (reactive networks), and microphysical response (radiative cooling and heating). See, for example, [67, 118, 166], which represent some of the earliest comprehensive simulations. For the most part, all these calculations have been able to fit the observations reasonably well, including the column density and Doppler width distributions, the size of absorbers [62], and the line Number evolution. Despite the fact that the cosmological models and parameters are different in each case, the simulations give roughly similar results provided that the proper ionization bias is used, b_ion ≡ (Ω_bh²)²/Γ, where Ω_b is the baryonic density parameter, h is the Hubble parameter and Γ is the photoionization rate at the hydrogen Lyman edge. (However, see [50] for a discussion of the sensitivity of statistical properties on numerical resolution.) A theoretical paradigm has thus emerged from these calculations in which Lya absorption lines originate from the relatively smaller 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 11), including fluctuations in underdense regions, spheroidal minihalos, and filaments extending over scales of a few Mpc.

Machacek et al. [115] expanded on earlier work to compare several Lya statistical measures from five different background cosmological models, including standard critical density Cold Dark Matter (CDM), open CDM, flat CDM with a cosmological constant, standard CDM with a tilted density spectrum, and a flat model with mixed hot and cold dark matter. All models were chosen to match local or low redshift observations, and most were also consistent with LOBE measurements of the CMBR. The calculations were designed to establish which statistics are sensitive to different cosmological models. In particular, they find that the line number count above a given column density threshold is relatively insensitive to background models. On the other hand, the shape of the optical depth probability distribution function is strongly correlated to the amount of small scale power in density fluctuations, and is thus a good discriminator among models on scales of a few hundred kpc.

Meiksin et al. [117] followed up with more detailed comparisons of Lyα systems in several cosmologies with observed high resolution QSO spectra. Although all models are consistent with previous studies in that they give reasonably good statistical agreement with observed Lya properties, under closer scrutiny none of the numerical models they considered passed all the tests, which included spectral flux, wavelet decomposed amplitude, and absorption line profile distributions. They suggest that comparisons might be improved, particularly in optically thin systems, by more energy injection into the IGM from late He⁺ reionization or supernovae-driven winds, or by a larger baryon fraction.

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 10¹⁴–10¹⁵ M_⊙, and the environment in these structures is composed of shock heated gas with temperatures of order 10⁷mbox — −10⁸ K which emits thermal bremsstrahlung and line radiation at X-ray energies. Also, because of their spatial size of ∼ 1 h^-1 Mpc and separations of order 50 h^-1 Mpc, they provide a measure of nonlinearity on scales close to the perturbation normalization scale 8 h^-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 observable redshift range 0 ≤ z ≤ 1.

4.5.1 Internal structure

Thomas et al. [155] 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 generally follow the universal density profile of Navarro et al. [125]. 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 [65]. 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 could distinguish between open and critically closed Universes.

Figure 11

Distribution of the gas density at redshift z = 3 from a numerical hydrodynamics simulation of the Lya forest with a CDM spectrum normalized to second year COBS observations, 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 × 10⁴ K, light blue = 3 × 10⁵ 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.2 Number density evolution

The evolution of the number density of rich clusters of galaxies can be used to compute Ω₀ and σ₈ (the power spectrum normalization on scales of 8 h^-1 Mpc) when numerical simulation results are combined with the constraint σ₈Ω ^0.5₀ ≈ 0.5, derived from observed present-day abundances of rich clusters. Bahcall et al. [24] computed the evolution of the cluster mass function in five different cosmological model simulations and find that the number of high mass (Coma-like) clusters in flat, low as models (i.e., the standard CDM model with σ₈ ≈ 0.5) decreases dramatically by a factor of approximately 10³ from z = 0 to z ≈ 0.5. For low Ω₀, high σ₈ models, the data result 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.3 ± 0.1 and σ₈ = 0.85 ± 0.5), and flat low density models with a cosmological constant (Ω₀ = 0.34 ± 0.13 and Ω₀ + Λ = 1).

4.5.3 X-ray luminosity function

The evolution of the X-ray luminosity function, as well as the number, size and temperature distribution of galaxy clusters are all potentially important discriminants of cosmological models and the underlying initial density power spectrum that gives rise to these structures. Because the X-ray luminosity (principally due to thermal bremsstrahlung emission from electron/ion interactions in the hot fully ionized cluster medium) is proportional to the square of the gas density, the contrast between cluster and background intensities is large enough to provide a window of observations that is especially sensitive to cluster structure. Comparisons of simulated and observed X-ray functions may be used to deduce the amplitude and shape of the fluctuation spectrum, the mean density of the Universe, the mass fraction of baryons, the structure formation model, and the background cosmological model.

Several groups [49, 56] have examined the properties of X-ray clusters in high resolution numerical simulations of a standard CDM model normalized to LOBE. Although the results are very sensitive to grid resolution (see [17] 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. On the other hand, similar calculations with different cosmological models [56, 52] suggest reasonable agreement of observed data with Cold Dark Matter + cosmological constant (ACDM), Cold + Hot Dark Matter (CHDM), and Open or low density CDM (OCDM) evolutions due to different universal expansions and density power spectra.

4.5.4 SZ effect

The Sunyaev-Zel’dovich (SZ) effect is the change in energy that CMB photons undergo when they scatter in hot gas typically found in cores of galaxy clusters. There are two main effects: thermal and kinetic. Thermal SZ is the dominant mechanism which arises from thermal motion of gas in the temperature range 10⁷–10⁸ K, and is described by the Compton y parameter

$$y = {\sigma _T}\int {\frac{{{n_e}{k_B}{T_e}}}{{{m_e}{c^2}}}} dl,$$

((7))

where σ_T = 6.65 × 10^-25 cm² is the Thomson cross-section, m_e, n_e and T_e are the electron rest mass, density and temperature, c is the speed of light, k_B is Boltzmann’s constant, and the integration is performed over the photon path. Photon temperature anisotropies are related to the y parameter by ΔT/T ≈ −2y in the Rayleigh-Jeans limit. The kinetic SZ effect is a less influential Doppler shift resulting from the bulk motion of ionized gas relative to the rest frame of the CMB.

Springel et al. [150] used a Tree/SPH code to study the SZ effects in a CDM cosmology with a cosmological constant. They find a mean amplitude for thermal SZ (y = 3.8 × 10^-6) just below current observed upper limits, and a kinetic SZ about 30 times smaller in power. Da Silva et al. [66] compared thermal SZ maps in three different cosmologies (low density + Λ, critical density, and low density open model). Their results are also below current limits: y ≈ 4 × 10^-6 for low density models with contributions from over a broad redshift range z ≤ 5, and y ≈ 1 × 10^-6 for the critical density model with contributions mostly from z < 1. However, further simulations are needed to explore the dependence of the SZ effect on microphysics, i.e., cooling, star formation, supernovae feedback.

4.6 Cosmological sheets

Cosmological sheets, or pancakes, form as overdense regions collapse preferentially along one axis. Originally studied by Zel’dovich [165] 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 [47, 145, 22] have been performed of these systems which include baryonic fluid with hydrodynamical shock heating, ionization, recombination, dark matter, thermal conductivity, and radiative cooling (Compton, bremsstrahlung, and atomic line cooling), in both one and two spatial dimensions to assert the significance of each physical process and to compute the fragmentation scale. See also [16] where fully general relativistic numerical calculations of cosmological sheets are presented in plane symmetry, including relativistic hydrodynamical shock heating and consistent coupling to spacetime curvature.

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 H₂. 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 H₂ cooled cosmological sheets may be associated with massive stars or star clusters. Anninos and Norman [18] 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 n_H2/n_H ∼ 3 × 10^-3, which cools the gas to 4 × 10^-3 eV and results in a fragmentation scale of 9 × 10³ M_⊙. This represents reductions of 50 and 10³ in temperature and Jeans mass respectively when compared, as in Figure 12, to the equivalent case in which hydrogen molecules were neglected.

Figure 12

Two different model simulations of cosmological sheets are presented: a six species model including only atomic line cooling (left, and a nine species model including also hydrogen molecules (right. 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 10⁷ (10⁵) solar masses.

However, the above calculations neglected important interactions arising from self-consistent treatments of radiation fields with ionizing and photo-dissociating photons and self-shielding effects. Susa and Umemura [153] studied the thermal history and hydrodynamical collapse of pancakes in a UV background radiation field. They solve the radiative transfer of photons together with the hydrodynamics and chemistry of atomic and molecular hydrogen species. Although their simulations were restricted to one-dimensional plane parallel symmetry, they suggest a classification scheme distinguishing different dynamical behavior and galaxy formation scenarios based on the UV background radiation level and a critical mass corresponding to 1 – 2σ density fluctuations in a standard CDM cosmology. These level parameters distinguish galaxy formation scenarios as they determine the local thermodynamics, the rate of H₂ line emissions and cooling, the amount of starburst activity, and the rate and mechanism of cloud collapse.

5 Summary

This review is intended to provide a flavor of the variety of numerical cosmological calculations performed of different events occurring 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. For the most part, the nature of the calculations dealing with the early or late Universe can be distinguished by their basic motivations. For example, calculations of early Universe phenomena are designed to explore alternative cosmological models or topologies, and in some cases, different theories of gravity. They also tend to study the nature of topological singularities, geometric effects, and the problem of initial conditions or the origin of matter distributions. Calculations of the late Universe are generally focused to establish bounds on cosmological parameters in the context of the standard model, to resolve the correct structure formation scenario, to model the complex multi-physics interactions operating at vastly different scales, and to systematically compare invariant measures against observed data for both model validation and interpreting observations.

Although a complete, self-consistent, and accurate description of our Universe is impractical considering the complex multi-scale and multi-physics requirements, a number of enlightening results have been demonstrated through computations. For example, both monotonic AVTD and chaotic oscillatory BLK behavior have been found in the asymptotic approach to the initial singularity in a number of inhomogeneous cosmological models, though some issues remain concerning the generic nature of the singularity, including the effect of nonlinear mode coupling of spatial gradients to the oscillatory history, and the behavior in non-vacuum spacetimes with arbitrary global topology. Numerical calculations suggest that scalar fields play an important complicated role in the nonlinear or chaotic evolution of cosmological models with consequences for the triggering (or not) of inflation and the subsequent dynamics of structure formation. It is possible, for example, that these fields can influence the details of inflation and have observable ramifications as fractal patterns in the density spectrum, gravitational waves, galaxy distribution, and cosmic microwave background anisotropies. All these effects require further studies. Numerical simulations have also been used to place limits on curvature anisotropies and cosmological parameters at early times by considering primordial nucleosynthesis reactions in anisotropic and inhomogeneous cosmologies.

Finally, the large collection of calculations performed of the post-recombination epoch related to large scale structure formation (for example, cosmic microwave background, gravitational lensing, Lyα absorption, and galaxy cluster simulations) have placed strong constraints on the standard model parameters and structure formation scenarios when compared to observations. Considering the range of models consistent with inflation, the preponderance of observational, theoretical and computational data suggest a best fit model of the late structure-forming Universe that is spatially flat with a cosmological constant and a small tilt in the power spectrum. These best fit model parameters, and in particular the introduction of a cosmological constant, are generally consistent with recent evidence of dark energy from supernovae and high precision CMBR observations.

Obviously many fundamental issues remain unresolved, including even the overall shape or topology of the cosmological model which best describes our Universe throughout its entire history. However, the field of numerical cosmology has matured in the development of computational techniques, the modeling of microphysics, and in taking advantage of current trends in computing technologies, to the point that it is now possible to perform high resolution multiphysics simulations and carry out reliable comparisons of numerical models with observational data.

6 Appendix: Basic Equations and Numerical Methods

Some basic equations relevant for fully relativistic as well as perturbative cosmological calculations are summarized in this section, including the complete Einstein equations, choices of kinematical conditions, initial data constraints, stress-energy-momentum tensors, dynamical equations for various matter sources, and the Newtonian counterparts on background isotropic models. References to numerical methods are also provided.

6.1 The Einstein equations

6.1.1 ADM formalism

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

$$d{s^2} = ( - {\alpha ^2} + {\beta _i}{\beta ^i})d{t^2} + 2{\beta _i}d{x^i}dt + {\gamma _{ij}}d{x^i}d{x^j},$$

((8))

where α and βⁱ are the lapse function and shift vector respectively, and γ_ij is the spatial 3-metric. The lapse defines the proper time between consecutive layers of spatial hypersurfaces, the shift propagates the coordinate system from 3-surface to 3-surface, and the induced 3-metric is related to the 4-metric via γ_μv = g_μv + n_μn_v. The Einstein equations are written as four constraint equations,

$$^{(3)}R + {K^2} - {K_{ij}}{K^{ij}} = 16\pi G{\rho _H},$$

((9))

$${\nabla _i}({K^{ij}} - {\gamma ^{ij}}K) = 8\pi G{s^j},$$

((10))

twelve evolution equations,

$$\begin{array}{*{20}{c}} {{\partial _t}{\gamma _{ij}} = {\mathcal{L}_\beta }{\gamma _{ij}} - 2\alpha {K_{ij}},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;} \\ {{\partial _t}{K_{ij}} = {\mathcal{L}_\beta }{K_{ij}} - {\nabla _i}{\nabla _j}\alpha + \alpha \left[ {{}^{(3)}{R_{ij}} - 2{K_{ik}}K_j^k + K\;{K_{ij}} - 8\pi G\left( {{S_{ij}} - \frac{1}{2}S{\gamma _{ij}} + \frac{1}{2}{\rho _H}{\gamma _{ij}}} \right)} \right]} \end{array},$$

((11))

and four kinematical or coordinate conditions for the lapse function and shift vector that can be specified arbitrarily (see Section 6.1.2). Here K_ij is the extrinsic curvature describing how the 3-metric evolves along a time-like normal vector. It is formally defined as the Lie derivative (${\mathcal{L}_n}$) of the spatial metric along the vector ${n^\mu },{K_{ij}} \equiv {\mathcal{L}_n}{\gamma _{ij}}/2$. Also,

$$\begin{array}{*{20}{c}} {{\mathcal{L}_\beta }{\gamma _{ij}} = {\nabla _i}{\beta _j} + {\nabla _j}{\beta _i},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;} \\ {{\mathcal{L}_\beta }{K_{ij}} = {\beta ^k}{\nabla _k}{K_{ij}} + {K_{ik}}{\nabla _j}{\beta ^k} + {K_{kj}}{\nabla _i}{\beta ^k},} \end{array}$$

((12))

where ∇_i is the spatial covariant derivative with respect to γ_ij, ⁽³⁾R_ij is the spatial Ricci tensor, K is the trace of the extrinsic curvature K_ij, and G is the gravitational constant. The matter source terms ρ_H, s^j, s_ij and s = s ⁱ _i as seen by the observers at rest in the time slices are obtained from the appropriate projections

$${\rho _H} = {n^\mu }{n^\nu }{T_{\mu \nu }},$$

((13))

$${s_i} = - \gamma _i^\mu {n^\nu }{T_{\mu \nu }},$$

((14))

$${s_{ij}} = - \gamma _i^\mu \gamma _j^\nu {T_{\mu \nu }},$$

((15))

for the energy density, momentum density and spatial stresses, respectively. Here c = 1, and Greek (Latin) indices refer to 4(3)-dimensional quantities.

It is worth noting that several alternative formulations of Einstein’s equations have been suggested, including hyperbolic systems [136] which have nice mathematical properties, and conformal traceless systems [147, 31] which make use of a conformal decomposition of the 3-metric and trace free part of the extrinsic curvature A_ij = K_ij - γ_ij K/3. Introducing ${\tilde \gamma _{ij}} = {e^{ - 4\psi }}{\gamma _{ij}}$ with e^4ψ γ^1/3 so that the determinant of ${\tilde \gamma _{ij}}$ is unity, and ${\tilde A_{ij}} = {e^{ - 4\psi }}{A_{ij}}$, evolution equations can be written in the conformal traceless system for $\psi ,{\tilde \gamma _j},K,{\tilde A_{ij}}$ and the conformal connection functions, though not all of these variables are independent. However, it is not yet entirely clear which of these methods is best suited for generic problems. For example, hyperbolic forms are easier to characterize mathematically than ADM and may potentially be more stable, but can suffer from greater inaccuracies by introducing additional equations or higher order derivatives. Conformal treatments are considered to be generally more stable [31], but can be less accurate than traditional ADM for short term evolutions [6].

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

6.1.2 Kinematic conditions

For cosmological simulations, one typically takes the shift vector to be zero, hence ${\mathcal{L}_\beta }{\gamma _{ij}} = {\mathcal{L}_\beta }{K_{ij}} = 0$. However, the shift can be used advantageously in deriving conditions to maintain the 3-metric in a particular form, and to simplify the resulting differential equations [59, 60]. See also [146] describing an approximate minimum distortion gauge condition used to help stabilize simulations of general relativistic binary clusters and neutron stars.

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

$$\alpha = {F_1}({x^\mu }){F_2}(\gamma ),$$

((16))

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

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

$${\nabla ^i}{\nabla _i}\alpha = \alpha [{K_{ij}}{K^{ij}} + 4\pi G({\rho _H} + s)] + {\beta ^i}{\nabla _i}K - {\partial _t}K,$$

((17))

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

6.1.3 Symplectic formalism

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

$$d{s^2} = {e^{ - \lambda /2}}{e^{\tau /2}}( - {e^{ - 2\tau }}d{\tau ^2} + d{\theta ^2}) + {e^{ - \tau }}[{e^P}d{\sigma ^2} + 2{e^P}Q\;d\sigma \;d\delta + ({e^P}{Q^2} + {e^{ - P}})d{\delta ^2}],$$

((18))

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

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

$$H = {H_1} + {H_2} = \frac{1}{2}\int_0^{2\pi } {d\theta (\pi _P^2 + {e^{ - 2P}}\pi _Q^2) + \frac{1}{2}\int_0^{2\pi } {d\theta \;{e^{ - 2\tau }}(P_{,\theta }^2 + {e^{2P}}Q_{,\theta }^2)} .} $$

((19))

The symplectic method approximates the evolution operator by

$${e^{H\Delta \tau }} = {e^{{H_2}\Delta \tau /2}}{e^{{H_1}\Delta \tau }}{e^{{H_2}\Delta \tau /2}} + \mathcal{O}{(\Delta \tau )^3},$$

((20))

although higher order representations are possible. If the two Hamiltonian components H₁ and H₂ are each integrable, their solutions can be substituted directly into the numerical evolution to provide potentially more accurate solutions with fewer time steps [36]. This approach is well-suited for studies of singularities if the asymptotic behavior is determined primarily by the kinetic subhamiltonian, a behavior referred to as Asymptotically Velocity Term Dominated (see Section 3.1.2 and [37]).

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

6.1.4 Regge calculus model

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

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

6.2 Sources of matter

6.2.1 Cosmological constant

A cosmological constant is implemented in the 3 + 1 framework simply by introducing the quantity - Λ/(8πG) as an effective isotropic pressure in the stress-energy tensor

$${T_{\mu \nu }} = - \frac{\Lambda }{{8\pi G}}{g_{\mu \nu }}.$$

((21))

The matter source terms can then be written as

$${\rho _H} = \frac{\Lambda }{{8\pi G}},$$

((22))

$${s_{ij}} = - \frac{\Lambda }{{8\pi G}}{\gamma _{ij}},$$

((23))

with sⁱ - 0.

6.2.2 Scalar field

The dynamics of scalar fields is governed by the Lagrangian density

$$\mathcal{L} = - \frac{1}{2}\sqrt { - g} [{g^{\mu \nu }}{\nabla _\mu }\phi \nabla \phi + \xi Rf(\phi ) + 2V(\phi )],$$

((24))

where R is the scalar Riemann curvature, V (φ) is the interaction potential, f(φ) is typically assumed to be f(φ) = φ², and ξ is the field-curvature coupling constant (ξ = 0 for minimally coupled fields and ξ = 1/6 for conformally coupled fields). Varying the action yields the Klein-Gordon equation

$${g^{\mu \nu }}{\nabla _\mu }{\nabla _\nu }\phi - \xi R\phi - {\partial _\phi }V(\phi ) = 0,$$

((25))

for the scalar field and

$$\begin{array}{*{20}{c}} {{T_{\mu \nu }} = (1 - 2\xi ){\nabla _\mu }\phi {\nabla _\nu }\phi + \left( {2\xi - \frac{1}{2}} \right){g_{\mu \nu }}{g^{\sigma \lambda }}{\nabla _\sigma }\phi {\nabla _\lambda }\phi \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;} \\ { - 2\xi \phi {\nabla _\mu }{\nabla _\nu }\phi + 2\xi \phi {g_{\mu \nu }}{g^{\sigma \lambda }}{\nabla _\sigma }{\nabla _\lambda }\phi + \xi {G_{\mu \nu }}{\phi ^2} - {g_{\mu \nu }}V(\phi ),} \end{array}$$

((26))

for the stress-energy tensor, where G_μv - R_gmv - g_μvR/2.

For a massive, minimally coupled scalar field [46],

$${T_{\mu \nu }} = {\partial _\mu }\phi \;{\partial _\nu }\phi - \frac{1}{2}{g_{\mu \nu }}{g^{\rho \sigma }}{\partial _\rho }\phi \;{\partial _\sigma }\phi - {g_{\mu \nu }}V(\phi ),$$

((27))

and

$${\rho _H} = \frac{1}{2}{\gamma ^{ij}}{\partial _i}\phi {\partial _j}\phi + \frac{1}{2}{\eta ^2} + V(\phi ),$$

((28))

$${s_i} = - \eta {\partial _i}\phi ,$$

((29))

$${s_{ij}} = {\gamma _{ij}}\left( { - \frac{1}{2}{\gamma ^{kl}}{\partial _k}\phi {\partial _l}\phi + \frac{1}{2}{\eta ^2} - V(\phi )} \right) + {\partial _i}\phi {\partial _j}\phi ,$$

((30))

where

$$\eta = {n^\mu }{\partial _\mu }\phi = \frac{1}{\alpha }({\partial _t} - {\beta ^k}{\partial _k})\phi ,$$

((31))

n^μ = (1, - βⁱ)/α, and V (φ) is a general potential which, for example, can be set to V = λφ⁴ in the chaotic inflation model. The covariant form of the scalar field equation (25) can be expanded as in [107] to yield

$$\frac{1}{\alpha }({\partial _t} - {\beta ^k}{\partial _k})\eta = \frac{1}{{\sqrt \gamma }}{\partial _i}(\sqrt \gamma {\gamma ^{ij}}{\partial _j}\phi ) + \frac{1}{\alpha }{\gamma ^{ij}}{\partial _i}\alpha {\partial _j}\phi + K\eta - {\partial _\phi }V(\phi ),$$

((32))

which, when coupled to Equation (31), determines the evolution of the scalar field.

6.2.3 Collisionless dust

The stress-energy tensor for a fluid composed of collisionless particles (or dark matter) can be written simply as the sum of the stress-energy tensors for each particle [161],

$${T_{\mu \nu }} = \sum {mn{u_\mu }{u_\nu },} $$

((33))

where m is the rest mass of the particles, n is the number density in the comoving frame, and u^μ is the 4-velocity of each particle. The matter source terms are

$${\rho _H} = \sum {mn{{(\alpha {u^0})}^2},} $$

((34))

$${s_i} = \sum {mn{u_i}(\alpha {u^0}),} $$

((35))

$${s_{ij}} = \sum {mn{u_i}{u_j}.} $$

((36))

There are two conservation laws: the conservation of particles ∇_μ(nu^μ) = 0, and the conservation of energy-momentum ∇_μT^μv = 0, where ∇_μ is the covariant derivative in the full 4-dimensional spacetime. Together these conservation laws lead to u^μ∇_μu^v = 0, the geodesic equations of motion for the particles, which can be written out more explicitly in the computationally convenient form

$$\frac{{d{x^i}}}{{dt}} = \frac{{{g^{i\alpha }}{u_\alpha }}}{{{u^0}}},$$

((37))

$$\frac{{d{u_i}}}{{dt}} = - \frac{{{u_\alpha }{u_\beta }{\partial _i}{g^{\alpha \beta }}}}{{2{u^0}}},$$

((38))

where xⁱ is the coordinate position of each particle, u⁰ is determined by the normalization u^μu_μ = - 1,

$$\frac{d}{{dt}} = {v^\mu }{\partial _\mu } = {\partial _t} + {v^i}{\partial _i}$$

((39))

is the Lagrangian derivative, and v^μ = u^μ/u⁰ is the “transport” velocity of the particles as measured by observers at rest with respect to the coordinate grid.

6.2.4 Ideal gas

The stress-energy tensor for a perfect fluid is

$${T_{\mu \nu }} = \rho h{u_\mu }{u_\nu } + P{g_{\mu \nu }},$$

((40))

where g_μv is the 4-metric, h - 1 + ε + P/ρ is the relativistic enthalpy, and ε, P, ρ and u_μ, are the specific internal energy (per unit mass), pressure, rest mass density and four-velocity of the fluid. Defining Vⁱ - uⁱ/u⁰ and

$$u = - {n_\mu }{u^\mu } = \alpha {u^0} = {(1 + {u^i}{u_i})^{1/2}} = {\left( {1 - \frac{{{V_i}{V^i}}}{{{\alpha ^2}}}} \right)^{ - 1/2}},$$

((41))

as the generalization of the special relativistic boost factor, the matter source terms become

$${\rho _H} = \rho h{u^2} - P,$$

((42))

$${s_i} = \rho hu{u_i},$$

((43))

$${s_{ij}} = P{\gamma _{ij}} + \rho h{u_i}{u_j}.$$

((44))

The hydrodynamics equations are derived from the normalization of the 4-velocity, u^μu_μ = −1, the conservation of baryon number, ∇_μ(ρu^μ) - 0, and the conservation of energy-momentum, ∇_μT^μv - 0. The resulting equations can be written in flux conservative form as [162]

$$\frac{{\partial D}}{{\partial t}} + \frac{{\partial (D{V^i})}}{{\partial {x^i}}} = 0,$$

((45))

$$\frac{{\partial E}}{{\partial t}} + \frac{{\partial (E{V^i})}}{{\partial {x^i}}} + P\frac{{\partial W}}{{\partial t}} + P\frac{{\partial (W{V^i})}}{{\partial {x^i}}} = 0,$$

((46))

$$\frac{{\partial {S_i}}}{{\partial t}} + \frac{{\partial ({S_i}{V^j})}}{{\partial {x^j}}} - \frac{{{S^\mu }{S^\nu }}}{{2{S^0}}}\frac{{\partial {g_{\mu \nu }}}}{{\partial {x^i}}} + \sqrt { - g} \frac{{\partial P}}{{\partial {x^i}}} = 0,$$

((47))

where $W = \sqrt { - g} {u^0},D = W\rho ,E = W\rho\epsilon ,\;{S_i} = W\rho h{u_i}$, and g is the determinant of the 4-metric satisfying $\sqrt { - g} = \alpha \sqrt \gamma$. A prescription for specifying an equation of state (e.g., P - (Γ-1)E/W = (Γ - 1)ρε for an ideal gas) completes the above equations.

When solving Equations (45, 46, 47), an artificial viscosity method is needed to handle the formation and propagation of shock fronts [162, 85, 84]. These methods are computationally cheap, easy to implement, and easily adaptable to multi-physics applications. However, it has been demonstrated that problems involving very high Lorentz factors are somewhat sensitive to different implementations of the viscosity terms, and can result in substantial numerical errors if solved using time explicit methods [126].

On the other hand, a number of different formulations [75] of these equations have been developed to take advantage of the hyperbolic and conservative nature of the equations in using high resolution and non-oscillatory shock capturing schemes (although strict conservation is only possible in flat spacetimes — curved spacetimes exhibit source terms due to geometry). These techniques potentially provide more accurate and stable treatments in highly relativistic regimes. A particular formulation used together with high resolution Godunov techniques and approximate Riemann solvers is the following [139, 26]:

$$\frac{{\partial \sqrt \gamma U(\vec w)}}{{\partial t}} + \frac{{\partial \sqrt { - g} {F^i}(\vec w)}}{{\partial {x^i}}} = \sqrt { - g} S(\vec w),$$

((48))

where

$$S(\vec w) = \left[ {0,{T^{\mu \nu }}\left( {\frac{{\partial {g_{vj}}}}{{\partial {x^\mu }}} - \Gamma _{\nu \mu }^\delta {g_{\delta j}}} \right),\alpha \left( {{T^{\mu 0}}\frac{{\partial \ln \alpha }}{{\partial {x^\mu }}} - {T^{\mu \nu }}\Gamma _{\nu \mu }^0} \right)} \right],$$

((49))

$${F^i}(\vec w) = \left[ {D\left( {{v^i} - \frac{{{\beta ^i}}}{\alpha }} \right),{S_j}\left( {{v^i} - \frac{{{\beta ^i}}}{\alpha }} \right) + P\delta _j^i,(E - D)\left( {{v^i} - \frac{{{\beta ^i}}}{\alpha }} \right) + P{v^i}} \right],$$

((50))

and $\vec w = (\rho ,{\nu _i},\epsilon), U (\vec w) = (D, {S_i}, E - D), E = \rho h{\tilde W^2} - P,{S_j} = \rho h{\tilde W^2}{v_j},D = \rho \tilde W,{v^i} = {\gamma ^{ij}}{v_j} = {u^i}/(\alpha {u^0}) + {\beta ^i}/\alpha$, and $\tilde W = \alpha {u^0} = {(1 - {\gamma _{ij}}{\nu ^i}{\nu ^j})^{ - 1/2}}$.

Although Godunov-type schemes are accepted as more accurate alternatives to AV methods, especially in the limit of high Lorentz factors, they are not immune to problems and should generally be used with caution. They may produce unexpected results in certain cases that can be overcome only with problem-specific fixes or by adding additional dissipation. A few known examples include the admittance of expansion shocks, negative internal energies in kinematically dominated flows, the ‘carbuncle’ effect in high Mach number bow shocks, kinked Mach stems, and odd/even decoupling in mesh-aligned shocks [135]. Godunov methods, whether they solve the Riemann problem exactly or approximately, are also computationally much more expensive than their simpler AV counterparts, and it is more difficult to incorporate additional physics.

A third class of computational fluid dynamics methods reviewed here is also based on a conservative hyperbolic formulation of the hydrodynamics equations. However, in this case the equations are derived directly from the conservation of stress-energy,

$${\nabla _\mu }{T^{\mu \nu }} = \frac{1}{{\sqrt { - g} }}{\partial _\mu }\left( {\sqrt { - g} {T^{\mu \nu }}} \right) + \Gamma _{\alpha \mu }^\nu {T^{\mu \alpha }} = 0,$$

((51))

to give

$$\frac{{\partial D}}{{\partial t}} + \frac{{\partial (D{V^i})}}{{\partial {x^i}}} = 0,$$

((52))

$$\frac{{\partial \varepsilon }}{{\partial t}} + \frac{{\partial (\varepsilon {V^i})}}{{\partial {x^i}}} + \frac{{\partial [\sqrt { - g} ({g^{0i}} - {g^{00}}{V^i})P]}}{{\partial {x^i}}} = {\Sigma ^0},$$

((53))

$$\frac{{\partial {S^j}}}{{\partial t}} + \frac{{\partial ({S^j}{V^i})}}{{\partial {x^i}}} + \frac{{\partial [\sqrt { - g} ({g^{ij}} - {g^{0j}}{V^i})P]}}{{\partial {x^i}}} = {\Sigma ^j},$$

((54))

with curvature source terms ${\Sigma ^\nu } = - \sqrt { - g} {T^{\beta \gamma }}\;\Gamma _{\beta \gamma }^\nu$. The variables D and Vⁱ are the same as those defined in the internal energy formulation, but now

$$\varepsilon = W\rho h{u^0} + \sqrt { - g} {g^{00}}P,$$

((55))

$${\mathcal{S}^i} = W\rho h{u^i} + \sqrt { - g} {g^{0i}}P$$

((56))

are different expressions for energy and momenta. An alternative approach of using high resolution, non-oscillatory, central difference (NOCD) methods [99, 100] has been applied by Anninos and Fragile [12] to solve the relativistic hydrodynamics equations in the above form. These schemes combine the speed, efficiency, and flexibility of AV methods with the advantages of the potentially more accurate conservative formulation approach of Godunov methods, but without the cost and complication of Riemann solvers and flux splitting.

NOCD and artificial viscosity methods have been discussed at length in [12] and compared also with other published Godunov methods on their abilities to model shock tube, wall shock and black hole accretion problems. They find that for shock tube problems at moderate to high boost factors, with velocities up to V ∼ 0.99, internal energy formulations using artificial viscosity methods compare quite favorably with total energy schemes, including NOCD methods and Godunov methods using either approximate or exact Riemann solvers. However, AV methods can be somewhat sensitive to parameters (e.g., viscosity coefficients, Courant factor, etc.) and generally suspect in wall shock problems at high boost factors (V > 0.95). On the other hand, NOCD methods can easily be extended to ultra-relativistic velocities (1-V < 10–11) for the same wall shock tests, and are comparable in accuracy to the more standard but complicated Riemann solver codes. NOCD schemes thus provide a reasonable alternative for relativistic hydrodynamics, though it should be noted that low order versions of these methods can be significantly more diffusive than either the AV or Godunov methods.

6.2.5 Imperfect fluid

The perfect fluid equations discussed in Section 6.2.4 can be generalized to include viscous stress in the stress-energy tensor,

$${T^{\mu \nu }} = \left[ {T_{ideal}^{\mu \nu }} \right] + \left[ {T_{viscous}^{\mu \nu }} \right]$$

((57))

$$ = \left[ {\rho h{u^\mu }{u^\nu } + P{g^{\mu \nu }}} \right] - \left[ {2\eta {\sigma ^{\mu \nu }} + \xi \theta ({g^{\mu \nu }} + {u^\mu }{u^\nu })} \right],$$

((58))

where η ≥ 0 and ξ ≥ 0 are the dynamic shear and bulk viscosity coefficients, respectively. Also, θ = ∇_μu^μ is the expansion of fluid world lines, σ^μv is the trace-free spatial shear tensor with

$${\sigma ^{\mu \nu }} = \frac{1}{2}\left( {{h^{\lambda \nu }}{\nabla _\lambda }{u^\mu } + {h^{\lambda \mu }}{\nabla _\lambda }{u^\nu }} \right) - \frac{1}{3}\theta {h^{\mu \nu }},$$

((59))

and h^μv = g^μv + u^μu^v is the projection tensor.

The corresponding energy and momentum conservation equations for the internal energy formulation of Section 6.2.4 become

$$\frac{{\partial E}}{{\partial t}} + P\frac{{\partial W}}{{\partial t}} + \frac{{\partial (E{v^i})}}{{\partial {x^i}}} + P\frac{{\partial (W{v^i})}}{{\partial {x^i}}} = \sqrt { - g} {u_\nu }{\nabla _\mu }T_{viscous}^{\mu \nu },$$

((60))

$$\frac{{\partial {S_j}}}{{\partial t}} + \frac{{\partial ({S_j}{v^i})}}{{\partial {x^i}}} + \sqrt { - g} \frac{{\partial P}}{{\partial {x^j}}} - \frac{{{S^\mu }{S^\nu }}}{{2{S^0}}} - \frac{{\partial {g_{\mu \nu }}}}{{\partial {x^j}}} = - \sqrt { - g} {g_{jv}}{\nabla _\mu }T_{viscous}^{\mu \nu }.$$

((61))

For the NOCD formulation discussed in Section 6.2.4 it is sufficient to replace the source terms in the energy and momentum equations (53, 53, 54) by

$${\Sigma ^\alpha } = - \sqrt { - g} {T^{\beta \gamma }}\Gamma _{\beta \gamma }^\alpha - \frac{\partial }{{\partial t}}(\sqrt { - g} T_{viscous}^{0\alpha }) - \frac{\partial }{{\partial {x^j}}}(\sqrt { - g} T_{viscous}^{j\alpha }).$$

((62))

6.3 Constrained nonlinear initial data

One cannot take arbitrary data to initiate an evolution of the Einstein equations. The data must satisfy the constraint equations (9) and (10). York [163] developed a procedure to generate proper initial data by introducing conformal transformations of the 3-metric ${\gamma _{ij}} = {\psi ^4}{\tilde \gamma _{ij}}$, the trace-free momentum components ${A^{ij}} = {K^{ij}} - {\gamma ^{ij}}K/3 = {\psi ^{ - 10}}{\tilde A^{ij}}$, and matter source terms ${s^i} = {\psi ^{ - 10}}{\hat s^i}$ and ${\rho _H} = {\psi ^{ - n}}{\hat \rho _H}$, where n > 5 for uniqueness of solutions to the elliptic equation (63) below. In this procedure, the conformal (or “hatted”) variables are freely specifiable. Further decomposing the free momentum variables into transverse and longitudinal components ${\hat A^{ij}} = \hat A_*^{ij} + {(\hat lw)^{ij}}$, the Hamiltonian and momentum constraints are written as

$${\hat \nabla _i}{\hat \nabla ^i}\psi - \frac{{\hat R}}{8}\psi + \frac{1}{8}{\hat A_{ij}}{\hat A^{ij}}{\psi ^{ - 7}} - \frac{1}{{12}}{K^2}{\psi ^5} + 2\pi G\hat \rho {\psi ^{5 - n}} = 0,$$

((63))

$${({\hat \nabla _i}{\hat \nabla ^i}w)^i} + \frac{1}{3}\hat \nabla \left( {{{\hat \nabla }_j}{w^j}} \right) + \hat R_j^i{w^j} - \frac{2}{3}{\psi ^6}{\hat \nabla ^i}K - 8\pi G{\hat s^i} = 0,$$

((64))

where the longitudinal part of ${\hat A^{ij}}$ is reconstructed from the solutions by

$${\left( {\hat lw} \right)^{ij}} = {\hat \nabla ^i}{w^j} + {\hat \nabla ^j}{w^i} - \frac{2}{3}{\hat \gamma ^{ij}}{\hat \nabla _k}{w^k}.$$

((65))

The transverse part of ${\hat A^{ij}}$ is constrained to satisfy ${\hat \nabla _j}\hat A_*^{ij} = \hat A_{*j}^j = 0$.

Equations (63) and (64) form a coupled nonlinear set of elliptic equations which must be solved iteratively, in general. The two equations can, however, be decoupled if a mean curvature slicing (K = K(t)) is assumed. Given the free data $K,{\hat \gamma _{ij}},{\hat s^i} and \hat \rho$, the constraints are solved for $\hat A_*^{ij},{(\hat lw)^{ij}}$ and ψ. The actual metric γ_ij and curvature K_ij are then reconstructed by the corresponding conformal transformations to provide the complete initial data. Anninos [7] describes a procedure using York’s formalism to construct parametrized inhomogeneous initial data in freely specifiable background spacetimes with matter sources. The procedure is general enough to allow gravitational wave and Coulomb nonlinearities in the metric, longitudinal momentum fluctuations, isotropic and anisotropic background spacetimes, and can accommodate the conformal-Newtonian gauge to set up gauge invariant cosmological perturbation solutions as free data.

6.4 Newtonian limit

The Newtonian limit is defined by 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 [124].

6.4.1 Dark and baryonic matter equations

The appropriate perturbation equations in this limit are easily derived for a background FLRW expanding model, assuming a metric of the form

$$d{s^2} = - (1 + 2\Phi )d{t^2} + a{(t)^2}(1 - 2\Phi ){\gamma _{ij}}d{x^i}d{x^j},$$

((66))

where

$${\gamma _{ij}} = \delta _j^i{\left( {1 + \frac{{k{r^2}}}{4}} \right)^{ - 2}},$$

((67))

and k = - 1, 0, + 1 for open, flat and closed Universes. Also, a ≡ 1/(1 + z) is the cosmological scale factor, z is the redshift, and Φ is the comoving inhomogeneous gravitational potential.

The governing equations in the Newtonian limit are the hydrodynamic conservation equations,

$$\frac{{\partial {{\tilde \rho }_b}}}{{\partial t}} + \frac{\partial }{{\partial {x^i}}}({\tilde \rho _b}v_b^i) + 3\frac{{\dot a}}{a}{\tilde \rho _b} = 0,$$

((68))

$$\frac{{\partial ({{\tilde \rho }_b}v_b^j)}}{{\partial t}} + \frac{\partial }{{\partial {x^i}}}({\tilde \rho _b}v_b^iv_b^j) + 5\frac{{\dot a}}{a}{\tilde \rho _b}v_b^j + \frac{1}{{{a^2}}}\frac{{\partial \tilde \rho }}{{\partial {x^j}}} + \frac{{{{\tilde \rho }_b}}}{{{a^2}}}\frac{{\partial \Phi }}{{\partial {x^j}}} = 0,$$

((69))

$$\frac{{\partial \tilde e}}{{\partial t}} + \frac{\partial }{{\partial {x^i}}}(\tilde ev_b^i) + \tilde p\frac{{\partial v_b^i}}{{\partial {x^i}}} + 3\frac{{\dot a}}{a}(\tilde e + \tilde p) = {\tilde S_{cool}},$$

((70))

the geodesic equations for collisionless dust or dark matter (in comoving coordinates),

$$\frac{{dx_d^i}}{{dt}} = v_d^i,$$

((71))

$$\frac{{dv_d^i}}{{dt}} = - 2\frac{{\dot a}}{a}v_d^i - \frac{1}{{{a^2}}}\frac{{\partial \tilde \Phi }}{{\partial {x^i}}},$$

((72))

Poisson’s equation for the gravitational potential,

$${\nabla ^2}\tilde \Phi = 4\pi G{a^2}({\tilde \rho _b} + {\tilde \rho _d} - {\tilde \rho _0}),$$

((73))

and the Friedman equation for the cosmological scale factor,

$$\frac{{da}}{{dt}} = {H_0}{\left[ {{\Omega _m}(\frac{1}{a} - 1) + {\Omega _\Lambda }({a^2} - 1) + 1} \right]^{1/2}}.$$

((74))

Here ${\tilde \rho _d},{\tilde \rho _b},\tilde p and \tilde e$ are the dark matter density, baryonic density, pressure and internal energy density in the proper reference frame, xⁱ and v ⁱ_b are the baryonic comoving coordinates and peculiar velocities, x ⁱ_d and v ⁱ_d are the dark matter comoving coordinates and peculiar velocities, ${\tilde \rho _0} = 3H_0^2{\Omega _0}/(8\pi G{a^3})$ is the proper background density of the Universe, Ω₀ is the total density parameter, Ω_m = Ω_b + Ω_d is the density parameter including both baryonic and dark matter contributions, Ω_Λ = Λ/(3H ²₀ ) is the density parameter attributed to the cosmological constant Λ, H₀ = 100 h km s^-1 Mpc^-1 is the present Hubble constant with 0.5 < h < 1, and ${\tilde S_{cool}}$ represents microphysical radiative cooling and heating rates which can include Compton cooling (or heating) due to interactions of free electrons with the CMBR, bremsstrahlung, and atomic and molecular line cooling. Notice that ‘tilded’ (‘non-tilded’) variables refer to proper (comoving) reference frame attributes.

An alternative total energy conservative form of the hydrodynamics equations that allows high resolution Godunov-type shock capturing techniques is

$$\frac{{\partial {\rho _b}}}{{\partial t}} + \frac{1}{a}\frac{\partial }{{\partial {x^i}}}({\rho _b}\tilde v_b^i) = 0,$$

((75))

$$\frac{{\partial ({\rho _b}\tilde v_b^j)}}{{\partial t}} + \frac{1}{a}\frac{\partial }{{\partial {x^i}}}({\rho _b}\tilde v_b^i\tilde v_b^j + p{\delta _{ij}}) + \frac{{\dot a}}{a}{\rho _b}\tilde v_b^j + \frac{{{\rho _b}}}{a}\frac{{\partial \tilde \Phi }}{{\partial {x^j}}} = 0,$$

((76))

$$\frac{{\partial E}}{{\partial t}} + \frac{1}{a}\frac{\partial }{{\partial {x^i}}}(E\tilde v_b^i + p\tilde v_b^i) + \frac{{2\dot a}}{a}E + \frac{{{\rho _b}\tilde v_b^i}}{a}\frac{{\partial \tilde \Phi }}{{\partial {x^j}}} = {S_{cool}},$$

((77))

with the corresponding particle and gravity equations

$$\frac{{dx_d^i}}{{dt}} = \frac{{\tilde v_d^i}}{a},$$

((78))

$$\frac{{d\tilde v_d^i}}{{dt}} = - \frac{{\dot a}}{a}\tilde v_d^i - \frac{1}{a}\frac{{\partial \tilde \Phi }}{{\partial {x^i}}},$$

((79))

$${\nabla ^2}\tilde \Phi = \frac{{4\pi G}}{a}({\rho _b} + {\rho _d} - {\rho _0}),$$

((80))

where ρ_b is the comoving density, ${\rho _0} = {a^3}{\tilde \rho _0},\tilde v_b^i$ is the proper frame peculiar velocity, p is the comoving pressure, $E = {\rho _b}\tilde v_b^2/2 + p/(\gamma - 1)$ is the total peculiar energy per comoving volume, and Φ the gravitational potential.

6.4.2 Primordial chemistry

The baryonic equations from the previous section are easily extended to include reactive chemistry of both atomic and molecular species (e.g., H, H⁺, He, He⁺, He⁺⁺, H^-, H ⁺₂ , H₂, and e^-) by assuming a common flow field, supplementing the total mass conservation equation (68) with

$$\frac{{\partial {{\tilde \rho }_j}}}{{\partial t}} + \frac{\partial }{{\partial {x^i}}}({\tilde \rho _j}v_b^i) + 3\frac{{\dot a}}{a}{\tilde \rho _j} = \sum\limits_i {\sum\limits_l {{k_{il}}(T)} {{\tilde \rho }_i}{{\tilde \rho }_l}} + \sum\limits_i {{I_i}{{\tilde \rho }_i}} $$

((81))

for each of the species, and including the effects of non-equilibrium radiative cooling and consistent coupling to the hydrodynamics equations. The k_il (T) are rate coefficients for the two body reactions and are generally incorporated in numerical codes as tabulated functions of the gas temperature T. The I_i are integrals evaluating the photoionization and photodissociation of the different species.

A fairly complete chemical network system useful for primordial gas phase compositions, including hydrogen molecules, consists of the following collisional, photoionization, and photodissociation chains

• Collisional reactions (primordial chain):

  • (1) H+e → H⁺ + 2e

  • (2) H+e → H⁺ + 2e

  • (3) H+e → H⁺ + 2e

  • (4) H + e → H⁺ + 2 e

  • (5) He⁺ + e → He⁺⁺ + 2e

  • (6) He⁺⁺ + e → He⁺ + γ

• Collisional reactions (Hz molecular chain):

  • (7) H + e → H^- + γ

  • (8) H + e → H^- + γ

  • (9) H + e → H^- + γ

  • (10) H + e → H^- + γ

  • (11) H + e → H^- + γ

  • (12) H + e → H^- + γ

  • (13) H + e → H^- + γ

  • (14) H + e → H^- + γ

  • (15) H + e → H^- + γ

  • (16) H + e → H^- + γ

  • (17) H + e → H^- + γ

  • (18) H ⁺₂ + H^- → H₂ + H

  • (19) H ⁺₂ + H^- → H₂ + H

• Photoionization reactions (primordial chain):

  • (20) H + γ → H⁺ + e

  • (21) He + γ → He⁺ + e

  • (22) He⁺ + γ → He⁺⁺ + e

• Photodissociation/ionization reactions (molecular chain):

  • (23) H^- + γ → H + e

  • (24) H^- + γ → H + e

  • (25) H^- + γ → H + e

  • (26) H^- + γ → H + e

  • (27) H₂ + γ → 2H

For a comprehensive description of the chemistry and explicit formulas modeling the kinetic and cooling rates relevant for cosmological calculations, the reader is referred to [92, 144, 54, 1, 21]. This reactive network is by no means complete, and in fact, ignores important coolants and contaminants (e.g., HD, LIH, and their intermediary products [151, 78, 48]) expected to form through non-equilibrium reactions at low temperatures and high densities. Although it is certainly possible to include even in three dimensional simulations, the inclusion of more complex reactants demands either more computational resources (to resolve both the chemistry and cooling scales) or an increasing reliance on equilibrium assumptions which can be very inaccurate.

6.4.3 Numerical methods

Many numerical techniques have been developed to solve the hydrodynamic and collisionless particle equations. For the hydrodynamic equations, the methods range from Lagrangian SPH algorithms with artificial viscosity [72, 88], to high resolution shock capturing Eulerian techniques on single static meshes [142, 134], nested grids [19], moving meshes [82], and adaptive mesh refinement [51]. For the dark matter equations, the canonical choices are treecodes [159] or PM and P3M methods [90, 68], although many variants have been developed to optimize computational performance and accuracy, including adaptive mesh, particle, and smoothing kernel refinement methods [45, 77, 130]. An efficient method for solving non-equilibrium, multi-species chemical reactive flows together with the hydrodynamic equations in a background FLRW model is described in [1, 21].

It is beyond the scope of this review to discuss algorithmic details of the different methods and their strengths and weaknesses. Instead, the reader is referred to [103, 77] for thorough comparisons of various numerical methods applied to problems of structure formation. Kang et al. [103] compare (by binning data at different resolutions) the statistical performance of five codes (three Eulerian and two SPH) on the problem of an evolving CDM Universe on large scales using the same initial data. The results indicate that global averages of physical attributes converge in rebinned data, but that some uncertainties remain at small levels. Frenk et al. [77] compare twelve Lagrangian and Eulerian hydrodynamics codes to resolve the formation of a single X-ray cluster in a CDM Universe. The study finds generally good agreement for both dynamical and thermodynamical quantities, but also shows significant differences in the X-ray luminosity, a quantity that is especially sensitive to resolution [17].

6.4.4 Linear initial data

The standard Zel’dovich solution [165, 68] is commonly used to generate initial conditions satisfying observed or theoretical power spectra of matter density fluctuations. Comoving physical displacements and velocities of the collisionless dark matter particles are set according to the power spectrum realization

$${\left| {\frac{{\delta \rho }}{\rho }(k)} \right|^2} \propto {k^n}{T^2}(k),$$

((82))

where the complex phases are chosen from a gaussian random field, T(k) is a transfer function [27] appropriate to a particular structure formation scenario (e.g., CDM), and n = 1 corresponds to the Harrison-Zel’dovich power spectrum. The fluctuations are normalized with top hat smoothing using

$$\sigma _8^2 = \frac{1}{{{b^2}}} = \int_0^\infty {4\pi {k^2}P(k){W^2}(k)dk,} $$

((83))

where b is the bias factor chosen to match present observations of rms density fluctuations in a spherical window of radius R_h = 8 ^h-i. Mpc. Also, P(k) is the Fourier transform of the square of the density fluctuations in Equation (82), and

$$W(k) = \frac{3}{{{{(k{R_h})}^3}}}(\sin (k{R_h}) - (k{R_h})\cos (k{R_h}))$$

((84))

is the Fourier transform of a spherical window of radius R_h.

Overdensity peaks can be filtered on specified spatial or mass scales by Gaussian smoothing the random density field [27]

$$\sigma ({r_o}) = \frac{1}{{{{(2\pi R_f^2)}^{3/2}}}}\int {\frac{{\delta \rho }}{\rho }(r')} \exp \left( { - \frac{{|{r_o} - r'{|^2}}}{{2R_f^2}}} \right){d^3}r'$$

((85))

on a comoving scale R_f centered at r = r_o (for example, R_f = 5 ^h-1 Mpc with a filtered mass of M_f ∼ 10¹⁵M_bd over cluster scales). Nσ peaks are generated by sampling different random field realizations to satisfy the condition v = σ(r₀)/σ-(R_f) = N, where σ(Rf) is the rms of Gaussian filtered density fluctuations over a spherical volume of radius Rf. Bertschinger [44] has provided a useful and publicly available package of programs called COS-MICS for computing transfer functions, CMB anisotropies, and gaussian random initial conditions for numerical structure formation calculations. The package solves the coupled linearized Einstein, Boltzman, and fluid equations for scalar metric perturbations, photons, neutrinos, baryons, and collisionless dark matter in a background isotropic Universe. It also generates constrained or unconstrained matter distributions over arbitrarily specifiable spatial or mass scales.

Bertschinger [44] has provided a useful and publicly available package of programs called COSMICS for computing transfer functions, CMB anisotropies, and gaussian random initial conditions for numerical structure formation calculations. The package solves the coupled linearized Einstein, Boltzman, and fluid equations for scalar metric perturbations, photons, neutrinos, baryons, and collisionless dark matter in a background isotropic Universe. It also generates constrained or unconstrained matter distributions over arbitrarily specifiable spatial or mass scales.