# A Database of COBE-Normalized CDM Simulations (Abbreviated Version)

## 1 INTRODUCTION

### 1.1 Importance of Numerical Simulations in Cosmology

Observations of the nearby universe reveal the existence of the large-scale structure. The visible matter is clumped into galaxies, and these galaxies are not distributed uniformly into space, but instead grouped into structures such as clusters, filaments, and walls, separated by deep voids. Velocity structures (deviations from Hubble flow) are observed as well. By contrast, observations of the Cosmic Microwave Background (CMB) reveal that the universe was extremely uniform near the epoch of recombination. Hence, the present large-scale structure must result from an evolutionary process that took place between recombination and the present. The most widely accepted scenario assumes that the present large-scale structure originates from the growth, by gravitational instability, of primordial density fluctuations present in the early universe. Any fluctuation larger than the Jeans length can grow by gravitational instability once the universe becomes matter-dominated. If the primordial density fluctuations originates from a Gaussian random process (the usual assumption), then the primordial density field is entirely described in terms of its power spectrum $`P(k)`$. The particular form of the power spectrum essentially depends upon the amount and nature of the various components (baryonic matter, dark matter, cosmological constant, and so on) present in the universe. If we assume a certain power spectrum, we can describe the primordial density field, and the formation and evolution of large scale-structure in the universe becomes an initial value problem: starting from the primordial density field, we can compute its evolution using the laws of general relativity. Unfortunately, this initial value problem is far too complex to be solved analytically. We can simplify the problem by noticing that the largest structures observed in the universe are significantly smaller than the horizon. This enables us to describe the evolution of the large-scale structure using Newtonian mechanics instead of general relativity (Peebles 1980, Chapter 2). Even so, the general problem cannot be solved analytically. This leaves two possible approaches: analytical approximations, or numerical simulations.

Two different kinds of analytical approximations have been considered. The first one is based on the fact that the initial fluctuations are small. We can expand the equations describing the evolution of these fluctuations in powers of the fluctuations, and solve them using perturbation theory. This approach is extremely useful in describing the early evolution of the fluctuations, and has led to very important results. However, it becomes inapplicable as soon as the fluctuations become nonlinear. Such fluctuations still have to grow by a factor of $`10^2`$ to reach the density of a cluster of galaxies, and $`10^5`$ to reach the density of a galaxy. Clearly, perturbation theory cannot be used to describe the late stages of large-scale structure formation. The second analytical approach consists of considering systems with a particular geometry (see, for instance, Zel’dovich 1970; Peebles 1980, §§19–21; Fillmore & Goldreich 1984a, 1984b; Bertschinger 1985a, 1985b). The most popular analytical models for large-scale structure formation are the Spherical Model, which assumes spherical symmetry, and the Pancake Model, which assumes planar symmetry. An important assumption of some of these models is that the system considered is isolated. For instance, the spherical model can describe the evolution of a self-gravitating spherical overdensity, but we must assume that any tidal influence from nearby structures can be neglected, an assumption that might be valid at late times but certainly not at early times.

These analytical approximations can therefore describe the universe at early times or at late times, but not both. This problem can be solved by using mixed schemes, that combine various analytical approximations in a way that allows an analytical description of the evolution of large-scale structure at all epochs. The most important ones are the Press-Schechter Approximation (Press & Schechter 1979), which combines perturbation theory with the spherical model, and the Zel’dovich Approximation (Zel’dovich 1970), which combines perturbation theory with the pancake model.

The alternative consists of using numerical simulation. Unlike analytical models, numerical simulations suffer from problems such as limited resolution and numerical noise. Also, simulations provide very little insight into the physical processes taking place, compared to analytical models. However, numerical simulations can describe the evolution of the large-scale structure entirely, from the initial conditions all the way to the present, without making any approximation or imposing any restriction on the geometry of the system. Cosmological N-body simulations have played a central role in the study of the formation and evolution of large-scale structure in the universe during most of the last two decades. These simulations have contributed the improve our understanding of the physical process of gravitational instability that leads to structure formation, have enable us to conceive and test various cosmological scenarios, and have produced simulated universes that can directly be compared with observations (Efstathiou & Eastwood 1981; Centrella & Melott 1983; Klypin & Shandarin 1983; Miller 1983; Shapiro, Struck-Marcell, & Melott 1983; White, Frenk, & Davis 1983; Davis et al. 1985; Efstathiou et al. 1985; Barnes & Hut 1986, 1989; Evrard 1986, 1987; Melott 1986; White et al. 1987a, 1987b; Frenk et al. 1988; Gramann 1988; Carlberg & Couchman 1989; Villumsen 1989; West, Oemler, & Dekel 1989; Couchman 1991; Fukushige et al. 1991; Hernquist, Bouchet, & Suto 1991; Martel 1991a; Moutarde et al. 1991; West, Villumsen, & Dekel 1991; Bouchet & Hernquist 1992; Fry, Melott, & Shandarin 1992; Park et al. 1992; Bahcall, Cen, & Gramann, 1993; Gramann, Cen, & Bahcall 1993; Melott & Shandarin 1993; Babul et al 1994; Pen 1995; Colombi, Bouchet, & Hernquist 1996; Moore, Katz, & Lake 1996: Yess & Shandarin 1996; Klypin, Nolthenius, & Primack 1997; Kravtsov, Klypin, & Khokhlov 1997; Navarro, Frenk, & White 1997; Gross et al. 1998; Thomas et al. 1998).<sup>1</sup><sup>1</sup>1Cosmological numerical simulations have been used intensively since 1981, and their results have appeared in hundreds of publications, so this list is necessarily incomplete. We decided to include only the key publications by each research group. We also excluded one- and two-dimensional simulations (for brevity), and simulations with hydrodynamics, which involve the next generation of numerical algorithms.

### 1.2 The Standard Model

Theoretical developments in particle theory and early-universe physics, combined with numerical simulations and observations of the large-scale structure of the universe, led to the emergence during the 1980’s of what became known as the Standard Cosmological Model. The inflationary scenario (Guth 1981; Linde 1982; Albretch & Steinhardt 1982) requires that the universe is spatially flat to a great accuracy. In a matter-dominated universe, in the absence of any exotic components such as a nonzero cosmological constant $`\mathrm{\Lambda }`$, this requires that the mean density of the universe is equal to its critical density, or, alternatively,

