# Halo mass function and scale-dependent bias from N-body simulations with non-Gaussian initial conditions

###### Abstract

We perform a series of high-resolution N-body simulations of cosmological structure formation starting from Gaussian and non-Gaussian initial conditions. We adopt the best-fitting cosmological parameters from the third- and fifth-year data releases of the Wilkinson Microwave Anisotropy Probe and we consider non-Gaussianity of the local type parameterised by eight different values of the non-linearity parameter . Building upon previous work based on the Gaussian case, we show that, when expressed in terms of suitable variables, the mass function of friends-of-friends haloes is approximately universal (i.e. independent of redshift, cosmology, and matter transfer function) to good precision (nearly 10 per cent) also in non-Gaussian scenarios. We provide fitting formulae for the high-mass end () of the universal mass function in terms of , and we also present a non-universal fit in terms of both and to be used for applications requiring higher accuracy. For Gaussian initial conditions, we extend our fit to a wider range of halo masses () and we also provide a consistent fit of the linear halo bias. We show that, for realistic values of , the matter power-spectrum in non-Gaussian cosmologies departs from the Gaussian one by up to two per cent on the scales where the baryonic-oscillation features are imprinted on the two-point statistics. Finally, using both the halo power spectrum and the halo-matter cross spectrum, we confirm the strong -dependence of the halo bias on large scales ( Mpc) which was already detected in previous studies. However, we find that commonly used parameterisations based on the peak-background split do not provide an accurate description of our simulations which present extra dependencies on the wavenumber, the non-linearity parameter and, possibly, the clustering strength. We provide an accurate fit of the simulation data that can be used as a benchmark for future determinations of with galaxy surveys.

###### keywords:

cosmology: theory, dark matter, large-scale structure – methods: N-body simulations – galaxies: haloes, clusters.## 1 Introduction

The detection of temperature anisotropies in the Cosmic Microwave Background (CMB) provided evidence that large-scale structure formation in the universe was seeded by small density fluctuations generated at early times. The statistical properties of these seeds are usually modelled with a Gaussian random field. Historically the Gaussian approximation was introduced for mathematical convenience. In the absence of a solid model for the generation of density fluctuations the Gaussian hypothesis was accepted on the basis of the central limit theorem (e.g. Bardeen et al. 1986 and references therein). The advent of inflationary models provided further support for Gaussianity. Small-amplitude curvature perturbations generated during a standard inflationary phase (single field, slow roll) are very nearly Gaussian distributed (e.g. Bartolo et al. 2004 and references therein).

However, many variants of the inflationary scenario predict appreciable levels of primordial non-Gaussianity. In terms of Bardeen’s gauge-invariant potential, , most of these models (but not all, see e.g. Creminelli et al., 2007) can be reduced to the form:

(1) |

where is an auxiliary Gaussian random field and quantifies the amount of primordial non-Gaussianity. On subhorizon scales, where denotes the usual peculiar gravitational potential related to density fluctuations via Poisson’s equation. The parameter thus has the same sign as the skewness of the density probability distribution function. This local form of non-Gaussianity (note that equation (1) applies in configuration space) can be obtained from a truncated expansion of the effective inflaton potential (Salopek & Bond, 1990; Falk et al., 1993; Gangui et al., 1994). The parameter thus encodes information about the inflaton physics. Standard inflation gives (Salopek & Bond, 1990; Maldacena, 2003). However, even in this case, the non-linear evolution of perturbations on superhorizon scales yields an observable of order unity (which, in reality, should be scale and redshift dependent; Bartolo et al. 2005, see also Pyne & Carroll 1996). Large values of naturally arise in multi-field inflation models (e.g. Linde & Mukhanov 1997; for an extensive review see Bartolo et al. 2004) and even in cyclic or ekpyrotic models of the universe with no inflation (Creminelli & Senatore, 2007; Buchbinder et al., 2008; Lehners & Steinhardt, 2008).

Observational constraints on have been derived studying three-point statistics of temperature fluctuations in the CMB (Komatsu & Spergel, 2001). The recent 5-year data from the Wilkinson Microwave Anisotropy Probe (WMAP) give at the 95 per cent confidence level (Komatsu et al., 2008). Parallel studies on the same dataset give using Minkowski functionals (Komatsu et al. 2008) and from wavelet decomposition (Curto et al., 2008). Some recent reanalyses of earlier 3-year WMAP data claim substantial evidence for positive : from the bispectrum of temperature fluctuations (Yadav & Wandelt, 2008) and from their one-point distribution function (Jeong & Smoot, 2007). On the other hand, a study of Minkowski functionals on the 3-year data gives (Hikage et al., 2008). Higher quality data are needed to improve these constraints. The upcoming Planck satellite should be able to reduce the uncertainty in to (Komatsu & Spergel, 2001).

Alternatively, one might use observational signatures of primordial non-Gaussianity imprinted in the large-scale structure (LSS) of the universe (e.g. Moscardini et al. 1991). Ideally, one would like to use high-redshift probes as the non-linear growth of density fluctuations quickly superimposes a strong non-Gaussian signal onto the primordial one so that the latter might then be difficult to recover. For instance, the large-scale distribution of neutral hydrogen in the era between hydrogen recombination and reionisation encodes information on (Pillepich et al., 2007). This could be probed by detecting the redshifted hyperfine 21-cm transition with very low-frequency radio arrays from space. In principle, an experiment of this kind can limit to (Pillepich et al. 2007, see also Cooray 2006). However, it is not clear yet whether such an experiment will ever be possible due to technical complexity and problematic foreground subtraction. At lower redshifts, can be constrained probing the statistics of rare events, as like as the mass function of galaxy groups and clusters (Matarrese et al., 1986, 2000; Koyama et al., 1999; Robinson & Baker, 2000; Robinson et al., 2000; LoVerde et al, 2008). Early attempts of using cluster counts to constrain have been rather inconclusive due to low-number statistics (see e.g. Willick 2000; Amara & Refregier 2004 and references therein). Even though cluster-mass estimates are still rather uncertain and massive objects are very rare, the observational perspectives look very promising. A number of galaxy surveys encompassing large fractions of the observable universe are being planned (e.g. ground-based surveys as DES, PanSTARRS, and LSST, and the satellite missions EUCLID and ADEPT) and could potentially lead to solid measurements of (e.g. Dalal et al., 2008; Carbone et al., 2008).

Name | Cosmology | ||||||

1.750 | +750 | 1200 | 20 | 50 | WMAP5 | ||

1.500 | +500 | 1200 | 20 | 50 | WMAP5 | ||

1.250 | +250 | 1200 | 20 | 50 | WMAP5 | ||

1.80 | +80 | 1200 | 20 | 50 | WMAP5 | ||

1.27 | +27 | 1200 | 20 | 50 | WMAP5 | ||

1.0 | 0 | 1200 | 20 | 50 | WMAP5 | ||

1.-27 | -27 | 1200 | 20 | 50 | WMAP5 | ||

1.-80 | -80 | 1200 | 20 | 50 | WMAP5 | ||

2.0 | 0 | 1200 | 20 | 50 | WMAP3 | ||

