Higher order moments of the density field in a parametrized sequence of non-Gaussian theories

Abstract

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, the dominant constituent of which 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 cosmic microwave background (CMB) (Kogut et al. 1996; Heavens 1998) and large-scale structure (Bouchet et al. 1993; Feldman, Kaiser & Peacock 1994; Gaztanaga 1994; Nusser, Dekel & Yahil 1994; Stirling & Peacock 1996; Colley 1997). 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 two-point correlation function, ξ(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 (1995) 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),

(1)

as

(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 (Fry 1984; Bernardeau 1992) for a scale-free spectrum. For a top-hat filter of the density field and a pure power-law spectrum, P(k)∝kⁿ, it can be shown (Peebles 1980; Goroff et al. 1986; Juszkiewicz, Bouchet & Colombi 1993; Bernardeau 1994) that to lowest non-trivial order in perturbation theory ,

(3)

(4)

with a very weak dependence on the matter density, Ω_m (Bouchet et al. 1995). These results have been confirmed by N-body simulations (Coles & Frenk 1991; Bouchet et al. 1992; Weinberg & Cole 1992; Lahav et al. 1993; Juszkiewicz et al. 1993; Lucchin et al. 1994; Bernardeau 1994; Lokas et al. 1995; Baugh, Gaztanaga & Efsatathiou 1995; Juszkiewicz et al. 1995; Colombi, Bouchet & Hernquist 1996) and are known to be insensitive to redshift-space distortions (Lahav et al. 1993; Colombi et al. 1996; Hivon et al. 1995) on the mildly non-linear scales where we will be working (see Section 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 (1998) or Kim & Strauss (1998). We expect that the Anglo-Australian Two Degree Field survey (2dF)¹ and the Sloan Digital Sky Survey (SDSS)² will provide excellent catalogues for estimation of the S_N for different subsamples 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 that 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. Of course there are also other, non-inflationary, models of structure formation that 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₃ and S₄ for arbitrary initial conditions in perturbation theory have been presented by Fry & Scherrer (1994) and Chodorowski & Bouchet (1996). Approximations to defect models have been studied by Jaffe (1994), Gaztanaga & Mahonen (1996) and Gaztanaga & Fosalba (1998). The latter authors also discussed an approximate calculation of the low-order moments for the model of Peebles (1999a; 1999b), which we shall discuss further below. Weinberg & Cole (1992) studied models obtained from non-linear mappings of initially Gaussian fields. Perhaps the closest ancestor to this work, however, is that of Coles et al. (1993), who used N-body simulations to calculate moments for a variety of non-Gaussian models, including a model similar to that developed below.

In this paper we investigate the predictions for the lowest moments, S₃ and S₄, in a simple, parametrized, non-Gaussian ‘model’. The model has useful interpolating properties. It has one parameter, ν. The limit ν→∞ recovers the Gaussian result, with lower values of ν being progressively more non-Gaussian. Low values of ν 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, and in Section 3 we discuss our technique for calculating the moments S₃ and S₄ using an N-body code, plus the tests we have run. In Section 4 we summarize the results, and in Section 5 we discuss the implications.

2 The non-Gaussian model

The COBE team (Kogut et al. 1996) first introduced a non-Gaussian model, where the real and imaginary parts of the Fourier- (or in their case ℓ-) space perturbations are drawn from independent distributions with ν 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 distribution becomes more Gaussian as ν→∞. It was shown in Kogut et al. (1996) that the COBE data prefer ν→∞, with the Gaussian model being five times more likely than any other model tested.

We have modified this model to specify distributions in real rather than k-space. We proceeded as follows: Let φ_i(x) be independent, Gaussian random fields of zero mean. Define, for integer ν,

(5)

Then ψ will be distributed, in real space, with zero mean. For ν→∞ the fluctuations will be Gaussian. As ν is lowered the model becomes progressively more non-Gaussian, with more ‘rare’ peaks. The limiting case ν=1 has the perturbations quadratic in a Gaussian field, much as in recently proposed non-Gaussian inflation models (Antoniadis, Mazur & Mottola 1997; Linde & Mukhanov 1997; Peebles 1999a; Peebles 1999b). The reduced moments of the initial conditions are easily calculated for a distribution. The mean is ν, and

〈(χ²-ν)²〉=2ν,

(6)

〈(χ²-ν)³〉=8ν,

(7)

〈(χ²-ν)⁴〉=48ν+3(2ν)².

(8)

A final feature of this model, not shared by other proposals in the literature that 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: 〈ψⁿ〉 scales as 〈ψ²〉^n/2. We show this feature explicitly for the case ν=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.

Figure 1.

The distribution of the initial density fluctuations for ν=5, smoothed with a Gaussian filter of R=1, 3 and 5 grid cells in a 64³ simulation. Note that the distribution is skew on all scales.

Open in new tabDownload slide

It is straightforward to generate initial conditions for this model with a power-law spectrum. For i=1,…,ν we generate independent realizations of a Gaussian field φ_i(x), with correlation function ξ_i(r), on a grid in the usual way. We then square φ_i and sum over i to obtain ψ, 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 ξ_i(r) so that the final correlation function ξ(r) has a given form. Writing and using the fact that φ_i is a Gaussian, one has

(9)

Note that since the φ_i are independent, the correlation functions for the components of ψ simply add. For simplicity we choose them to be equal.

It is important to notice that ξ0 for all separations, so the ‘integral constraint’ [i.e. the integral of ξ(r) over all space must vanish] cannot hold for this spectrum. However, we are interested in power-law spectra for which ξ0 anyway, so this technical point will not affect us. We shall choose the power-law index of the φ field n=−2 so that n=2n+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=−1 by choosing

(10)

where is the power spectrum of the fields φ_i. This gives , or\

(11)

We set k_*=0.2 h Mpc⁻¹ and adjust Δ_* so that , or σ₈≃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=−1 form we expect. As 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 and low k. This is to be expected: is given by an autocorrelation of , and in the simulations 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 φ field to ensure an initial power-law spectrum for . While 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 ψ field here — it is the convolution integral in the φ field that requires this numerical correction.

To sum up: the procedure outlined above provides an algorithm for generating realizations of a random field, the one-point function of which is a (as has been explicitly checked in the simulations, see Section 4) and the two-point function of which has a power-law spectrum to within 1 per cent over more than an order of magnitude in scale. We show this explicitly for ν=1 in Fig. 2. This initial realization of the ψ 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.

Figure 2.

The power spectrum of the initial conditions for ν=1 generated by the method described in the text. The solid line is the n=−1 target spectrum and the squares indicate the average power spectrum from 50 realizations of the ICs. The deviation is less than 1 per cent over more than an order of magnitude in scale.

Open in new tabDownload slide

3 Calculating the moments

We use N-body simulations to calculate the moments S₃ and S₄ in the quasi-linear regime. A discussion of several of the relevant numerical issues can be found in Baugh et al. (1995), Juszkiewicz et al. (1995), Colombi et al. (1996), Szapudi & Colombi (1996), Kim & Strauss (1998) and Jain & Bertschinger (1998). We use a particle—mesh (PM) code, described in detail in Meiksin, White & Peacock (1999). We have simulated critical-density universes (Ω_m=1) with power-law spectra in boxes of size 150 h⁻¹ Mpc≤L_box≤250 h⁻¹ 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_box, where (see below).

We have run simulations with either 64³ or 128³ particles and a 64³ or 128³ force ‘mesh’. All the simulations were started at 1+z=20 to 30 and run to the present (z=0). The evolution was carried out in log of the scalefactor a=(1+z)⁻¹. The time-step was dynamically chosen as a small fraction of the inverse square root of the maximum acceleration, with an upper limit of Δaa=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. (1999), the code reproduces the non-linear power spectrum excellently, compared with semi-analytic fitting functions (e.g. Peacock & Dodds (1996)) or other N-body codes. The initial conditions were generated from an initial density field, with either an N_part or N_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 cells to a side, as described in Peacock & Dodds (1996). Tests indicate that 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. (1999).

At the end of each simulation 8 192 sphere centres were thrown down at random within the volume and the number of particles in each of multiple concentric top-hat spheres was computed. The radii of the cells were constrained to be larger than ≃2.5 (mesh) cells and smaller than 0.2R_box. As the smallest sphere is more than 5 grid cells across, the effects caused by 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 (Bouchet, Adam & Pellat 1985). 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 that the larger radius cells, with radii approaching 0.2R_box, were affected by the finite box size, with S₄ being more affected than S₃. For this reason we restricted the largest cell in the analysis to be less than 0.15R_box. Analytic arguments (Hui & Gaztanaga 1998) suggest that our boxes are large enough to avoid finite-volume biases at the level of a few per cent 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

(12)

where and 〈…〉 represents an ensemble average (which we approximate by an average over the simulations), we have (Peebles 1980)

(13)

14

15

Our estimates of ξ_N are used to calculate S₃ and S₄, with the errors (including correlations) estimated from the run-to-run scatter and propagated in the usual manner. Note that we average the μ_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 〈xy〉≠〈x〉/〈y〉 (Hui & Gaztanaga 1998). From the work of Hui & Gaztanaga (1998), we estimate that their so called ‘estimation-biases’ are O(1 per cent) for the cases of interest here. We show in Fig. 3 the estimates of S₄ 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₄ 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 realizations used.

Figure 3.

Estimates of S₄(R=10 h⁻¹ Mpc) from the cells thrown for each of a set of n=−1/2 simulations as described in the text. The average of these values is the dashed line, and the S₄ computed from the entire suite of simulations before averaging is the rightmost (solid) point with error bar. This value agrees well with the perturbation theory expectation, shown as the solid line.

Open in new tabDownload slide

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 (Kim & Strauss 1998). This method minimizes errors resulting from 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 (1998), is a valid expansion only near the peak of the distribution. To remove unphysical oscillations it must be regularized. Kim & Strauss (1998) 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₄, making implementation of this method difficult. For these reasons we have chosen to use the more traditional moments method. Another alternative to throwing cells (Lokas et al. 1995) 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 (1998) 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 equations 3 and 4) when . 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 n≃0 the perturbation theory results are obtained only when ξ₂ is quite small, in contrast to the case where n≃-1, for which the perturbation theory results are a good approximation even for . This agrees with the results of Colombi et al. (1996); 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₃ and S₄, when . This was the case for either Gaussian initial conditions or runs with the non-Gaussian model of Section 2 with ν≫1. This is an important result, because it argues that the simulations reported in Meiksin et al. (1999) can be safely used to estimate correlations between power-spectrum bins in the mildly non-linear regime (Meiksin & White 1999).

Figure 4.

The simulation results for S₃ and S₄ for three power-law spectra with n=−1 (solid squares), n=−1/2 (open circles) and n=0 (open triangles) as a function of ξ₂. The perturbation theory results, valid for ξ₂<1, are shown as the horizontal dashed lines. Error bars (shown for S₄) are calculated from the run-to-run scatter in the simulations. For S₃ the 1σ error bars are approximately the same size as the points themselves. The error bars on the points are correlated.

Open in new tabDownload slide

Having established that the code performs as expected, we simulated models with n=−1 as described in Section 2 for varying values of ν. The results are described in the next section. Before leaving these tests, however, we show in Fig. 5 the values of S₃ and S₄ for the n=−1 spectra with Gaussian initial conditions, but varying Ω_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 Ω_m have a negligible impact upon the values of S₃ and S₄ that we calculate (Lahav et al. 1993; Hivon et al. 1995; Colombi et al. 1996).

Figure 5.

The value of S₃ and S₄ from the n=−1 simulations in real space for a critical-density universe (solid squares), an Ω_m=1-Ω_Λ=0.3 universe (open circles) and in redshift space in the Λ model (open triangles).

Open in new tabDownload slide

4 Results

Our main results are shown in Fig. 6. Here we plot S₃(10 h⁻¹ Mpc) and S₄(10 h⁻¹ Mpc) as a function of ν. We have chosen this scale because, for our normalization, ≃0.1 at 10 h⁻¹ Mpc. For n=−1 and Gaussian initial conditions we would expect S₃≃2.857 and S₄≃13.89 in the weakly non-linear regime. These values are shown as the dotted lines in Fig. 6. Note that S₃ and S₄ increase rapidly as the model is made more non-Gaussian.

Figure 6.

The moments, S₃ and S₄, as a function of ν for the non-Gaussian model described in Section 2. These values are evaluated at 10 h⁻¹ Mpc where ≃0.1. The limit ν=1 has the perturbations quadratic in a Gaussian field, as in recently proposed isocurvature models. The limit ν→∞ (the rightmost point) recovers the Gaussian result. The predictions of perturbation theory assuming Gaussian initial conditions, S₃≃2.857 and S₄≃13.89, are shown as the horizontal dashed lines.

Open in new tabDownload slide

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 equations 13–15) would suggest that and . This is to be compared with the Gaussian prediction that and . Defining S′_N such that S′_N would be independent of ν and ξ₂ for a field, we have