$$\mathrm{\Omega }_0=1,$$
(1)   
where $`\mathrm{\Omega }_08\pi G\overline{\rho }_0/3H_0^2`$ is the density parameter, $`\overline{\rho }_0`$ is the mean density of the universe, and $`H_0`$ is the Hubble constant (throughout this paper, we use subscripts 0 to designate present values). In this scenario, the large-scale structure of the universe that we observe today results from the growth, by gravitational instability, of small density perturbations present at recombination, which originate from quantum processes in the early universe. There are two difficulties with this scenario. First, there is strong evidence that the amount of “ordinary matter” in the universe is insufficient to satisfy equation (1). Primordial nucleosynthesis provides stringent upper limit to the baryonic content of the universe, and shows that the present baryon contribution to the density parameter, $`\mathrm{\Omega }_{\mathrm{B0}}`$, is less than $`0.026h^2`$, where $`h`$ is the Hubble constant in units of $`100\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ (Krauss & Kernan 1995; Copi, Schramm, & Turner 1995; Krauss 1998). Furthermore, dynamical studies of rich clusters of galaxies show that the contribution to the density parameter of the matter that clusters at that scale is $`\mathrm{\Omega }_{\mathrm{clusters}}=0.2\pm 0.1`$ (Gott et al. 1974; Carlberg et al. 1996; Lin et al. 1996). Second, observations of the temperature fluctuations in the Cosmic Microwave Background (CMB) provide an upper limit to the density fluctuations at recombination, of order $`\delta \rho /\rho 10^5`$. Such fluctuations would grow by gravitational instability, to reach an amplitude of order $`10^2`$ by the present, clearly insufficient to explain the origin of large-scale structure and galaxies.

These two difficulties were solved by postulating the existence of a component known as dark matter. This dark matter, which can be detected only through its gravitational influence, makes up the difference between the amount of matter required to satisfy equation (1) and the amount of matter that is observed or indirectly measured. We can then reconcile the dynamical estimates of $`\mathrm{\Omega }_0`$ with equation (1) by assuming that the dark matter is distributed more smoothly than the luminous matter that makes up galaxies and clusters, an idea known as biasing (Kaiser 1984). However, since the dynamical estimates of $`\mathrm{\Omega }_0`$ exceed the limit imposed by primordial nucleosynthesis, some amount of dark matter must be clustered on galactic and cluster scales, even though the bulk of the dark matter is smoothly distributed into space. Also, density fluctuations in the dark matter start growing when the universe becomes matter-dominated, and by the time recombination occurs, the dark matter fluctuations have already grown by a factor of $`20h^2/\mathrm{\Omega }_{\mathrm{B0}}`$ (Kolb & Turner 1990, §9.5). These fluctuations provide potential wells into which the baryons fall soon after recombination. This enables us to reconcile the CMB measurements, — which are sensitive to the fluctuations in the baryonic matter, but not the dark matter — with the existence of galaxies and large-scale structure, Many candidates for dark matter have been suggested, several of them emerging from particle theory. These various forms of dark matter are usually classified as Hot Dark Matter (HDM) and Cold Dark Matter (CDM).

Cosmological numerical simulations performed in the 1980’s have shown that an $`\mathrm{\Omega }_0=1`$ universe in which most matter is in the form of Cold Dark Matter, and the distribution of luminous matter is biased relative to the distribution of dark matter, can successfully reproduce all observations of large-scale structure that were then known (see, e.g. Davis et al. 1985), while satisfying the constraints imposed by the inflationary scenario, primordial nucleosynthesis, and the CMB (unlike HDM, which as difficulties explaining structure formation at galactic scales). These simulations played a central role in establishing the $`\mathrm{\Omega }_0=1`$, biased CDM model as the Standard Cosmological Model.

### 1.3 Non-Standard Models

The Standard Model, which has been hailed by many as the “final answer” to the problem of structure formation and evolution in the universe, ran into serious problems during the 1990’s. In this subsection, we briefly review these problems.

#### 1.3.1 The Age Problem

In a flat, matter-dominated universe, the age of the universe $`t_0`$ is 2/3 of the Hubble time, that is $`t_0=2/3H_0=6.52\times 10^9h^1\mathrm{years}`$. For $`h`$ in the range $`0.51`$, this corresponds to an age in the range $`6.5213.04\times 10^9\mathrm{years}`$. Measurements of the ages of globular clusters indicate that the oldest clusters are certainly older than $`9.5\times 10^9\mathrm{years}`$, and most likely in the range $`1113\times 10^9\mathrm{years}`$ (Jimenez et al. 1996; Chaboyer 1998; Chaboyer et al. 1998; Jimenez 1998). These measurements are only marginally consistent with the standard model, in the sense that they require a Hubble constant near its smallest possible value, $`h0.5`$. However, recent observations have significantly reduced the range of plausible values for the Hubble constant, showing that $`h`$ is likely to be in the range $`0.650.75`$ (Freedman, W. L. 1998, and references therein). In the standard model, this corresponds to an age in the range $`8.6910.03\times 10^9\mathrm{years}`$. The upper end of this range is still compatible with the measured ages of globular clusters, but just barely.

#### 1.3.2 Large-Scale Structure

Until 1992, the amplitude of the primordial density fluctuation power spectrum was unknown. We were free to tune this amplitude in order to reproduce the correct amount of galaxy clustering observed today. The latter is usually characterized by the rms density fluctuation $`\sigma _8`$ at a scale of $`8h^1\mathrm{Mpc}`$. Observations of clusters of galaxies show that $`\sigma _80.6`$ for the standard model (Viana & Liddle 1996, see eq.  below). The discovery by the COBE DMR experiment of degree-scale fluctuations in the CMB temperature (Smoot et al. 1992) has eliminated this freedom, by fixing the amplitude power spectrum. A COBE-normalized Standard Model produces too much structure at cluster scales (Barlett & Silk 1993). The resulting value of $`\sigma _8`$ is $`1.22`$ for $`h=0.5`$ (Bunn & White 1997), too large by a factor of 2, and becomes even larger for larger $`h`$. By combining the COBE result with observations of the present large-scale structure, we obtain a constraint on the quantity $`\mathrm{\Omega }_0h`$, which is

$$0.2\mathrm{\Omega }_0h0.3$$
(2)   
(Peacock & Dodds 1994). Since $`h`$ is certainly larger than 0.5, this implies $`\mathrm{\Omega }_0<0.6`$.

#### 1.3.3 The Baryon Catastrophe

In the Standard Model, the baryon fraction of the universe is small. Primordial nucleosynthesis imposes the constraint that $`\mathrm{\Omega }_{\mathrm{B0}}h^2<0.026`$ (Krauss & Kernan 1995; Copi, Schramm, & Turner 1995; Krauss 1998). For $`h=0.65`$, this corresponds to $`\mathrm{\Omega }_{\mathrm{B0}}=0.061`$. Hence, if $`\mathrm{\Omega }_0=1`$, at most 6% of the matter if the universe is composed of baryons, the rest being dark matter. However, observations of X-ray clusters reveal that the baryon fraction in these clusters is $`0.1h^{1.5}`$, or 19% for $`h=0.65`$ (Briel, Henry, & Boringer 1992). Hence, X-ray clusters contain a large excess of baryons relative to dark matter compared with the average values in the universe, a situation often referred to as the baryon catastrophe. This problem could be solved if we can think of a physical process that would concentrate the baryons inside clusters, creating a bias relative to the dark matter. However, no such physical process is known. It is much simplier to assume that the density parameter is less than unity. In this case, the universal baryon fraction is not $`\mathrm{\Omega }_{\mathrm{B0}}`$, but $`\mathrm{\Omega }_{\mathrm{B0}}/\mathrm{\Omega }_0`$. This baryon fraction is $`0.1h^{1.5}`$ according to observations of X-ray clusters, and smaller than $`0.026/\mathrm{\Omega }_0h^2`$ according to primordial nucleosynthesis. Combining these two results, we get

$$\mathrm{\Omega }_0h^{1/2}<0.26.$$
(3)   
This rules out the Standard Model (unless $`h<0.07`$ !). For $`h>0.5`$, this limit becomes $`\mathrm{\Omega }_0<0.37`$.

#### 1.3.4 Evolution of Cluster Abundance

In the Standard Model, the density parameter $`\mathrm{\Omega }`$ is unity at all times, and density fluctuations can grow by gravitational instability all the way to the present. In other models, density fluctuations can grow at early times, when $`\mathrm{\Omega }`$ is near unity. But eventually, $`\mathrm{\Omega }`$ drops significantly below unity, and the density fluctuations “freezes-out.” Hence, in a model with $`\mathrm{\Omega }_0<1`$, the present abundance of clusters should be comparable to the abundance immediately after freeze-out, since not much growth has taken place since then. Conversely, in the Standard Model, the present abundance of clusters should be larger than the past one, since the growth of density fluctuations never freezes out. Bahcall, Fan, & Cen (1997) and Bahcall & Fan (1998) have determined the mass of three massive distant clusters, located at redshift $`z>0.5`$, and showed that in a $`\mathrm{\Omega }_0`$ universe, there should be only $`10^3`$ such clusters at $`z>0.5`$. They conclude that the density parameter is in the range

$$0.1<\mathrm{\Omega }_0<0.35$$
(4)   
(Bahcall 1999).

#### 1.3.5 Distant Type I Supernovae

The relationship between the luminosity distance $`D_L`$ and the redshift is model-dependent. If standard candles can be observed at cosmological distances, then the $`D_L(z)`$ relationship can be inferred, and limits can be placed on the value of the cosmological parameters. This method was recently applied to samples of distant (“High-$`z`$”) Type I supernovae, by two independent research teams (Garnevich et al. 1998, and references therein; Perlmutter et al. 1998, and references therein). Applied to models with a nonzero cosmological constant $`\mathrm{\Lambda }`$, their observations provide severe constraints in the $`\mathrm{\Omega }_0\lambda _0`$ phase space (where $`\lambda _0\mathrm{\Lambda }/3H_0^2`$). Not only the Standard Model is excluded with a high degree of confidence, but open models ($`\mathrm{\Omega }_0<1`$), without a cosmological constant are also excluded, unless $`\mathrm{\Omega }_0`$ is very small. Observations of the CMB (White 1998; Tegmark et al. 1998) provide a different constraint, that does not rule out the Standard Model. However, combining the CMB and Type I supernovae observations leads to separate determinations of $`\mathrm{\Omega }_0`$ and $`\lambda _0`$. The preferred values are $`\mathrm{\Omega }_00.3`$ and $`\lambda _00.7`$.

#### 1.3.6 Anthropic Considerations

In models such as chaotic inflation, in which the observed big bang is just one of an infinite number of expanding regions in each of which the fundamental property takes a different value (Linde 1986, 1987, 1988), and models in which a state vector is derived for the universe which is a superposition of terms with different values of the fundamental property (e.g. Hawking 1983, 1984; Coleman 1988), the probability of observing any particular values of the cosmological parameters is conditioned by the existence of observers in those “subuniverses” in which the parameters take these values (Efstathiou 1995; Vilenkin 1995; Weinberg 1996; Martel, Shapiro, & Weinberg 1998). This probability is proportional to the fraction of matter which is destined to condense out of the background into mass concentrations large enough to form observers. Using this approach, Martel et al. (1998) calculated the relative likelihood of observing any given value of the cosmological constant $`\mathrm{\Lambda }`$ within the context of the flat CDM model normalized to COBE, and found that small but finite value of the cosmological constant, in the range suggested by observations, are favored over the value $`\mathrm{\Lambda }=0`$. Garriga, Tanaka, & Vilenkin (1998) have performed a similar analysis, but applied the the density parameter, and found that intermediate values of $`\mathrm{\Omega }_0`$ are more likely to be observed than values near 0 or near 1. Anthropic arguments do not favor the values that the parameters take in the Standard Model, $`\mathrm{\Omega }_0=1`$ and $`\lambda _0=0`$.

#### 1.3.7 Alternatives to the Standard Model

The problems listed above strongly argue against the Standard Model, and forces us to consider alternatives. The age problem, large-scale structure problem, and the baryon catastrophe can all be solved by considering Open CDM, or OCDM models, in which the density parameter $`\mathrm{\Omega }_0<1`$. However, such models do not satisfy the flatness requirement of inflation. CDM models with a nonzero cosmological constant $`\lambda _0`$ equal to $`1\mathrm{\Omega }_0`$, known as $`\mathrm{\Lambda }`$CDM models, satisfy this flatness requirement, and the addition of the cosmological constant improves the age and large-scale structure problems, while providing a better agreement to the supernovae data. Recently, several authors have shown that it is possible to reconcile the inflationary scenario with an open universe, thus eliminating the flatness requirement (Ratra & Peebles 1994; Bucher, Goldhaber, & Turok 1995; Yamamoto, Sasaki, & Tanaka 1995; Linde 1995; Linde & Mezhlumian 1995). This not only supports open, matter-dominated models, but also allows for the possibility of an open universe with $`\lambda _00`$ and $`\mathrm{\Omega }_0+\lambda _0<1`$.

The large-scale structure problem can also be solved by introducing a “tilt” in the primordial power spectrum. In this Tilted CDM, or TCDM model, the primordial power spectrum $`P(k)`$ at large scales does no have the Harrison-Zel’dovich form $`P(k)k`$, but instead varies as $`P(k)k^n`$, where the primordial exponent $`n`$ can differ from unity. The universe might also contain a mixture of two different forms of dark matter, one cold and one hot, a model known as CHDM. Finally, the universe might contain a smooth component whose pressure $`p`$ and density $`\rho `$ are related by an equation of state $`p=w\rho `$, a concept known as “quintessence” (Caldwell, Dave, & Steinhardt 1998; see also Fry 1985; Charlton & Turner 1987; Silveira & Waga 1994; Martel 1995; Martel & Shapiro 1998). The cosmological constant is a particular form of quintessence, corresponding to $`w=1`$; other forms have been suggested, such as domain walls, textures, or strings.

The Standard Model had no free parameters. The values of $`\mathrm{\Omega }_0`$ and $`\lambda _0`$ were fixed at 1 and 0, respectively. The value of $`h`$ had to be close to $`0.5`$ in order to avoid conflicts with ages of globular clusters, and the primordial power spectrum was assumed to be a CDM spectrum with no tilt. With the emergence of alternative models, there are now many free parameters. The density parameters $`\mathrm{\Omega }_0`$ no longer has to be unity. The cosmological constant $`\lambda _0`$ can be nonzero and, if open inflation is correct, a nonzero $`\lambda _0`$ does not have to be equal to $`1\mathrm{\Omega }_0`$. The Hubble constant can vary over a certain range without conflicting with observations, and the slope $`n`$ of the primordial power spectrum does not have to be 1. In models such as CHDM and quintessence models, there are additional parameters: the contribution of each component to the density parameter, and in the case of quintessence models, the coefficient $`w`$ appearing in the equation of state.

During the 1980’s, N-body simulations have played a central role in establishing the Standard Model, and then went more or less into hibernation, as efforts were invested into adding more physics to the original algorithms (hydrodynamics in particular). The emergence of alternative cosmological models has lead to a renewal of interest in N-body simulations. Such numerical simulations are essential for testing cosmological models against observations. Furthermore, they are useful from a theoretical viewpoint, since they can reveal how each cosmological parameter affects the process of large-scale structure formation.

### 1.4 The Need for a Database

Numerical methods such as the Particle-Mesh algorithm (PM) and the Particle-Particle/Particle-Mesh algorithm ($`\mathrm{P}^3\mathrm{M}`$) are well documented. Details of the algorithms can be found in textbooks (e.g. Hockney & Eastwood 1981) and papers (e.g. Efstathiou et al. 1985). Hence, any researcher can easily access all the information and knowledge necessary to develop such algorithms. However, the effort required to develop, test, and optimize a PM or $`\mathrm{P}^3\mathrm{M}`$ algorithm from scratch can be quite substantial, and can be regarded as a waste of effort, since it essentially amounts to “reinventing the wheel.” Also, performing simulations with large number of particles can demand a substantial investment in resources such as computer time, which is also wasteful if these simulations, or similar ones, have already been performed by other researchers. Consequently, it is a common practice among researchers to share either their programs or the results of their simulations.

Klypin & Holtzman (1997) have combined into a single package their version of the PM algorithm and programs for generating initial conditions and analyzing the results. This package has been made available to the astronomical community, and can can be downloaded from a world-wide-web site. This allows other researchers interested in performing cosmological numerical simulations to “get started” immediately, without having to develop and test any computer program. However, installing and running these programs might pose some difficulties depending upon the kind of computer resources available to the user.

We use a different approach. Instead of making our programs available to the astronomical community (something we might do eventually), it is the results of the simulations themselves that we are making available. We performed a very large number of numerical simulations, a total of 160, for 68 different cosmological models. This constitutes by far the largest database of cosmological simulations ever assembled, and it is still growing as more simulations are being performed. We are making this database available to the astronomical community (see §3.4 below). This approach is complementary to the one used by Klypin & Holtzman. By providing the results of the simulations, we eliminate the need for researchers to perform themselves these simulations, and the same simulations can be used by many different researchers. However, someone might be interested in simulating a cosmological model which is not included in the database, in which case the algorithm of Klypin & Holtzman can be used. Alternatively, we can, upon request, perform additional simulations and include them in the database. An interesting question is whether the results of simulations from the database can be analyzed using Klypin & Holtzman programs. In principle, this should be possible. The output files in the database are not written in the same format as the ones produced by Klypin & Holtzman’s PM code, but it is fairly trivial to write a program that “translate” files from one format to another.

There is an important difference that must be pointed out. The program of Klypin & Holtzman is based on the PM algorithm, while the simulations in the database were performed using a $`\mathrm{P}^3\mathrm{M}`$ algorithm. For a same number of particles, the $`\mathrm{P}^3\mathrm{M}`$ algorithm has a length resolution superior to the one of the PM algorithm by a factor of order 6 (depending upon the particular choice of smoothing length). However, since the PM algorithm is significantly faster than the $`\mathrm{P}^3\mathrm{M}`$ algorithm, it is possible to make up for the lack of resolution of the PM code by simply using more particles. We used $`64^3`$ particles in all simulations.<sup>2</sup><sup>2</sup>2We intend to add simulations with $`128^3`$ particles to the database in a near future. These have the same length resolutions as PM calculations with $`384^3`$ particles (such as the ones performed by Gross et al. ).

As mentioned above, there are numerous alternatives to the Standard Model. In this paper we consider CDM models in which the only components are ordinary matter (dark and baryonic) and possibly a nonzero cosmological constant (thus excluding CHDM models, and generic quintessence models). We consider the three cases $`\mathrm{\Omega }_0=1`$, $`\lambda _0=0`$ (the Einstein-de Sitter model), $`\mathrm{\Omega }_0<1`$, $`\lambda _0=0`$, and $`\mathrm{\Omega }_0+\lambda _0=1`$. We also allow the primordial power spectrum to have a tilt. These models are usually referred to as Tilted CDM (TCDM), Tilted, Open CDM (TOCDM), and Tilted, Lambda CDM (T$`\mathrm{\Lambda }`$CDM).

## 2 THE NUMERICAL SIMULATIONS

### 2.1 The Algorithm

All simulations presented in this paper were done using the P<sup>3</sup>M algorithm (Hockney & Eastwood 1981; Efstathiou et al. 1985). The computational volume is a cubic box of comoving size $`L_{\mathrm{box}}`$ and comoving volume $`V_{\mathrm{box}}=L_{\mathrm{box}}^3`$ with triply periodic boundary conditions, expanding with Hubble flow. The matter distribution inside the computational volume is represented by $`N`$ equal-mass particles. The forces on particles are computed by solving Poisson’s equation on a cubic grid using a Fast Fourier Transform method. The resulting force field represents the Newtonian interaction between particles down to a separation of a few mesh spacings. At shorter distances the computed force is significantly smaller than the physical force. To increase the dynamical range of the code, the force at short distance is corrected by direct summation over pairs of particles separated by less than some cutoff distance $`r_e`$. With the addition of this so-called short-range correction, the code accurately reproduces the Newtonian interaction down to the softening length $`ϵ`$, which is a fraction of the grid spacing. The system is evolved forward in time using a second order Runge-Kutta time-integration scheme with a variable time step.

Our particular version of the $`\mathrm{P}^3\mathrm{M}`$ algorithm uses supercomoving variables (Martel & Shapiro 1998; see also Shandarin 1980). In these variables, the position $`\stackrel{~}{𝐫}`$, peculiar velocity $`\stackrel{~}{𝐯}`$, time $`\stackrel{~}{t}`$, density $`\stackrel{~}{\rho }`$, and peculiar gravitational potential $`\stackrel{~}{\varphi }`$ are related to their Eulerian counterparts by

$`\stackrel{~}{𝐫}`$ $`=`$ $`{\displaystyle \frac{𝐫}{ar_{}}},`$ (5)   
$`\stackrel{~}{𝐯}`$ $`=`$ $`{\displaystyle \frac{a𝐯t_{}}{r_{}}},`$ (6)   
$`d\stackrel{~}{t}`$ $`=`$ $`{\displaystyle \frac{dt}{a^2t_{}}},`$ (7)   
$`\stackrel{~}{\rho }`$ $`=`$ $`{\displaystyle \frac{a^3\rho }{\rho _{}}},`$ (8)   
$`\stackrel{~}{\varphi }`$ $`=`$ $`{\displaystyle \frac{a^2\varphi t_{}^2}{r_{}^2}},`$ (9)   
where

$`\rho _{}`$ $`=`$ $`\overline{\rho }_0={\displaystyle \frac{3H_0^2\mathrm{\Omega }_0}{8\pi G}},`$ (10)   
$`t_{}`$ $`=`$ $`{\displaystyle \frac{2}{H_0(\mathrm{\Omega }_0a_0^3)^{1/2}}}.`$ (11)   
In these equations, $`a(t)`$ is the Friedmann-Robertson-Walker scale factor, $`a_0`$ is its present value, and $`r_{}`$, is a free parameter whose value is chosen according to the characteristic length scale of the problem. These variables are similar to the more standard comoving variables in many respects. In particular, equations (5) and (8) imply that a volume expanding with Hubble flow remains fixed in supercomoving variables, and that the mean density inside that volume remains constant. The main difference is in the change of time variable given by equation (7). In supercomoving coordinates, the time $`\stackrel{~}{t}`$ is negative, and equal to $`\mathrm{}`$ at the big bang. In an Einstein-de Sitter model, $`\stackrel{~}{t}=1`$ at present. This change of time variable has the virtue of eliminating the cosmological drag term in the momentum equation.

In all simulations, we set $`r_{}=L_{\mathrm{box}}/a_0`$. Equation (5) then implies that the box size in supercomoving variables is unity at all times.

The time-evolution of the scale factor $`a(t)`$ is governed by the Friedmann equation. For universes composed of ordinary, nonrelativistic matter and a nonzero cosmological constant $`\lambda _0`$, the Friedmann equation takes the form

$$\left(\frac{1}{a}\frac{da}{dt}\right)^2=H(t)^2=H_0^2\left[(1\mathrm{\Omega }_0\lambda _0)\left(\frac{a}{a_0}\right)^2+\mathrm{\Omega }_0\left(\frac{a}{a_0}\right)^3+\lambda _0\right].$$
(12)   
In supercomoving variables, there is a precise normalization for the scale factor, which depends upon the particular cosmological model. For the models considered in this paper, the solution of the Friedmann equation and the present value of the scale factor are the following:

(a) Einstein-de Sitter model ($`\mathrm{\Omega }_0=1`$, $`\lambda _0=0`$)

$$a=\stackrel{~}{t}^2,a_0=1.$$
(13)   
(b) Open models ($`\mathrm{\Omega }_0<1`$, $`\lambda _0=0`$)

$$a=(\stackrel{~}{t}^{\mathrm{\hspace{0.17em}2}}1)^1,a_0=(1\mathrm{\Omega }_0)/\mathrm{\Omega }_0.$$
(14)   
(c) Flat models with nonzero cosmological constant ($`\mathrm{\Omega }_0+\lambda _0=1`$)

$$\stackrel{~}{t}=\frac{1}{2}_1^a\frac{dy}{y^{3/2}(1+y^3)^{1/2}},a_0=\left(\frac{\lambda _0}{\mathrm{\Omega }_0}\right)^{1/3}.$$
(15)   
Notice that the solutions for $`a(\stackrel{~}{t})`$ do not depend explicitly upon the cosmological parameters, which are absorbed in the definition of $`a_0`$. Hence, for all models included in the database, there are only 3 different solutions of the Friedmann equation. This is one of the most useful properties of supercomoving variables. For simplicity, we shall drop the tilde notation for supercomoving variables in the remainder of this paper, except in §3.3.

### 2.2 The Power Spectrum

For all simulations presented in this paper, we use the Cold Dark Matter (CDM) power spectrum of Bardeen et al. (1986), with the normalization of Bunn & White (1997). The power spectrum at redshift $`z`$ is given by

$$P(k,z)=2\pi ^2\left(\frac{c}{H_0}\right)^{3+n}\delta _H^2^2(z,0)k^nT_{\mathrm{CDM}}^2(k),$$
(16)   
where $`c`$ is the speed of light, $`(z,0)\delta _+(0)/\delta _+(z)`$ is the linear growth factor between redshift $`z`$ and the present, and $`\delta _+`$ is the linear growing mode (see eqs. – below), $`n`$ is the tilt, and $`T_{\mathrm{CDM}}`$ is the transfer function, given by

$$T_{\mathrm{CDM}}(q)=\frac{\mathrm{ln}(1+2.34q)}{2.34q}\left[1+3.89q+(16.1q)^2+(5.46q)^3+(6.71q)^4\right]^{1/4}$$
(17)   
(Bardeen et al. 1986), with $`q`$ is defined by

$`q`$ $`=`$ $`\left({\displaystyle \frac{k}{\mathrm{Mpc}^1}}\right)\alpha ^{1/2}(\mathrm{\Omega }_0h^2)^1\mathrm{\Theta }_{2.7}^2,`$ (18)   
$`\alpha `$ $`=`$ $`a_1^{\mathrm{\Omega }_{\mathrm{B0}}/\mathrm{\Omega }_0}a_2^{(\mathrm{\Omega }_{\mathrm{B0}}/\mathrm{\Omega }_0)^3},`$ (19)   
$`a_1`$ $`=`$ $`(46.9\mathrm{\Omega }_0h^2)^{0.670}\left[1+(32.1\mathrm{\Omega }_0h^2)^{0.532}\right],`$ (20)   
$`a_2`$ $`=`$ $`(12.0\mathrm{\Omega }_0h^2)^{0.424}\left[1+(45.0\mathrm{\Omega }_0h^2)^{0.582}\right]`$ (21)   
(Hu & Sugiyama 1996, eqs. \[D-28\] and \[E-12\]), where $`\mathrm{\Theta }_{2.7}`$ is the temperature of the cosmic microwave background in units of 2.7K, and $`\delta _H`$ is the density perturbation at horizon crossing (Liddle & Lyth 1993). Fits for $`\delta _H`$ are given by Bunn & White (1997), as follows,

$$10^5\delta _H=\{\begin{array}{cc}1.95\mathrm{\Omega }_0^{0.350.19\mathrm{ln}\mathrm{\Omega }_00.17\stackrel{~}{n}}e^{(\stackrel{~}{n}+0.14\stackrel{~}{n}^2)},\hfill & \lambda _0=0\text{;}\hfill \\ 1.94\mathrm{\Omega }_0^{0.7850.05\mathrm{ln}\mathrm{\Omega }_0}e^{(0.95\stackrel{~}{n}+0.169\stackrel{~}{n}^2)},\hfill & \lambda _0=1\mathrm{\Omega }_0\text{;}\hfill \end{array}$$
(22)   
where $`\stackrel{~}{n}n1`$.

### 2.3 Setting up Initial Conditions

We assume that the initial fluctuations originate from a Gaussian random process. The initial density contrast can then be represented as a superposition of plane waves with random phases, and amplitudes related to the power spectrum $`P(k)`$, where $`𝐤`$ is the wavenumber, and $`k=|𝐤|`$. In an infinite universe, all values of $`𝐤`$ are allowed. The power spectrum is therefore continuous, and the number of modes (that is, plane waves) present in the initial density contrast is infinite. The simulations, however, are performed inside a finite comoving cubic volume $`V_{\mathrm{box}}=L_{\mathrm{box}}^3`$ with periodic boundary conditions. This periodicity implies that only modes with wavenumbers $`𝐤=(k_x,k_y,k_z)=(n_x,n_y,n_z)k_0`$, where $`n_x`$, $`n_y`$, $`n_z`$ are integers and $`k_02\pi /L_{\mathrm{box}}`$ is the fundamental wavenumber, can be present in the simulated initial conditions. Furthermore, since the initial conditions are represented by particles, the components $`k_x`$, $`k_y`$, $`k_z`$ of the wavenumber cannot exceed the nyquist frequency $`k_{\mathrm{nyq}}=N^{1/3}k_0/2`$, where $`N`$ is the number of particles in the computational volume, and $`N^{1/3}`$ is the number of particles along one dimension. Modes with higher wavenumber cannot be represented because of undersampling. Hence we are faced with the task of representing continuous initial conditions using a discrete sample of plane waves. This key aspect of any numerical cosmological simulation is, surprisingly, seldom discussed in the literature. Here we present a detailed description.

In a periodic universe with comoving cubic volume $`V_{\mathrm{box}}=L_{\mathrm{box}}^3`$, the density contrast $`\delta `$ can be decomposed into a sum of plane waves,

$$\delta (𝐫)=\underset{𝐤}{}\delta _𝐤^{\mathrm{disc}}e^{i𝐤𝐫},$$
(23)   
where $`𝐫`$ is the comoving, or supercomoving, position, and $`\delta _𝐤^{\mathrm{disc}}`$ is the amplitude of the mode with wavenumber $`𝐤`$. The superscript “disc” stands for “discrete.” The real universe is of course not periodic, in which case all values of $`𝐤`$ are allowed. To convert equation (23) from the discrete limit to the continuous limit, consider first any function $`f(𝐤)`$ that is summed over all possible values of $`𝐤`$. In the discrete limit, we have

$$\underset{𝐤}{}f_𝐤^{\mathrm{disc}}=\underset{\mathrm{all}\mathrm{V}.\mathrm{E}.}{}f_𝐤^{\mathrm{disc}}=\frac{1}{k_0^3}\underset{\mathrm{all}\mathrm{V}.\mathrm{E}.}{}f_𝐤^{\mathrm{disc}}k_0^3=\frac{V_{\mathrm{box}}}{(2\pi )^3}\underset{\mathrm{all}\mathrm{V}.\mathrm{E}.}{}f_𝐤^{\mathrm{disc}}_{\mathrm{V}.\mathrm{E}.}d^3k,$$
(24)   
where “V.E.” represents a volume element in $`𝐤`$-space, which is a cube of volume $`k_0^3=V_{\mathrm{box}}/(2\pi )^3`$ centered at $`𝐤`$. Assuming that the function $`f`$ does not vary significantly over one volume element, we can pull it inside the integral,

$$\underset{𝐤}{}f_𝐤^{\mathrm{disc}}\frac{V_{\mathrm{box}}}{(2\pi )^3}\underset{\mathrm{all}\mathrm{V}.\mathrm{E}.}{}_{\mathrm{V}.\mathrm{E}.}f_𝐤^{\mathrm{disc}}d^3k.$$
(25)   
Of course, integrating over the volume element, and then summing over all volume elements, is effectively like integrating over all $`𝐤`$-space, so equation (25) reduces to

$$\underset{𝐤}{}f_𝐤^{\mathrm{disc}}\frac{V_{\mathrm{box}}}{(2\pi )^3}f_𝐤^{\mathrm{disc}}d^3k=f_𝐤^{\mathrm{cont}}d^3k,$$
(26)   
where the superscript “cont” stands for “continuous.” The continuous and discrete functions are related by

$$f_𝐤^{\mathrm{cont}}=\frac{V_{\mathrm{box}}}{(2\pi )^3}f_𝐤^{\mathrm{disc}}.$$
(27)   
Using these formulae, we can rewrite equation (23) as

$$\delta (𝐫)=d^3k\delta _𝐤^{\mathrm{cont}}e^{i𝐤𝐫},$$
(28)   
where

$$\delta _𝐤^{\mathrm{cont}}=\frac{V_{\mathrm{box}}}{(2\pi )^3}\delta _𝐤^{\mathrm{disc}}.$$
(29)   
To find the relationships between $`\delta _𝐤^{\mathrm{disc}}`$, $`\delta _𝐤^{\mathrm{cont}}`$, and the power spectrum, consider the rms density fluctuation $`\sigma _x`$ at some particular scale $`x`$. This quantity is given by

$$\sigma _x^2=\frac{V_{\mathrm{box}}}{(2\pi )^3}d^3k|\delta _𝐤^{\mathrm{disc}}|^2W(kx).$$
(30)   
where $`W`$ is a window function. We present the derivation of this result in Appendix A. In the continuous limit, $`\sigma _x`$ is related to the power spectrum, by

$$\sigma _x^2=\frac{1}{(2\pi )^3}d^3kP(k)W(kx)$$
(31)   
(see, e.g. Bunn & White , eqs.  and ). By combining equations (29), (30), and (31), we get

$$P(k)=V_{\mathrm{box}}|\delta _𝐤^{\mathrm{disc}}|^2=\frac{(2\pi )^6}{V_{\mathrm{box}}}|\delta _𝐤^{\mathrm{cont}}|^2.$$
(32)   
Both $`P(k)`$ and $`\delta _𝐤^{\mathrm{cont}}`$ have dimensions of a volume while $`\delta _𝐤^{\mathrm{disc}}`$ is dimensionless. Notice that the form of these expressions depends upon the actual definition of the Fourier Transform, which tends to vary among authors.

To set up initial conditions, we lay down the particles on a cubic lattice, and displace each particle by an amount $`\mathrm{\Delta }𝐫`$ given by

$$\mathrm{\Delta }𝐫=i\underset{𝐤}{}\frac{G_𝐤\delta _𝐤^{\mathrm{disc}}𝐤}{k^2}e^{i𝐤𝐫},$$
(33)   
where $`𝐫`$ is the unperturbed position, $`\delta _𝐤^{\mathrm{disc}}=|\delta _𝐤^{\mathrm{disc}}|e^{i\varphi _𝐤}`$ is a complex number with amplitude $`|\delta _𝐤^{\mathrm{disc}}|=[P(k)/V_{\mathrm{box}}]^{1/2}`$ and phase $`\varphi _𝐤`$ chosen randomly between 0 and $`2\pi `$ with uniform probability, and the sum extends over all modes included in the initial conditions (see §2.4). As Efstathiou et al. (1985) point out, assuming random phases would be sufficient to ensure that the initial conditions are Gaussian, in the continuous limit (that is, in an infinite universe). However, this assumption is insufficient in the discrete limit (that is, in a finite universe with periodic boundary conditions). To ensure the Gaussianity of the initial conditions, it is necessary, and sufficient, to include the Gaussian factor $`G_𝐤`$, a random number chosen from a Gaussian distribution with mean 0 and dispersion 1. This guarantees the initial conditions are Gaussian, even though there might be a lack of resolution at some scales; $`G_𝐤`$ does not change the spectral amplitude of the fluctuations.

To compute the initial peculiar velocity field, we assume that the initial time of the calculation is early enough for the perturbation to be in the linear regime, but late enough so that the linear decaying mode can be neglected. The initial peculiar velocity of the particles are then related to their displacements by

$$𝐯_i=\left(\frac{1}{\delta _+}\frac{d\delta _+}{dt}\right)_{z_i}\mathrm{\Delta }𝐫,$$
(34)   
where $`z_i`$ is the initial redshift of the simulations, $`\mathrm{\Delta }𝐫`$ is computed using equation (33), and $`\delta _+`$ is the linear growing mode of the perturbation, which depends upon the cosmological model. For the Einstein-de Sitter model ($`\mathrm{\Omega }_0=1`$, $`\lambda _0=0`$), the growing mode is

$$\delta _+(z)=(1+z)^1.$$
(35)   
For open models ($`\mathrm{\Omega }_0<1`$, $`\lambda _0=0`$), the growing mode is

$$\delta _+(z)=1+\frac{3}{x}+3\left(\frac{1+x}{x^3}\right)^{1/2}\mathrm{ln}\left[(1+x)^{1/2}x^{1/2}\right]$$
(36)   
(Peebles 1980), where

$$x=(\mathrm{\Omega }_0^11)(1+z)^1.$$
(37)   
Finally, for flat models with a cosmological constant ($`\mathrm{\Omega }_0+\lambda _0=1`$), the growing mode is given by

$$\delta _+(z)=\left(\frac{1}{y}+1\right)^{1/2}_0^y\frac{dw}{w^{1/6}(1+w)^{3/2}}$$
(38)   
(Martel 1991b), where

$$y=\frac{\lambda _0}{\mathrm{\Omega }_0}(1+z)^3.$$
(39)   
### 2.4 The Simulations

We set the comoving length of the computational volume $`L_{\mathrm{box}}`$ equal to $`128\mathrm{Mpc}`$ (present length units). The total mass of the system is $`M_{\mathrm{sys}}=3H_0^2\mathrm{\Omega }_0L_{\mathrm{box}}^3/8\pi G=5.821\times 10^{17}\mathrm{\Omega }_0h^2\mathrm{M}_{}`$. We use $`N=64^3=262,144`$ particles of mass $`M_{\mathrm{part}}=M_{\mathrm{sys}}/N=2.220\times 10^{12}\mathrm{\Omega }_0h^2\mathrm{M}_{}`$. We solve Poisson’s equation on a $`128^3`$ grid. In all simulations, $`ϵ`$ and $`r_e`$ were set equal to 0.3 and 2.7 mesh spacings, respectively. This corresponds, in physical units, to a comoving softening length $`ϵ=300\mathrm{kpc}`$. This is a reasonable value for gravity-only cosmological simulations. At smaller scales, hydrodynamical effects become important and cannot be ignored. The dynamical range in length of the algorithm is $`L_{\mathrm{box}}/ϵ=467`$.

The ratio of the nyquist wavenumber $`k_{\mathrm{nyq}}`$ to the fundamental wavenumber $`k_0`$ is $`N^{1/3}/2=32`$. Hence each component $`k_i`$, $`i=x`$, $`y`$, $`z`$, of the wavenumber can take 65 values; $`k_i=n_ik_0`$, with $`32n_i32`$. The initial conditions can therefore represent $`65^3=\mathrm{274\hspace{0.17em}625}`$ modes. However, the reality condition requires that the amplitudes of modes with equal and opposite wavenumbers are related by $`\delta _𝐤^{\mathrm{disc}}=(\delta _𝐤^{\mathrm{disc}})^{}`$ \[in order for $`\delta (𝐫)`$ to be real\]. Furthermore, we exclude modes with $`|𝐤|=(k_x^2+k_y^2+k_z^2)^{1/2}>k_{\mathrm{nyq}}`$. This reduces the actual number of modes represented in the initial conditions to $`\mathrm{68\hspace{0.17em}532}`$.

All simulations start at an initial redshift $`z_i=24`$. The algorithm produces “dumps” (snapshots of the system) at numerous intermediate redshifts, up to the present. These redshifts where chosen by imposing that the dumps are equally space in conformal time $`\eta `$, defined by $`d\eta a_0dt/a(t)`$. We set the difference $`\mathrm{\Delta }\eta `$ between consecutive dumps equal to $`L_{\mathrm{box}}/c`$. Thus, if $`t`$ and $`t^{}`$ are the times corresponding to 2 consecutive dumps, they are related by

$$\frac{L_{\mathrm{box}}}{c}=_t^{}^t[1+z(t)]𝑑t.$$
(40)   
This particular choice results in most dumps being concentrated near the present. Typically, about half of the dumps are between redshifts $`z=1`$ and $`z=0`$. Since the relationship between time and redshift, $`z(t)`$, is model-dependent, the redshifts where dumps are made depend upon the cosmological parameters $`\mathrm{\Omega }_0`$, $`\lambda _0`$, and $`H_0`$ (but not $`\sigma _8`$). Every simulation also produces a dump at $`z=z_i=24`$, and one at $`z=0`$. The number of dumps per simulation varies between 44 and 128.

## 3 THE DATABASE

### 3.1 The Cosmological Models

The power spectrum described in §2.2 is characterized by 6 independent parameters: (1) the density parameter $`\mathrm{\Omega }_0`$, (2) the contribution $`\mathrm{\Omega }_{\mathrm{B0}}`$ of the baryonic matter to the density parameter, (3) the cosmological constant $`\lambda _0`$, (4) the Hubble constant $`H_0`$, (5) the temperature $`T_{\mathrm{CMB}}`$ of the Cosmic Microwave Background, and (6) the tilt $`n`$ of the power spectrum. In order to keep the size of the parameter space at a manageable level, we set $`T_{\mathrm{CMB}}=2.7K`$ and $`\mathrm{\Omega }_{\mathrm{B0}}=0.015h^2`$, thus reducing the dimensionality of the parameter space to 4. Also, the normalization of the power spectrum is often described in terms of the rms density fluctuation $`\sigma _8`$ at a scale of $`8h^1\mathrm{Mpc}`$. The value of $`\sigma _8`$ is a function of the 6 aforementioned parameters. We invert this relation, treating $`\sigma _8`$ as an independent parameter, and the tilt $`n`$ as a dependent one. The independent parameters in the database are therefore $`\mathrm{\Omega }_0`$, $`\lambda _0`$, $`H_0`$, and $`\sigma _8`$. For each model, we performed up to 3 different simulations, with different realizations of the initial conditions (this amounts to choosing a different set of random numbers for the phases $`\varphi _𝐤`$ of the complex numbers $`\delta _𝐤^{\mathrm{disc}}`$, and the Gaussian factors $`G_𝐤`$).

An important question was to decide which models should be included in the database. Our goal here is not to find the “ultimate model,” which provides the best match to current observations. This would defeat the purpose of having a database, and furthermore, as former supporters of the Standard Model can appreciate, the “best model” can eventually be proven incorrect by new observations. Our intention is to provide an adequate coverage of the parameter phase-space. However, we do not want to invest much effort into simulating models that are considered “unlikely,” because some of the parameters have extreme values. With this in mind, we performed 160 simulations, which provide a broad coverage of the parameter phase-space, but we favored “likely” regions of the parameter phase-space over “unlikely” ones, by performing more simulations in these regions. For instance, we consider models with Hubble constant varying in the range $`H_0=5085\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, but 139 of the calculations (87%) have a Hubble constant in the more plausible range $`H_0=6575\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$.

The value of the parameters are given in Table 1 for the entire database (with $`H_0`$ in units of $`\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$). The first 4 columns contain the values of the 4 independent parameters $`\mathrm{\Omega }_0`$, $`\lambda _0`$, $`H_0`$, and $`\sigma _8`$. The dependent parameter $`n`$ is in the fifth column. The sixth and seventh columns contain the number of dumps per simulation and the codes of the simulations respectively (see §3.2). The parameter phase-space coverage of the database is illustrated in Figure 1. The top left panel shows a projection of the 4-dimensional parameter phase-space onto the $`\mathrm{\Omega }_0\lambda _0`$ plane. The dots indicate the cases for which there are simulations in the database. The number next to each dot indicates the number of simulations for that particular combination of $`\mathrm{\Omega }_0`$ and $`\lambda _0`$. This panel includes all simulations in the database. The top right panel shows the same projection, but for a subset of the simulations, all simulations with $`H_0=65\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. The remaining 4 panels show different projections and different subsets. As we see, the coverage of the parameter phase space is quite dense. The biggest “hole” is seen in the $`\mathrm{\Omega }_0\lambda _0`$ projection (top panels). There are currently no simulations for open models with a nonzero cosmological constant ($`\lambda _00`$ and $`\mathrm{\Omega }_0+\lambda _0<1`$) in the database. As we pointed out in §1.3, these models are certainly worth considering, and we intend to include such models in the database in the near future.

Several interesting quantities can be computed directly from the parameters. One of them is the age of the universe. For $`\lambda _00`$ models, $`t_0`$ is given by

$$t_0=\frac{1}{H_0}_0^1\left[\frac{x}{\lambda _0x^3+(1\mathrm{\Omega }_0\lambda _0)x+\mathrm{\Omega }_0x}\right]^{1/2}𝑑x$$
(41)   
(see, e.g., Martel 1990). It is of course independent of $`\sigma _8`$. Table 2 gives the ages is Gigayears for the various models included in the database.

Another interesting quantity is $`\sigma _8^{\mathrm{clus}}`$, the value of $`\sigma _8`$ inferreded from observations of clusters of galaxies. Using the X-ray temperature distribution function of clusters, Viana & Liddle (1996) have produced an empirical formula for $`\sigma _8^{\mathrm{clus}}`$,

$$\sigma _8^{\mathrm{clus}}=0.6\mathrm{\Omega }_0^{C(\mathrm{\Omega }_0)},$$
(42)   
where

$$C(\mathrm{\Omega }_0)=\{\begin{array}{cc}0.36+0.31\mathrm{\Omega }_00.28\mathrm{\Omega }_0^2,\hfill & \lambda _0=0\text{ ;}\hfill \\ 0.590.16\mathrm{\Omega }_0+0.06\mathrm{\Omega }_0^2,\hfill & \lambda _0=1\mathrm{\Omega }_0\text{ .}\hfill \end{array}$$
(43)   
Table 3 gives the values of $`\sigma _8^{\mathrm{clus}}`$ for the various models included in the database. Notice that these values do not always match the actual values of $`\sigma _8`$ used for the calculations (fourth column of Table 1).

### 3.2 Nomenclature

Each cosmological model, that is, each combination of the four parameters $`\mathrm{\Omega }_0`$, $`\lambda _0`$, $`H_0`$, and $`\sigma _8`$, is identified by a two-character code composed of an uppercase letter and a lowercase letter. For instance, the Einstein-de Sitter model with $`H_0=65\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ and $`\sigma _8=1.0`$ is identified by the code Xa. The letters were chosen for practical reasons, and the reader should not try to find some logic in these choices. Each simulation is identified by a 3-character code, composed of the two-character code for the model, plus a digit to identified the simulation. For instance, the three simulations for the Xa model are identified by the codes Xa1, Xa2, and Xa3. The codes for the entire database are given in the last column of Table 1. Each simulation produces many output files, or dumps, which are snapshots of the system at various redshift. Each file is identified by a 7-character code, which consists of the 3-character code of the simulation, an underscore, and a 3-digit number which identifies the file. For instance, the first output file created by the simulation Xa1 is called Xa1\_001, and contains a snapshot of the system at the initial redshift $`z_i=24`$. The next file created by that simulation is called Xa1\_002, and contains a snapshot at $`z=22.079`$, and so on. The last file is called Xa1\_058, and contains a snapshot at the present ($`z=0`$). The lists of redshifts where dumps are available can be obtained from the authors. There is one such list for every combination of $`\mathrm{\Omega }_0`$, $`\lambda _0`$, and $`H_0`$ included in the database.

