# Nonlinear Cosmological Structure Formation

## 1 Introduction

The power spectrum $`P(k)`$ of matter fluctuations is a fundamental quantity in cosmology. It provides the most basic statistical measure of gravitational clustering, and an accurate determination of the power spectrum is among the most important goals of every major galaxy survey. Since each cosmological model has its distinct power spectrum, one can hope to obtain crucial information such as the matter content of the Universe and the nature of the primordial fluctuations by comparing the observed power spectrum with theoretical predictions for various plausible models.

Measurements of the power spectrum require extensive galaxy surveys covering a large volume of space. It is only recently that such surveys have been completed and are starting to reveal intriguing results. For instance, much interest has been generated by the compilation of the power spectrum determined from the APM and IRAS galaxy catalogs and the effort to reconstruct the underlying linear mass power spectrum (Peacock & Dodds 1994; Peacock 1997). The recently completed Las Campanas survey of over 23,000 galaxies has yielded a power spectrum incompatible with the standard cold dark matter (CDM) model but still consistent with a class of open CDM, cold+hot dark matter (C+HDM), and tilted CDM models (Lin et al. 1996). Larger ongoing projects such as the CfA Century, the 2dF, and the Sloan surveys promise to provide measurements with reduced error bars and data points on larger scales for better constraints.

On the theoretical side, calculations of the power spectrum falls into two regimes: the linear and the nonlinear, which are characterized by the amplitude of the density fluctuations $`\delta =\delta \rho /\rho 1`$ and $`>1`$, respectively. Since the observed power spectrum spans both linear and nonlinear regimes, theoretical calculations in both regimes must be performed before all data points can be fully utilized as constraints. In the linear regime, the computation of the power spectrum has become standard practice. It is obtained by integrating the coupled, linearized Einstein and Boltzmann equations that describe the evolution of the metric and density perturbations throughout the history of the Universe. Several groups have written numerical codes for such calculations, and a publicly-available version is described in Ma & Bertschinger (1995) and placed at http://arcturus.mit.edu. This code computes, in both synchronous and conformal Newtonian gauges, the evolution of the phase-space distributions of photons, baryons, cold dark matter, and both massless and massive neutrinos. It is therefore applicable to most CDM, C+HDM, and CDM with a cosmological constant (LCDM) models. To gain deeper insight and for convenience, analytical approximations for the linear power spectra in these models have also been published. Details about the linear power spectrum are discussed in Section 2.

Gravitational clustering, however, eventually becomes a nonlinear process. The existence of galaxies and clusters is a manifestation of the nonlinear nature of gravity. Determining the nonlinear power spectrum $`P_{\mathrm{nl}}`$ is therefore an important task. Not surprisingly, $`P_{\mathrm{nl}}`$ is more difficult to obtain than the linear power spectrum. Higher-order perturbation theories help to extend the range of validity to the quasi-linear regime, but the fully nonlinear power spectrum can be calculated only from numerical simulations. Thus far, theoretical predictions of $`P_{\mathrm{nl}}`$ have been carried out for only a few models at limited epochs because it is difficult and time-consuming to perform numerical simulations with sufficient dynamic range to allow calculations of $`P_{\mathrm{nl}}`$ over a wide range of scales. Section 3 summarizes some recent work and presents new analytic and numerical results for nonlinear clustering in C+HDM models.

## 2 Linear Power Spectrum

The power spectrum $`P(k)`$ of matter fluctuations is related to the density field $`\delta `$ in $`k`$-space by

$$\delta (\stackrel{}{k})\delta ^{}(\stackrel{}{k}^{})=P(k)\delta _D(\stackrel{}{k}\stackrel{}{k}^{}),$$
(1)   
where $`\delta _D`$ is the Dirac-delta function. The Fourier transform of the power spectrum is the two-point correlation function $`\xi (r)=\delta (\stackrel{}{x_1})\delta (\stackrel{}{x_2})`$, where $`r|\stackrel{}{x_1}\stackrel{}{x_2}|`$. For a Gaussian density field (as predicted in most inflationary theories), its statistical property is entirely determined by $`P(k)`$.