2.750 | +750 | 1200 | 20 | 50 | WMAP3 | ||

3.0 | 0 | 150 | 3 | 70 | WMAP5 | ||

3.250 | +250 | 150 | 3 | 70 | WMAP5 |

Primordial non-Gaussianity is also expected to modify the clustering properties of massive cosmic structures forming out of rare density fluctuations (Grinstein & Wise, 1986; Matarrese et al., 1986; Lucchin et al., 1988; Koyama et al., 1999). Also in this case, however, the non-linear evolution of the mass density generally superimposes a stronger signal than that generated by primordial non-Gaussianity onto the galaxy three-point statistics. The galaxy bispectrum is thus sensitive to only at high redshift (Verde et al., 2000; Scoccimarro et al., 2004; Sefusatti & Komatsu, 2007).

Recently, Dalal et al. (2008) have shown analytically that primordial non-Gaussianity of the local type is expected to generate a scale-dependent large-scale bias in the clustering properties of massive dark-matter haloes. This is a consequence of the fact that large and small-scale density fluctuations are not independent when . Similar calculations have been presented by Matarrese & Verde (2008), Slosar et al. (2008), Afshordi & Tolley (2008), and McDonald (2008). Numerical simulations by Dalal et al. (2008) are in qualitative agreement with the analytical predictions confirming the presence of a scale-dependent bias. Using these analytical models for halo biasing to describe the clustering amplitude of luminous red galaxies and quasars from the Sloan Digital Sky Survey, Slosar et al. (2008) obtained at the 95 per cent confidence level. This shows that LSS studies are competitive with CMB experiments to constrain primordial non-Gaussianity but also calls for more accurate parameterisations of the mass function and clustering statistics of dark-matter haloes arising from non-Gaussian initial conditions.

Most of the analytic derivations of the non-Gaussian halo mass function (Matarrese et al., 2000; LoVerde et al, 2008, e.g.) are based on the extended Press-Schechter model (Press & Schechter, 1974; Bond et al., 1991) which, in the Gaussian case, is known to produce inaccurate estimates of halo abundance (Sheth & Tormen, 1999; Jenkins et al., 2001). Similarly, the scale dependent bias is obtained either using the peak-background split model (Slosar et al., 2008) or assuming that haloes form from the highest linear density peaks (Matarrese & Verde, 2008). Both techniques have limited validity in the Gaussian case (Jing, 1998; Porciani et al., 1999; Sheth & Tormen, 1999). In this paper we test the accuracy of the excursion-set model and the peak-background split in the non-Gaussian case. This extends the previous studies of Kang et al. (2007), Grossi et al. (2007) and Dalal et al. (2008) for the halo mass function and of Dalal et al. (2008) for the halo bias by exploring more realistic values for with simulations of better quality. In practice, we run a series of high-resolution N-body simulations where we follow the process of structure formation starting from Gaussian and non-Gaussian initial conditions. The halo mass function and bias extracted from the simulations are then compared with the existing analytical models and used to build accurate fitting formulae. These will provide a benchmark for future determinations of non-Gaussianity with galaxy surveys.

The paper is organised as follows. In Section 2 we describe our N-body simulations. In Sections 3, 4, and 5 we present our results for the halo mass function, the matter power spectrum and the halo bias, respectively. In Section 6 we discuss the implications of our results for the analysis by Slosar et al. (2008). Our conclusions are summarised in Section 7.

Name | ||||||
---|---|---|---|---|---|---|

WMAP3 | 0.73 | 0.76 | 0.95 | 0.24 | 0.042 | 0.76 |

WMAP5 | 0.701 | 0.817 | 0.96 | 0.279 | 0.0462 | 0.721 |

## 2 N-body simulations

### 2.1 Specifics of the simulations

We use the lean version of the tree-PM code Gadget-2 (Springel, 2005) kindly made available by Volker Springel to follow the formation of cosmic structure in a flat cosmology. We run three different series of simulations (each containing collisionless particles) that differ in the adopted cosmology, box size (and thus force softening length, ), and initial redshift (details are summarised in Table 1). The assumed cosmological parameters are listed in Table 2. For our series and they coincide with the 5-yr WMAP best estimates (Komatsu et al., 2008). The combined 3-yr WMAP+LSS results by Spergel et al. (2007) are instead used for series .

We produce non-Gaussian initial conditions directly applying equation (1) after having generated the Gaussian random field with standard Fourier techniques. We consider eight values for the parameter : . The first five are within the current constraints from CMB data (Komatsu et al., 2008), while the three largest values are useful to compare with previous work. Within each series of simulations, we use the same set of random phases to generate the Gaussian potential . This facilitates the comparison between different runs by minimising sample variance.

The linear matter transfer function, , is computed using the Linger code (Bertschinger, 2001) and is applied after creating the non-Gaussian potential in equation (1). Particle displacements and velocities at are generated using the Zel’dovich approximation (Zel’dovich, 1970). A critical discussion of this choice is presented in the Appendix.

Particle positions and velocities are saved for 30 time steps logarithmically spaced in between and . Dark-matter haloes are identified using the standard friends-of-friends (FOF) algorithm with a linking length equal to 0.2 times the mean interparticle distance. We only considered haloes containing at least 100 particles.

Our first two series of simulations only include large periodic boxes covering a volume of where we can study haloes with masses ranging from up to . These simulations will be used to analyse both the mass function and the bias of dark-matter haloes. On the other hand, the third series includes simulations covering a volume of . They will be used to study the mass function and the bias of low-mass haloes with .

### 2.2 A note on the definition of

The definition of given in equation (1) depends on the cosmic epoch at which it is applied. The reason for this time dependence is that both potentials and decay with time proportionally to with the linear growth factor of density fluctuations and the Robertson-Walker scale factor.

In this paper, we define by applying equation (1) at early times, namely at . Other authors have adopted different conventions. Grossi et al. (2007) use the linearly-extrapolated fields at to define . Therefore, their values of the parameter need to be divided by the factor to match ours. In the WMAP5 cosmology, . On the other hand, Dalal et al. (2008) apply equation (1) at , the redshift at which they generate the initial conditions for the simulations. This agrees with our definition to better than 0.01 per cent.

The sign convention for the non-linearity parameter might possibly generate further ambiguity. In our simulations, positive values correspond to positive skewness of the mass-density probability distribution function. The same convention has been adopted by Grossi et al. (2007), Kang et al. (2007) and Dalal et al. (2008).

Acronym | Reference | Functional form | Parameters |
---|---|---|---|

PS | Press & Schechter (1974) | ||

ST | Sheth & Tormen (1999) | ||

J | Jenkins et al. (2001) | ||

R | Reed et al. (2003) | ||

W | Warren et al. (2006) | A=0.7234, a=1.625, b=0.2538, c=1.1982 | |

T | Tinker et al. (2008) | vary with halo overdensity |

## 3 The Halo Mass Function