### 3.3 Conversion to Physical Units

The positions and velocities stored in the dumps are expressed in supercomoving variables. They can be converted to physical units using equations (5)–(11). The positions $`𝐫`$ and velocities $`𝐮=H(z)𝐫+𝐯`$ in physical units are given by

$`𝐫`$ $`=`$ $`{\displaystyle \frac{L_{\mathrm{box}}\stackrel{~}{𝐫}}{1+z}},`$ (44)   
$`𝐮`$ $`=`$ $`L_{\mathrm{box}}\left[{\displaystyle \frac{H(z)\stackrel{~}{𝐫}}{1+z}}+{\displaystyle \frac{\mathrm{\Omega }_0^{1/2}H_0(1+z)\stackrel{~}{𝐯}}{2a_0^{1/2}}}\right].`$ (45)   
In these expressions, we have reintroduced the tilde notation for the supercomoving variables. The expression of $`𝐫`$ is the same for all models, but the one for $`𝐯`$ is model dependent. After eliminating $`H(z)`$ using equation (12) and $`a_0`$ using equations (13)–(15), we obtain the following expressions:

(a) Einstein-de Sitter model ($`\mathrm{\Omega }_0=1`$, $`\lambda _0=0`$)

$$𝐮=H_0L_{\mathrm{box}}\left[(1+z)^{1/2}\stackrel{~}{𝐫}+\frac{(1+z)\stackrel{~}{𝐯}}{2}\right].$$
(46)   
(b) Open models ($`\mathrm{\Omega }_0<1`$, $`\lambda _0=0`$)

