# Analytical Approximation to the Nonlinear Power Spectrum of Gravitational Clustering

## 1 Introduction

The power spectrum of matter fluctuations $`P(k)`$ provides the most fundamental statistical measure of gravitational clustering. When the amplitude of the density fluctuations is small, the power spectrum can be calculated easily from the linear perturbation theory of gravitational collapse. In the nonlinear regime, however, determination of the fully evolved power spectrum for a given cosmological model requires numerical simulations. Since nearly all observable astronomical systems have experienced some nonlinear collapse, it would provide much physical insight and practical convenience to devise a general analytical approximation (based on simulation results) for the nonlinear power spectrum for a wide range of cosmologically interesting models.

By fitting to $`N`$-body results, Hamilton et al. (1991) studied scale-free models with a power-law spectrum and presented a simple analytical formula that relates the spatially averaged nonlinear and linear two-point correlation function, $`\overline{\xi }_{\mathrm{nl}}(r)`$ and $`\overline{\xi }_\mathrm{l}(r_0)`$, where $`r`$ is related to its pre-collapsed linear scale $`r_0`$ by $`r_0=r(1+\overline{\xi }_{\mathrm{nl}})^{1/3}`$. This transformation then appeared to be magically independent of the spectral index $`n`$ assumed in the model. Further tests against numerical simulations, however, found significant errors when the Hamilton et al. function was applied to models with $`n<1`$ (Jain, Mo, & White 1995; Padmanabhan 1996). Jain et al. (1995) instead proposed $`n`$-dependent formulas to relate $`\overline{\xi }`$ and $`P(k)`$ in the linear and nonlinear regimes in both scale-free models and the standard CDM model. For the more realistic CDM model, for which the spectral index changes from the primordial value $`n1`$ on large scales to nearly $`3`$ on small scales, they used an effective index given by $`n_{\mathrm{eff}}=d\mathrm{ln}P(k)/d\mathrm{ln}k|_{k_c}`$, where $`k_c`$ is the scale at which the rms mass fluctuation $`\sigma `$ is unity. The index $`n_{\mathrm{eff}}`$ therefore reflects the slope of the power spectrum at the length scale where nonlinearity becomes important. Peacock and Dodds (1996) extended this work to allow for a low $`\mathrm{\Omega }_\mathrm{m}`$ and a non-zero cosmological constant.

No previous work has investigated in detail the subject of linear to nonlinear mapping in cold+hot dark matter (C+HDM) models that assume massive neutrinos are a component of the dark matter. This is perhaps because the physics in C+HDM models is generally more complicated than in CDM or LCDM (CDM with a cosmological constant) models due to the additional length scale associated with the free streaming of the neutrinos (Ma 1996). Nevertheless, massive neutrinos remain a prime dark matter candidate, and the recent evidence for neutrino masses from the Super-Kamiokande experiment has made this possibility particularly intriguing (Fukuda et al. 1998). Although neither Jain et al. or Peacock & Dodds has tested these models, one may surmise that their formulas can be naturally extended to C+HDM models as long as the spectral index in the formula is calculated from the C+HDM power spectrum. This unfortunately does not work. Both fitting functions underestimate the nonlinear density variance,<sup>1</sup><sup>1</sup>1The notations $`\mathrm{\Delta }`$ and $`P(k)`$ here are the same as in Jain et al., and are equivalent to $`\mathrm{\Delta }^2`$ and $`P(k)/(2\pi )^3`$ in Peacock & Dodds . $`\mathrm{\Delta }_{\mathrm{nl}}=4\pi k^3P_{\mathrm{nl}}(k)`$, at $`k>2h`$ Mpc<sup>-1</sup> in C+HDM models, and the errors in Peacock & Dodds , for example, reach $`50`$% at $`k10h`$ Mpc<sup>-1</sup>. The linear to nonlinear transformation is therefore regulated by more than simply $`n_{\mathrm{eff}}`$. (Smith et al. (1998) recently reported agreement between the Peacock-Dodds formula and results from two C+HDM simulations. Their simulation resolution of $`0.3h^1`$ Mpc, however, limited their test to only the mildly nonlinear regime, and could not probe the nonlinear regime where the large discrepancies reside.)