The shape of the linear power spectrum depends on the cosmological parameters assumed in a given theory of structure formation. It is governed by the linear perturbation theory of gravitational clustering, and can be computed by integrating the coupled, linearized Einstein, Boltzmann, and fluid equations for the metric and density perturbations. Some important parameters that affect the shape of the power spectrum are: (1) Primordial spectral index $`n`$: $`P(k)k^n`$. An example is the Harrison-Zeldovich (Harrison 1970; Zeldovich 1972) spectrum which takes $`n=1`$. Most versions of inflationary models also predict nearly $`n=1`$ power spectra. (2) Matter-radiation equality time: The equality time $`t_{\mathrm{eq}}`$ is defined to be the epoch in the thermal history of the Universe when the energy density in matter equals that in radiation. The equality redshift scales as $`1+z_{\mathrm{eq}}\mathrm{\Omega }_\mathrm{m}h^2`$, where $`\mathrm{\Omega }_\mathrm{m}`$ is the density parameter in matter and $`h`$ is the Hubble constant in units of 100 km/s/Mpc. This parameter controls the location of the peak of $`P(k)`$ in CDM-type models since the density perturbations with wavenumbers $`k>k_{\mathrm{eq}}`$ enter the horizon in the radiation-dominated era and cannot grow appreciably, whereas perturbations with $`k<k_{\mathrm{eq}}`$ enter the horizon in the matter-dominated era and can grow unimpeded. (3) Nature of dark matter. The pure CDM model, for example, exhibits a characteristic $`P(k)k^3`$ slope at high $`k`$, while the power spectrum for the pure hot dark matter (HDM) model is cut off exponentially below the free-streaming scale of the massive neutrinos due to the phase mixing in the neutrino phase-space distribution. The hybrid C+HDM models, which are parameterized by the neutrino mass density $`\mathrm{\Omega }_\nu `$ and CDM mass density $`\mathrm{\Omega }_c=1\mathrm{\Omega }_\nu \mathrm{\Omega }_b`$ (where $`\mathrm{\Omega }_b`$ is the mass density in baryons), exhibit intermediate behavior. The suppression in the power in both the cold and hot components due to neutrino free-streaming generally increases with increasing $`\mathrm{\Omega }_\nu `$.

The solid curves in Figure 1 show the linear power spectra for the pure CDM and three C+HDM models with $`\mathrm{\Omega }_\nu =0.1,0.2`$, and 0.3. The corresponding neutrino masses in the four models are 0, 2.3, 4.6, and 7 eV, respectively. (Only one of the three species of neutrinos is assumed massive.) The models have a total matter density of $`\mathrm{\Omega }_m=1`$ and $`h=0.5`$, and all are normalized to the COBE rms quadrupole $`Q_{\mathrm{rms}\mathrm{PS}}=18\mu K`$ (Gorski et al. 1996), as evidenced by the convergence of the curves at small $`k`$. The differences at large $`k`$ reflect the distinct properties of cold and hot dark matter: the larger $`\mathrm{\Omega }_\nu `$ is, the weaker is the gravitational clustering on small length scales. For comparison, a COBE-normalized low-density model with $`\mathrm{\Omega }_m=0.3`$, a cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, and $`h=0.75`$ is also shown (dashed curve).

Analytical approximations with less than $`10`$% error are very convenient as an input for many calculations. The pure CDM model is well-approximated by

$`P_c(q,a,\mathrm{\Omega }_\nu =0)`$ $`=`$ $`Aa^2k^n`$ (2)   
$`\times `$ $`\left[{\displaystyle \frac{\mathrm{ln}(1+\alpha _1q)}{\alpha _1q}}\right]^2{\displaystyle \frac{1}{[1+\alpha _2q+(\alpha _3q)^2+(\alpha _4q)^3+(\alpha _5q)^4]^{1/2}}},q={\displaystyle \frac{k}{\mathrm{\Gamma }h}},`$   
where $`A`$ is the overall normalization factor, $`a`$ is the expansion factor, $`n`$ is the primordial spectral index, $`k`$ is the wavenumber in units of Mpc<sup>-1</sup>, and $`\alpha _1=2.34,\alpha _2=3.89,\alpha _3=16.1,\alpha _4=5.46`$, and $`\alpha _5=6.71`$ (Bardeen et al. 1986). The normalization $`A`$ can be determined from the 4-year data from the COBE satellite, and $`A=2420h^4`$ Mpc<sup>4</sup> for $`n=1`$ flat models (Bunn & White 1997). The shape parameter $`\mathrm{\Gamma }`$ characterizes the dependence on cosmological parameters, and a good approximation is found to be given by $`\mathrm{\Gamma }=\mathrm{\Omega }_\mathrm{m}h/\mathrm{exp}[\mathrm{\Omega }_\mathrm{b}(1+1/\mathrm{\Omega }_\mathrm{m})]`$ (Efstathiou et al. 1992; Sugiyama 1995). The error in the approximation for the standard CDM model with $`\mathrm{\Omega }_\mathrm{b}=0.05`$, for example, is $`<10`$%. A higher-accuracy fit for this much-studied model can be achieved by setting $`\mathrm{\Omega }_\mathrm{b}=0`$ in equation (2) (i.e. setting $`\mathrm{\Gamma }=\mathrm{\Omega }_\mathrm{m}h`$) and modifying the coefficients to $`\alpha _1=2.205,\alpha _2=4.05,\alpha _3=18.3,\alpha _4=8.725`$, and $`\alpha _5=8.0`$. The fractional error relative to the direct numerical result is then reduced to smaller than 1% for $`k<40h`$ Mpc<sup>-1</sup> (Ma 1996).