$$𝐮=H_0L_{\mathrm{box}}\left[(1+\mathrm{\Omega }_0z)^{1/2}\stackrel{~}{𝐫}+\frac{\mathrm{\Omega }_0(1+z)\stackrel{~}{𝐯}}{2(1\mathrm{\Omega }_0)^{1/2}}\right].$$
(47)   
(c) Flat models with nonzero cosmological constant ($`\mathrm{\Omega }_0+\lambda _0=1`$)

$$𝐮=H_0L_{\mathrm{box}}\left\{\left[\mathrm{\Omega }_0(1+z)^3+\lambda _0\right]^{1/2}\frac{\stackrel{~}{𝐫}}{1+z}+\frac{\mathrm{\Omega }_0^{2/3}(1+z)\stackrel{~}{𝐯}}{2\lambda _0^{1/6}}\right\}.$$
(48)   
### 3.4 Technical Considerations

The database contains 160 simulations for 68 cosmological models. For each simulation, there is a dump at $`z=z_i=24`$ and one at $`z=0`$, plus numerous dumps at intermediate redshifts. There is a total of $`\mathrm{11\hspace{0.17em}973}`$ dumps in the database. Each dump contains $`6N=\mathrm{1\hspace{0.17em}572\hspace{0.17em}864}`$ numbers, the coordinates of the position and velocity for each particle in the simulation. These numbers are stored in single precision (32 bits), although the simulations themselves were performed in double precision. Each dump is a binary file (IEEE 754 standard) of size $`6.3\mathrm{Mb}`$, which contains first the $`x`$-coordinates of all particles, followed by the $`y`$-coordinates, the $`z`$-coordinates, the $`v_x`$-coordinates, the $`v_y`$-coordinates, and finally the $`v_z`$-coordinates. Figure 2 shows a sample FORTRAN program that reads a file from the database. The size of the entire database is $`75.3\mathrm{Gb}`$. The database currently resides on archival tapes at the High Performance Computing Facility, University of Texas, where the simulations were performed.