This Letter differs from previous work in two ways. First, the simple analytical formula presented here closely approximates the fully nonlinear power spectrum of mass fluctuations at $`z<5`$ in the previously unexplored C+HDM models as well as LCDM models with varying $`\mathrm{\Omega }_\mathrm{m}`$. Numerical simulations of four COBE-normalized flat C+HDM and flat LCDM models are performed to calibrate the coefficients in the analytical formula. Second, the formula introduced here depends directly on $`\sigma _8`$ (the rms linear mass fluctuation on $`8h^1`$ Mpc scale) instead of a spectral index as in previous work. This is achieved by recognizing a scaling rule (see § 3), which also greatly simplifies the construction of the analytical formula. This work also extends into the nonlinear regime a previous investigation of the effects of neutrino free-streaming on the linear C+HDM power spectrum (Ma 1996).

## 2 Input Linear Power Spectrum

For a wide range of CDM and LCDM models that assume neutrinos are massless, a good approximation to the linear power spectrum is given by

$$P(k,a,\mathrm{\Omega }_\nu =0)=\frac{Ak^n[D(a)/D_0]^2\left[\mathrm{ln}(1+\alpha _1q)/\alpha _1q\right]^2}{[1+\alpha _2q+(\alpha _3q)^2+(\alpha _4q)^3+(\alpha _5q)^4]^{1/2}},$$
(1)   
where $`k`$ is the wavenumber in units of Mpc<sup>-1</sup>, $`q=k/\mathrm{\Gamma }h`$, $`\mathrm{\Gamma }`$ is a shape parameter, 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 shape parameter $`\mathrm{\Gamma }`$ characterizes the dependence on cosmological parameters and is well approximated 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; see also Bunn & White 1997). The function $`D(a)`$ is the linear growth factor, whose present value is $`D_0=D(a=1)`$, and it can be expressed as $`D(a)=ag`$, where the relative growth factor $`g`$ is well approximated by $`g(\mathrm{\Omega }_\mathrm{m}(a),\mathrm{\Omega }_\mathrm{\Lambda }(a))=2.5\mathrm{\Omega }_\mathrm{m}(a)[\mathrm{\Omega }_\mathrm{m}(a)^{4/7}\mathrm{\Omega }_\mathrm{\Lambda }(a)+\left(1+\mathrm{\Omega }_\mathrm{m}(a)/2\right)\left(1+\mathrm{\Omega }_\mathrm{\Lambda }(a)/70\right)]^1`$ (Lahav et al. 1991; Carroll et al. 1992). In LCDM models, $`g1`$ until the universe becomes $`\mathrm{\Lambda }`$-dominated at $`1+z\mathrm{\Omega }_\mathrm{m}^{1/3}`$; the value of $`g`$ then decreases with increasing $`a`$. The normalization factor $`A`$ can be chosen by fixing the value of $`\sigma _8`$; if instead the COBE normalization is desired, it is $`A=\delta _H^2(c/H_0)^{n+3}/(4\pi )`$, where (for flat models) $`\delta _H=`$ $`1.94\times 10^5\mathrm{\Omega }_\mathrm{m}^{0.7850.05\mathrm{ln}\mathrm{\Omega }_\mathrm{m}}\mathrm{exp}(0.95\stackrel{~}{n}0.169\stackrel{~}{n}^2)`$ with $`\stackrel{~}{n}=n1`$ (Bunn & White 1997).

The linear power spectra for the C+HDM models require additional treatment since the effect of massive neutrinos on the shape of the power spectrum is both time and scale dependent. It is found that by introducing a second shape parameter, $`\mathrm{\Gamma }_\nu =a^{1/2}\mathrm{\Omega }_\nu h^2`$, to characterize the neutrino free-streaming distance, one can obtain a good approximation to the linear power spectra (density averaged over the cold and hot components) in flat C+HDM models at $`z<5`$ when neutrinos are adequately nonrelativistic (Ma 1996):