The linear power spectra for the C+HDM models require additional treatment since the effect of neutrino free-streaming on the shape of the power spectrum is both time- and scale-dependent. It is found (Ma 1996) that by introducing a second shape parameter

$$\mathrm{\Gamma }_\nu =a^{1/2}\mathrm{\Omega }_\nu h^2$$
(3)   
to characterize the neutrino free-streaming distance, accurate approximations to the linear power spectra in C+HDM models can be obtained. The cold and hot components clearly have different power spectra due to their different thermal properties. For the CDM component in C+HDM models, a good approximation is given by

$$P_c(k,a,\mathrm{\Omega }_\nu )=P_c(k,a,\mathrm{\Omega }_\nu =0)\left(\frac{1+b_1x^{b_4/2}+b_2x^{b_4}}{1+b_3x_0^{b_4}}\right)^{\mathrm{\Omega }_\nu ^{1.05}},x=\frac{k}{\mathrm{\Gamma }_\nu },x_0=x(a=1),$$
(4)   
where $`P_c(k,a,\mathrm{\Omega }_\nu =0)`$ for the pure CDM model is given by equation (2), and $`b_1=0.01647,b_2=2.803\times 10^5,b_3=10.90`$, and $`b_4=3.259`$ for $`k`$ in units of Mpc<sup>-1</sup>. The functional form for the ratio $`P_c(k,a,\mathrm{\Omega }_\nu )/P_c(k,a,0)`$ is chosen to have the asymptotic behavior $`a^{2(f_{\mathrm{}}1)}`$, which can be derived analytically, and $`f_{\mathrm{}}=(5\sqrt{124\mathrm{\Omega }_\nu /25}1)/4`$ is the asymptotic growth rate for $`k\mathrm{}`$, where $`1f_{\mathrm{}}\mathrm{\Omega }_\nu ^{1.05}`$ is a good approximation. For the HDM component in C+HDM models, an accurate approximation is given by

$$P_\nu (k,a,\mathrm{\Omega }_\nu )=P_c(k,a,\mathrm{\Omega }_\nu )\left(\frac{e^{c_1x}}{1+c_2x^{1/2}+c_3x+c_4x^{3/2}+c_5x^2}\right),x=\frac{k}{\mathrm{\Gamma }_\nu h},$$
(5)   
where $`c_1=0.0015,c_2=0.1207,c_3=0.1015,c_4=0.01618`$, and $`c_5=0.001711`$ for $`k`$ in units of Mpc<sup>-1</sup>.

It is often useful to have an analytic approximation for the density-weighted power spectrum $`P(k)=\{\mathrm{\Omega }_\nu P_\nu ^{1/2}+(1\mathrm{\Omega }_\nu )P_c^{1/2}\}^2`$ that measures the total gravitational fluctuations contributed by the two components. Here the CDM and baryons have been assumed to have the same power (i.e., $`P_c=P_b`$), which is a good approximation for the range of redshifts and $`\mathrm{\Omega }_b`$ studied in this paper. The functional form used for the CDM spectrum $`P_c`$ in equation (4) works well here, and a good approximation for the density-weighted power spectrum in C+HDM models is given by

