Large-scale non-Gaussian mass function and halo bias: tests on N-body simulations

Constraining primordial non-Gaussianity offers a powerful test of the generation mechanism of cosmological perturbations in the early universe. While standard single-field models of slow-roll inflation lead to small departures from Gaussianity, non-standard scenarios allow for a larger level of non-Gaussianity (Bartolo et al. (2004) and references therein). The standard observables to constrain non-Gaussianity are the Cosmic Microwave Background (CMB) and the Large Scale Structure (LSS) of the Universe. A powerful technique is based on the abundance (Matarrese et al., 2000; Verde et al., 2001; LoVerde et al., 2008; Koyama et al., 1999; Robinson & Baker, 2000; Robinson et al., 2000) and clustering (Grinstein & Wise, 1986; Matarrese et al., 1986; Lucchin et al., 1988) of rare events, such as dark matter density peaks, as they trace the tail of the underlying matter distribution. Theoretical predictions on various observational aspects of non-Gaussianity have been extensively tested against N-body simulations, leading to different and sometimes conflicting results (Kang et al., 2007; Grossi et al., 2007; Dalal et al., 2007; Desjacques et al., 2008; Pillepich et al., 2008).

Dalal et al. (2007) and Matarrese & Verde (2008) showed that primordial non-Gaussianity affects the clustering of dark matter halos inducing a scale-dependent bias on large scales. Not only this effect has been already exploited to place stringent constraints on non-Gaussianity (Slosar et al., 2008; Afshordi & Tolley, 2008), but also it is particularly promising for constraining non-Gaussianity from future surveys, which will provide a large sample of galaxy clusters over a volume comparable to the horizon size (e.g., DES, PanSTARRS, BOSS, LSST, ADEPT, EUCLID) (Dalal et al., 2007; Carbone et al., 2008; Afshordi & Tolley, 2008; Seljak, 2008). Bartolo et al. (2005) showed that even for small primordial non-Gaussianities, the evolution of perturbations on super-Hubble scales yields extra contributions. The amplitude of these contributions is comparable to the forecasted errors of some planned surveys, opening up the possibility of measuring them.

In light of this, it is important to use N-body experiments to test the validity of theoretical predictions for halo-bias in non-Gaussian framework. Indeed, all proposed analytic biasing expressions have been derived in the extended Press-Schechter framework which assumes spherical collapse dynamics, sharp k-space filtering and Gaussian initial conditions. The validity of the extrapolation of the extended Press-Schechter approach to the non-Gaussian case can be tested independently by considering also the halo mass function. It is thus also important to test and calibrate on N-body simulations the predictions of the non-Gaussian halo mass function (Kang et al., 2007; Grossi et al., 2007; Dalal et al., 2007) and of the non-Gaussian halo bias simultaneously. This is what we set out to do here.

In this paper we start by reviewing the analytic predictions for the Gaussian and non-Gaussian halo abundance and clustering (§2). In §3 we describe the numerical simulations with Gaussian and non-Gaussian initial conditions. In §4 we present the test for the non-Gaussian mass function. In §5 and 6 we test the analytic predictions of Gaussian and non-Gaussian large scale bias against N-body simulations. In §7 we compare our results with the literature. Finally, we conclude in §8.

Deviations from Gaussian initial conditions are commonly parameterized in terms of the dimensionless $f_{\rm NL}$ parameter (Salopek & Bond, 1990; Gangui et al., 1994; Verde et al., 2000; Komatsu & Spergel, 2001):

where $\Phi$ denotes the gravitational potential and $\phi$ is a Gaussian random field. As noted by e.g., LoVerde et al. (2008), Afshordi & Tolley (2008) and Pillepich et al. (2008), different authors use different conventions. Here $\Phi$ denotes Bardeen’s gauge-invariant potential which, on sub-Hubble scales, reduces to the usual Newtonian peculiar gravitational potential but with a negative sign. In addition, there are two conventions for normalizing Eq. (1): the LSS and the CMB one. In the LSS convention $\Phi$ is linearly extrapolated at $z=0$. In the present paper we use this convention. In the CMB convention $\Phi$ is instead primordial: thus $f_{\rm NL}=[g(z=\infty)/g(0)]f_{\rm NL}^{CMB}\sim 1.3f^{CMB}_{NL}$, where $g(z)$ denotes the linear growth suppression factor in non Einstein-de-Sitter Universes.

