# A Simple Method for Computing the Non-Linear Mass Correlation Function with Implications for Stable Clustering

## Abstract

We propose a simple and accurate method for computing analytically the mass correlation function for cold dark matter and scale-free models that fits N-body simulations over a range that extends from the linear to the strongly non-linear regime. The method, based on the dynamical evolution of the pair conservation equation, relies on a universal relation between the pair-wise velocity and the smoothed correlation function valid for high and low density models, as derived empirically from N-body simulations. An intriguing alternative relation, based on the stable-clustering hypothesis, predicts a power-law behavior of the mass correlation function that disagrees with N-body simulations but conforms well to the observed galaxy correlation function if negligible bias is assumed. The method is a useful tool for rapidly exploring a wide span of models and, at the same time, raises new questions about large scale structure formation.

<sup>2</sup><sup>2</sup>affiliationtext: Département de Physique Théorique, Université de Genève, CH-1211 Genève, Switzerland<sup>3</sup><sup>3</sup>affiliationtext: Copernicus Astronomical Center, 00-716 Warsaw, Poland   

Understanding the origin and evolution of the clustering pattern of galaxies is one of the most important goals of cosmology. Until now, this problem has been investigated using a four-fold path: (1) perturbation theory (for a review of recent advances, see Juszkiewicz & Bouchet (1996) and references therein); (2) a kinetic description, adapted from the BBGKY hierarchy, used in plasma physics (Peebles (1980), §IV); (3) N-body simulations (e.g., Jenkins et al. (1998), hereafter VIRGO); (4) and semi-analytic fits to N-body results, based on the so-called universal scaling hypothesis (see Hamilton et al. (1991); Jain et al. (1995); Peacock & Dodds (1996); Ma (1998)). The advantages and limitations of these methods are often complementary. For example, applying perturbation theory often leads to analytic results for a wide class of models while the N-body simulations allow a study of only one model at a time. On the other hand, perturbation theory works only in the weakly non-linear regime while N-body experiments describe the fully non-linear dynamics, albeit over a limited dynamical range. The subject of the present study is an analytic ansatz for the evolution of the two-point correlation function of density fluctuations spanning the linear and non-linear regime, which builds on all four methods described above.

Our approach is based on the pair conservation equation, which relates the mean (pair-weighted) relative velocity of a pair of particles to the time evolution of the correlation function in a self-gravitating gas:

$$\frac{a}{3[1+\xi (x,a)]}\frac{\overline{\xi }(x,a)}{a}=\frac{v_{12}(x,a)}{Hr},$$
(1)   
where $`a(t)`$ is the expansion factor with $`a=1`$ at present, $`r=ax`$ is the proper separation and $`H(a)`$ is the Hubble parameter (see Davis & Peebles (1977); Peebles (1980)). Here $`\overline{\xi }(x,a)`$ represents the two-point correlation function averaged over a ball of comoving radius $`x`$:

$$\overline{\xi }(x,a)=\frac{3}{x^3}_0^x\xi (y,a)y^2𝑑y.$$
(2)   
The approximate solution of (1) is known in the large separation limit, where $`|\xi |1`$ (linear regime); the stable clustering hypothesis is often invoked to describe the small separation limit, where $`\xi 1`$ (non-linear regime). Hence, equation (1) is “a guide to speculation on the behavior of the correlation function” (Peebles (1980), p.268) since an assumed $`v_{12}`$ implies a function $`\xi `$ that should agree with the weak and strong field limits, and interpolates between these limits in a reasonable way.

