The appropriate perturbation equations in this limit are easily derived for a background FLRW expanding model, assuming a metric of the form
where and for open, flat and closed Universes. Also, is the cosmological scale factor, is the redshift, and is the comoving inhomogeneous gravitational potential.The governing equations in the Newtonian limit are the hydrodynamic conservation equations,
the geodesic equations for collisionless dust or dark matter (in comoving coordinates), Poisson’s equation for the gravitational potential, and the Friedman equation for the cosmological scale factor, Here , , and are the dark matter density, baryonic density, pressure and internal energy density in the proper reference frame, and are the baryonic comoving coordinates and peculiar velocities, and are the dark matter comoving coordinates and peculiar velocities, is the proper background density of the Universe, is the total density parameter, is the density parameter including both baryonic and dark matter contributions, is the density parameter attributed to the cosmological constant , is the present Hubble constant with , and 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’ (‘nontilded’) variables refer to proper (comoving) reference frame attributes.An alternative total energy conservative form of the hydrodynamics equations that allows high resolution Godunovtype shock capturing techniques is
with the corresponding particle and gravity equations where is the comoving density, , is the proper frame peculiar velocity, is the comoving pressure, is the total peculiar energy per comoving volume, and is the gravitational potential.
The baryonic equations from the previous section are easily extended to include reactive chemistry of both atomic and molecular species (e.g., , , , , , , , , and ) by assuming a common flow field, supplementing the total mass conservation equation (68) with
for each of the species, and including the effects of nonequilibrium radiative cooling and consistent coupling to the hydrodynamics equations. The are rate coefficients for the two body reactions and are generally incorporated in numerical codes as tabulated functions of the gas temperature . The 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)  
(2)  
(3)  
(4)  
(5)  
(6)  
Collisional reactions ( molecular chain):


(7)  
(8)  
(9)  
(10)  
(11)  
(12)  
(13)  
(14)  
(15)  
(16)  
(17)  
(18)  
(19)  
Photoionization reactions (primordial chain):


(20)  
(21)  
(22)  
Photodissociation/ionization reactions (molecular chain):


(23)  
(24)  
(25)  
(26)  
(27)  
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 nonequilibrium 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.
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 P^{3}M 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 nonequilibrium, multispecies 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 Xray cluster in a CDM Universe. The study finds generally good agreement for both dynamical and thermodynamical quantities, but also shows significant differences in the Xray luminosity, a quantity that is especially sensitive to resolution [17].
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
where the complex phases are chosen from a gaussian random field, is a transfer function [27] appropriate to a particular structure formation scenario (e.g., CDM), and corresponds to the Harrison–Zel’dovich power spectrum. The fluctuations are normalized with top hat smoothing using where is the bias factor chosen to match present observations of rms density fluctuations in a spherical window of radius . Also, is the Fourier transform of the square of the density fluctuations in Equation (82), and is the Fourier transform of a spherical window of radius .Overdensity peaks can be filtered on specified spatial or mass scales by Gaussian smoothing the random density field [27]
on a comoving scale centered at (for example, with a filtered mass of over cluster scales). peaks are generated by sampling different random field realizations to satisfy the condition , where is the rms of Gaussian filtered density fluctuations over a spherical volume of radius .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.
http://www.livingreviews.org/lrr20012 
© Max Planck Society and the author(s)
Problems/comments to 