The deviations from Gaussianity we are after become important on very large scales $k\la 0.03$h/Mpc and for massive halos. Therefore, one needs to perform N-body simulations on very large boxes, yet with enough resolution to identify massive virialized structures at different redshifts.

Suitable initial conditions have been set up following the method described in more detail in Grossi et al. (2008) (see also Grossi et al. (2007); Viel et al. (2009)). In brief, a random realization of a Gaussian gravitational potential, $\Phi_{\rm L}$, normalized to be the one linearly extrapolated at $z=0$, is generated in Fourier space, then it is inverse-Fourier transformed back to real space and added to the non-Gaussian term, $\Phi_{\rm NL}=f_{\rm{NL}}\left(\Phi_{\rm L}^{2}-\langle\Phi_{\rm L}^{2}\rangle\right)$. The resulting field $\Phi_{L}+\Phi_{\rm NL}$ that is linear and at $z=0$, is transformed back in Fourier space. We eventually modulate the power-law spectrum using the transfer function and compute the corresponding density field, which we then scale back to the initial conditions redshift ($z=60$). The corresponding gravitational potential is then used to displace particles according to the Zel’dovich approximation. This method allows one to simulate non-Gaussian models having power spectra which are all consistent with that of the Gaussian case and was already used by Viel et al. (2009).

In order to check the reliability of the initial conditions generation, we have performed a specific test: using $256^{3}$ particles in a box of size 1000 Mpc$/h$, primordial density fields (extrapolated linearly at $z=0$) were generated and smoothed using spherical top-hat filters of different radii $r_{s}=4$, $6$, $8$, $12$ Mpc/$h$. The smoothed skewness was then extracted from the fields and compared to the analytical prediction for $f_{\rm NL}=100,200,500,1000$, as shown in Figure 1.

The set of simulations used in this work assumes the ‘concordance’ $\Lambda$CDM model. We fix the relevant parameters consistently with those derived from the analysis of the WMAP 5-year data (Komatsu et al., 2008): $\Omega_{m0}=0.26$ for the matter density parameter, $\Omega_{\Lambda 0}=0.74$ for the $\Lambda$ contribution to the density parameter, $h=0.72$ for the Hubble parameter (in units of $100$ km s^-1 Mpc^-1). The initial power spectrum adopts the Cold Dark Matter (CDM) transfer function suggested by Eisenstein & Hu (1999), has a spectral index $n=0.96$ and is normalized in such a way that $\sigma_{8}=0.8$. In all experiments, performed using the GADGET-2 numerical code (Springel, 2005), switching off the hydrodynamical part, we consider a box of (1200 Mpc/$h)^{3}$ with $960^{3}$ particles: the corresponding particle mass is then $m\approx 1.4\times 10^{11}h^{-1}M_{\odot}$. The gravitational force has a Plummer-equivalent softening length of $\epsilon_{l}=25h^{-1}$ kpc. The runs produced 15 outputs from the initial redshift ($z=60$) to the present time. The 5 simulations consider different amounts of primordial non-Gaussianity, parametrized by the $f_{\rm NL}$ parameter: $f_{\rm NL}=0$ (i.e. the reference Gaussian case) and $f_{\rm NL}=\pm 100,\pm 200$. The catalogues of dark matter haloes are extracted from the simulations using the standard friends-of-friends algorithm adopting a linking length of $0.2$ times the mean interparticle distance; only objects with at least $32$ particles are considered.

We thus measure the halo bias in the simulations as

where $P_{hm}(k,z,M)$ denotes the cross-power spectrum of dark matter with halos of mass $M$ at scale $k$, and for the simulation snapshot at redshift $z$. Similarly $P_{mm}(k,z)$ denotes the dark matter power spectrum. Here and hereafter the subscript $s$ denotes quantities measured from the simulation.

