# Higher order moments of the density field in a parameterized sequence of non-gaussian theories

## 1 Introduction

The standard paradigm for the formation of large scale structure is that quantum fluctuations during an inflationary epoch seeded initially small, gaussian fluctuations in density, which grew through the action of gravitational instability in a universe whose dominant constituent is cold dark matter (CDM). This theory has proved very predictive and agrees with a wide range of observational data, however several of the assumptions are difficult to test with high precision. The assumption that the initial fluctuations have a gaussian distribution holds more generally than just in the inflationary CDM model, as the central limit theorem implies that fluctuations emerging from many uncorrelated, random processes will be nearly gaussian. Though it is indeed a plausible assumption, and a prediction in the simplest inflationary models, it has received only limited observational support through measurements of the CMB  and large-scale structure . In the latter case the situation is made more difficult by the action of gravity which turns an initially gaussian random field into a non-gaussian field once the modes become non-linear.

For gaussian fluctuations the only non-trivial moment is the 2-point correlation function, $`\xi (r)`$, or its Fourier Transform the power spectrum $`P(k)`$. If the fluctuations are non-gaussian the higher-order moments of the field carry additional information. The evolution of the higher order moments induced by gravitational instability in an initially gaussian random field has been studied extensively (see Strauss & Willick  for a review). Of particular interest are the $`S_N`$, defined in terms of the volume averaged correlation functions (cumulants of the probability distribution function)

$$\overline{\xi }_N\frac{1}{V^N}_Vd^3r_1d^3r_2\mathrm{}d^3r_N\xi (𝐫_1,\mathrm{},𝐫_N)$$
(1)   
as

$$\overline{\xi }_N=S_N\overline{\xi }_2^{N1}.$$
(2)   
In the mildly non-linear regime the growth of gaussian initial conditions by gravitational instability predicts that the $`S_N`$ are independent of scale  for a scale-free spectrum. For a top-hat filter of the density field and a pure power-law spectrum, $`P(k)k^n`$, it can be shown  that to lowest non-trivial order in perturbation theory ($`\overline{\xi }_21`$),

$`S_3`$ $`=`$ $`{\displaystyle \frac{34}{7}}(n+3)`$ (3)   
$`S_4`$ $`=`$ $`{\displaystyle \frac{60712}{1323}}{\displaystyle \frac{62}{3}}(n+3)+{\displaystyle \frac{7}{3}}(n+3)^2`$ (4)   
with a very weak dependence on the matter density, $`\mathrm{\Omega }_0`$ . These results have been confirmed by N-body simulations  and are known to be insensitive to redshift space distortions  on the mildly non-linear scales where we will be working (see §3).