An approximate universal relation between the pair-wise velocity and the smoothed correlation function has been conjectured by Hamilton et al. (1991) on the basis of N-body simulation results and further explored in Nityananda & Padmanabhan (1994) and Padmanabhan & Engineer (1998). In the past, the relation has been used in attempts to derive a general functional that converts directly from a linear to a non-linear mass correlation. In this paper, we present a simple extension of the relation that applies to both high- and low-density models, but take a different approach to obtaining the non-linear correlation function. Namely, we use the universal relation to close (1); we then evolve the resulting partial differential equation to compute the non-linear correlation function. This turns out to be a fast and surprisingly accurate method that matches N-body results for a wide variety of cold dark matter (CDM) models. As a stand-alone computer program, the algorithm can be adopted on a programmable calculator; or, it can easily be incorporated in more sophisticated programs that predict other cosmological properties, such as CMBFAST (Seljak & Zaldarriaga (1996)).

Figure 1 clearly illustrates the nearly model-independent relationship between the pair-wise velocity and the smoothed correlation function, observed in N-body simulations for a wide range of perturbation spectra. We define this relation as

$$V[f\overline{\xi }]\frac{v_{12}}{Hr}.$$
(3)   
Compared to Hamilton et al. (1991), a novel feature is plotting the relation in terms of $`f(\mathrm{\Omega })\overline{\xi }`$, where $`fd\mathrm{ln}D/d\mathrm{ln}a`$ is the standard linear density growth factor, rather than $`\overline{\xi }`$ alone — a difference that is essential for extending the relation to low-density models. The evolved, non-linear clustering of scale-free spectra with $`n=1,2`$ as well as the CDM family of models produces a very similar relation between $`v_{12}/Hr`$ and $`\overline{\xi }`$. An excellent fit to the functional relation $`V[f\overline{\xi }(x,a)]`$ in Figure 1, based on the $`n=1`$ curve, is given by