In principle, the quantity one is interested in would be the bias of the halo power spectrum $b_{hh}=\sqrt{P_{hh}/P_{mm}}$, but $b_{hm}$ is a less noisy quantity (the shot noise of the finite number of halos is greatly suppressed in the estimate of $P_{hm}(k)$). The quantity $b_{hm}$ is not guaranteed to be identical to $b_{hh}$ if bias has a stochastic component that does not correlate with the matter density field. In Fig. 2 we show that this is not the case and that there is good agreement on large scales between $b_{hh}$ and $b_{hm}$, justifying using the less noisy $b_{hm}$ as an estimator for $b_{hh}$.

We compare the halo mass function of the non-Gaussian simulations with the theoretical predictions of Eqs. (5), (6) and (7) that is, including our ansatz for the the non-spherical collapse correction: $\delta_{c}\longrightarrow\sqrt{q}\delta_{c}$. For clarity we show here the non-Gaussian to Gaussian mass function ratio, i.e. the factor $R_{NG}(M,z)$. The comparison between theory and simulations results is shown in Fig. 6 for a few redshift snapshots and for $f_{NL}=\pm 100$, and in Fig. 7 for $f_{NL}=\pm 200$ for the same redshifts. Dashed lines are the mass function of Matarrese et al. (2000)–Eq. (6) – and dot-dashed lines are that of LoVerde et al. (2008)–Eq. (7)–.

Contrary to Kang et al. (2007) and Dalal et al. (2007), we conclude that both Matarrese et al. (2000) and LoVerde et al. (2008) are good descriptions of the non-Gaussian correction to the mass function, once the correction for non-spherical collapse is included.

Figs. 6 and 7 seem to indicate that LoVerde et al. (2008) may be a better fit for small masses and Matarrese et al. (2000) at high masses. This is not surprising: the Edgeworth expansion works well away from the extreme tails of the distribution (i.e. for moderate $\delta_{c}/\sigma_{M}$), while the saddle-point-approximation used in Matarrese et al. (2000), is expected to work better at the very tails of the distribution (very high $\delta_{c}/\sigma_{M}$). We expect that the mass function of Matarrese et al. (2000) will be a better fit at very high masses or larger $f_{NL}$. This will be further explored in future work.

The large-scale, linear halo Eulerian bias for the Gaussian case is (Mo et al., 1997; Scoccimarro et al., 2001)

where $q=0.75$ and $p=0.3$, account for non-spherical collapse and are a fit to numerical simulations. Here, $\sigma_{M}$ denotes the rms value of the dark matter fluctuation field smoothed on a scale $R$ corresponding to the Lagrangian radius of the halos of mass $M$; $z_{f}$ denotes the halo formation redshift and $z_{o}$ denotes the halo observation redshift. As we are interested in massive halos, we expect that $z_{f}\simeq z_{o}$. As the non-Gaussian halo bias correction is proportional to $b^{G}-1$, the dependence of $b^{G}$ on whether the selected halos underwent a recent merger (i.e. $z_{f}\sim z_{o}$) or are old halos (i.e. $z_{f}\gg z_{o}$) affects the amplitude of the non-Gaussian correction (Slosar et al., 2008; Carbone et al., 2008). Before we trust our simulation to accurately describe the non-Gaussian halo bias we check whether we recover the Gaussian one and whether the linear halo bias approximation is a good description for the scales, redshifts and mass-ranges we are interested in. Gao et al. (2005) show that analytical predictions for the Gaussian halo bias are in reasonable agreement with simulations and that the bias for low-mass halos shows strong dependence on formation time but high mass halos (the ones we are interested in) do not. The halo bias for the Gaussian simulation and the comparison with the theory prediction is shown in Fig. 8. Except for the Gaussian halo bias $b_{0}^{G}\equiv 1+\delta_{c}(z_{o})/(\sigma_{M}^{2}D(z_{o}))$ defined in Efstathiou et al. (1988) and Kaiser (1984) indicated by the dotted (blue) line, the simulated data agree with the theoretical expectations at different redshifts. In particular, in Fig. 8, the black solid line represents the total Gaussian bias of Eq. (11), the dashed (red) line represents the contribution from the first line of Eq. (11), and, finally, the dot-dashed (green) line is $1+q(b_{0}^{G}-1)$. The small difference when using $z_{f}\sim z_{o}$ implies that, for the Gaussian halo bias of very massive halos ($M\mathrel{\raise 1.16pt\hbox{$>$}\kern-7.0pt\lower 3.06pt\hbox{{$\scriptstyle\sim$}}}10^{13}$M_⊙), it is reliable to assume that the correction from the “non-spherical collapse” can be encapsulated in the factor $q$ in front of $\delta_{c}(z_{o})/(\sigma_{M}^{2}D(z_{o}))$.