One of the long standing efforts in cosmology is to determine the mass function of dark matter haloes – i.e. the number of haloes per unit volume per unit mass at redshift – from the statistical properties of the linear density field. Analytical work has suggested that, when expressed in terms of suitable variables, the functional form of should be universal to changes in redshift and cosmology (Press & Schechter, 1974; Bond et al., 1991; Sheth & Tormen, 1999). N-body simulations have shown that this is approximately true when structure formation is seeded by Gaussian perturbations (Jenkins et al., 2001; Evrard et al., 2002; White, 2002; Warren et al., 2006; Tinker et al., 2008).

Following these studies, we describe the halo abundance in our simulations through the following functional form

(2) |

where is the mean background matter density today, and is the variance of the linear density field

(3) |

with the corresponding power spectrum and some window function with mass resolution (here top-hat in real space). The validity of equation (2) has been widely tested against numerical simulations and useful parameterisations for have been provided starting from Gaussian initial conditions (Sheth & Tormen, 1999; Jenkins et al., 2001; Warren et al., 2006). These fitting functions have an accuracy ranging from 5 to 20 per cent depending on redshift, cosmology, and the exact definition of halo masses. Recently, Tinker et al. (2008) have detected deviations from universality in : redshift-dependent corrections are needed to match the mass function in simulations with an accuracy of 5 per cent. This result is based on haloes identified with the spherical overdensity algorithm. It is well known that the mass function of FOF haloes shows a more universal scaling even though other halo finders might be more directly linked to actual observables (Jenkins et al., 2001; Tinker et al., 2008). Deviations from universality for FOF haloes will be further discussed in Section 3.3. One should anyway keep in mind that baryonic physics can cause 30 per cent deviations in with respect to the pure dark-matter case (Stanek et al., 2008).

### 3.1 Halo mass function from Gaussian initial conditions

The halo mass functions extracted from our Gaussian simulations – Run 1.0 (triangles), Run 2.0 (squares), and Run 3.0 (circles) – are presented in Figure 1. The combination of different box sizes allows us to cover the very wide range which roughly corresponds to the mass interval at . Figure 1 has been obtained by combining data from snapshots at redshifts . Note that, at a fixed redshift, larger values of correspond to higher masses. On the other hand, with increasing the redshift, larger values of are associated with a given halo mass. Even though datapoints correspond to different redshifts and cosmologies, they all form a well defined sequence. This indicates that the function is universal to good approximation. For a given , outputs at a fixed redshift scatter around the universal sequence by 10-15 per cent. A number of fitting formulae have been proposed in the literature to parameterise this sequence. In Figure 1, we compare some of them (summarised in Table 3) with our datapoints. Fractional deviations between models and data are shown in the bottom panel. Barring the classical Press-Schechter result, all the fitting formulae describe our data to better than 20 per cent. The best agreement is found all over the mass range with Warren et al. (2006) followed by Jenkins et al. (2001) which both show deviations from our data at the 10 per cent level. The Sheth & Tormen (1999) model also provides an accurate description of the data for small halo masses but tends to overestimate the abundance of the most massive objects. On the other hand, the fit by Reed et al. (2003) tends to underestimate the high-mass tail of the mass function. Overall our findings are in good agreement with Heitmann et al. (2006) and Lukić et al. (2007).

Following Warren et al. (2006) and Tinker et al. (2008), we fit the outcome of the simulations with the function

(4) |

The best-fitting parameters have been determined through minimisation using the Markov Chain Monte Carlo method, and read:

(5) | |||||

In terms of the parameterisation given in Warren et al. (2006) and reported in Table 3, this corresponds to . The fit in equation (3.1) describes our dataset up to deviations of a few per cent over the entire mass and redshift ranges for Run 1.0 and Run 2.0, while it shows larger deviations (up to nearly 10 per cent) towards the high-mass end of Run 3.0, (see Figure 1). It is important to remember, however, that Run 3.0 covers a much smaller volume than the others and thus is more severely affected by sample variance.

### 3.2 The universal halo mass function from non-Gaussian initial conditions

Is the function universal also in the non-Gaussian case? This question is addressed in Figure 2 where we show the output of our main series of simulations at four redshifts () to test the scaling of the mass function in terms of . Only bins containing at least 20 haloes are considered. Within a certain tolerance, the halo mass functions at different masses and redshifts all lie on the same curve for a given . The scatter of the points at a fixed redshift around this curve roughly amounts to 10 per cent, and it becomes smaller towards our largest values of .

We thus generalise equation (2) to non-Gaussian initial conditions by assuming that

(6) |

and we provide a fitting formula for . Given the similarity to the Gaussian case, we still adopt the functional form given in equation (4) but let the parameters vary with . The best-fitting values have been determined in two steps. We first used a Markov Chain Monte Carlo method to determine at fixed through minimisation. The results suggest that the dependence for each parameter of the mass function can be accurately described by polynomials of different orders. Eventually, we used the data to derive the coefficients of these polynomials.

Parameter | ||
---|---|---|

A | 1.694 | -0.00199 |

B | 0.566 | -0.00029 |

C | 1.151 | -0.00071 |

D | 0.287 | -0.00030 |

Parameter | |||||
---|---|---|---|---|---|

A | 1.708 | -2.07 | 3.1 | 0 | 0 |

B | 0.560 | -0.46 | +12.46 | -2.36 | +2.65 |

C | 1.150 | -0.76 | +2.7 | 0 | 0 |

D | 0.293 | -0.16 | -14.07 | +3.88 | -3.94 |

.

The degree of complexity required to fit the simulation data grows considerably with increasing . For (a range that fully encloses the values currently allowed by CMB studies), the mass-function parameters in equation (4) are well approximated by the linear relation

(7) |

Table 4 lists the corresponding best-fitting parameters.
The quality of this fitting formula is assessed in
the left panel of Figure 3, where the mass function for the
simulations with is compared with
the corresponding fit. Residuals (shown in the bottom panel)
are smaller than 5 per cent all over the range
corresponding to the mass interval
at .

On the other hand, equation (7)
is not suitable to account for values of substantially larger than 250.
To obtain an accurate fit of the universal halo mass function over the range
we had to consider
polynomials up to order in :

(8) |

and

The best-fitting values of the parameters above are listed in Table
5 while the corresponding functions are compared with
the simulation data in the right panel of
Figure 3. Also in this case residuals are smaller than
5 per cent for .

The universality of the fitting formula in equation (6) has been further tested against our non-Gaussian simulation of the WMAP3 cosmology, Run2.750, which has not been used to determine the best-fitting parameters. This blind check shows that, in the range (roughly corresponding to at ), the provided fit reproduces the mass function with an accuracy of 5 per cent.

We warn the readers against extending our fitting formulae beyond their range of validity, in particular at low halo masses. The simulations of our main series resolve haloes with 100 particles. For , our analytical formulae for the mass function have been derived using only haloes that are more massive than this limit. Moreover, since the high-mass tail of the mass function is enhanced (suppressed) for positive (negative) values of with respect to the Gaussian case, mass conservation requires that the opposite effect is seen at lower masses. We have directly tested the goodness of our fit towards the smaller masses using Run3.250 (which has a boxsize 8 times smaller than for the simulations in the main series but the same number of particles) and indeed found that the fitting formulae in equations (8) and (3.2) systematically overestimate the abundance of small mass haloes by 10-30 per cent. We will address the low-mass tail of the mass function for in future work.