$$\frac{P(k,a,\mathrm{\Omega }_\nu )}{P_c(k,a,\mathrm{\Omega }_\nu =0)}=g(x,\mathrm{\Omega }_\nu )=\left(\frac{1+d_1x^{d_4/2}+d_2x^{d_4}}{1+d_3x_0^{d_4}}\right)^{\mathrm{\Omega }_\nu ^{1.05}},x=\frac{k}{\mathrm{\Gamma }_\nu },x_0=x(a=1),$$
(6)   
where $`P_c(k,a,\mathrm{\Omega }_\nu =0)`$ again is given by equation (2), and the coefficients are $`d_1=0.004321,d_2=2.217\times 10^6,d_3=11.63`$, and $`d_4=3.317`$ for $`k`$ in units of Mpc<sup>-1</sup>.

It should be noted that the C+HDM power spectra do not obey the simple evolution $`P(k)a^2`$ in equation (2) for the flat pure CDM model. The free-streaming of massive neutrinos slows down with time, allowing the neutrinos to cluster gravitationally and the neutrino density perturbations to grow on increasingly smaller scales. As a result, the growth of the C+HDM power spectra is both scale- and time-dependent. This effect is taken care of by the time-dependent parameter $`\mathrm{\Gamma }_\nu `$ in equation (3), which is built in in equations (4)-(6) via the scaled variable $`x`$.

## 3 Nonlinear Gravitational Clustering

### 3.1 Numerical Simulations

The linear theory of gravitational clustering is an elegant and powerful theory. It describes accurately the growth of density perturbations from the early Universe until a redshift of $`100`$, and the solution to the theory involves straightforward time-integration of a set (albeit a large set) of coupled ordinary differential equations. Gravitational clustering, however, eventually becomes a nonlinear process, and the study of the fully nonlinear process ultimately relies on numerical simulations.

Cosmological simulations can be roughly divided into two categories: (1) $`N`$-body simulations that deal only with dissipationless gravitational interactions of dark matter; and (2) hydrodynamical simulations that model gaseous dissipation via cooling and heating processes in addition to the gravitational interactions among the dark matter and gas. Depending on the models, the simulations are started at a redshift between 20 and 100 when the rms density fluctuations are well below unity, and the initial conditions of the simulations are generated from the linear power spectrum at that epoch.

As an illustration of gravitational clustering in collisionless $`N`$-body simulations, Figure 2 shows the projection of the smoothed matter distribution at four redshifts, $`z=5,3,1.5`$, and 0, from an $`N`$-body simulation of the $`\mathrm{\Omega }_\nu =0.1`$ C+HDM model. The model assumes $`\mathrm{\Omega }_c=0.85`$, $`\mathrm{\Omega }_\mathrm{b}=0.05`$, and $`h=0.5`$, and is normalized to the COBE quadrupole $`Q_{\mathrm{rms}\mathrm{PS}}=18\mu K`$ (Gorski et al. 1996). The simulation box is 100 Mpc comoving on a side, and the force softening is 50 kpc. A total of $`128^3`$ cold and $`128^3`$ hot particles are used. Each panel in Figure 2 is 100$`\times `$100 Mpc comoving, and shows the projection along one axis of the entire simulation box. Periodic boundary conditions are adopted in the simulation, so for example, the dense clumps near the central left and right edges of the box in the $`z=0`$ panel (lower right) are part of the same cluster of dark matter halos. The darkest several halos at $`z=0`$ all have masses above $`10^{14}M_{}`$.

The power spectrum of the density fluctuations offers the lowest-order measurement of the gravitational clustering exhibited in Figure 2. In Figure 3, we plot the corresponding nonlinear power spectra computed from the particle spatial distributions at $`z=3,1.5`$, and 0 shown in Figure 2. The linear power spectra given by equation (6) are also shown for comparison. The hierarchical, or “bottom-up”, nature of gravitational collapse in these models is evident: the high-$`k`$ modes (i.e. small length scales) have become strongly nonlinear, while the low-$`k`$ modes are still following the linear power spectrum. The fact that the three lowest $`k`$ modes are still linear at $`z=0`$ ensures that our choice of the simulation box size (100 Mpc) is large enough to include all waves that have gone nonlinear at the present epoch. It is also interesting to note that the point of departure from linearity moves to the left of the figure as $`z`$ decreases, indicating that objects become nonlinear on increasingly larger length scales as the Universe evolves.