$$P(k,a,\mathrm{\Omega }_\nu )=P(k,a,\mathrm{\Omega }_\nu =0)\left(\frac{1+d_1x^{d_4/2}+d_2x^{d_4}}{1+d_3x_0^{d_4}}\right)^{\mathrm{\Omega }_\nu ^{1.05}},$$
(2)   
where $`x=k/\mathrm{\Gamma }_\nu `$, $`x_0=x(a=1)`$, $`P(k,a,\mathrm{\Omega }_\nu =0)`$ for the pure CDM model is given by eq. (1), and $`d_1=0.004321,d_2=2.217\times 10^6,d_3=11.63`$, and $`d_4=3.317`$ for $`k`$ in Mpc<sup>-1</sup>. (The scale factor in $`\mathrm{\Gamma }_\nu `$ assumes COBE normalization at $`a=1`$. If other normalization is used, replace $`a`$ by $`a\sigma _8/\sigma _8^{\mathrm{cobe}}`$.)

## 3 Nonlinear Power Spectrum

Numerical simulations of structure formation in two flat C+HDM models with neutrino fraction $`\mathrm{\Omega }_\nu =0.1`$ and 0.2 and two flat LCDM models with matter density $`\mathrm{\Omega }_\mathrm{m}=0.3`$ and 0.5 are performed in order to obtain the nonlinear power spectra of matter fluctuations. All four simulations are performed in a (100 Mpc)<sup>3</sup> comoving box. The gravitational forces are computed with a particle-particle particle-mesh (P<sup>3</sup>M) algorithm (Bertschinger & Gelb 1991; Ma et al. 1997) with a comoving Plummer force softening length of 50 kpc. An identical set of random phases is used in the initial conditions for all four runs. The primordial power spectrum has an index of $`n=1`$, with density fluctuations drawn from a random Gaussian field. A total of $`128^3`$ simulation particles are used to represent the cold dark matter. For the C+HDM models, $`128^3`$ and $`10\times 128^3`$ particles are used to represent the hot component in the $`\mathrm{\Omega }_\nu =0.1`$ and 0.2 models, respectively. Although the large particle number is needed to finely sample the velocity phase space (Ma & Bertschinger 1994), tests performed for the $`\mathrm{\Omega }_\nu =0.2`$ model show that the power spectrum itself is little affected when $`128^3`$ hot particles are used. Since structure forms too late in flat C+HDM models with $`\mathrm{\Omega }_\nu >0.2`$ (Ma et al. 1997 and references therein), only the $`\mathrm{\Omega }_\nu =0.1`$ and 0.2 models are studied here. Both models assume $`\mathrm{\Omega }_\mathrm{b}=0.05`$ and $`h=0.5`$, while $`\mathrm{\Omega }_{\mathrm{cdm}}=0.85`$ and 0.75 for the two models, respectively. The two LCDM models chosen for the simulations have $`(\mathrm{\Omega }_\mathrm{m},\mathrm{\Omega }_\mathrm{\Lambda },h)=(0.3,0.7,0.75)`$ and $`(0.5,0.5,0.7)`$, and $`\mathrm{\Omega }_\mathrm{b}=0`$. All four models are normalized to the 4-year COBE results (Bennett et al. 1996; Gorski et al. 1996).