On the other hand, for Gaussian initial conditions, we combined simulations
with different box sizes to derive the fitting function in equations
(4) and
(3.1). This allowed us to extend the validity
of our fit to the much wider mass range
.

Our fitting formulae give three different approximations for
the universal mass function in the Gaussian case. In general, the fit given
in equations (4) and
(3.1) has to be preferred
as it has been obtained from
a richer dataset spanning a much wider range of halo masses.
However for masses
above at , the fit in equations
(4) and (7) and Table 4 provides
the most accurate representation of our
data. In any case, the different
fitting functions never deviate by more than 3-4 per cent.
Also note that our two fitting functions for the non-Gaussian simulations
agree by better than 1 per cent for and by
a few per cent for and .

### 3.3 The limit of universality: redshift dependence

Regardless of the value of ,
we have found that the halo mass function is universal, when written in
terms of , with an accuracy of roughly 10 per cent.
If one is interested in giving
analytical approximations for the halo mass function which are more accurate
than the universal fit,
it is necessary to introduce redshift-dependent corrections
(see also Tinker et al. 2008 for the Gaussian case).
In the left panel of Figure 4, we show
how well the universal fit
(whose parameters are listed in Table 5) describes
the simulation outputs at .
At and
for masses
the fitting formula deviates for the data by more than
10 per cent. The smaller the redshift, the worse is the agreement between the
data points and the universal fit.
The bigger the , the less critical is the comparison.

In this Section we provide a non-universal fit which is very accurate at
low redshift. In particular, we write:

(10) |

where is the rms deviation of the linear density field at . We approximate with the functional form given in equation (4) but now let the parameters vary with both and . Markov Chain Monte Carlo fitting suggests that each parameter of the mass function can be accurately described as follows:

(11) |

Parameter | ||||
---|---|---|---|---|

A | 1.82 | 2.85 | 4.53 | -5.92 |

B | 0.578 | 5.30 | 7.53 | -7.73 |

C | 1.15 | 9.52 | 9.08 | -4.42 |

D | 0.294 | 4.92 | 4.67 | -3.80 |

The best fitting parameters for and are listed in Table 6, while the quality of the fitting formula is assessed in the right panel of Figure 4. Residuals are smaller than 5 per cent all over the mass range, indicating that for and the fit of equations 10, 4, and 11 has to be preferred to the universal fit given in the previous Section. On the other hand, for higher values of and for higher redshifts, the universal fit gives a better and more economic (in terms of parameters) description of the data.

### 3.4 Comparison with theoretical models

The halo mass function arising from mildly non-Gaussian
initial conditions can be modelled by generalizing the Press-Schechter
formalism.
Using the saddle-point approximation to evaluate the probability for the
linear density field to be above a given threshold value,
Matarrese et al. (2000) have derived a model for .
More recently, LoVerde et al (2008) presented another expression for the mass
function by using
the Edgeworth asymptotic expansion for the probability density function
of the linear density field.
In both cases, only leading-order corrections in have been
accounted for.
In absolute terms, these models are not expected to be accurate as
they should suffer from the same shortcomings as the Press-Schechter model
in the Gaussian case.
However, they can be used to compute the fractional non-Gaussian correction
(Verde et al., 2000; LoVerde et al, 2008).
In Figure 5, we use this quantity
to test the models against our simulations.
Datapoints with errorbars show
the N-body output at and 1, while
the dotted lines indicate the models of Matarrese et al. (2000) and LoVerde et al (2008)
as indicated in equations (B.6)^{1}^{1}1This fixes a typo in equation
(68) of Matarrese et al. (2000) and (4.19) of LoVerde et al (2008),
respectively.
It is evident that the models overestimate the non-Gaussian correction.
Following the indications in Carbone et al. (2008),
we also show a modified version of the models which is obtained by lowering
the critical threshold for halo collapse as
(solid lines in Figure 5).
Such a correction vastly improves the agreement with the simulations.

Dalal et al. (2008) proposed to fit the halo mass function in terms
of the convolution between and a Gaussian kernel in
with a -dependent mean and variance. Figure 5 shows that their
fit tends to overestimate the non-Gaussian corrections especially for
large, positive values of and masses .
On the other hand, for it has a similar accuracy as the
formulae derived from the Press-Schechter formalism corrected with the
reduced threshold.

The good agreement between the fractional non-Gaussian corrections derived from the modified PS models and from the simulations is not enough to derive from future observations of galaxy clusters. In fact the ratio is not an observable: the only quantity that we can hope to compare with observations is the mass function. In order to make predictions for , the models for the fractional non-Gaussian correction need to be multiplied by a Gaussian mass function. This step might introduce relatively large systematic errors (see Figure 1) which could degrade any measurement of based on the cluster mass function. We address this issue in Figure 6 where we plot the fractional deviation of some model predictions for the function with respect to the simulation output (results are very similar for different values of ). We consider the model by LoVerde et al (2008) corrected with the factor and multiplied by three different Gaussian models: Sheth & Tormen (1999), Warren et al. (2006), and our fit with . Note that some of the final outcomes systematically differ by 10-20 per cent over the entire mass range covered by the simulations. This clearly shows that a careful measurement of the Gaussian mass function is necessary to avoid a biased estimation of the non-linearity parameter. Note that, for , the models by Matarrese et al. (2000) and LoVerde et al (2008) (both with the reduced collapse threshold) combined with our Gaussian fit are in rather good agreement with the numerical mass functions (similar results are obtained using the Gaussian fit by Warren et al. (2006) for masses below a few ). Perhaps not surprisingly, no model describes the simulation data for all the values of as well as our fitting formulae for the non-Gaussian mass function given in Sections 3.2 and 3.3.

### 3.5 Summary of accuracy and range of validity of the mass function fits

In order to facilitate the use of our fitting formulae for the halo mass function we summarize here their accuracy and range of validity.

## 4 Matter Power Spectrum

In this section we study how non-Gaussian initial conditions
influence the power spectrum of the mass density field.
At tree level, the power spectrum does not depend on in
Eulerian perturbation theory. However, one-loop corrections
make the power spectrum -dependent.
Qualitatively, theoretical expectations are that
positive (negative) values of
tend to enhance (suppress) the amplitude of the power spectrum on non-linear
scales.
In Figure
7 we plot the ratio of power spectra
extracted from the simulations
of our main series at redshifts and 1.
The matter power spectrum of non-Gaussian models appears to deviate
already by a few per cent at Mpc. As expected, deviations
become more severe with increasing the wavenumber .
Our results are in agreement with the perturbative
calculations
by Taruya et al. (2008).
We note, however, that Grossi et al. (2008) found
smaller deviations between the non-Gaussian and Gaussian power spectra
at larger values of and .