It is also instructive to compare the nonlinear structures at the same epoch in different cosmological models. Figure 4 shows the smoothed matter distribution at $`z=1.5`$ from simulations of the four models in Figure 1: a LCDM, and C+HDM with $`\mathrm{\Omega }_\nu =0.1,0.2`$, and 0.3 (clockwise starting from upper left). The LCDM model has $`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, and $`h=0.75`$, and the three C+HDM models have $`h=0.5`$, $`\mathrm{\Omega }_\mathrm{b}=0.05`$, and $`\mathrm{\Omega }_c+\mathrm{\Omega }_\nu +\mathrm{\Omega }_\mathrm{b}=1`$. All models are COBE-normalized, and all simulations are performed with the same box size and force resolution as the $`\mathrm{\Omega }_\nu =0.1`$ run described earlier. (The only exception is that the $`\mathrm{\Omega }_\nu =0.2`$ and 0.3 simulations used $`10\times 128^3`$ instead of $`128^3`$ particles to represent the HDM component.) The initial conditions for all four simulations are generated with the same random phases; structures therefore appear at similar locations in all panels. The effect of $`\mathrm{\Omega }_\nu `$ on structure formation, however, is striking: the higher $`\mathrm{\Omega }_\nu `$ is, the fewer collapsed objects there are at a given epoch. This trend simply reflects the decrease in the small-scale linear power shown in Figure 1, which results from neutrino free-streaming. Note that the model with $`\mathrm{\Omega }_\nu =0.3`$ (lower left panel), which corresponds to a neutrino mass of 7 eV, has very few structures even at $`z=1.5`$. It is also interesting to note that the upper two panels of Figure 4 for the LCDM and $`\mathrm{\Omega }_\nu =0.1`$ C+HDM models look very similar. This is because the linear power spectra for the two models are in fact very similar at $`k>0.1h`$ Mpc<sup>-1</sup>, as shown in Figure 1. The corresponding linear and nonlinear power spectra at $`z=1.5`$ for the three C+HDM models in Figure 4 are shown in Figure 5.

The power spectra shown thus far for the C+HDM models are the density-averaged $`P(k)=\{\mathrm{\Omega }_\nu P_\nu ^{1/2}+(1\mathrm{\Omega }_\nu )P_c^{1/2}\}^2`$, where $`P_\nu `$ and $`P_c`$ are the individual power spectra for the hot and cold components. (It has been assumed $`P_c=P_{\mathrm{baryon}}`$, which is a good approximation for the redshift range of interest here.) The two components evolve distinctly due to their different thermal velocities, and the details of the shape and the growth of the linear power spectra are discussed in Ma (1996). Here, the nonlinear power spectra for the separate components are presented in Figure 6. For clarity, only two epochs, $`z=3`$ and 0, are plotted. Notice how at $`z=3`$, the HDM spectrum remains linear to almost $`k=2h`$Mpc<sup>-1</sup> but the CDM remains linear to only $`k0.5h`$Mpc<sup>-1</sup>, while at $`z=0`$, both components become nonlinear at $`k0.3h`$Mpc<sup>-1</sup>. This is largely due to the slowing down of the neutrinos, which makes it easier for the neutrino particles to fall into the CDM potential wells at later times, and therefore become nonlinear at smaller $`k`$. More precisely, the median thermal velocity of neutrinos is

$$v3(1+z)\frac{k_BT_{0,\nu }}{m_\nu c}=15(1+z)\frac{10\mathrm{eV}}{m_\nu }\mathrm{km}/\mathrm{s},$$
(7)   
so the neutrinos in the $`\mathrm{\Omega }_\nu =0.1`$ model ($`m_\nu =2.3`$ eV) have slowed down by a factor of 4 since $`z=3`$ to about 60 km/s today.

### 3.2 Analytical Approximations

Just as it was useful to cast the linear power spectra in simple functional forms for a wide range of models (see equations (1)-(6)), analytical approximations to the nonlinear power spectrum can also provide valuable insight to the process of gravitational clustering. It is, however, more difficult to compute the nonlinear $`P_{\mathrm{nl}}`$ because the numerical simulations that have sufficient dynamic range to allow calculations of $`P_{\mathrm{nl}}`$ over a wide range of scales are generally much more CPU-intensive than the integration of the linearized Boltzmann equations. The behavior of the power spectrum in the nonlinear regime therefore is less well-understood.

An early attempt to relate linear and nonlinear quantities was focused on the spatially-averaged two-point correlation function $`\overline{\xi }(r)`$ in models with $`\mathrm{\Omega }_\mathrm{m}=1`$ and a power-law power spectrum (Hamilton et al. 1991). It was found that if a given nonlinear scale $`r_{\mathrm{nl}}`$ is identified with its pre-collapsed linear scale $`r_\mathrm{l}`$ by $`r_\mathrm{l}^3=r_{\mathrm{nl}}^3(1+\overline{\xi }_{\mathrm{nl}})`$, then there exists a simple, universal function relating the linear and nonlinear two-point correlation function:

$$\overline{\xi }_{\mathrm{nl}}(r_{\mathrm{nl}})=F[\overline{\xi }_\mathrm{l}(r_\mathrm{l})].$$
(8)   
At that time, this transformation appeared to be magically independent of the spectral index $`n`$ of the linear power spectrum assumed in the model, and the generality of this formula rendered the task of reconstructing the primordial spectrum from the observed nonlinear clustering of galaxies less complicated. However, further tests against numerical simulations showed that this formula worked well only for spectral index $`n>1`$, and it erred by factors up to 3 and 10 for $`n=1.5`$ and $`2`$, respectively (Jain et al. 1995). In particular, the nonlinear $`\overline{\xi }_{\mathrm{nl}}`$ was found to rise more sharply with increasing $`\overline{\xi }_\mathrm{l}`$ for models with more negative $`n`$, so Jain et al. introduced an $`n`$-dependent formula to accommodate this feature. For the more realistic cosmological models such as the CDM where the spectral index is a function of scale, they proposed to use an effective index $`n_{\mathrm{eff}}`$, defined to be $`n_{\mathrm{eff}}=d\mathrm{ln}P/d\mathrm{ln}k_0`$, where $`k_0`$ is the scale at which the rms mass fluctuation $`\sigma (R=1/k_0)`$ is unity. The index $`n_{\mathrm{eff}}`$ therefore represents the slope of the power spectrum at a scale where nonlinear effects are becoming important. Later work extended the analysis to models with varying matter density $`\mathrm{\Omega }_\mathrm{m}`$ and cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }`$ (Peacock & Dodds 1996), where instead of using an effective index for all scales, the local slope $`n(k)=d\mathrm{ln}P/d\mathrm{ln}k`$ was adopted.