Eq. (8) shows that the redshift and scale dependence of the non-Gaussian correction can be factorized as a term that depends only on redshift and one that depends only on $k$ and $M$. The $M$-dependence is expected to be very weak at large scales ($k<0.03$ h/Mpc). Here we will test the mass, scale and redshift dependence of the non-Gaussian halo bias and we calibrate its normalization on the simulations.

In Fig. 9 we show the dependence on halo mass of $\Delta b/b_{L}$. We define the quantity

where $(\Delta b/b_{L})_{theory}$ is given by Eq. (8). To study the mass dependence, we evaluate the theory at fixed mass $\widehat{M}=10^{14}M_{\odot}$. We compute the bias from the simulations taking halos in six different mass bins. Fig. 9 includes only scales $k<0.03$ $h$/Mpc, different lines correspond to different redshift snapshots between $z=0$ and $z=1.5$. As expected, there is no noticeable dependence on halo mass.

Having confirmed the expected weak dependence on halo mass for masses $M>10^{13}M_{\odot}/h$ and on scales $k<0.03$ $h$/Mpc, we can study the redshift and scale dependence of $\Delta b/b_{L}$, considering halos of different masses above $10^{13}M_{\odot}/h$.

The redshift dependence of $\Delta b/b_{L}$, $(\Delta b/b^{G}_{L})_{s}[2f_{NL}\alpha_{M}(k)q]^{-1}$ is shown in Fig. 10 where $M>10^{13}M_{\odot}/h$ and scales $k<0.026$ $h$/Mpc were used.

In applying the correction $\delta_{c}/\sigma_{M}^{2}\longrightarrow q\delta_{c}/\sigma_{M}^{2}$ to $\Delta b/(b_{G}-1)$ we have actually corrected $b^{G}_{0}$, i.e. we have employed the same approximation used for the green dot-dashed line of Fig. 8, giving Eq. (9). Eq.(9) in fact is only the consequence of our correction to the Gaussian halo bias. Note that the approximation $z_{f}\sim z_{0}$ we employed here is expected to hold for rare–massive–halos and Fig. 8 shows that this is a good approximation. A detailed study of the dependence of the non-Gaussian halo bias correction on the formation redshift of the halos will be presented elsewhere.

There seems to be an indication that the $q$-correction factor for the large-scale bias correction may slightly depend on the value of $f_{NL}$: in particular the figure shows that it could be slightly smaller than $q$ for $f_{NL}$ large and negative and smaller for $f_{NL}$ large and positive. This is not unexpected: the presence of non-Gaussianity may alter the dynamics of non-spherical collapse (e.g., through tidal forces – see e.g., Desjacques (2008) – or by significantly changing the redshift for collapse with respect to the Gaussian case). At this stage, however, this trend is not highly significant and further study will be left to future work.

Finally, we show the scale dependence of Eq. (9), $(\Delta b/b^{G}_{L})_{s}[2f_{NL}\delta_{c}(z)q]^{-1}$, in Fig. 11. The thin lines correspond to different redshifts and the thick black line to their average. The dotted line is the theory prediction with $q=0.75$. Note that there is an excellent agreement on the scales of interest, e.g., $k<0.03$ $h$/Mpc. On smaller scales the effect of non-Gaussianity is very small and the measurement become extremely noisy. These results are in qualitative agreement with the findings of Pillepich et al. (2008).