Because of the size of the database, it would be impractical (if not impossible) to install it on a web site or an anonymous ftp site where it could be easily retrieved by the user. This might change in the future, but currently the only way to access the database is to contact the authors, preferably by E-mail, and send a list of the dumps requested. Then, the authors and the user can choose the best strategy for transferring file, according to the computer resources and needs of the user. Requests should be sent to database@galileo.as.utexas.edu.

## 4 ANALYSIS OF THE SIMULATIONS

The goal of this paper is to present the database, and describe its content. We supplement this description by analyzing the final state of each simulation, which corresponds to the present. We focus on four particular aspects of the present large-scale structure: the rms density fluctuation, the two-point correlation function, the moments of the peculiar velocity field, and the properties of clusters. In this abbreviated version of the paper, we only present the analysis of the rms density fluctuation and the two-point correlation function. The full version of the paper can be obtained by contacting the authors.

### 4.1 RMS Density Fluctuation

The present rms density fluctuation $`\sigma _8`$ at scale $`8h^1\mathrm{Mpc}`$ is treated as an independent parameter in our simulations. However, we do not set up the state of the system at present. Instead, we set up initial conditions at high redshift ($`z_i=24`$), and evolve the system numerically all the way to the present. We adjust the initial conditions in such a way that the density fluctuation at presents ends up being equal to the desired value of $`\sigma _8`$. To achieve this, we assume that the power spectrum evolves with time according to linear perturbation theory (hence the presence of the factor $`^2`$ in equation ). Actually, we do not expect the actual value of $`\sigma _8`$ to be precisely equal to the desired one, for several reasons. Let us designate by $`\sigma _8^{\mathrm{cont}}`$ the desired value of $`\sigma _8`$ for each simulation, that is, the quantity appearing in the fourth column of Table 1. The superscript “cont” indicates that this the value in the real universe, where the wavenumber $`𝐤`$ varies continuously. We are representing the initial conditions using a finite number of discrete modes, which is clearly an approximation. We designate by $`\sigma _8^{\mathrm{disc}}`$ the value of $`\sigma _8`$ resulting from this approximation. This value is given by

