The appropriate perturbation equations in this limit are easily derived for a background FLRW expanding model, assuming a metric of the form
The governing equations in the Newtonian limit are the hydrodynamic conservation equations,
An alternative total energy conservative form of the hydrodynamics equations that allows high resolution Godunov-type shock capturing techniques is
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
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):
Collisional reactions ( molecular chain):
Photoionization reactions (primordial chain):
Photodissociation/ionization reactions (molecular chain):
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.
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 , moving meshes , and adaptive mesh refinement . For the dark matter equations, the canonical choices are treecodes  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.  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.  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 .
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 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
Overdensity peaks can be filtered on specified spatial or mass scales by Gaussian smoothing the random density field 
Bertschinger  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.
© Max Planck Society and the author(s)