Our results have two important practical implications. First, the widespread habit of using the Gaussian matter power spectrum to determine non-Gaussian bias parameters leads to scale-dependent systematic errors that might become severe when high-precision is required. Second, primordial non-Gaussianity modifies the power-spectrum on the scales where baryonic oscillations (BAOs) are present. Reversing the argument, two-point statistics could be also used to constrain the value of . Note however, that all probes based on galaxy clustering will suffer from uncertainties in the bias parameter (and its scale dependence) that might hinder a measure of based on the study of BAOs. On the other hand, weak lensing studies will directly measure the matter power spectrum. The target of many future wide-field missions is to provide estimates at the per cent level. For parameter estimation, a comparable accuracy will be required on model spectra within a wide range of wavenumbers centred around Mpc (Huterer & Takada, 2005). Therefore, even values of within the current CMB constraints could imprint detectable effects in the matter power spectrum at the scales of interest. The key question is whether one can discern the effect of and, consequently, how much primordial non-Gaussianity will affect the estimate of the other cosmological parameters. We will get back to this in future work.

## 5 Halo Clustering

The clustering of dark-matter haloes is biased relative to that of the underlying mass distribution by an amount which depends on halo mass, redshift, and the scale at which the clustering is considered (see e.g. Mo & White, 1996; Catelan et al., 1998; Smith et al., 2007). For Gaussian initial conditions, this has been widely tested against numerical simulations (e.g. Sheth et al., 2001; Seljak & Warren, 2004; Tinker et al., 2005).

In general, the halo bias can be quantified using either the power spectrum of the halo density field, , or the cross-spectrum between the halo and the underlying matter density field, . In the two cases the bias reads

(12) |

or

(13) |

where is the matter power spectrum. If the bias due to halo formation is local and deterministic then apart from measurement errors. However, in the presence of a stochastic component that does not correlate with the density field . In practice, however, the measurement of all power spectra is affected to some level by shot noise due to the discrete nature of dark-matter haloes and N-body particles. If the distribution of the tracers can be approximated as the Poisson sampling of an ideal density field, then the measured power spectrum corresponds to that of the underlying field plus the mean volume per particle (Peebles, 1980). Discreteness effects are thus expected to be negligible for and due to the large number density of particles in the simulations. On the other hand, massive haloes are rare and, being extended objects, cannot be modelled as the Poisson sampling of a continuous distribution (Mo & White 1996, Magliocchetti & Porciani 2003, Porciani in preparation). It is not clear then how to correct for the discreteness effect in their power spectrum (Smith et al., 2007). For these reasons we use in our analysis and we adopt (without performing any discreteness correction) only to verify the results (see Figure 8).

### 5.1 Halo bias from Gaussian initial conditions

It is well known that the halo bias factor from Gaussian initial conditions is approximately scale-independent for small values of the wavenumber . We will refer to this asymptotic value on large scales as the “linear bias” and denote it by . Similarly to the halo mass function, when expressed in terms of , the linear bias assumes a universal form which, within a given accuracy, is independent of redshift and just weakly dependent on cosmology (e.g. Sheth & Tormen, 1999; Seljak & Warren, 2004).

We measure the linear bias for the haloes in our simulations as follows.
We first determine the functions and
by directly applying equations (12) and (13).
Within the statistical uncertainties,
both functions approach asymptotically to a constant on large scales (
Mpc). We use the average of the bias function
measured in the range Mpc (4 -bins)
as our estimate of the linear bias. The standard error of the mean is
used to quantify the corresponding statistical
uncertainty.^{2}^{2}2Consistently, in
what follows, we use the rms value as the error on .
For , we assume that the relative error
is the same as in the Gaussian case.

In Figure 9 we show the linear bias obtained from Run 1.0 (triangles), Run 2.0 (squares) and Run 3.0 (circles) as a function of . Simulation data from snapshots between and are compared with the commonly used parameterisations listed in Table 7. Our results are in good agreement with the fit by Sheth et al. (2001) for large masses and with that by Tinker et al. (2005) for smaller masses. Note that by combining simulation boxes we are able to explore a larger interval of than previous studies.

Acronym | Reference | Functional form | Parameters and Variables |
---|---|---|---|

MW | Mo & White (1996) | ||

ST | Sheth et al. (2001) | ||

SW | Seljak & Warren (2004) | with | |

T | Tinker et al. (2005) |

Given that no existing model for the linear bias accurately reproduces our results over the entire mass range spanned by the simulations, we decided to derive a new fitting formula. In particular, we parameterised the outcome of our simulations as

(14) |

and used minimisation to find

(15) | |||||

This fit (which reproduces the numerical data with great accuracy in the range ) should be considered as the linear bias naturally associated with the mass function given in equations (4) and (3.1).

### 5.2 Halo bias from non-Gaussian initial conditions

Recent analytical models, have suggested that the halo bias arising from non-Gaussian initial conditions of the local type does not tend to a constant on large scales. Rather, the deviation from the Gaussian case should follow

(16) | |||||

where , Mpc is the Hubble radius,
is the matter transfer function,
and is the linear growth factor of matter perturbations
normalised to unity at
(Dalal et al., 2008; Matarrese & Verde, 2008; Slosar et al., 2008; Afshordi & Tolley, 2008; McDonald, 2008).^{3}^{3}3The factor
is needed since
Dalal et al. (2008) and Slosar et al. (2008)
normalise the growth factor to be during matter domination.
The numerical simulations by
Dalal et al. (2008) have indeed shown that the halo bias is scale dependent
even for small values of
in non-Gaussian cosmologies (with )
and found qualitative agreement with equation (16).
In Figure 10, we show how the bias depends on scale
in our simulations which also consider smaller values of .
Our results confirm the presence of a strongly scale-dependent bias.
Larger values of correspond to a more marked
scale dependence. Note, however, that for Mpc
the non-Gaussian deviation changes sign. On these scales,
the halo-matter and halo-halo spectra emerging from non-Gaussian perturbations
has actually
less power than in the Gaussian case. The opposite happens with the
matter power spectrum (even to a larger degree) and
the net effect
is a negative .
This result implies that equation (16)
can only hold asymptotically on very large scales.
This is not suprising if interpreted within the peak-background-split
formalism where the large-scale bias is linked to the first derivative
of the mass function with respect to . In the non-Gaussian case
the bias is composed of two parts, a scale-independent term and the correction
given in equation (16).
Since the halo mass function changes shape when
is varied, also the constant bias should depend on for a fixed
halo mass.
Increasing corresponds to a larger abundance of massive haloes
and to a slightly smaller constant bias with respect to the Gaussian case.
Likely, this is what makes the in the simulations negative for
positive .
To proceed with a detailed analysis of our simulations,
we find it convenient to rewrite equation (16) as

(17) |

where and . In Figure 11 we test the scaling of with redshift, linear bias and wavenumber for (where we have the best signal-to-noise ratio at high halo mass). Similar results are obtained with different values of . The quantity shown is which should correspond to if the analytical model provides a good description of the data. This quantity is indicated by a dashed line. The following two trends clearly emerge from the data. For small values of , the model overestimates the data by 20-70 per cent increasing with and independently of . On smaller scales, discrepancies become more and more severe. At Mpc, the model is systematically a factor of 5 higher than the data. The -dependence of is therefore different than in equation (16).