Figure 1 contrasts the linear (dotted) and nonlinear (solid) power spectra at various redshifts for the four simulated models. The hierarchical nature of gravitational collapse in these models is evident: the high-$`k`$ modes have become strongly nonlinear, whereas the low-$`k`$ modes still follow the linear power spectrum. The dashed curves are from the analytical approximation described below. Figure 2 illustrates the dependence of the linear to nonlinear transformation on both time and cosmological parameters by plotting at various redshifts the ratio of the nonlinear and linear density variance, $`\mathrm{\Delta }_{\mathrm{nl}}(k)/\mathrm{\Delta }_\mathrm{l}(k_0)`$, against the linear $`\mathrm{\Delta }_\mathrm{l}(k_0)`$ (where $`\mathrm{\Delta }=4\pi k^3P`$). Note that $`\mathrm{\Delta }_{\mathrm{nl}}`$ and $`\mathrm{\Delta }_\mathrm{l}`$ are evaluated at different wavenumbers, where $`k_0=k(1+\mathrm{\Delta }_{\mathrm{nl}})^{1/3}`$ corresponds to the pre-collapsed scale of $`k`$, as introduced in Hamilton et al. (1991). The dependence of $`\mathrm{\Delta }_{\mathrm{nl}}(k)/\mathrm{\Delta }_\mathrm{l}(k_0)`$ on $`\mathrm{\Omega }_\mathrm{m}`$, $`\mathrm{\Omega }_\nu `$, and time, however, demonstrates that the universality seen by Hamilton et al. is limited to the power-law models with $`n>1`$ studied there. The linear to nonlinear mapping is clearly more complicated in more realistic cosmological models (such as LCDM and C+HDM) whose spectral index changes continuously with length scale.

Despite the apparent complication from the time dependence in Figure 2, I find the scaling rule given below very helpful in simplifying the construction of an analytical approximation. As Figure 3 shows, the time dependence is largely removed when $`\mathrm{\Delta }_{\mathrm{nl}}/\mathrm{\Delta }_\mathrm{l}`$ is plotted against the scaled quantity $`\stackrel{~}{\mathrm{\Delta }_\mathrm{l}}`$,

$$\stackrel{~}{\mathrm{\Delta }_\mathrm{l}}=\frac{\mathrm{\Delta }_\mathrm{l}}{\sigma _8^\beta },\beta =0.7+10\mathrm{\Omega }_\nu ^2$$
(3)   
instead of $`\mathrm{\Delta }_\mathrm{l}`$. For example, $`\beta =0.7`$ for the CDM and LCDM models, while $`\beta =0.8`$ and 1.1 for the $`\mathrm{\Omega }_\nu =0.1`$ and 0.2 C+HDM models. The deviations in the two LCDM models at late times ($`a>0.7`$) result from the retardation in the relative growth factor $`g`$ in $`\mathrm{\Omega }_\mathrm{m}<1`$ models. As shown in eq. (4) below, they can be easily accounted for by including a factor of $`g^3`$ in the highly nonlinear regime in the fitting formula.

The scaling behavior in eq. (3) can be understood in terms of the time dependence of the local slope of the linear power spectrum near the scale where nonlinear effects become important. An example is the effective spectral index $`n_{\mathrm{eff}}+3=d\mathrm{ln}\mathrm{\Delta }_\mathrm{l}/d\mathrm{ln}k|_{k_c}`$ considered by Jain et al. (1995), where $`k_c`$ is the scale at which the rms linear mass fluctuation $`\sigma `$ equals unity. Figure 4 illustrates the dependence of $`n_{\mathrm{eff}}`$ on time and cosmological parameters. The value of $`n_{\mathrm{eff}}`$ generally increases with time because as $`\sigma `$ grows, the wavenumber $`k_c`$ at which $`\sigma =1`$ decreases, and the spectral index at $`k_c`$ (i.e. $`n_{\mathrm{eff}}`$) becomes larger since the slope of the power spectrum for the models studied here always increases with decreasing $`k`$. In Figure 4, $`n_{\mathrm{eff}}`$ exhibits the fastest growth in the $`\mathrm{\Omega }_\nu =0.2`$ C+HDM model at $`\sigma _8>0.3`$ because the neutrino free streaming effect is more prominent in higher $`\mathrm{\Omega }_\nu `$ models, which acts to suppress structure growth below the free streaming scale, causing $`\mathrm{\Delta }_\mathrm{l}`$ to bend more at $`0.1<k<1h`$ Mpc<sup>-1</sup> (see the dotted curves in Fig. 1). At $`\sigma _8<0.3`$ in the $`\mathrm{\Omega }_\nu =0.2`$ model, on the other hand, $`n_{\mathrm{eff}}`$ stays nearly constant because it is probing the nearly flat, $`k>5h`$ Mpc<sup>-1</sup> part of $`\mathrm{\Delta }_\mathrm{l}`$. Despite this interesting behavior, Figure 4 shows that for all models, the dependence of $`n_{\mathrm{eff}}+3`$ on $`\sigma _8`$ is well approximated by a power law at $`\sigma _8>0.3`$: $`d\mathrm{ln}(n_{\mathrm{eff}}+3)/d\mathrm{ln}\sigma _8\beta `$, where $`\beta =0.7+10\mathrm{\Omega }_\nu ^2`$ as given in eq. (3). This therefore explains why replacing the factor $`n_{\mathrm{eff}}+3`$ in the earlier work with $`\sigma _8^\beta `$ works well here.