Here we examine the nonlinear mapping of the power spectrum in the C+HDM models, which has not been explored yet. Figure 7 shows how the nonlinear density variance $`\mathrm{\Delta }_{\mathrm{nl}}=4\pi k_{\mathrm{nl}}^3P_{\mathrm{nl}}(k_{\mathrm{nl}})`$ diverges from the linear $`\mathrm{\Delta }_\mathrm{l}=4\pi k_\mathrm{l}^3P_\mathrm{l}(k_\mathrm{l})`$, where the linear and nonlinear wavenumbers are related by $`k_\mathrm{l}^3=k_{\mathrm{nl}}^3/(1+\mathrm{\Delta }_{\mathrm{nl}})`$. The squares are obtained from simulations of the $`\mathrm{\Omega }_\nu =0.2`$ (left curve) and 0.1 (right) C+HDM models. The curves generally obey the asymptotic condition $`\mathrm{\Delta }_{\mathrm{nl}}=\mathrm{\Delta }_\mathrm{l}`$ for $`\mathrm{\Delta }_\mathrm{l}1`$, and $`\mathrm{\Delta }_{\mathrm{nl}}\mathrm{\Delta }_\mathrm{l}^{1.5}`$ in the highly-nonlinear, stable clustering regime. The mapping, however, is clearly not independent of cosmology: the larger $`\mathrm{\Omega }_\nu `$ is in a model, the faster $`\mathrm{\Delta }_{\mathrm{nl}}/\mathrm{\Delta }_\mathrm{l}`$ increases at $`\mathrm{\Delta }_\mathrm{l}1`$. This trend reflects the $`n`$-dependence pointed out by Jain et al. (1995) and can be explained by the different shapes of $`P_\mathrm{l}`$ shown in Figure 1, where the models with higher $`\mathrm{\Omega }_\nu `$ have less power and hence more negative slope at high $`k`$. Since all models are normalized to COBE and have the same amplitude at low $`k`$, the critical $`k_0`$ (where $`\sigma (k_0)=1`$) increases for larger $`\mathrm{\Omega }_\nu `$. The effective index $`n_{\mathrm{eff}}`$ defined at $`k_0`$ is therefore more negative for larger $`\mathrm{\Omega }_\nu `$, resulting in the steeper rise of $`\mathrm{\Delta }_{\mathrm{nl}}`$ in Figure 7.