Measurements of the $`S_N`$ have been performed on several galaxy catalogues in two and three dimensions; see e.g. Table 1 in Hui & Gaztanaga  or Kim & Strauss . We expect that the Anglo-Australian Two Degree Field survey (2dF<sup>1</sup><sup>1</sup>1http://meteor.anu.edu.au/$``$colless/2dF) and the Sloan Digital Sky Survey (SDSS<sup>2</sup><sup>2</sup>2http://www.sdss.org) will provide excellent catalogues for estimation of the $`S_N`$ for different sub-samples to high order. Assuming growth through gravitational instability from initially small gaussian fluctuations one can use the additional information contained in the $`S_N`$ to relate the properties of the observed galaxy distribution to the fluctuations in the underlying dark matter. The observations of CMB anisotropies have enhanced our faith in gravity as the engine of growth. The extra assumption required in this step is that the initial conditions were gaussian.

As with other assumptions in the standard paradigm, our assumption of gaussianity for the initial fluctuations should be tested against observations. One difficulty with testing for non-gaussianity has been the lack of a simple, predictive theory which describes what form, out of the infinite possibilities, the non-gaussianity should take. While inflationary models generically predict gaussian fluctuations, not all inflationary models are “generic” in this sense. And of course there are other, non-inflationary, models of structure formation which predict non-gaussian perturbations. The best known of these non-gaussian theories are models based on topological defects. However they have problems fitting the observational data and are difficult to model properly requiring expensive simulations even for the initial conditions.

Moments have been predicted for a variety of non-gaussian models under various approximations. Results for the unsmoothed moments $`S_3`$ and $`S_4`$ for arbitrary initial conditions in perturbation theory have been presented by . Approximations to defect models have been studied by . The latter authors also discussed an approximate calculation of the low-order moments for the model of Peebles \[1999a, 1999b\], which shall discuss further below. Weinberg & Cole  studied models obtained from non-linear mappings of initially Gaussian fields. Perhaps the closest ancestor to this work, however, is that of , who used N-body simulations to calculate moments for a variety of non-gaussian models, including a $`\chi _1^2`$ model similar to that developed below.

In this paper we investigate the predictions for the lowest moments, $`S_3`$ and $`S_4`$, in a simple, parameterized, non-gaussian “model”. The model has useful interpolating properties. It has one parameter, $`\nu `$. The limit $`\nu \mathrm{}`$ recovers the gaussian result, with lower values of $`\nu `$ being progressively more non-gaussian. Low values of $`\nu `$ mimic cosmological models with rare, high-amplitude peaks.

The outline of the paper is as follows: in the next section we discuss the non-gaussian model in more detail, in §3 we discuss our technique for calculating the moments $`S_3`$ and $`S_4`$ using an N-body code, plus the tests we have run. In §4 we summarize the results, and in §5 we discuss the implications.

## 2 The Non-Gaussian Model

The COBE team  first introduced a non-gaussian model where the real and imaginary parts of the Fourier- (or in their case $`\mathrm{}`$-) space perturbations, are drawn from independent $`\chi _\nu ^2`$ distributions with $`\nu `$ degrees of freedom, adjusted to have zero mean and scaled to give the right power spectrum. While the precise physical significance of this model is unclear, it serves as a well defined reference with useful interpolating properties, as the $`\chi _\nu ^2`$ distribution becomes more gaussian as $`\nu \mathrm{}`$. It was shown in  that the COBE data prefer $`\nu \mathrm{}`$, with the gaussian model being 5 times more likely than any other model tested.

We have modified this model to specify $`\chi _\nu ^2`$ distributions in real rather than $`k`$-space. We proceeded as follows: Let $`\varphi _i(x)`$ be independent, gaussian random fields of zero mean. Define, for integer $`\nu `$,

$$\psi \frac{1}{\nu }\underset{i=1}{\overset{\nu }{}}\left(\varphi _i^2\varphi _i^2\right).$$
(5)   
Then $`\psi `$ will be $`\chi _\nu ^2`$ distributed, in real-space, with zero mean. For $`\nu \mathrm{}`$ the fluctuations will be gaussian. As $`\nu `$ is lowered the model becomes progressively more non-gaussian, with more “rare” peaks. The limiting case $`\nu =1`$ has the perturbations quadratic in a gaussian field, much as in recently proposed non-gaussian inflation models \[1997, 1997, 1999a, 1999b\]. The reduced moments of the initial conditions are easily calculated for a $`\chi _\nu ^2`$ distribution. The mean is $`\nu `$, and

$`(\chi ^2\nu )^2`$ $`=`$ $`2\nu ,`$ (6)   
$`(\chi ^2\nu )^3`$ $`=`$ $`8\nu ,`$ (7)   
$`(\chi ^2\nu )^4`$ $`=`$ $`48\nu +3(2\nu )^2.`$ (8)   
A final feature of this model, not shared by other proposals in the literature which involve non-linear mappings of Gaussian density fields, is that the distribution of fluctuations will remain approximately the same when smoothed on a variety of length scales. This is related to the scaling of this model: $`\psi ^n`$ scales as $`\psi ^2^{n/2}`$. We show this feature explicitly for the case $`\nu =5`$ in Fig. 1, where the distribution of initial fluctuations from one of our simulations is given for different smoothing scales, spanning a factor of 5 in scale.

It is straightforward to generate initial conditions for this model with a power-law spectrum. For $`i=1,\mathrm{},\nu `$ we generate independent realizations of a gaussian field $`\varphi _i(x)`$, with correlation function $`\xi _i(r)`$, on a grid in the usual way. We then square $`\varphi _i`$ and sum over $`i`$ to obtain $`\psi `$, which we use as an initial field to generate displacements using the Zel’dovich approximation as described below.

The only remaining step is to choose the correlation functions $`\xi _i(r)`$ so that the final correlation function $`\xi _\psi (r)`$ has a given form. Writing $`\delta \varphi _i^2\varphi _i^2\varphi _i^2`$ and using the fact that $`\varphi _i`$ is a gaussian one has

$$\delta \varphi _i^2(𝐱_a)\delta \varphi _i^2(𝐱_b)=2\xi _\varphi ^2(𝐱_a𝐱_b).$$
(9)   
Note that since the $`\varphi _i`$ are independent, the correlation functions for the components of $`\psi `$ simply add. For simplicity we choose them to be equal.

It is important to notice that $`\xi _\psi 0`$ for all separations, so the “integral constraint” (i.e. the integral of $`\xi (r)`$ over all space must vanish) cannot hold for this spectrum. However we shall be interested in power-law spectra for which $`\xi 0`$ anyway, so this technical point will not affect us. We shall choose the power law index of the $`\varphi `$ field $`n_\varphi =2`$ so that $`n_\psi =2n_\varphi +3=1`$, close to the slope observed for CDM models on the scales of interest. In the absence of finite box size corrections, we get $`n_\psi =1`$ by choosing

$$\mathrm{\Delta }_\varphi ^2(k)=\mathrm{\Delta }_{}^2\left(\frac{k}{k_{}}\right),$$
(10)   
where $`\mathrm{\Delta }_\varphi ^2k^3P_\varphi (k)/(2\pi ^2)`$ is the power spectrum of the fields $`\varphi _i`$. This gives $`\xi _\varphi (r)=(\pi /2)\mathrm{\Delta }_{}^2(k_{}r)^1`$ or

$$\mathrm{\Delta }_\psi ^2(k)=\frac{2}{\nu }\left(\frac{\pi }{2}\right)^2\mathrm{\Delta }_{}^4\left(\frac{k}{k_{}}\right)^2.$$
(11)   
We set $`k_{}=0.2h\mathrm{Mpc}^1`$ and adjust $`\mathrm{\Delta }_{}`$ so that $`\mathrm{\Delta }_\psi ^2(k_{})=0.25`$, or $`\sigma _8\frac{1}{2}`$. The same normalization is used in the test simulations described below.

We have checked whether the resulting power spectra for our non-gaussian models match the $`n_\psi =1`$ form we expect. Since the evolution is a non-linear function of the initial conditions, it is important to do this comparison on the initial conditions rather than the evolved spectrum. We find that there is a departure from a pure power-law at both high-$`k`$ and low-$`k`$. This is to be expected: $`\mathrm{\Delta }_\psi ^2`$ is given by an auto-correlation of $`\mathrm{\Delta }_\varphi ^2`$, and in the simulations $`\mathrm{\Delta }_\varphi ^2`$ is missing modes with wavelengths longer than the box or smaller than the mesh scale. We have corrected for this effect by modifying the input spectrum of the $`\varphi `$ field to ensure an initial power-law spectrum for $`\mathrm{\Delta }_\psi ^2`$. While $`\mathrm{\Delta }_\psi ^2`$ changes by more than 3 orders of magnitude over the range of scales we simulate, the correction factor we apply always lies between 1 and 3 (in power), with the largest correction at low-$`k`$. Note that we are not directly modifying the initial power spectrum of the $`\psi `$ field here – it is the convolution integral in the $`\varphi `$ field that requires this numerical correction.

To sum up: The procedure outlined above provides an algorithm for generating realizations of a random field whose 1-point function is a $`\chi _\nu ^2`$ (as has been explicitly checked in the simulations, see §4) and whose 2-point function has a power-law spectrum to within 1% over more than an order of magnitude in scale. We show this explicitly for $`\nu =1`$ in Fig. 2. This initial realization of the $`\psi `$ field can then be evolved using N-body simulations to study the effect of gravitational growth on these initial conditions, as described in the next section.

## 3 Calculating the moments

We use $`N`$-body simulations to calculate the moments $`S_3`$ and $`S_4`$ in the quasi-linear regime. A discussion of several of the relevant numerical issues can be found in . We use a particle-mesh (PM) code, described in detail in Meiksin, White & Peacock . We have simulated critical density universes ($`\mathrm{\Omega }_0=1`$) with power-law spectra in boxes of size $`150h^1\mathrm{Mpc}L_{\mathrm{box}}250h^1`$Mpc on a side, so that the fundamental mode was always well in the linear regime and the box represented a “fair” sample of the universe. To minimize finite volume effects, we will work on scales less than about $`0.15R_{\mathrm{box}}`$, where $`L_{\mathrm{box}}^3=(4\pi /3)R_{\mathrm{box}}^3`$ (see below).

We have run simulations with either $`64^3`$ or $`128^3`$ particles and a $`64^3`$ or $`128^3`$ force “mesh”. All the simulations were started at $`1+z=20`$ to 30 and run to the present ($`z=0`$). The evolution was done in log of the scale factor $`a=(1+z)^1`$. The time step was dynamically chosen as a small fraction of the inverse square root of the maximum acceleration, with an upper limit of $`\mathrm{\Delta }a/a=4`$ per cent per step. This resulted in a final particle position error of less than 0.1 per cent of the box size. As described in Meiksin et al. , the code reproduces the non-linear power spectrum excellently as compared with semi-analytic fitting functions (e.g. Peacock & Dodds ) or other N-body codes. The initial conditions were generated from an initial density field, with either an $`N_{\mathrm{part}}`$ or $`N_{\mathrm{mesh}}`$ FFT, using the Zel’dovich approximation. The particles were initially placed either on a uniform grid, or at random within cells in a structure with $`N_{\mathrm{part}}^{1/3}`$ cells to a side, as described in Peacock & Dodds . Tests indicate the most reliable conditions are when the resolution imposed by finite particle number matches the force resolution from the mesh and when the particles were displaced from random positions within cells, i.e. the conditions used in the simulations described in Meiksin et al. .

At the end of each simulation $`8,192`$ sphere centers were thrown down at random within the volume and the number of particles in each of multiple concentric top-hat spheres computed. The radii of the cells was constrained to be larger than $`2.5`$ (mesh) cells and smaller than $`0.2R_{\mathrm{box}}`$. Since the smallest sphere is more than $`5`$ grid cells across, the effects due to the finite resolution of the PM code are minimal. An explicit check with the higher resolution simulations verified this expectation, though to achieve per cent level accuracy we need to work on scales larger than 4 mesh cells (8 in diameter) for the larger boxes where the grid scale is less non-linear. This is also to be expected: in the non-linear regime the structure on the grid scale is determined mainly by the collapse of larger wavelength modes, which should be well evolved by the PM code at all times . Thus as the grid scale becomes more non-linear, the effects of the finite force grid on perturbations a few grid cells across are lessened. We have been careful to ensure that the non-linear scale at $`z=0`$ has a wavelength of several mesh zones. Comparison with different box sizes and resolutions suggested the larger radius cells, with radii approaching $`0.2R_{\mathrm{box}}`$, were affected by the finite box size, with $`S_4`$ being more affected than $`S_3`$. For this reason we restricted the largest cell in the analysis to be less than $`0.15R_{\mathrm{box}}`$. Analytic arguments  suggest that our boxes are large enough to avoid finite volume biases at the few per cent level on these scales.

The moments of the counts-in-cells distribution were calculated and averaged over many ($`>100`$) simulations with different realizations of the same initial spectrum and statistics. Writing

$$\mu _M=\left(\frac{N\overline{N}}{\overline{N}}\right)^M$$
(12)   
where $`\overline{N}=N`$ and $`\mathrm{}`$ represents an ensemble average (which we approximate by an average over the simulations) we have 

$`\overline{\xi }_2`$ $`=`$ $`\mu _2\overline{N}^1`$ (13)   
$`\overline{\xi }_3`$ $`=`$ $`\mu _33\mu _2\overline{N}^1+2\overline{N}^2`$ (14)   
$`\overline{\xi }_4`$ $`=`$ $`\mu _46\mu _3\overline{N}^13\mu _2^2+11\mu _2\overline{N}^26\overline{N}^3.`$ (15)   
Our estimates of $`\overline{\xi }_N`$ are used to calculate $`S_3`$ and $`S_4`$, with the errors (including correlations) estimated from the run-to-run scatter and propagated in the usual manner. Note we average the $`\mu _M`$ over the simulations and then compute the $`S_N`$ rather than averaging the $`S_N`$ computed from each simulation. This avoids biases introduced because $`x/yx/y`$ . From the work of Hui & Gaztanaga  we estimate that their so called “estimation-biases” are $`𝒪(1\%)`$ for the cases of interest here. We show in Fig. 3 the estimates of $`S_4`$ from the cells thrown for a set of simulations with gaussian initial conditions and power-law index $`n=1/2`$. The average of these values is shown as the dashed line and the $`S_4`$ computed from the entire suite of simulations as described above is the final (solid) point with error bar. This value agrees well with the perturbation theory expectation – the solid line – whereas the naive average does not. In this Figure the error bars are overestimated for each realization by the small number of cells thrown, but this error is reduced by the large number of realization used.

An alternative to the counts-in-cells method described above, is to perform a maximum likelihood fit to the entire distribution function of the counts . This method minimizes errors due to finite volume effects and shot-noise. However it requires one to know a priori a functional form for the distribution function. The Edgeworth expansion, used by Kim & Strauss , is a valid expansion only near the peak of the distribution. To remove unphysical oscillations it must be regularized. Kim & Strauss  remove the oscillations by convolving the Edgeworth expansion with a Poisson distribution, to account for the very sparse sampling they performed. However once the number of particles in a cell becomes appreciably greater than unity this convolution no longer regularizes the distribution and a different expansion is needed. Typically models for the distribution of counts involve many higher order moments, beyond $`S_4`$, making implementation of this method difficult. For these reasons we have chosen to use the more traditional moments method. Another alternative to throwing cells  is to interpolate the final density field on a grid, and compute moments of the smoothed field by summing over the grid points (the smoothing is done using FFT methods). For our implementation – using top-hat cells and small grids – this method, while significantly faster, was not as accurate. For fine enough grids and smooth interpolation (e.g. CIC or TSC) we would expect to obtain the same results as the direct counts-in-cells method, but we did not investigate this in detail. Finally the method of Szapudi  for throwing an effectively infinite number of cells was too CPU intensive. For the low order moments we are considering, the traditional method allows sufficient cells to be thrown that Szapudi’s method is not required.

As a test we first simulated gaussian initial conditions with pure power-law spectra with $`n=1`$, $`n=0.5`$ and $`n=0`$, for which the results are known analytically (see Eqs. 3, 4) when $`\overline{\xi }_21`$. The agreement between the numerical and analytic results, shown in Fig. 4, was very good, showing that numerical effects were under control. It is interesting to note that for $`n0`$ the perturbation theory results are obtained only when $`\overline{\xi }_2`$ is quite small, in contrast to the case where $`n1`$, for which the perturbation theory results are a good approximation even for $`\overline{\xi }_20.1`$. This agrees with the results of , see their Fig. 5. For $`n=1`$ (our fiducial model) the numerical results and analytic predictions agreed to better than 5 per cent, for both $`S_3`$ and $`S_4`$, when $`\overline{\xi }_21`$. This was the case for either gaussian initial conditions, or runs with the non-gaussian model of §2 with $`\nu 1`$. This is an important result, since it argues that the simulations reported in  can be safely used to estimate correlations between power spectrum bins in the mildly non-linear regime .

Having established that the code performs as expected, we simulated models with $`n_\psi =1`$ as described in §2 for varying values of $`\nu `$. The results are described in the next section. Before leaving these tests however, we show in Fig. 5 the values of $`S_3`$ and $`S_4`$ for the $`n=1`$ spectra with gaussian initial conditions, but varying $`\mathrm{\Omega }_m`$ and in redshift rather than real space. Our redshift space results use the distant observer approximation, i.e. for each of the simulations we added the $`z`$-component of the velocity, in units of the Hubble constant, to the $`z`$-component of the position before throwing the cells. These results confirm that redshift space distortions and variations in $`\mathrm{\Omega }_m`$ have a negligible impact upon the values of $`S_3`$ and $`S_4`$ that we calculate .

## 4 Results

Our main results are shown in Fig. 6. Here we plot $`S_3(10h^1\mathrm{Mpc})`$ and $`S_4(10h^1\mathrm{Mpc})`$ as a function of $`\nu `$. We have chosen this scale because, for our normalization, $`\overline{\xi }_20.1`$ at $`10h^1\mathrm{Mpc}`$. For $`n=1`$ and gaussian initial conditions we would expect $`S_32.857`$ and $`S_413.89`$ in the weakly non-linear regime. These values are shown as the dotted lines in Fig. 6. Note that $`S_3`$ and $`S_4`$ increase rapidly as the model is made more non-gaussian.

Once the model is non-gaussian we no longer have reason to expect that the $`S_N`$ will be independent of scale. However non-linear effects process the initial non-gaussianity in a non-trivial way, as seen in Fig. 6. Simple scaling arguments from the initial conditions, see Eqs. (13-15), would suggest that $`\overline{\xi }_3^{\mathrm{init}}=\sqrt{8/\nu }\overline{\xi }_2^{3/2}`$ and $`\overline{\xi }_4^{\mathrm{init}}=(12/\nu )\overline{\xi }_2^2`$. This is to be compared to the gaussian prediction that $`\overline{\xi }_3^{\mathrm{gauss}}\overline{\xi }_2^2`$ and $`\overline{\xi }_4^{\mathrm{gauss}}\overline{\xi }_2^3`$. Defining $`S_N^{}`$ such that $`S_N^{}`$ would be independent of $`\nu `$ and $`\overline{\xi }_2`$ for a $`\chi _\nu ^2`$ field we have

$`S_3^{}`$ $``$ $`S_3(\overline{\xi }_2\nu )^{1/2}`$ (16)   
$`S_4^{}`$ $``$ $`S_4(\overline{\xi }_2\nu )`$ (17)   
for unsmoothed fields. For a pure $`\chi _\nu ^2`$ field, $`S_3^{}=\sqrt{8}2.82`$ and $`S_4^{}=12`$. By calculating the moments of the initial density field from which the Zel’dovich displacements are generated, we have explicitly checked that these values are in fact attained in the initial conditions on the grid (i.e. without any smoothing). Any deviation is at the 1 per cent level for $`S_3^{}`$ and $`S_4^{}`$ over a wide range of $`\nu `$. We show the evolved values of $`S_3^{}`$ and $`S_4^{}`$ vs $`\nu `$ in Fig. 7, with the gaussian point at $`\nu =\mathrm{}`$ omitted. Notice that gravity has modified the initial conditions so that the $`S_N^{}`$ are no longer independent of $`\nu `$. The curves rise to larger $`\nu `$, indicating that the gravitationally induced skewness and kurtosis are a larger fraction of the initial conditions as the model becomes more gaussian.

In addition, the effects of non-linearity change the scaling with $`\overline{\xi }_2`$. For example, for the $`\nu =1`$ model the $`S_N`$ are much less dependent on $`\overline{\xi }_2`$ than the $`S_N^{}`$, varying by $`𝒪(10\%)`$ for $`0.03\overline{\xi }_20.3`$. This behaviour can be explained as the gravitational contribution becoming more important than the initial conditions as clustering evolves, as anticipated by Fry & Scherrer . This result suggests that one must use care in inferring from the scale-independence of the $`S_N`$ that the initial fluctuations are gaussian. To properly make this inference it may be necessary to use very large scale or high redshift measurements.

## 5 Conclusions

We have used N-body simulations to calculate the low order moments of the mass density field for a sequence of non-gaussian models with power-law spectra. We obtain good agreement with the predictions of perturbation theory for gaussian initial conditions. As the real-space perturbations are made progressively more non-gaussian the moments depart significantly from the values of Eqs. (3,4) as shown in Fig. 6. Moreover, gravitational evolution modifies the scaling of the $`S_N`$ with $`\overline{\xi }_2`$ and $`\nu `$. For example the scale-dependence of the $`S_N`$ does not follow that of the initial conditions over the range the simulations probe: $`0.03\overline{\xi }_20.3`$.

## Acknowledgments

I would like to thank Avery Meiksin and John Peacock for their help in developing the PM code used in this paper and Lam Hui and Michael Strauss for useful conversations on measuring $`S_3`$ and $`S_4`$ from the simulations. I am grateful to Pablo Fosalba, Josh Frieman, Bob Scherrer and David Weinberg for discussions on non-gaussianity, Andrew Liddle for enlightening conversations on inflationary models and Joanne Cohn for useful comments on the manuscript. M.W. is supported by the NSF.