Annu. Rev. Astron. Astrophys. 1998. 36: 599-654 Copyright © 1998 by Annual Reviews. All rights reserved

2.5. Numerical and Physical Limitations

Cosmological simulations suffer from limited dynamic range, numerical errors, and neglected physics. These effects have recently begun to receive greater attention as increasing computer power enables more rigorous testing. The introduction of gas dynamics into simulation codes has also presented a new set of numerical issues, stimulating practitioners to assess more carefully the limitations of their results.

In purely gravitational N-body simulations of structure formation, the numerical issues include dynamic range, force accuracy, and time integration accuracy. Most codes use the leapfrog scheme (Section 2.1); although this is much less accurate than the high-order schemes used for globular cluster or Solar System integration, tests suggest that with a timestep small enough so that all particles travel less than a softening length each timestep, the statistical results of cosmological interest converge (Quinn et al 1997). Most workers tune their force evaluation algorithm to have pairwise errors of at most a few percent (e.g. Pearce & Couchman 1997), which Hernquist et al (1993) concluded is adequate given the masking effect of two-body relaxation. This leaves dynamic range as the biggest concern. Three distinct kinds of dynamic range are needed for a faithful simulation: mass resolution (number of particles), initial power spectrum sampling (range of wavenumbers present in the initial conditions), and spatial resolution (force-softening length compared with box size).

Mass resolution was perhaps the most severe problem when simulations had fewer than 10⁵ equal-mass particles, but it is giving way to the relentless advance of technology (N = 10⁹ will become the standard on supercomputers within a very few years) and to multimass methods (e.g. Katz et al 1994, Splinter 1996). Particle masses of 10⁹ M or less are probably satisfactory for simulating the gross dynamics of halos of galaxies like the Milky Way, and this resolution is readily available today.

Resolution of the initial power spectrum is related to the mass resolution because the initial density field is sampled at the particle positions. Little et al (1991) suggested that initial high-frequency power is unimportant on scales that become nonlinear, although Splinter et al (1998) argued to the contrary. Tormen & Bertschinger (1996), Cole (1997) suggested ways to extend the sampling of the initial power spectrum by adding long waves to an evolved simulation.

Force resolution is probably the most significant dynamic range issue for N-body simulations [although Melott & Shandarin (1989) argued that mass resolution is equally important]. Resolving the dynamics of galaxies in simulations sufficiently large to sample large-scale tidal fields requires 1-kpc resolution in boxes of size 100 Mpc, a dynamic range of 10⁵ in length scale, about a factor of 10 beyond the current state of the art of high-resolution codes. Smaller force softening requires more timesteps (Quinn et al 1997) and leads to increased two-body relaxation (Huang et al 1993). The practical figure of merit is the ratio of force softening distance to mean interparticle spacing. With P³M and tree codes, this ratio is typically set to about 0.1, while with PM codes, it is about unity. Most practitioners have accepted the benefits of higher-resolution enabling codes to follow gravitational collapse to high density. Low force resolution leads to excessive merging of halos (e.g. Carlberg 1994, Gelb & Bertschinger 1994a, Zurek et al 1994, Moore et al 1996). However, Melott et al (1997) showed that high resolution leads to spurious heating and fragmentation for the collapse of planar Zel'dovich pancakes, and they advocated a force softening no smaller than the mean interparticle spacing. Splinter et al (1998) found less dramatic effects with realistic three-dimensional initial conditions, yet they also concluded that cosmological N-body simulations cannot be trusted on scales smaller than the mean interparticle spacing. While not all would agree, the issues they raise merit further investigation.

When gas dynamics is added, several new types of errors can arise. The most obvious (and also often the most subtle) are numerical problems with the hydrodynamics algorithm, which can best be diagnosed by testing against exact solutions and other codes (Kang et al 1994; CS Frenk, SDM White, P Bode, JR Bond, GL Bryan, et al, unpublished manuscript). Two-body relaxation of dark matter can also cause spurious heating of the gas solely by the gravitational coupling; Steinmetz & White (1997) studied this effect and advised on the resolution requirements to avoid it.

Self-gravity and cooling lead to potential problems for gas dynamics codes. Gnedin & Bertschinger (1996) demonstrated the need for consistent treatment of gravity and hydrodynamics in order to ensure local energy conservation. Bate & Burkert (1997) compared SPH with an Eulerian code and showed that SPH may give incorrect results unless the minimum resolvable SPH mass is less than the Jeans mass. Truelove et al (1997) found a similar requirement for Eulerian hydrodynamics by using adaptive mesh refinement, as did Owen & Villumsen (1997) in two-dimensional simulations using SPH for gas coupled with PM for dark matter. Evrard et al (1994) had previously noted the strong resolution dependence of gas dynamical results for cold self-gravitating gas. Even in a cluster in hydrostatic equilibrium, Anninos & Norman (1996) found that the X-ray luminosity depends strongly on numerical resolution. Clearly these resolution effects need further study.

With gravity and cooling but no heating or star formation, the baryons in high-resolution simulations collapse to high density, impeding further progress (e.g. Katz & White 1993). Recent work has focused on the role of photoionization heating in suppressing galaxy formation (Efstathiou 1992;, Navarro & Steinmetz 1997 and references therein). In the course of investigating these effects, Weinberg et al (1997a) discovered an interplay between mass resolution, photoionization, and radiative cooling that led to incorrect results on galaxy formation when photoionization was included at low resolution. Their experience is a cautionary reminder that when physics is added or codes are run in a new regime, careful testing must precede astrophysical conclusions.