The data also drop a hint that, for Mpc,
the scaling with might only persist up to a maximum value of ,
. For it appears that the values
of are always smaller than expected from the extrapolation of the
trend determined at smaller .
The value of seems to depend both on redshift and
wavenumber and roughly corresponds to constant halo mass for a given .
However, uncertainties in at these high masses become very large and
it is difficult to judge how robust the existence of really
is.
We note anyway that when we tried to fit data at different redshifts (for
a given and Mpc)
by adding a variable normalisation constant
in front of equation (16),
we systematically obtained significantly different fits
(at a confidence level of
2.5 ) at different redshifts. This trend disappears when only the lowest
values of are considered at each redshift for the fit.

Data from simulations with all the considered values of are shown with different symbols and colors in Figure 12. Each panel refers to a particular wavenumber bin (indicated by the label in units of Mpc). The model in equation (16) is again indicated by a dashed line. Note that, in most cases, it substantially deviates from the simulation data. In particular, measured from the simulations shows a much stronger -dependence than the analytical formula, as already seen in Figure 11. In general, the overall amplitude of drops by an extra factor of with respect to when moving from Mpc to Mpc independently of and . Also, does not seem to scale linearly with while its linear dependence on appears to be solid, at least for . We thus introduce a correcting factor defined by

(18) |

and we measure it by fitting the simulation data for and at constant values of and . We use an effective variance weighted least squares method to simultaneously account for errorbars on both bias parameters. The best-fitting values are reported in Table 8 and can be used to compute the function by interpolation. The final expression for , corrected with the factor, is shown in Figure 12 with solid lines.

0.0117 | |||||||
---|---|---|---|---|---|---|---|

0.0191 | |||||||

0.0303 | |||||||

0.0494 | |||||||

0.0804 |

Data in Table 8 have an amazing regularity. Apart from a normalisation constant, each column (row) shows the same linear trend with (). This suggests that, within the explored parameter range ( Mpc and ),

(19) |

We thus use this equation to fit the original data for the halo bias from Gaussian and non-Gaussian initial conditions and find

(20) | |||||

at the 68.3 per cent confidence level. Note that we computed the power spectra in finite-sized bins of the wavenumber, so that there is some degree of ambiguity in associating the results with a given value of . Unfortunately the choice plays a role in determining as is a steep function of on the scales of interest. In Table 8 and in equation (20), we have used the arithmetic mean of the wavenumbers contributing to a given bin. If one instead uses the logarithmic center of the bin, is slightly reduced with a best-fitting value of . The parameters and are unaltered. Therefore, a systematic contribution should be added to the error budget of .

equations (16) and (18) assume that the Gaussian bias is constant with but this is only approximately true in the simulations. The fit in equation (20), the Table 8 and the Figures (11) and (12) have been obtained by identifying with the actual bias measured in the Gaussian simulation at each wavenumber. If, instead, the estimate for introduced in Section 5.1 is used, one gets , and , slightly improving the goodness of the fit.

The quadratic dependence of on is rather surprising as it cannot be straightforwardly derived from the simple models listed above. It might possibly arise from higher-order terms which have been neglected in the expansion that leads to equation (16). Anyway, it is clearly present in the simulations as it can be seen by looking at the variation of along a given row in Table 8. Within the range of of physical interest, the effect is rather small: the coefficient only corresponds to a few percent correction. Note that a quadratic term breaks the symmetry in the amplitude of between non-linearity parameters with opposite sign and identical absolute value. It is hard to directly test this against our simulations as we just have two runs with and both of them correspond to rather small where the uncertainties in are large. An alternative explanation for a non-vanishing could be that it artificially derives from imposing a linear relation in to data that do not scale linearly for . Indeed, just using datapoints with small values of we derive bigger values of for large (more or less in line with ). Therefore, what is robust is that at least one of the scalings with or with is incorrect in equation (16). We found that a scaling proportional to (with and two adjustable parameters) does slightly better (in terms of reduced ) than , at least for Mpc. However, since the scaling with has a sound theoretical basis (Mo & White, 1996; Catelan et al., 1998) we preferred to quote our results as in equation (19). From the statistical point of view, the parameters (20) provide an acceptable description of the simulation data to high confidence for all values of . However, they are particularly accurate for , while (with the same and ) has to preferred for smaller values of .

The linear correction in should be thought of as the first-order term
of a series expansion in the wavenumber.
We attempted to determine the corresponding quadratic term by
considering larger values of in the fit (one bin more, up to
Mpc).
However, values of become small compared with the numerical errors
and we found that the quadratic parameter is badly constrained
by the data () while the other
parameters remain nearly unchanged (and get larger uncertainties).
Also note that the Gaussian bias starts to depart
from at Mpc and it is not clear whether
equation (16) should still be expected to hold in this regime.

Dalal et al. (2008) derived an expression for which coincides
with equation (16) but does not include the
linear transfer function.
Theoretically, this is hard to understand, as
non-Gaussianity is generated well before matter-radiation equality
and one should account for the linear evolution of density perturbations.
Anyway, due to the different -dependence, their expression for
provides a better fit to the simulation data than equation (16)
when both models are allowed to vary in amplitude
with a tunable free parameter.^{4}^{4}4The best-fitting value for
this coefficient reads for the model of Dalal et al. (2008) and
for the
model of
Slosar et al. (2008) None of them, however, provides such an accurate fit to the data
as our
equations (19) and (20), which improve the by at
least a factor of .

## 6 Discussion

Slosar et al. (2008) have used equation (16) to constrain by considering measures of the clustering amplitude of luminous red galaxies (LRGs) and quasars from the Sloan Digital Sky Survey. Combining all datasets, they found at 95 per cent confidence. How would this result change based on our simulations? Disentangling the different contributions, the strongest constraints to in Slosar et al. (2008) come from the angular power spectrum of quasars with photometric redshifts in the range and a mean bias of . Weaker limits are also contributed from the power spectrum of spectroscopic LRGs and the angular spectrum of photometric LRGs. (with a bias of at ). Figure (12) and Table 8 suggest that at the scales of interest ( Mpc) the model given in equation (16) tends to overestimate the scale-dependent bias seen in the simulations. Therefore, to match an observed , a larger value of is required than predicted by the analytic model. When applied to the data by Slosar et al. (2008), our correction should thus give somewhat looser limits on . Because of the strong -dependence of the function it is impossible to give more precise estimates without fitting the power-spectrum data. Note, however, that a steeper -dependence potentially makes determinations of even more competitive with respect to studies of CMB anisotropies.

## 7 Summary

We use a series of high-resolution N-body simulations to study the mass function and the clustering properties of dark-matter haloes arising from Gaussian and non-Gaussian initial conditions. In particular, we follow the formation of structure in a universe characterised by the best-fitting parameters from the third- and fifth-year WMAP data releases. We consider non-Gaussianity of the local type and we use eight different values of () enclosing the parameter space currently allowed by studies of the cosmic microwave background. Our main results can be summarised as follows.

(i) The mass function of dark-matter haloes varies with . Larger values of the non-linearity parameter correspond to higher abundances of the most massive haloes. Analytical models based on the Press-Schechter method (Matarrese et al., 2000; LoVerde et al, 2008) are compatible with our simulated results for the ratio of the Gaussian and non-Gaussian mass functions only if the critical threshold for halo collapse is lowered to . An accurate fit of the Gaussian is necessary to derive the non-Gaussian mass function from the aforementioned ratio.