$$(\sigma _8^{\mathrm{disc}})^2=^2(z_i,0)\underset{𝐤}{}|\delta _{𝐤,i}^{\mathrm{disc}}|^2W(k\mathrm{}),$$
(49)   
where $`\mathrm{}8h^1\mathrm{Mpc}`$. Hence, $`\sigma _8^{\mathrm{disc}}`$ is computed by summing over all modes present in the initial conditions, ignoring the Gaussian factor, and then extrapolating to the present using linear perturbation theory.

The introduction of the Gaussian factor in equation (33) further modifies the value of $`\sigma _8`$. We designate by $`\sigma _8^{\mathrm{gauss}}`$ the value of $`\sigma _8`$ resulting from the presence of this factor,

$$(\sigma _8^{\mathrm{gauss}})^2=^2(z_i,0)\underset{𝐤}{}G_𝐤^2|\delta _{𝐤,i}^{\mathrm{disc}}|^2W(k\mathrm{}).$$
(50)   
Finally, we designate by $`\sigma _8^{\mathrm{num}}`$ the “numerical” value of $`\sigma _8`$, which is the actual rms density fluctuation inside the computational volume at present, obtained from the numerical simulation. This value should differ from $`\sigma _8^{\mathrm{gauss}}`$ from several reasons. First, the numerical algorithm has a finite accuracy, owing to the fact that the time step is finite and that the gravitational force is softened at short distances. Second, the evolution of a single mode would never follow precisely the exact solution when the system is represented by a finite number of particles. Third, and more importantly, equations (49) and (50) use linear perturbation theory to extrapolate from the initial conditions to the present, but this can only be approximate, as mode coupling introduces nonlinear effects at small scale.

We investigated the importance of these various effects by computing the various values of $`\sigma _8`$ for all simulations. The value of $`\sigma _8^{\mathrm{cont}}`$ is imposed. The values of $`\sigma _8^{\mathrm{disc}}`$ and $`\sigma _8^{\mathrm{gauss}}`$ can be computed directly using equations (49) and (50) (these values are provided automatically by the code that generates initial conditions). To evaluate $`\sigma _8^{\mathrm{num}}`$, we used a direct, somewhat brute-force approach. For each simulation, we located one million (!) spheres of radius $`8h^1\mathrm{Mpc}`$ at random locations inside the computational volume at present, and computed the density contrast $`\delta _{\mathrm{sph}}`$ inside each sphere, using

$$\delta _{\mathrm{sph}}=\frac{N_{\mathrm{sph}}\overline{N}_{\mathrm{sph}}}{\overline{N}_{\mathrm{sph}}},$$
(51)   
where $`N_{\mathrm{sph}}`$ is the number of particles inside the sphere, and $`\overline{N}_{\mathrm{sph}}`$ is the “mean number,” given by