$$V[x]=\{\begin{array}{cc}\frac{2}{3}x\hfill & x<0.15\hfill \\ 0.7x\mathrm{exp}(0.31x^{0.61})\hfill & 0.15x<20\hfill \\ 3.3x^{0.17}\hfill & x20\hfill \end{array}$$
(4)   
valid for $`x\stackrel{<}{}10^3`$. In this paper, we use this fitting formula, designed for the $`n=1`$ curve, as the expression for $`V[x]`$ in (1) to be applied to all models. We find this to be sufficient to reproduce N-body results to within 10% accuracy over the range of models and scales shown in Figure 2, which extends deep into the non-linear regime. To push to lower density models and improve the accuracy further, it would be simple to modify the algorithm for $`V[x]`$ to include the $`\mathrm{\Omega }`$-dependent rise near $`x\stackrel{>}{}20`$.

Starting with the linear correlation function, and armed only with information about the background cosmological evolution, we propose to dynamically obtain the non-linear $`\xi `$ as a function of separation and time. Here then is our idealized procedure in three steps.

1. Reformulate: We first rewrite the partial differential equation as

$$\frac{\mathrm{ln}\overline{\xi }}{\mathrm{ln}a}=\mathrm{\hspace{0.33em}3}\frac{(1+\xi )}{\overline{\xi }}V[f\overline{\xi }]$$
(5)   
where $`V[x]`$ is given in equation (4).

2. Initialize: The initial conditions are set by the linear correlation function at a red shift $`z=1+1/a_i`$ such that $`\xi (x,a_i)`$ is less than unity for all $`x`$ of interest. Our procedure assumes that only the amplitude, and not the shape of the linear correlation function has changed over this interval, as occurs for cold dark matter models with a primordial spectrum of adiabatic density perturbations. Hence, this procedure will not apply to cosmological models with a late-time decaying neutrino, but will apply to hot dark matter models wherein the shape of the linear power spectrum is set by red shift $`z100`$.

3. Evolve: We numerically solve the partial differential equation and evolve $`\overline{\xi }`$. At each step in the evolution, we use $`\xi =\overline{\xi }\times (1\overline{\gamma }/3)`$ with $`\overline{\gamma }d\mathrm{ln}\overline{\xi }/d\mathrm{ln}x`$ to determine the correlation function. The value of $`f`$ is updated at each step in $`a`$ as appropriate for the cosmology.

The remarkable results are shown in Figure 2. Here we see that our simple procedure gives excellent agreement with N-body simulations. Based on the span of behavior in the cosmic time evolution and the shape of correlation function, we expect this procedure should be valid for a wide range of cosmological models, including quintessence (Caldwell, Dave, & Steinhardt (1998)) and models with tilted spectra.

Figure 2 demonstrates that we have obtained a simple and powerful new tool for rapidly and accurately obtaining the non-linear power spectrum for a wide range of models. However, the physical origin of the nearly model-independent relation $`V[f\overline{\xi }]`$ is not understood in detail. In the linear regime, perturbation theory predicts $`v_{12}/Hr=(2/3)f\overline{\xi }`$. In the non-linear regime, Padmanabhan et al. (1996) have suggested that insight may be obtained by comparison to the case of the gravitational collapse of a spherical top hat mass distribution. Using their solution (eqs. (16-19) in their paper), we find that $`v_{12}/Hr`$ is linearly proportional to $`f\overline{\xi }`$ times a slowly decreasing function of $`\overline{\xi }`$ for a surprisingly wide range of $`\overline{\xi }1`$, including $`\overline{\xi }=10`$, the turnover point in Figure 1. In particular, for $`\overline{\xi }=10`$, $`v_{12}/Hr=1.77f(\mathrm{\Omega })`$, which is similar to $`v_{12}/Hr2f(\mathrm{\Omega })`$ in Figure 1.

In the strongly non-linear regime, $`\overline{\xi }10`$, Figure 1 shows a visible difference in the shape of $`V[f\overline{\xi }]`$ between the high density and low density models. This may be due to the suppression of linear growth, which occurs at late times in low density models and leads to the enhanced clustering on small scales relative to large scales. However, this has a negligible effect on the computed non-linear correlation function. For example, using the curve for $`\mathrm{\Lambda }`$CDM shown in Figure 1 as the basis of $`V`$ in our procedure, we find the amplitude of the non-linear correlation function differs by only $`10\%`$ at $`r0.1`$ Mpc/h. For the models shown, this accuracy is comparable to what is obtained by Hamilton et al. (1991) and Peacock & Dodds (1996). The advantage here is that our method can be immediately applied to new types of CDM models (e.g., quintessence cosmologies) without having to run new N-body simulations to fix fitting parameters.

An important issue raised by our ansatz is the validity of the stable clustering hypothesis. The stable clustering regime corresponds to the limit where particle pairs detach from the Hubble flow and $`v_{12}/Hr1`$. Figure 1 shows that $`v_{12}/Hr`$ first overshoots unity by a factor of two and then rebounds towards unity. However, it is not clear whether it converges to unity at $`\overline{\xi }1000`$ or possibly oscillates if the simulations are extended to higher values of $`\overline{\xi }`$.

It is interesting to compare the predictions of our ansatz if the relation between $`v_{12}/Hr`$ and $`\overline{\xi }`$ is modified to enforce more rapid convergence to stable clustering. A ready example is an alternative ansatz based on the pair conservation equation recently proposed by Juszkiewicz et al. (1999) (hereafter JSD). Their ansatz for $`v_{12}`$, based on an interpolation between the behavior predicted by perturbation theory in the weakly non-linear regime and stable clustering in the strongly non-linear regime is given by

$$v_{12}(x,a):=\frac{2}{3}Hrf\overline{\overline{\xi }}(x,a)\left[\mathrm{\hspace{0.17em}1}+\alpha \overline{\overline{\xi }}(x,a)\right]$$
(6)   
where $`\overline{\overline{\xi }}\overline{\xi }/(1+\xi )`$ and $`\alpha `$ is a function which controls the strength of the non-linear feedback. Here we use $`\alpha =1.81.1\gamma `$, based on perturbation theory (see Scoccimarro & Frieman (1996)), where $`\gamma `$ is the slope of the correlation function at $`\xi =1`$. The key point, as shown in Figure 3, is that the pair-wise velocity rapidly approaches the stable clustering limit by $`\overline{\xi }10`$, and remains there to within $`20\%`$ on smaller scales in the more strongly non-linear regime. This means particle pairs separated by $`1`$ Mpc/h have the rough behavior of virialized objects, such as clusters and galaxies.

The correlation function $`\xi `$ obtained by closing the pair conservation equation with (6), as shown in Figure 3, displays a power-law behavior in the non-linear regime with index $`1.7`$, in disagreement with N-body simulations of CDM but curiously similar to the galaxy correlation function observed in the APM survey (Maddox et al. (1996)). This result is not unique to the JSD ansatz; substituting any shape similar to that shown in the top panel of Figure 3 for $`V[f\overline{\xi }]`$ into our ansatz would produce a similar effect on the mass correlation function. In the past, the disagreement between the dark matter power spectrum observed in simulations, which shows no evidence of power-law behavior, and the simple power-law observed in the galaxy correlation function has been attributed to a scale-dependent bias. The conventional picture is that the bias function, $`b(r)`$, is just so as to cause a non-power-law behavior in the dark matter to be translated into a power-law behavior of the galaxy correlation function, $`\xi _g(r)=\xi (r)b^2(r)`$.

Our present findings concerning the dependence of the mass correlation function on the model independent $`v_{12}/Hr`$ vs. $`f\overline{\xi }`$ relation suggests a radical possibility. Perhaps the bias factor is completely negligible and, instead, CDM models are missing some important physical feature (e.g. mechanics of galaxy formation, or some property of the dark matter) which causes rapid convergence to stable clustering in the non-linear regime ($`v_{12}/Hr1`$ without substantial overshoot) and to power-law behavior of the correlation function. The ansatz illustrated in Figure 3, which enforces stable clustering by fiat, may be implicitly describing a modified CDM model of galaxy formation which incorporates the new physical feature. The difference between N-body simulations and observations, whether due to bias in the conventional picture or to something more radical, such as a modification to CDM, can perhaps be determined empirically by studying red shift distortion on the scales where N-body simulation suggests overshoot in $`v_{12}/Hr`$ and the ansatz of Figure 3 does not.

In sum, our studies have produced a simple recipe for computing the non-linear power spectrum for a wide range of models. (Upon publication of the paper, we will make a program available at feynman.princeton.edu/$``$steinh.) A key feature of the method is the universal function relating the pair velocity to the mass correlation function which does not converge rapidly to the stable clustering limit, but, rather, overshoots by a factor of two. This feature is responsible for the fact that the mass correlation function does not approach a power-law. Our studies have also raised several interesting issues in structure formation: why is $`v_{12}/Hr`$ vs. $`f\overline{\xi }`$ so similar for a wide range of models? can the universal relation be derived from theory? for what range of models is the relation model-independent? how does the universal relation ultimately approach the stable clustering limit at small scales (if it does at all)? and, does the success of a universal relation based on the stable clustering hypothesis, as in Figure 3, suggest a viable, alternative explanation for the power-law behavior of the galaxy correlation function?

###### Acknowledgements.

We would like to thank Dick Bond, Marc Davis, Josh Frieman, Andrew Hamilton, Chung-Pei Ma, Jim Peebles, David Spergel, and Roman Scoccimarro for useful conversations. We also thank Bhuvnesh Jain and Volker Springel and the VIRGO collaboration for providing N-body simulation results. This work was carried out, in part, at the Isaac Newton Institute for Mathematical Sciences during their program on Structure Formation in the Universe. We would like to thank the organizers, N. Turok and V. Rubakov, and the staff of the Institute for their support and kind hospitality. The work of RRC and PJS was supported by the US Department of Energy grant DE-FG02-91ER40671. The work of RJ was supported by grants from the Polish Government (KBN grants No. 2.P03D.008.13 and 2.P03D.004.13), the Tomalla Foundation of Switzerland, by the Poland-US M. Skłodowska-Curie Fund. FRB and RJ were supported by the Franco-Polish collaboration grant (Jumelage).