We conclude that Eq. (9), with $q\sim 0.75$, provides a good fit to non-Gaussian simulations.

After the discussion of §2, it should be clear that if $f_{NL}^{CMB}$ were used in the theoretical predictions or $S_{3}$ were not linearly extrapolated to $z=0$, then any constraints on non-Gaussianity so obtained would have to be re-scaled by a factor $\sim 1.3$. This seems to be the case of some work in the literature. On the other hand the q-correction factor effectively introduces a re-scaling of a factor $\sim 0.75^{3/2}=0.65$ for the mass function case and $0.75$ for the bias case. It is a coincidence that, for the halo bias, $1.3\times 0.75\sim 1$, thus the $f_{NL}$ normalization mistake cancels out with the spherical collapse approximation error. This fortuitous cancellation does not happen to the same level in the mass function $0.75^{3/2}\times 1.3\sim 0.8$, explaining perhaps some of the claimed discrepancy of the simulations with the analytic mass function predictions, and the claimed agreement with the halo bias predictions. Another possible source of inaccuracy would be an inconsistent treatment of the redshft evolution of $\delta_{c}$ and $S_{3}$ (see discussion is §2).

In Fig. 12 we compare our theoretical predictions with the results presented in Pillepich et al. (2008) and Desjacques et al. (2008). The left panel shows our simulation results at $z=0$ for $f_{NL}=\pm 100$ and our theoretical predictions. Additionally we show the fit presented by Pillepich et al. (2008), Eqs. (8) and (9), evaluated for the suitable values of $f_{NL}$ accounting for the different notations for $f_{NL}$. We also adopt our cosmological parameters when converting $\sigma_{M}$ to $M$. The right panel shows the simulation results presented in Pillepich et al. (2008) at z=0 and their fitting formula at $z=0$ and $z=1$. We over plot our theoretical models evaluated for their cosmological parameters and for the corresponding values of $f_{NL}$. Moreover we add the data points from Desjacques et al. (2008) for $z=0$ and $z=1$, suitably rescaled by the differences of the $f_{NL}$ value used. The mass function fits of Fig. 12 differ for large masses, in the regime where simulations errors become large; the fits are however consistent given the individual points error-bars.

Our theoretical formulae for the non-Gaussian mass function (Eqs. 5, 6, 7) and for the non-Gaussian halo bias (Eq. 9) are physically motivated expressions that have been tested on N-body simulations. They have the advantage over fitting formulae that they can be more robustly interpolated and extrapolated to cosmologies and parameters that have not been directly simulated and they are more robust over parameters ranges where the simulations have low signal-to-noise. Compared to simple fitting formulae, Eq. 6,7 and 9 have the disadvantage that they require the calculation of some numerical integrals. To overcome this, we supply tabulated values for $S_{3,M}^{(1)}$, $\sigma_{M}$ and $\alpha_{M}$ for a WMAP5 cosmology in the range of interest at www.ice.csic.es/personal/verde/nongaussian.html.

The $q$-correction we find here has implications for previously reported and forecasted constraints on non-Gaussianity. In Table 1 we report present and forecasted constraints on $f_{\rm NL}$ from the literature rescaled to $f_{\rm NL}^{CMB}$ and corrected for our factor $q$.

This confirms that constraints on $f_{NL}$ achievable using the non-Gaussian halo bias are competitive with CMB constraints ($f_{NL}\sim 5$ for Planck and $f_{NL}\sim 3$ for a CMBPol-type mission, Babich & Zaldarriaga (2004); Yadav, Komatsu & Wandelt (2007)).