$$\overline{N}_{\mathrm{sph}}=\frac{4\pi \mathrm{}^3N}{3V_{\mathrm{box}}}$$
(52)   
(notice that $`\overline{N}_{\mathrm{sph}}`$ is not an integer). The value of $`\sigma _8^{\mathrm{num}}`$ is then given by

$$(\sigma _8^{\mathrm{num}})^2=10^6\underset{\mathrm{all}\mathrm{spheres}}{}(\delta _{\mathrm{sph}})^2.$$
(53)   
We plot these various values of $`\sigma _8`$ against each others in Figure 3. The top-left panel shows the effect of discreteness. All dots are located below the dashed line, indicating that $`\sigma _8^{\mathrm{disc}}<\sigma _8^{\mathrm{cont}}`$. This is caused not as much by the discreteness itself as by the fact that modes outside the range $`k_0kk_{\mathrm{nyq}}`$ are missing in the initial conditions. Still, the effect is very small, 7% in the worst case (which happens to be model Uc). As we see on the top right panel, the effect of introducing of the Gaussian factor $`G_𝐤`$ is quite important, and causes a spread of order 10% in the value of $`\sigma _8`$. This is primarily an effect of undersampling. The modes with wavenumbers comparable to $`k_0`$ are very few, but contribute significantly to $`\sigma _8`$. This is consequence of the fact that our computational volume is actually too small to constitute a “fair sample” of the universe. With a larger volume, there would still be very few modes with wavenumbers of order $`k_0`$, but these modes would be farther from the peak of the power spectrum, and therefore would contribute less to $`\sigma _8`$. The bottom left panel shows the effect of actually performing the simulation. The values of $`\sigma _8^{\mathrm{gauss}}`$ and $`\sigma _8^{\mathrm{num}}`$ are comparable for small $`\sigma _8`$, but differ at large $`\sigma _8`$ where nonlinear effects become important. This panel shows that (1) the onset of nonlinearity occurs at $`\sigma _80.6`$, (2) The effect of nonlinearity on the value of $`\sigma _8`$ is small, of order 10%, and (3) the effect can go either way, it can either increase or decrease the value of $`\sigma _8`$, and the occurrences of these two cases are comparable.

The bottom right panel shows the combined effect of discreteness, Gaussian factor, and nonlinearity. The spread is quite large. The value $`\sigma _8^{\mathrm{num}}`$ that comes out of the simulation can differ by as much as 20% from the value $`\sigma _8^{\mathrm{cont}}`$ we intended to obtain, that is, the values given in the fourth column of Table 1. The reader should be aware of this fact when selecting a particular model from the database. Consequently, we listed in Table 4 the values of $`\sigma _8^{\mathrm{num}}`$ for all simulations.

### 4.2 2-point Correlation Function

The distribution of galaxies in the universe can be described statistically using N-point correlation functions. The first and most important of these functions is the 2-point correlation function $`\xi (r)`$, which measures the excess probability of finding two galaxies separated by a distance $`r`$. The 2-point correlation function can be estimated from the particle distributions by locating spherical shells around particles and counting the number of particles inside these shells. If $`N(r_1,r_2)`$ is the number of particles inside a spherical shell of inner radius $`r_1`$ and outer radius $`r_2`$ centered on a given particles, and $`N(r_1,r_2)`$ is the average of $`N(r_1,r_2)`$ over all particles, then by definition

$$N(r_1,r_2)=\frac{4\pi (r_2^3r_1^3)n}{3}+4\pi n_{r_1}^{r_2}\xi (r)r^2𝑑r.$$
(54)   
where $`n`$ is the number density of particles. The first term in equation (54) gives the correct answer in the case of a uniform distribution. The second term represents the effect of the correlation. We can solve this equation for $`\xi `$ (see Martel 1991a). After some algebra, we get

$$\xi (x)=\frac{1}{r^3(x)\mathrm{ln}10}\frac{d}{dx}\left[\frac{N(r(x),r_1)}{4\pi n}\frac{r^3(x)r_1^3}{3}\right],$$
(55)   
where $`x\mathrm{log}_{10}r`$. To compute $`\xi `$, we first evaluate $`N(r(x),r_1)`$ by computing the spacing between all pairs of particles, and counting pairs in bins equally spaced in $`x`$. We then compute the derivative in equation (55) using a standard five-point finite difference operator. We computed $`\xi (r)`$ at present for all simulations in the database. The results are plotted in Figures 4–7. For models with more than one simulation (most have three), we averaged the curves. They were actually so similar, in all cases, that error bars in Figures 4–7 would be too small to be seen. This shows that the 2-point correlation function depends essentially upon the cosmological model, with very little dependence upon the particular realization of the initial conditions. Each panel in Figures 4–7 corresponds to a particular combination of $`\mathrm{\Omega }_0`$, $`\lambda _0`$, and $`H_0`$, with different values of $`\sigma _8`$ represented by different curves. The dashed lines show the observed galaxy 2-point correlation function,

$$\xi (r)=\left(\frac{r}{5.4h^1\mathrm{Mpc}}\right)^{1.77}$$
(56)   
(Peebles 1993, eq. \[7.32\]). Models with $`\sigma _8=0.3`$ fail to reproduce the observed correlation function on three counts. The slope of $`\xi (r)`$ is too shallow, the amplitude is too small, and there tend to be a “kink” in the correlation function at separations of order $`13\mathrm{Mpc}`$. As $`\sigma _8`$ increases, the kink goes away, and the amplitude and slope increase. Models with $`\sigma _8=0.8`$ provide the best fit to the observations, almost independently of the values of the other parameters! For larger values of $`\sigma _8`$, the slope and amplitude are too large.

All curves have a shoulder at small separations, where the slope drops significantly. Martel (1991a) argues that this is a consequence of the softening of the force at small distances. The effect of this softening is to “take” pairs of particles that would have a separation $`r`$ less than the softening length $`ϵ=300\mathrm{kpc}`$ in the absence of softening, and transfer them to separations $`rϵ`$, resulting in a flattening of the correlation function. To illustrate this, we indicated in all panels of Figures 4–7 the location of the softening length $`ϵ`$ by a thick line. This line is located right in the middle of the “shoulder” for all curves. With higher force resolution, pairs that are now located at separations $`rϵ`$ would be located instead at separations $`rϵ`$, and the fit to the observed slope would be improved.

## 5 SUMMARY AND PROSPECTS

Using a P<sup>3</sup>M algorithm with $`64^3`$ particles, we have performed 160 cosmological simulations, for 68 cosmological models. This constitutes the largest database of cosmological simulations ever assembled. We covered a four-dimensional parameter phase space by varying the density parameter $`\mathrm{\Omega }_0`$, the cosmological constant $`\lambda _0`$, the Hubble constant $`H_0`$, and the rms density fluctuation $`\sigma _8`$. We are making this database available to the astronomical community. We also performed a limited analysis of the simulations. Our results are the following:

(1) The present rms density fluctuation $`\sigma _8^{\mathrm{num}}`$ differs from the one expected by linearly extrapolating the initial power spectrum to the present, because of the combined effects of having a finite number of modes in the initial conditions, introducing a Gaussian factor in the initial conditions, and having nonlinear coupling between modes. The first effect is negligible. The second effect is of order 10% or less, but this probably depends upon the size of the computational volume \[the one used for the simulations, $`(128\mathrm{Mpc})^3`$, is a little too small to constitute a fair sample of the universe\]. The effect of nonlinearity is negligible for $`\sigma _8<0.6`$, and of order 10% or less for larger $`\sigma _8`$. This can go either way: nonlinearities can either increase of decrease $`\sigma _8`$ relative to what linear theory predicts.

(2) The observed two-point correlation function $`\xi (r)`$ is well reproduced by models with $`\sigma _8=0.8`$, nearly independently of the values of the parameters $`\mathrm{\Omega }_0`$, $`\lambda _0`$, and $`H_0`$. For models with $`\sigma _8<0.8`$, the correlation function is too small, its slope is too shallow, and it has a kink at separations $`r=13\mathrm{Mpc}`$. For models with $`\sigma _8>0.8`$, the correlation function is too large and its slope is too steep.

(3) At small separations, $`r<1\mathrm{Mpc}`$, the velocity moments satisfy the relations $`|V_\mathrm{R}|H_0r`$ and $`V_{\mathrm{PP}}2^{1/2}V_{\mathrm{PL}}`$, indicating that small clusters have reached virial equilibrium. At larger separations, $`|V_\mathrm{R}|`$ increases above the Hubble velocity, indicating that clusters are accreting matter from the field. The velocity moments depend essentially upon $`\mathrm{\Omega }_0`$ and $`\sigma _8`$, and not $`\lambda _0`$ and $`H_0`$. The pairwise particle velocity dispersions are much larger than the observed pairwise galaxy velocity dispersion, except for models with $`\mathrm{\Omega }_0=0.2`$ and $`\sigma _80.4`$. But if the velocity dispersion of galaxies is biased relative to the velocity dispersion of dark matter, then models with larger values of $`\mathrm{\Omega }_0`$ or $`\sigma _8`$ can be reconciled with observations.

(4) The multiplicity functions are decreasing for small values of for models with $`\sigma _80.3`$. At larger values of $`\sigma _8`$, the multiplicity functions have a horizontal plateau, whose length increases with $`\sigma _8`$. For models with $`\sigma _8>0.9`$, the multiplicity functions have a $``$ shape which results from the merging of intermediate-size clusters. For all models, clusters have densities in the range $`100\overline{\rho }_01000\overline{\rho }_0`$. A simple analytical model suggest that clusters have a density $`\rho 178\overline{\rho }`$ when they reach virial equilibrium. Our results suggest that many clusters have reached that equilibrium in the past, when $`\overline{\rho }`$ was larger than $`\overline{\rho }_0`$ (this could be checked by performing a cluster analysis on earlier dumps). The spin parameters $`\lambda `$ are in the range $`0.0080.2`$, with the median near 0.05, and the distributions of elongations favors prolate shapes ($`e_2>e_1`$) over oblate shapes ($`e_1<e_2`$). These results indicate the absence of rotationally supported disks in these simulations.

The database is growing. We are currently adding new simulations to the original 160 simulations described in this paper. There are at least seven different motivations for performing additional simulations.