(ii) We find that, in perfect analogy with the Gaussian case (Jenkins et al., 2001), the halo mass function assumes an approximately universal form. This means that, when expressed in terms of suitable variables, its dependence on redshift and cosmology is erased to good precision (nearly 10 per cent). We parameterise the -dependence of the universal mass function and provide an accurate fit for its high-mass end. For and for masses , the best-fitting parameters for the non-Gaussian halo mass function in equation (4) are given in equation (7) and Table 4. This fit reproduces the mass function of friend-of-friends haloes with an accuracy of 5 per cent on top of a systematic contribution (up to 10 per cent) due to the non perfect universality. For applications requiring higher precision, an additional formula is provided: for and the fit in equations (10), (4), and (11) has to be preferred to the universal fit. On the other hand, for higher values of and for higher redshifts, the universal fit gives a better and more economic (in terms of parameters) description of the data. In the Gaussian case, we extend the fit to a larger interval of halo masses () by combining simulations with different box sizes: – see equations (4) and (3.1). Our fitting function provides a precious tool to forecast constraints on from future surveys and to analyze current datasets.

(iii) The matter power-spectrum in non-Gaussian cosmologies departs from the Gaussian one already on very large scales. For values of within the current CMB constraints these scale-dependent deviations can be as high as two per cent at and increase with wavenumber. The discrepancy is systematic: models with positive have more large-scale power than the Gaussian case and models with negative have less. This warns against the widespread habit of using the Gaussian matter power spectrum to determine non-Gaussian bias parameters when high-precision is required. It also suggests that primordial non-Gaussianity modifies the shape and amplitude of the baryonic-oscillation feature in the two-point statistics and the convergence power spectrum in weak-lensing studies.

(iv) We present an accurate fitting formula for the linear bias of dark matter haloes arising from Gaussian initial conditions extending previous work to larger mass intervals. This, together with the mass function fit mentioned above, can be used to constrain parameters of halo-occupation models from clustering data.

(v) Finally, using the halo-matter cross spectrum, we confirm the strong -dependence of the halo bias on large scales ( Mpc) which was already detected by Dalal et al. (2008). However, we show that commonly used parameterisations based on the peak-background split overestimate the effect for Mpc. The discrepancy increases with the wavenumber and at Mpc in the simulations changes sign with respect to the models. On top of this, the analytic model for the scale-dependent part of the bias requires corrections which depend on the non-linearity parameter, the wavenumber and, possibly, also on redshift and clustering strength. equations (18) and (19) with the best-fitting parameters listed in (20) provide a fitting formula which accurately reproduces the outcome of the simulations for Mpc and . This fit should be employed to constrain from future clustering data at low and high redshift.

## Acknowledgments

All simulations were performed at the Swiss National Supercomputing Center (CSCS) in Manno, Switzerland. AP and OH acknowledge support from the Swiss National Science Foundation. We thank Volker Springel for kindly making the lean version of the Gadget-2 code available to us. We acknowledge discussions with Robert E. Smith and Tommaso Giannantonio. While our paper was being submitted, a related work by Desjacques, Seljak & Iliev (2008) appeared as a preprint. Their results based on a different halo finder are in qualitative agreement with ours with the exception of the non-Gaussian bias correction for small halo masses () which they find to exceed the model by Slosar et al. (2008).

## References

