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^{5} equal-mass particles, but it is giving way
to the
relentless advance of technology (*N* = 10^{9} 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^{9}
*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^{5} 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^{3}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.