(1) Additional Simulations for the same Models: For the sake of providing a good coverage of the parameter phase-space, we have limited the number of simulations per model to 3 or less. We can perform additional simulations for models already included in the database, if there is a need for doing so. This could be the case if, for some reason, a particular model (that is, a particular combination of $`\mathrm{\Omega }_0`$, $`\lambda _0`$, $`H_0`$, and $`\sigma _8`$) becomes particularly interesting, and deserves more scrutiny. Also, having more simulations per model has the virtue of improving the statistics. For instance, the size of the error bars in Figures 8–12 would be reduced if we had more than 3 simulations per model. Finally, for gravitational lensing simulations, it is necessary to combine dumps generated by different simulations, and having more simulations can be desirable (see, e.g. Premadi, Martel, & Matzner 1998).

(2) Different Box Size: All simulations included in the database were performed using a computational box of size $`L_{\mathrm{box}}=128\mathrm{Mpc}`$. As described in §2.4, the softening length is comparable to the scale where nongravitational effects become important. Hence, there is no reason to consider smaller box sizes, unless we want to use fewer particles. There are, however, reasons for considering larger boxes. As we pointed out in §4.1, a box of size $`128\mathrm{Mpc}`$ is too small to constitute a “fair sample” of the universe. Using boxes of size $`256\mathrm{Mpc}`$ or even $`512\mathrm{Mpc}`$ would certainly provide a better, “fairer” description of the large-scale structure, even with the same number of particles.

(3) Larger Number of Particles: The is no point increasing the number of particles as long as we keep the box size at $`128\mathrm{Mpc}`$, since the resolution would be increased at scales where nongravitational effects are important. However, if larger boxes are used, the number of particles can be increased accordingly in order to maintain the resolution of the algorithm at small scale. If we continue to adopt $`300\mathrm{kpc}`$ as the resolution scale of the algorithm, simulations in $`(256\mathrm{Mpc})^3`$ and $`(512\mathrm{Mpc})^3`$ boxes could be performed with $`128^3`$ and $`256^3`$ particles, respectively.

(4) New Background Models: The 4-dimensional parameter phase-space considered in this paper is quite large, and the set of 68 cosmological models included in the database covers a small fraction of it. There are several “holes” in the projections shown in Figure 1. In particular, there are no simulations for open models with a nonzero cosmological constant ($`\lambda _00`$, $`\mathrm{\Omega }_0+\lambda _0<1`$). Simulations for additional background models could be added to the database, either to provide a better coverage of the parameter phase-space, or because there is a particular model we are interested in, “we” designating either the authors, or other researchers sending us a special request. Actually, the original database contained only 151 simulations for 65 cosmological models. Following a special request by Hamana (1998), we added 9 simulations to the database, for 3 new models: Ea, Pa, and Xg.

(5) Additional Parameters: The current database covers a 4-parameter phase space, because we held the CMB temperature $`T_{\mathrm{CMB}}`$ and the baryon density parameter $`\mathrm{\Omega }_{\mathrm{B0}}`$ at values of 2.7 and $`0.015h^2`$, respectively. The CMB temperature is known so accurately that treating it as a variable parameter would be pointless. This is not the case for the baryon density parameter. According to primordial nucleosynthesis, the quantity $`\mathrm{\Omega }_{\mathrm{B0}}h^2`$ has an allowed range from 0.01 to 0.026 (Krauss & Kernan 1995; Copi et al. 1995; Krauss 1998). Furthermore, X-ray observations of clusters of galaxies suggest that the ratio of gas mass to dark matter mass in these clusters exceeds the mean value in the universe, a phenomenon known as “the baryon catastrophe” (Briel et al. 1992; White et al. 1993; Martel et al. 1994). There is at present no definitive explanation for this phenomenon, but one possible explanation is that primordial nucleosynthesis is somehow incorrect, and predicts a value of $`\mathrm{\Omega }_{\mathrm{B0}}`$ which is too small.

(6) Different Components: All simulations in the database used a CDM power spectrum as initial conditions. There are, however, several other models that constitute interesting alternatives to the CDM model, which could be added to the database. One of them is the Hot Dark Matter model (HDM), though this model has fallen out of favor in recent years, due to its inability to form galaxies inside deep voids such as Boötes. A more interesting alternative is the mixed Cold + Hot Dark Matter model (CHDM), which contains both a cold dark matter component and a massive neutrino component. This model introduces one additional parameter, the contribution $`\mathrm{\Omega }_{\nu 0}`$ of the neutrinos to the mean energy density of the universe.

(7) Different Cosmologies: The cosmological models included in the database contain only non-relativistic matter and a nonzero cosmological constant. It would be very interesting to consider models with other components. Possible candidates include domain walls, cosmic strings, or relativistic particles (Fry 1985; Charlton & Turner 1987; Silveira & Waga 1994; Martel 1995; Martel & Shapiro 1998). Recently, these various candidates have been combined into a single concept called “quintessence” (Caldwell, Dave, & Steinhardt 1998). The effects of these various components is twofold: First, the presence of these components modifies the expansion rate of the universe and the growth rate of density perturbations, thus changing the history of large-scale structure formation. Second, they might affect the shape and normalization of the primordial power spectrum, in ways that remain to be determined (none of these models were considered by Bunn & White ).

###### Acknowledgements.

This work benefited from stimulating discussions with Paul Shapiro. We are pleased to acknowledge the support of NASA Grants NAG5-2785, NAG5-7363, and NAG5-7821, NSF Grants PHY93 10083, PHY98 00725 and ASC 9504046, the University of Texas High Performance Computing Facility through the office of the vice president for research. HM acknowledges the support of a fellowship provided by the Texas Institute for Computational and Applied Mathematics.    

## Appendix A Calculation of $`\sigma _x^{\mathrm{disc}}`$

The mass inside a sphere centered at $`𝐫_0`$ is given by

$$M(𝐫_0)=_{\mathrm{sph}(𝐫_0)}\overline{\rho }_{\mathrm{com}}(1+\delta )d^3r=\overline{\rho }_{\mathrm{com}}\left[V_{\mathrm{sph}}+_{\mathrm{sph}(𝐫_0)}d^3r\underset{𝐤}{}\delta _𝐤^{\mathrm{disc}}e^{i𝐤𝐫}\right],$$
(A1)   
where $`\overline{\rho }_{\mathrm{com}}`$ is the average comoving density, $`V_{\mathrm{sph}}`$ is the volume of the sphere, and the integral is computed over that volume. The relative mass excess in the sphere is given by

$$\frac{\mathrm{\Delta }M}{M}(𝐫_0)=\frac{1}{V_{\mathrm{sph}}}_{\mathrm{sph}(𝐫_0)}d^3r\underset{𝐤}{}\delta _𝐤^{\mathrm{disc}}e^{i𝐤𝐫}.$$
(A2)   
We introduce the following change of variables,

$$𝐫=𝐫_0+𝐲.$$
(A3)   
In $`𝐲`$-space, the sphere is now located at the origin, and equation (A3) becomes

$$\frac{\mathrm{\Delta }M}{M}(𝐫_0)=\frac{1}{V_{\mathrm{sph}}}_{\mathrm{sph}(0)}d^3y\underset{𝐤}{}\delta _𝐤^{\mathrm{disc}}e^{i𝐤𝐫_0}e^{i𝐤𝐲}.$$
(A4)   
We now square this expression, and get

$$\left(\frac{\mathrm{\Delta }M}{M}\right)^2(𝐫_0)=\frac{9}{16\pi ^2x^6}\left[_{\mathrm{sph}(0)}d^3y\underset{𝐤}{}\delta _𝐤^{\mathrm{disc}}e^{i𝐤𝐫_0}e^{i𝐤𝐲}\right]\left[_{\mathrm{sph}(0)}d^3z\underset{𝐤^{}}{}\delta _𝐤^{}^{\mathrm{disc}}e^{i𝐤^{}𝐫_0}e^{i𝐤^{}𝐳}\right],$$
(A5)   
where $`x`$ is the radius of the sphere. The rms density contrast at scale $`x`$ is obtained by averaging the above expression over all possible locations of the sphere inside the computational box,

$`\sigma _x^2`$ $``$ $`\left({\displaystyle \frac{\mathrm{\Delta }M}{M}}\right)^2_{V_{\mathrm{box}}}={\displaystyle \frac{1}{V_{\mathrm{box}}}}{\displaystyle _{V_{\mathrm{box}}}}d^3r_0\left({\displaystyle \frac{\mathrm{\Delta }M}{M}}\right)^2(𝐫_0)`$ (A6)   
$`=`$ $`{\displaystyle \frac{1}{V_{\mathrm{box}}}}{\displaystyle \frac{9}{16\pi ^2x^6}}{\displaystyle _{V_{\mathrm{box}}}}d^3r_0{\displaystyle _{\mathrm{sph}(0)}}d^3y{\displaystyle _{\mathrm{sph}(0)}}d^3z{\displaystyle \underset{𝐤}{}}{\displaystyle \underset{𝐤^{}}{}}\delta _𝐤^{\mathrm{disc}}\delta _𝐤^{}^{\mathrm{disc}}e^{i𝐤𝐲}e^{i𝐤^{}𝐳}e^{i(𝐤+𝐤^{})𝐫_0}.`$   
The integral over $`V_{\mathrm{box}}`$ reduces to

$$_{V_{\mathrm{box}}}d^3r_0e^{i(𝐤+𝐤^{})𝐫_0}=V_{\mathrm{box}}\delta _{𝐤,𝐤^{}}.$$
(A7)   
We substitute this expression in equation (A6), and use the Kronecker $`\delta `$ to eliminate the summation over $`𝐤^{}`$. Equation (A6) reduces to

$$\sigma _x^2=\frac{9}{16\pi ^2x^6}\underset{𝐤}{}|\delta _𝐤^{\mathrm{disc}}|^2\left[_{\mathrm{sph}(0)}d^3ye^{i𝐤𝐲}\right]^2.$$
(A8)   
The remaining integral can be evaluated easily. Equation (A8) reduces to

$$\sigma _x^2=\underset{𝐤}{}|\delta _𝐤^{\mathrm{disc}}|^2W(kx),$$
(A9)   
where

$$W(y)\frac{9}{y^6}(\mathrm{sin}yy\mathrm{cos}y)^2.$$
(A10)   
Using equation (26), we can rewrite this expression in an integral form,

$$\sigma _x^2=\frac{V_{\mathrm{box}}}{(2\pi )^3}d^3k|\delta _𝐤^{\mathrm{disc}}|^2W(kx).$$
(A11)   
Figure Captions