We have considered $(1.2$ Gpc/h)³ size and $960^{3}$ particles N-body simulations with non-Gaussian initial conditions, with non-Gaussianity parameter $f_{NL}=\pm 100$, $f_{NL}=\pm 200$ and a reference Gaussian simulation ($f_{NL}=0$). The clustering properties and the abundance of the simulation’s halos were then compared with independent simulations and theoretical predictions. We find good agreement between different simulations, indicating that the initial conditions set-up is under control. We find that the Press-Schechter-based description of the non-Gaussian correction to the Gaussian mass function of Matarrese et al. (2000) and LoVerde et al. (2008) is a good fit to the simulations, provided that:
a) The Press-Schechter-based description is used to compute the ratio between Gaussian and non-Gaussian mass function.
b) The critical density $\delta_{c}$ is corrected to account for non-spherical collapse dynamics.
This is summarized in our Eq. (5) and in Eqs. (6) and (7) for the non-Gaussian mass functions of Matarrese et al. (2000) and LoVerde et al. (2008), respectively. For large thresholds this correction is equivalent to a re-scaling of the spherical collapse threshold $\delta_{c}\times\sqrt{q}$ where $q=0.75$. The $q$-correction is thus equivalent to a reduction of $f_{NL}$ by a factor $\sim 1.5$ because in the mass function, to leading order $f_{NL}$ multiplies $\delta_{c}^{3}$.

We find that the non-Gaussian halo bias prescription of Dalal et al. (2007), Matarrese & Verde (2008), Slosar et al. (2008) and Afshordi & Tolley (2008) provides a good description of the scaling of the large-scale halo clustering of the simulations. In particular, we have tested separately the predicted redshift, scale and $f_{NL}$ dependence. The overall amplitude of the effect, however, should be corrected by a factor $\sim q$ which can also be understood in the context of ellipsoidal collapse or as a modification to the excursion set ansatz and the sharp-k space filtering (see Eq. (9)). There is an indication that this correction may be slightly dependent on $f_{\rm NL}$. This is not unexpected, but the signal-to-noise of the effect is too small in the current simulations to draw robust conclusions. We also find that on large ($k<0.03$ $h$/Mpc) scales, as expected, the fractional correction to the non-Gaussian halo bias is independent of mass. On smaller scales a dependence on mass is expected, but the simulations do not have sufficient signal-to-noise to verify it. The $q$–correction to the non-Gaussian halo bias modifies current and forecasted constraints reported in the literature as indicated in our Table 1.

The formulae we presented here for the non-Gaussian mass function (Eq. 5, 6,7) and non-Gaussian halo bias (Eq. 9) are physically motivated expressions which provide good fits to a suite of N-body simulations. As such, they can be more robustly interpolated and extrapolated than simple fitting functions (in www.ice.csic.es/personal/verde/nongaussian.html we provide useful quantities for ease of use of these equations). We confirm that the non-Gaussian halo bias offers a robust and highly competitive test of primordial non-Gaussianity.

Computations have been performed on the IBM-SP5 at CINECA (Consorzio Interuniversitario del Nord-Est per il Calcolo Automatico), Bologna, with CPU time assigned under an INAF-CINECA grant, and on the IBM-SP4 machine at the “Rechenzentrum der Max-Planck-Gesellschaft” at the Max-Planck Institut fuer Plasmaphysik with CPU time assigned to the “Max-Planck-Institut für Astrophysik” and at the “Leibniz-Rechenzentrum” with CPU time assinged to the Project “h0073”. The authors thank A. Pillepich and C. Porciani for dicussions and for making available the simulations outputs for comparison. We also thank N. Dalal, U. Seljak, A. Riotto and M. Maggiore for discussions. LV is supported by FP7-PEOPLE-2007-4-3-IRG n. 202182 and CSIC I3 grant n. 200750I034. CC is supported through a Beatriu de Pinos grant. This research was supported by the DFG cluster of excellence Origin and Structure of the Universe. SM and LM acknowledge partial support by ASI contract I/016/07/0 “COFIS”, ASI-INAF I/023/05/0, ASI-INAF I/088/06/0 and ASI contract Planck LFI Activity of Phase E2. LV thanks the Santa Fe cosmology workshop 2008 , LV and CC thank the Galileo Galilei Institute for theoretical physics in Florence, where part of this work was carried out, and INFN for partial support. LV thanks B. Wandelt for discussions.