- Afshordi & Tolley (2008) Afshordi N.,Tolley A.J., arXiv:0806.104
- Amara & Refregier (2004) Amara, A., & Refregier, A. 2004, MNRAS, 351, 375
- Bardeen et al. (1986) Bardeen J.M., Bond J.R., Kaiser N., Szalay A.S., 1986, ApJ, 304, 15
- Bartolo et al. (2004) Bartolo N., Komatsu E., Matarrese S., Riotto A., 2004, PhR, 402, 103
- Bartolo et al. (2005) Bartolo N., Matarrese S., Riotto A., 2005, JCAP, 10, 010
- Bertschinger (2001) Bertschinger, E. 2001, ApJS, 137, 1
- Bond et al. (1991) Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, Apj, 379, 440
- Buchbinder et al. (2008) Buchbinder, E. I., Khoury, J., & Ovrut, B. A. 2008, Physical Review Letters, 100, 171302
- Catelan et al. (1998) Catelan, P., Lucchin, F., Matarrese, S., & Porciani, C. 1998, MNRAS, 297, 692 Matarrese,
- Carbone et al. (2008) Carbone C., Verde L., Matarrese S., 2008, ApJ, 684, 1
- Cohn & White (2008) Cohn J.D., White M., 2008, MNRAS, 385, 2025
- Cooray (2006) Cooray, A. 2006, Physical Review Letters, 97, 261301
- Creminelli & Senatore (2007) Creminelli P., Senatore L., 2007, JCAP, 11, 010
- Creminelli et al. (2007) Creminelli P., Senatore L., Zaldarriaga M., 2007, JCAP, 03, 19
- Crocce et al. (2006) Crocce M., Pueblas S., Scoccimarro R., 2006, MNRAS, 373, 369
- Curto et al. (2008) Curto, A., Martinez-Gonzalez, E., Mukherjee, P., Barreiro, R. B., Hansen, F. K., Liguori, M., & Matarrese, S. 2008, arXiv:0807.0231
- Dalal et al. (2008) Dalal N., Dore’ O., Huterer D.,Shirokov A., 2008, Phys. Rev. D77, 123514
- Desjacques et al. (2008) Desjacques, V., Seljak, U., & Iliev, I. T. 2008, arXiv:0811.2748
- Evrard et al. (2002) Evrard, A. E., et al. 2002, ApJ, 573, 7
- Falk et al. (1993) Falk T., Rangarajan R., Srednicki M., 1993, ApJ, 403, 1
- Gangui et al. (1994) Gangui A., Lucchin F., Matarrese S., Mollerach S., 1994, ApJ, 430, 447
- Grinstein & Wise (1986) Grinstein, B., & Wise, M. B. 1986, ApJ, 310, 19
- Grossi et al. (2007) Grossi M., Dolag K., Branchini E., Matarrese S., Moscardini L., 2007, MNRAS, 382, 1261
- Grossi et al. (2008) Grossi M., Branchini E., Dolag K., Matarrese S., Moscardini L., 2008, arXiv:0805.0276
- Grossi et al. (2009) Grossi M., Verde L., Carbone C., Dolag K., Branchini E., Iannuzzi F., Matarrese S., Moscardini L., 2009, arXiv:0902.2013
- Heitmann et al. (2006) Heitmann K., Lukić Z., Habib S., Ricker P.M., 2006, ApJ, 642, 85
- Hikage et al. (2008) Hikage, C., Matsubara, T., Coles, P., Liguori, M., Hansen, F. K., & Matarrese, S. 2008, MNRAS, 389, 1439
- Huterer & Takada (2005) Huterer, D., & Takada, M. 2005, Astroparticle Physics, 23, 369
- Jenkins et al. (2001) Jenkins A., 2001, MNRAS, 321, 372
- Jeong & Smoot (2007) Jeong, E., & Smoot, G. F. 2007, arXiv:0710.2371
- Jing (1998) Jing, Y. P. 1998, ApJ, 503, L9
- Kang et al. (2007) Kang X., Norberg P., Silk J., 2007, MNRAS, 376, 343
- Knebe et al. (2009) Knebe A., Wagner C., Knollmann S., Diekershoff T., Krause F., 2009, arXiv:0904.0083
- Komatsu & Spergel (2001) Komatsu, E., & Spergel, D. N. 2001, PhRvD, 63, 063002
- Komatsu et al. (2008) Komatsu, E., et al. 2008, arXiv:0803.0547
- Koyama et al. (1999) Koyama, K., Soda, J., & Taruya, A. 1999, MNRAS, 310, 1111
- Lehners & Steinhardt (2008) Lehners, J.-L., & Steinhardt, P. J. 2008, PhRvD, 77, 063533
- Linde & Mukhanov (1997) Linde A., Mukhanov V., 1997, PhRvD, 56, 535
- LoVerde et al (2008) LoVerde M., Miller A., Shandera S., Verde L., 2008, JCAP, 04, 014
- Lucchin et al. (1988) Lucchin, F., Matarrese, S., & Vittorio, N. 1988, ApJl, 330, L21
- Lukić et al. (2007) Lukić Z., Heitmann K., Habib S., Bashinsky S., Ricker P.M., 2007, ApJ, 671, 1160
- Magliocchetti & Porciani (2003) Magliocchetti, M., & Porciani, C. 2003, MNRAS, 346, 186
- Maldacena (2003) Maldacena J., 2003, JHEP, 05, 013
- Matarrese et al. (1986) Matarrese, S., Lucchin, F., & Bonometto, S. A. 1986, ApJ., 310, L21
- Matarrese et al. (2000) Matarrese S., Verde L., Jimenez R., 2000, ApJ, 541, 10
- Matarrese & Verde (2008) Matarrese S., Verde L., 2008, ApJL, 677, 77
- McDonald (2008) McDonald P., arXiv:0806.1061
- Mo & White (1996) Mo, H. J., & White, S. D. M. 1996, MNRAS, 282, 347
- Moscardini et al. (1991) Moscardini, L., Matarrese, S., Lucchin, F., & Messina, A. 1991, MNRAS, 248, 424
- Peebles (1980) Peebles, P. J. E. 1980, Research supported by the National Science Foundation. Princeton, N.J., Princeton University Press, 1980. 435 p.,
- Pillepich et al. (2007) Pillepich, A., Porciani, C., & Matarrese, S. 2007, ApJ, 662, 1
- Porciani et al. (1999) Porciani, C., Catelan, P., & Lacey, C. 1999, ApJ, 513, L99
- Press & Schechter (1974) Press W.H., Schechter P. 1974, ApJ, 187, 425
- Pyne & Carroll (1996) Pyne T., Carroll S.M., 1996, PhRvD, 53, 2920
- Reed et al. (2003) Reed D., Gardner J., Quinn T., Stadel J., Fardal M., Lake G., Governato F., 2003, MNRAS, 346, 365
- Robinson & Baker (2000) Robinson, J., & Baker, J. E. 2000, MNRAS, 311, 781
- Robinson et al. (2000) Robinson, J., Gawiser, E., & Silk, J. 2000, ApJ, 532, 1
- Salopek & Bond (1990) Salopek D.S., Bond J.R., 1990, PhRvD, 42, 3936
- Scoccimarro et al. (2004) Scoccimarro, R., Sefusatti, E., & Zaldarriaga, M. 2004, PhRvD, 69, 103513
- Sefusatti & Komatsu (2007) Sefusatti, E., & Komatsu, E. 2007, PhRvD, 76, 083004
- Seljak & Warren (2004) Seljak U., Warren M.S., 2004, MNRAS, 355, 129
- Sheth et al. (2001) Sheth R.K., Mo H.J., Tormen G., 2001, MNRAS, 323, 1
- Sheth & Tormen (1999) Sheth R.K., Tormen G. 1999, MNRAS, 308, 119
- Slosar et al. (2008) Slosar A., Hirata C., Seljak U., Ho S., Padmanabhan N., 2008, JCAP, 08, 031
- Smith et al. (2007) Smith R.E., Scoccimarro R., Sheth R.K. 2007, PhRvD, 75, 3512
- Smith (2008) Smith R.E., 2008, arXiv0810.1960S
- Spergel et al. (2007) Spergel, D. N., et al. 2007, ApJS, 170, 377
- Springel (2005) Springel, V. 2005, MNRAS, 364, 1105
- Stanek et al. (2008) Stanek, R., Rudd, D.,& Evrard, A. E. 2008, arXiv:0809.2805
- Taruya et al. (2008) Taruya A., Koyama K., Matsubara T., 2008, arXiv:0808.4085
- Tinker et al. (2005) Tinker J.L., Weinberg D.H., Zheng Z., Zehavi I., 2005, ApJ, 631, 41
- Tinker et al. (2008) Tinker J.L., Kravtsov A.V., Klypin A., Abazajian K., Warren M.S., Yepes G., Gottlober S., Holz D.E., 2008, arXiv:0803.2706
- Verde et al. (2000) Verde, L., Wang, L., Heavens, A. F., & Kamionkowski, M. 2000, MNRAS, 313, 141
- Warren et al. (2006) Warren M.S., Abazajian K., Holz D.E., Teodoro L., 2006, ApJ, 646, 881
- White (2002) White, M. 2002, ApJS, 143, 241
- Willick (2000) Willick, J. A. 2000, ApJ, 530, 80
- Yadav & Wandelt (2008) Yadav, A. P. S., & Wandelt, B. D. 2008, Physical Review Letters, 100, 181301
- Zel’dovich (1970) Zel’dovich, Y. B. 1970, AAP, 5, 84

## Appendix A Initial conditions and Zel’dovich transients

The initial positions and velocities of the particles in our N-body simulations have been generated using the Zel’dovich approximation. This method introduces long-lasting artificial transients in the growth of perturbations which might alter the halo mass function at the high-mass end even at very late epochs (Crocce et al., 2006). It is therefore important to start the simulation at a sufficiently high redshift to ensure that all transients have decayed within the cosmic time at which the simulation output is used for science applications. Alternatively, less stringent requirements on the initial redshift are necessary if one uses second order Lagrangian perturbation theory to displace particles at the initial time (Crocce et al., 2006).

The simpler Zel’dovich approximation (which only requires the calculation of the gravitational potential) is much more widespread. In this case, a few authors have investigated how to compute the optimal starting redshift (see e.g. Lukić et al. 2007) as well as quantified the effects of the initial redshift on the halo mass function (Tinker et al., 2008) and on the internal properties of dark matter haloes (Knebe et al., 2009). These studies show that simulations of the concordance cosmology (and gaussian initial conditions) with initial redshifts of (depending on the simulations specifications) have converged to the correct solution by (at least for halo masses ). Even though our