The simple scaling behavior introduced by eq. (3) allows one to approximate the evolution of the nonlinear power spectrum directly in terms of $`\sigma _8`$ instead of $`n_{\mathrm{eff}}`$. Combining these factors, I find that a close approximation for the nonlinear power spectrum is given by

$$\frac{\mathrm{\Delta }_{\mathrm{nl}}(k)}{\mathrm{\Delta }_\mathrm{l}(k_0)}=G\left(\frac{\mathrm{\Delta }_\mathrm{l}}{g_0^{1.5}\sigma _8^\beta }\right),G(x)=[1+\mathrm{ln}(1+0.5x)]\frac{1+0.02x^4+c_1x^8/g^3}{1+c_2x^{7.5}},$$
(4)   
where $`\beta `$ is given by eq. (3), $`k_0=k(1+\mathrm{\Delta }_{\mathrm{nl}})^{1/3}`$, and $`g_0=g(\mathrm{\Omega }_\mathrm{m},\mathrm{\Omega }_\mathrm{\Lambda })`$ and $`g=g(\mathrm{\Omega }_\mathrm{m}(a),\mathrm{\Omega }_\mathrm{\Lambda }(a))`$ are, respectively, the relative growth factor<sup>2</sup><sup>2</sup>2 Even for C+HDM models, $`g`$ in eq. (4) is taken to be the familiar $`g(\mathrm{\Omega }_\mathrm{m}(a),\mathrm{\Omega }_\mathrm{\Lambda }(a))`$ for $`\mathrm{\Omega }_\nu =0`$ models. The true relative growth factor for C+HDM models is in fact given by eq. (2), but due to its complicated scale dependence, attempts thus far to incorporate this factor directly in eq. (4) have not led to approximations with high accuracies (see Ma 1998). at present day and at $`a`$ discussed in § 2. The time dependence is in factors $`\sigma _8^\beta `$ and $`g`$. For CDM and LCDM, a good fit is given by $`c_1=1.08\times 10^4`$ and $`c_2=2.10\times 10^5`$. For C+HDM, a good fit is given by $`c_1=3.16\times 10^3`$ and $`c_2=3.49\times 10^4`$ for $`\mathrm{\Omega }_\nu =0.1`$, and $`c_1=6.96\times 10^3`$ and $`c_2=4.39\times 10^4`$ for $`\mathrm{\Omega }_\nu =0.2`$. The dependence of $`c_1`$ and $`c_2`$ on $`\mathrm{\Omega }_\nu `$ can in principle be cast in a functional form (see Ma 1998), but since the allowed range of $`\mathrm{\Omega }_\nu `$ is narrow, separate coefficients are given here in order to obtain the highest possible fitting accuracy. The accuracy of eq. (4) is illustrated in Figure 1 (dashed curves), where the rms error for each model ranges from 3% to 10% for $`k<10h`$ Mpc<sup>-1</sup> at all times except $`z>4`$, when the errors are about 15%.