However, neither formula proposed by Jain et al. (1995) or Peacock & Dodds (1996) can be extended to the C+HDM models. The dotted curves in Figure 7 illustrate the large discrepancies in the Peacock-Dodds fitting function when it is applied to the $`\mathrm{\Omega }_\nu =0.1`$ and 0.2 models. In retrospect, it is not surprising their formulas do not apply: Although the power spectrum at a given epoch is determined entirely from the spatial distribution of the particles, the evolution of the power spectrum depends on the particle velocities as well. A general formula for the linear to nonlinear transformation therefore must depend on both the shape and the growth rate of $`P(k)`$. Since both formulas are designed for models without massive neutrinos and depend only on the shape of $`P(k)`$, they cannot be applied to C+HDM models.

Here we propose a new analytical approximation for the mapping of linear and nonlinear power spectrum in CDM as well as C+HDM models: <sup>1</sup><sup>1</sup>1An improved approximation with a higher accuracy has since appeared in Ma (1998).

$$\mathrm{\Delta }_{\mathrm{nl}}(k_{\mathrm{nl}})=\mathrm{\Delta }_\mathrm{l}(k_\mathrm{l})G\left[\mathrm{\Delta }_l(k_\mathrm{l})\right],$$
(9)   
where

$$G(x)=\frac{1+a_1x^4+a_2x^8/g^{2.5}}{1+a_3x^4+a_4x^{7.5}/g^2},$$
(10)   
and the coefficients are $`a_1=4756,a_2=384.6,a_3=3732`$, and $`a_4=24.20`$. The functional form of $`G`$ is chosen to give the appropriate asymptotic behavior $`\mathrm{\Delta }_{\mathrm{nl}}\mathrm{\Delta }_\mathrm{l}`$ in the linear regime ($`x1`$) and $`\mathrm{\Delta }_{\mathrm{nl}}\mathrm{\Delta }_\mathrm{l}^{3/2}`$ in the stable clustering regime ($`x1`$). The dependence of $`G`$ on $`\mathrm{\Omega }_\nu `$ comes from the function $`g`$ of equation (6), which gives the relative amplitudes of the power spectra in the C+HDM models and the pure CDM model. This function is analogous to the commonly-used growth factor $`g(a=1,\mathrm{\Omega }_\mathrm{m},\mathrm{\Omega }_\mathrm{\Lambda })=\frac{5}{2}\mathrm{\Omega }_\mathrm{m}[\mathrm{\Omega }_\mathrm{m}^{4/7}\mathrm{\Omega }_\mathrm{\Lambda }+(1+\mathrm{\Omega }_\mathrm{m}/2)(1+\mathrm{\Omega }_\mathrm{\Lambda }/70)]^1`$ for LCDM models (Lahav et al. 1991; Carroll, Press & Turner 1992). They differ, however, in that the growth factor is scale-independent in LCDM models but is a function of scale in C+HDM models since neutrinos only affect the growth below the free-streaming scale, as discussed in Section 2.

## 4 Summary

The power spectrum is a fundamental measure of gravitational clustering in cosmology. We have discussed the features and evolution of both the linear and nonlinear power spectra in various cosmological models.

The linear power spectrum $`P_l`$ is calculated from time integration of the coupled, linearized Einstein, Boltzmann, and fluid equations for the metric and density perturbations. The shape and growth of the power spectrum depend on cosmological parameters such as the total matter density $`\mathrm{\Omega }_\mathrm{m}`$, the neutrino fraction $`\mathrm{\Omega }_\nu `$, and the Hubble constant $`h`$. Simple analytical functions were presented that can approximate the linear $`P_l`$ in both CDM and C+HDM models with less than 10% error.

The fully nonlinear power spectrum $`P_{\mathrm{nl}}`$ has to be computed from numerical simulations. Simulation output for CDM, LCDM, and C+HDM models was presented to illustrate the sensitive dependence of structure formation and evolution on cosmological parameters. We also discussed the present understanding of the relation between the linear and nonlinear power spectra, and proposed a new approximation which applies to CDM as well as C+HDM models.

The supercomputing time for this work was generously provided by the National Scalable Cluster Project at the University of Pennsylvania and the National Center for Supercomputing Applications.