(16)

(17)

for unsmoothed fields. For a pure field, and S′₄=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′₃ and S′₄ over a wide range of ν. We show the evolved values of S′₃ and S′₄ vs ν in Fig. 7, with the Gaussian point at ν=∞ omitted. Notice that gravity has modified the initial conditions so that the S′_N are no longer independent of ν. The curves rise to larger ν, indicating that the gravitationally induced skewness and kurtosis constitute a larger fraction of the initial conditions as the model becomes more Gaussian.

Figure 7.

The modified moments, S′₃ and S′₄, defined in the text, as a function of ν as in Fig. 6. The dashed lines show the values S′₃≃2.82 and S′₄=12 of the initial conditions. Note that gravitational instability has induced a dependence on ν.

Open in new tabDownload slide

In addition, the effects of non-linearity change the scaling with ξ₂. For example, for the ν=1 model the S_N are much less dependent on ξ₂ than the S′_N, varying by O(10 per cent) for . This behaviour can be explained as the gravitational contribution becoming more important than the initial conditions as clustering evolves, as anticipated by Fry & Scherrer (1994). 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 make this inference properly 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 equations (3) and (4), as shown in Fig. 6. Moreover, gravitational evolution modifies the scaling of the S_N with ξ₂ and ν. For example, the scale-dependence of the S_N does not follow that of the initial conditions over the range the simulations probe: .

Acknowledgments

I 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₃ and S₄ 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. My work is supported by NSF—9802362.

References