The functional form of $`G(x)`$ in eq. (4) 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`$). In the mildly nonlinear regime, $`0.1<\mathrm{\Delta }_\mathrm{l}<1`$, the pre-factor $`[1+\mathrm{ln}(1+0.5x)]`$ is introduced to approximate the non-negligible positive slope of $`\mathrm{\Delta }_{\mathrm{nl}}/\mathrm{\Delta }_\mathrm{l}`$. This factor is needed because $`\mathrm{\Delta }_{\mathrm{nl}}`$ and $`\mathrm{\Delta }_\mathrm{l}`$ are evaluated at different wavenumbers $`k`$ and $`k_0`$, where the pre-collapsed $`k_0`$ is always smaller than $`k`$. Due to the steep positive slope of $`\mathrm{\Delta }_\mathrm{l}`$ in this region, $`\mathrm{\Delta }_\mathrm{l}`$ at $`k_0`$ is noticeably smaller than at $`k`$, and $`\mathrm{\Delta }_{\mathrm{nl}}(k)/\mathrm{\Delta }_\mathrm{l}(k_0)`$ is thus significantly above unity. Without the logarithmic pre-factor in eq. (4) to account for this elevation, the approximation to $`\mathrm{\Delta }_{\mathrm{nl}}`$ can be underestimated by up to 30% at $`0.1<\mathrm{\Delta }_\mathrm{l}<1`$.

## 4 Summary and Discussion

This paper presents a single formula, eq. (4), that accurately approximates the fully nonlinear power spectrum of matter fluctuations for redshift $`z<5`$ for flat CDM, LCDM, and C+HDM models with varying matter density $`\mathrm{\Omega }_\mathrm{m}`$ and neutrino fraction $`\mathrm{\Omega }_\nu `$. Eqs. (1), (2) and (4) together offer a complete description of the shape and time evolution of the matter power spectrum in both linear and nonlinear regime for a wide range of cosmologically interesting models. Figure 1 summarizes the analytical and simulation results for four representative models. Depending on the models and epochs, the rms errors are between 3% and 10% for $`k<10h`$ Mpc<sup>-1</sup> at $`z<4`$ in eq. (4). In comparison, the Peacock & Dodds formula has a rms error between 6% and 17% for the two LCDM simulations studied here, and the error reaches 50% at $`k10h`$ Mpc<sup>-1</sup> for the two C+HDM models.

In contrast to the scale-free models studied in Hamilton et al. (1991), the relation between the linear and nonlinear power spectrum is far from universal. The different panels of Figure 2 illustrate the dependence on cosmological parameters $`\mathrm{\Omega }_\mathrm{m}`$ and $`\mathrm{\Omega }_\nu `$. Moreover, within a given model, the fact that the curves for different times in Figure 2 do not overlap is important because it implies that the ratio $`\mathrm{\Delta }_{\mathrm{nl}}(k)/\mathrm{\Delta }_\mathrm{l}(k_0)`$ depends not only on $`\mathrm{\Delta }_\mathrm{l}(k_0)`$ but also on time, or equivalently, the overall amplitude of $`\mathrm{\Delta }_\mathrm{l}(k_0)`$. This amplitude (or time) dependence is present, albeit in a somewhat subtle form, in earlier work. The effective spectral index used by Jain et al. (1995) is clearly time-varying as shown in Figure 4 of this paper. The local spectral index $`n=d\mathrm{ln}P/d\mathrm{ln}k|_{k_0}`$ used by Peacock and Dodds (1996) depends on time implicitly through the factor $`\mathrm{\Delta }_{\mathrm{nl}}`$ in the relation $`k_0=k(1+\mathrm{\Delta }_{\mathrm{nl}})^{1/3}`$. In comparison, in this paper I have adopted the commonly used parameter $`\sigma _8`$ instead of a spectral index to characterize this time dependence, and have shown that the scaling behavior in eq. (3) and Figure 3 absorbs this dependence in C+HDM as well as LCDM models and results in the simple formula, eq. (4).

I thank the referee, Andrew Hamilton, for insightful comments that have helped to improve the manuscript, and Robert Caldwell for stimulating discussions. Supercomputing time for the numerical simulations is provided by the National Scalable Cluster Project at the University of Pennsylvania and the National Center for Supercomputing Applications. A Penn Research Foundation Award is acknowledged.