Higher moments of primordial non-Gaussianity and N-body simulations

Higher moments of primordial non-Gaussianity and N-body simulations

[    [    [

We perform cosmological N-body simulations with non-Gaussian initial conditions generated from two independent fields. The dominant contribution to the perturbations comes from a purely Gaussian field, but we allow the second field to have local non-Gaussianity that need not be weak. This scenario allows us to adjust the relative importance of non-Gaussian contributions beyond the skewness, producing a scaling of the higher moments different from (and stronger than) the scaling in the usual single field local ansatz. We compare semi-analytic prescriptions for the non-Gaussian mass function, large scale halo bias, and stochastic bias against the simulation results. We discuss applications of this work to large scale structure measurements that can test a wider range of models for the primordial fluctuations than is usually explored.

a]Saroj Adhikari, a]Sarah Shandera, b]Neal Dalal Prepared for submission to JCAP IGC-14/2-1

Higher moments of primordial non-Gaussianity and N-body simulations

  • Institute for Gravitation and the Cosmos, The Pennsylvania State University, University Park, PA 16802, USA

  • Department of Astronomy, University of Illinois, 1002 W. Green St., Urbana, IL 61801, USA

E-mail: sza5154@psu.edu, shandera@gravity.psu.edu, dalaln@illinois.edu



1 Introduction

Recently, the Planck satellite mission reported tight constraints on primordial non-Gaussianity from measurements of the Cosmic Microwave Background (CMB) [1]. They find the amplitude of the local, equilateral and orthogonal type bispectrum to be , , and respectively at level. While these results show that the primordial fluctuations were remarkably Gaussian, they still leave room for interesting signatures of primordial physics to be found in statistics beyond the power spectrum: minimal, single field models predict non-Gaussianity that is about two orders of magnitude below these constraints. Non-Gaussianity is such an informative tool for inflationary physics that it is crucial to push observational bounds as far as possible. A number of forecasts have shown that the use of clustering data from future large scale structure surveys can constrain with [2, 3, 4, 5, 6, 7]. In addition, the galaxy bispectrum is a promising way to probe many shapes to the level [8, 9, 10]. The constraints from large scale structure (LSS) —whether consistent with the CMB measurements or in tension — will provide interesting and useful complementary results.

In the inflationary scenario, non-Gaussian signatures in the primordial density fluctuations depend on the details of interactions of the inflaton or other fields relevant for generating the perturbations; therefore, any detection (or the lack) of primordial non-Gaussianity tells us about the dynamics of those fields. Furthermore, if non-Gaussianity is detected, we expect to see patterns in the correlations that are consistent with perturbation theory (in contrast to apparently independent statistics at each order, for example). The simplest such pattern is in the relative amplitudes of higher order statistics, or, how the amplitudes of the correlation functions of the gravitational potential, , scale with order . The relative scaling of higher moments falls into a fairly narrow range of behaviors in inflationary models [11].

Thanks to the non-linearity of structure formation, the statistics of objects in the late universe contain information about the entire series of higher order moments of the initial density fluctuations. This is especially straightforward to see in number counts of galaxy clusters [12]. While information about higher moments is clearly non-trivial to extract when non-Gaussianity is weak, some constraints with current data are already affected by assumptions about higher order moments. For example, [13] reports constraints on primordial non-Gaussianity using a sample of 237 X-ray clusters from the ROSAT All-Sky Survey [14]. Their analysis clearly suggests that it is important to take the scaling of moments into account when deriving constraints on from cluster number counts. The results presented in this paper will be useful for further studies in the same direction. Tighter constraints could likely be achieved with current data by combining the X-ray clusters with SZ detected clusters that have been separately used to constrain non-Gaussianity [15, 16, 17]. In addition, the eROSITA survey [6] and third generation surveys detecting clusters via the Sunyaev-Zel’dovich effect are expected to provide much larger samples of clusters in the next few years. It will be interesting to revisit the non-Gaussian analysis when that data becomes available. Other large scale structure statistics, including the power spectrum and bispectrum of galaxies, also contain contributions from higher order primordial statistics. Constraints on non-Gaussianity from those observables are still being developed, and will ultimately be very powerful (constraints obtained to date include [18, 19, 20, 21, 22, 23, 24, 25, 26]). To make full use of them it is important to understand the signatures of the full spectrum of non-Gaussian models predicted from inflation.

In this paper, we will use a slight variant of the popular local ansatz of primordial non-Gaussianity in order to generate two different scalings of the -point functions and study how the scalings can be distinguished. For the usual local ansatz with hierarchical scaling, many groups have already performed simulations [27, 28, 29, 30, 31, 32] and found that the semi-analytic prescription of [12] for the non-Gaussian mass function works well. Here we will also test the validity of the mass function proposed in [11] and used in [13] for another scaling that is well motivated by inflationary models. In addition to looking at the mass function, we will also look at a signature of the shape (momentum dependence) of local type non-Gaussianity using the scale dependent halo bias [33], and at the stochastic bias on large scales [34]. Subleading contributions to the bias and leading contributions to the stochastic bias are sensitive to non-Gaussianity beyond the skewness.

The rest of the paper is structured as follows: in Section 2 we introduce our two field model and the scaling of the higher moments. In Section 3 we present the analytical calculations for the mass function of halos, bias, and stochasticity (with details relegated to the appendices). In Section 4 we describe our simulations, and then compare the output from simulations to the analytical calculations in Section 5. Finally, we conclude in Section 6.

2 Model

In the usual local ansatz, the Bardeen potential has non-Gaussian statistics thanks to a contribution from a non-linear, local function of a Gaussian field :


Here parametrizes the size of the non-linear term and the level of non-Gaussianity. Non-Gaussianity of this type is usually thought of as produced by a light field that is not the inflaton [35, 36, 37, 38, 39, 40], and in fact cannot be generated by single field inflation proceeding along the attractor solution with modes in the Bunch-Davies vacuum [41, 42].

The two-point statistics (the power spectrum) of the fluctuations are well measured from CMB observations for about three decades in scale [43, 44, 45]. The amplitude of the three point function gives the skewness of the distribution of , and is also well-constrained by the CMB [1]. Higher order correlation functions are more difficult to measure in the data. We would like to see if we can get some handle on the structure of the scaling of these higher order statistics of . For this purpose, we define the dimensionless moments :


where the subscript indicates that we take the connected part of the -point function.

Two scalings:

For the local model given by Eq. (2.1), the moments scale approximately as


where comes out of combinatorics. The numerical coefficients in this scaling are not quite precise because of the difference in integrals over momenta at each order (see Appendix A), but the parametric dependence on the amplitude of the skewness and the total power is fixed. This is the behavior of the moments (which we label ‘hierarchical’ scaling) for not too large i.e. when , where . However, if , the moments scale differently:


where . In the single source case, this is far too non-Gaussian to be consistent with observations. Therefore, we consider the following two source model [46, 34, 11]:


where the two Gaussian fields and are uncorrelated i.e. .

To express the correlations in the gravitational potential in terms of the amplitude of fluctuations in just one of the source fields we define


the ratio of the contribution of the Gaussian part of the field to the total Gaussian power. is defined by


For simplicity we assume both Gaussian components have constant, identical spectral indices but different amplitudes: with , and similarly for (with the same index but a different amplitude ). With this choice, is a scale-independent constant.

The power spectrum, bispectrum, and trispectrum in our model are:


where we have defined as


Here, is the infrared cutoff for a boxsize of . In general, we do not know of the size of the universe beyond our observable universe, but for our purposes, the simulation box size is the natural choice. is the scale leaving the horizon at the initial epoch [47]. The above integrals converge for large values of . For the computations to compare with simulations we set , where is the number of particles in a simulation. To arrive at the expressions quoted above, we have used


In Equations (2.8), (2.9) and (2.10) we have included terms that are usually sub-dominant in the case of single field, weakly non-Gaussian local ansatz. In our model, these terms are important when the field is strongly non-Gaussian. To discuss the observational constraints on , let’s consider , and . Observational constraint from small non-Gaussianity can be satisfied by making and , in which case the Gaussian contributions () dominate the total power spectrum in Eq.(2.8) as well. Notice that the non-Gaussian contribution to the power can shift the spectral index slightly, so that when the field is strongly non-Gaussian the measured spectral index is close to, but not identical to, the spectral index of the Gaussian components111The new, integral terms have slightly different shapes than the usual terms. The difference, however, is small—approximately described by which has a weak dependence on . These terms are also infrared divergent. For the purpose of comparing with the results from N-body simulations, the box size of the simulation provides a natural cutoff [30]. We will only look at quantities well enough inside the volume that the arbitrary size doesn’t affect our results., .

If and , then one can generate the feeder scaling in Eq.(2.4) without being inconsistent with the current observations of power spectrum and bounds on non-Gaussianity. The feeder scaling dominates when the condition is satisfied. This scaling, or a hybrid between hierarchical and feeder scaling, arises in particle physics scenarios where a second, non-Gaussian field couples to the inflaton and provides an extra source for the fluctuations [11, 48]. However, those scenarios differ from the model here because they most often generate bispectra not of the local type222One can also argue on statistical grounds that our observable universe is unlikely to have local non-Gaussianity of the type written in Eq.(2.5) with large if inflation lasts much longer than 55 e-folds [49, 50]).. Here we are using the two field, local model primarily as a phenomenological tool, easy to implement in N-body simulations, to generate the feeder-type scaling of moments rather than as the output of a particular inflation model. The mass function is sensitive only to the integrated moments, not the shape, so this is a useful test of how different scalings affect the number of objects in a non-Gaussian cosmology.

3 Abundance and Clustering Statistics

In this section we present the analytic predictions for the effect of locally non-Gaussian, two source initial conditions on the abundance and clustering of dark matter halos.

At large scales, the evolution of the density contrast is well described by linear perturbation theory, and the density contrast in Fourier space at redshift is given by , where


Here is the transfer function, is the growth function, is the Hubble scale today, and is the energy density in matter today (compared to critical density). The smoothed density contrast, similarly, is , where


is the smoothing function (here the Fourier transform of the real space top-hat). Note that we will generally suppress the dependence in and in this paper, and usually write and only. We can now compute the connected n-point functions of the smoothed density contrast in real space:


and therefore the dimensionless moments of smoothed density fields: . Note that the dimensionless moments are redshift independent. Eq.(3.3) can be separated into two components: (i) of and (ii) of , corresponding to contributions to the hierarchical and feeder scalings respectively. Some of these integrals, with a brief summary of our method to evaluate them, are given in Appendix A.

3.1 Mass function

We follow previous studies of non-Gaussian mass functions in that we calculate the ratio of the number density in the presence of non-Gaussianity, , to the number density for Gaussian initial conditions, , for a particular halo mass at some redshift using the Press-Schechter formalism [51]. The fractional volume of dark matter in the collapsed structures (halos) is proportional to , where is the critical value of the smoothed density contrast above which the dark matter in a region collapses to form halos and is a probability density function (PDF). Here we have written the smoothing scale in terms of the mass rather than the smoothing radius ; they are simply related by , where is the mean matter density of the universe. Then, the number density (mass function) is given by:


For a Gaussian PDF, one can easily perform the integration to find a prediction for the mass function. However, the result is known to be only an approximation and in practice Gaussian mass functions are calibrated on simulations [52].

To apply the method above in case of non-Gaussian initial conditions we need a non-Gaussian PDF. The Petrov expansion [53] (which generalizes the Edgeworth expansion [54, 55]) expresses a non-Gaussian PDF as a series in the cumulants of the distribution. In terms of ’s, this is:


where , and ’s are the Hermite polynomials defined as . The second sum is over the non-negative integer members of the set that satisfy


For each set . This series can be integrated term by term to obtain , and with Eq.(3.4) gives ratio of non-Gaussian mass function to Gaussian mass function [13]




The two superscripts indicate that the set of terms that are of the same order depends on whether higher order cumulants have hierarchical () or feeder () scaling. Formally grouping terms assuming the scalings in Eq(2.3) or Eq.(2.4) are exact gives (for )


The prime stands for the derivative with respect to the halo mass . In the hierarchical case, the sets still satisfy Eq.(3.6), while for feeder scaling the are the sets of non-negative integer solutions to .

Eq.(3.7) assumes that the two point statistics of the smoothed linear density contrast (i.e ) are the same for the non-Gaussian and the Gaussian cases (not that the Gaussian contributions to are the same). However, this is difficult to maintain at all scales in our simulations. We require the two point clustering statistics to match at a particular scale , but do not correct for the shift to the spectral index coming from the non-Gaussian term in the power spectrum. As a result, on scales other than we need to make a distinction between for the non-Gaussian cosmology and for the Gaussian cosmology. In that case, there is an extra (mass or scale dependent) factor in the ratio , where


This factor is typically quite close to one. For example, in the mass range () at , for single field case. For our smallest feeder scaling simulation (with ), in this mass scale and redshift range, . The factor deviates away from unity more at larger redshifts and at mass scales far from .

In addition, the derivation above assumes the same constant of proportionality between and regardless of the level of non-Gaussianity. The standard Press-Schechter constant of proportionality is two, but for the non-Gaussian case it is reasonable to fix the constant by requiring . Gaussian and non-Gaussian cosmologies with identical will have slightly different normalization factors. This factor shifts further away from 2 as the level of non-Gaussianity in the initial conditions is increased. Integrating various truncations of the expanded PDFs indicates we expect a difference from 2 between about 0.5% and 2% for the amplitudes of non-Gaussianity we consider in this work. So, we will introduce an extra factor that multiplies our analytical mass function to fit with the simulation results.

3.2 Large scale bias

On large scales, where density fluctuations are in the linear regime, the clustering of halos is expected to trace the clustering of the underlying matter field. The proportionality constant relating halo clustering to matter clustering is called the halo bias. Local type non-Gaussianity can modify halo bias, compared to Gaussian universes, by coupling the amplitude of short wavelength modes to that of long wavelength modes. Since the coupling occurs in the gravitational potential field (with constant amplitude), rather than in the density field, local type non-Gaussianity introduces a new, scale-dependent term relating the power spectrum of halos to the power spectrum of the linear dark matter field. On large scales (small wave numbers) the non-Gaussian term can dominate and the analytic prediction for the bias is relatively simple.

The potential use of the halo (or galaxy) bias as a probe of primordial non-Gaussianity was first demonstrated in [33]. An analytic derivation capturing the first order effect of non-Gaussian initial conditions had been presented much earlier in [56] and clarified and improved following the Dalal et al result in [57, 58, 59]. The halo bias in models with two sources for the primordial fluctuations was considered in [34, 60], and in [61] which gives some theoretical predictions for the model we consider here. In addition, [62] previously performed N-body simulations for a model where the kurtosis was larger than in the single field case by a factor of two. That work corresponds to our scenario with . Here we are primarily interested in values of that are very small so that the scaling is feeder type, but we also consider cases with which have intermediate scaling. Measurements of the power spectra of several different galaxy populations have been used to place constraints on primordial non-Gaussianity, at roughly the level, depending on the population and treatment of systematic errors [18, 19, 20, 21, 22, 23, 24, 25, 26].

The general form of the large scale bias for our two-field model Eq.(2.5) can be calculated using the peak-background split formalism [19, 62]. Appendix C has the peak-background-split calculation that results in Eq.(C.10) as the expression of bias. In case of small, local non-Gaussianity ( for the case with a second, strongly non-Gaussian field), the expression for large scale bias reduces to:

Here we have used subscripts on the bias coefficients to emphasize that the scale-independent term depends on all sources of the fluctuations (), while the scale-dependent term depends only on those sources with a primordial non-Gaussian component (). Recall that in the more frequently quoted expression with a single, non-Gaussian source, these two bias coefficients are both equal to the Gaussian bias plus scale-independent corrections proportional to the level of non-Gaussianity. When multiple sources are present, the first bias coefficient can be split into terms which are include the Gaussian bias for each source, while the bias coefficient in the second, scale-dependent term is to lowest order the Gaussian bias for the non-Gaussian source. In addition, Eq.(3.2) demonstrates that the large scale halo bias is a probe of for the hierarchical scaling and a probe of for the feeder scaling. In other words, the dominant contribution to the non-Gaussian bias is proportional to the amplitude of the bispectrum as expected.

To compare the analytic expressions against simulation results, we will use the size of the simulation volume to truncate integrals in the infrared and will not fit the power spectrum on scales very close to the simulation box size ()333First, for i.e modes approaching the scale of the box size of our simulations, the sample variance is large. Second, the bias for feeder scaling depends on which has a sharply decreasing behavior near of the simulation; this effect runs with the box size. If we are interested in scales near the () for our feeder scalings, then we will either have to do simulations with larger box sizes, or generate initial conditions differently such that when averaged over many realizations of the initial conditions we don’t get this feature..

3.3 Large scale stochasticity

Using the peak-background-split method, we can also calculate the halo-halo power spectrum in terms of the underlying matter distribution . We can then define the large scale stochasticity, or cross-correlation coefficient as


The calculations and the corresponding expressions for and are given in Appendix C. Peak-background split calculations suggest that in the large scale limit (small ), gives the fraction of power in the initial conditions from contributions of the field with non-Gaussianity (i.e. . For small non-Gaussianity, in case of hierarchical scaling, this is , while in case of feeder scaling, this is . For the single field cases, there is no large scale stochasticity (). We will test these expectations against results from simulations. Large scale stochasticity for the Gaussian case and non-Gaussian case with hierarchical scaling was discussed in [62]; their comparison with simulations showed a small discrepancy between their analytic expression and the simulation results. Our results below indicate one possible source for at least some of that discrepancy.

4 Simulations

The simulations for this project were performed using the popular GADGET-2 code [63]. The initial conditions were generated using second order Lagrangian perturbation theory (2LPT) [64]. We used the code from [64] for generating local type (single field) non-Gaussian initial conditions using 2LPT and modified the code—discussed next—to generate initial conditions for our two field model.

Initial conditions:

First, two Gaussian random fields ( and ) were generated using the power spectrum of our fiducial Gaussian cosmology (). The amplitude of the power spectrum for the and fields were multiplied by the appropriate factors and respectively. Then, the field was squared and multiplied by . The total non-Gaussian field was obtained by adding the three components as in Eq. (2.5). We ensured of the generated non-Gaussian field was that of our specified cosmology for all of our parameter sets by renormalizing the field by a factor of . These are the only modifications done to the 2LPT code. The rest of the code generated the required displacement and velocity fields as usual. The redshift of all our initial conditions is .

N-body simulations:

The public version of the GADGET-2 code was used to perform cosmological collisionless dark matter only simulations. All simulations were done with particles in a cube of side . This gives the mass of a single particle to be . The force softening length was set to be of the inter-particle distance. Simulations with Gaussian initial conditions ( and ) were also performed with the same seeds as the field (as it has the dominant contribution for small ) to compare with our feeder models; another set of Gaussian simulations was performed with and to compare with the hierarchical simulations. For each set of parameters listed in Table 1, we ran four simulations with different seeds. All simulation results reported in this paper are average over the four simulations. Similarly, the errors reported are the 1 standard deviation of the four simulations. The AHF halo finder [65] was used to identify halos which were then used to get the mass function of dark matter halos and power spectra of halos (halo-matter cross power spectrum and halo-halo autospectrum). In all our analyses, we only use halos with number of dark matter particles .

Parameter space of simulations:

Since our method of generating the feeder scaling produces a slightly different bispectrum shape than the hierarchical case, the scale dependence of for the two scalings also differ. For comparison purposes, we will define at the scale of corresponding to a halo mass of , as the ratio:


In Table 1, we list the parameter sets of that we have simulated. Notice that these parameter sets are not consistent with the small reported by the Planck mission from bispectrum measurements. However, it is necessary to use parameter sets with larger values of in order to get useful results from N-body simulations.

Name (Eq:(4.1))
M993 993 0.1 50000 0.290 2.144
F677 677 0.00005 0.198 4287
H500 500 1 500 0.145 0.0027
M384 384 0.11925 20620 0.112 0.5089
F215 215 0.00003 0.063 2923
F122 122 0.00003 0.036 1933
H99 99 1 99 0.029 0.0001
F70 70 0.00003 0.020 1303
Table 1: Parameter space of our simulations. In the first column, we name our models following a simple naming convention. The first letter of the name stands for the type of scaling of the model: F stands for feeder scaling, H stands for hierarchical scaling, M stands for mixed scaling (when neither component is negligible). The number following is the approximate , also listed in the second column. For example: F70 means that the scaling of the model is feeder and has . The quantity in the third column is defined in Eq.(2.6), and gives the ratio of power in the Gaussian part of the contribution (), to the total Gaussian power (). The second last column is the dimensionless skewness computed at a halo mass scale . The last column is the ratio of the dimensionless skewness from the feeder contribution to that of the hierarchical contribution and indicates the relative importance of the feeder term.

5 Results and Discussion

Mass Functions:

The hierarchical scaling has been considered a number of times already with the prediction from the Edgeworth-series formalism [12] providing good fit to the outputs from simulations [27, 28, 29, 30, 31, 32]. Our focus will be on the feeder type scaling.

In the following sections, we will use the mass function truncated up to . From error analysis (Appendix B), we see that we gain little by adding higher terms for the feeder mass function. To this order, for the feeder case, the ratio of non-Gaussian to Gaussian mass function is:


where we have chosen to rewrite all the quantities in terms of the halo mass rather than and ; the dependence of the moments ’s and has been suppressed for clarity, as is done throughout the paper. We use this expression to fit to the results from N-body simulations.

Similarly, for hierarchical scaling of moments, the expression for the non-Gaussian mass function including terms upto is:


where the dependence of and has been suppressed for clarity. In the above mass function formulae for feeder scaling Eq.(5.1) and for hierarchical scaling Eq.(5.2), the expression for is calculated for a given model using Eq.(3.11) and the factors are fit for each of our simulation results. Both of these factors approach unity for small non-Gaussianity. To compute the cumulants necessary to calculate the dimensionless moments in the above formulae, we have used the Monte-Carlo method described in [30]. See Appendix A for details.

Let us now compare and discuss the results from simulations and calculations from the Edgeworth series formulation. First, as a check of our simulations, we did a purely hierarchical scaling parameter set—single field, simulation—which can be compared directly with previous works. Then, for the two source case, we study the parameter sets listed in Table 1 which allow feeder scaling as well as mixed scaling (i.e, neither term in Eqs.(2.9),(2.10) is negligible). In this section, we will present the simulation results and the corresponding Edgeworth mass functions. We analyzed our mass function results with a simple two parameter () chisquare minimization procedure. The errors reported in the best fit values increase the reduced chisquare of the fit by unity when added to the best fit values. We find that all simulations (H99, H500, F70, F122, F215) prefer a reduced and different values for , which are of the expected size and increase appropriately with . We will discuss the dependence of this factor on and the type of scaling later.

First, let us present our mass function fits. We use obtained by performing chisquare minimization together for H99, H500, F70, F122, F215 by forcing the same but allowing overall rescaling factors () for each case. This is a reasonable procedure if one interprets shifting as allowing departures from the assumption of spherical collapse, which should be relatively independent of the level of non-Gaussianity. However, the normalization requirement (which has not been analytically enforced) suggests that should depend on the level of non-Gaussianity. The top panel of Figure 1 has our mass function results for two simulations with hierarchical scalings.

123456Mass ()H500, single field
11.522.53Mass ()H99, single field
123456F215-0.5-0.3-0.10.1Mass ()
11.522.53F70-0.3-0.2-0.100.1Mass ()
Figure 1: Top: the simulation results and semi-analytic prediction Eq.(5.2) for hierarchical simulations. Left panel: () with . Right panel: () with . Bottom: The simulation results and semi-analytic prediction Eq.(5.1) for feeder simulations. Left panel: () and . Right panel: () with . In addition, we also plot the fractional difference between the simulation results and the prediction from the truncated mass function Eq.(5.1) but using (), and assuming that scale as Eq.(2.4). A negative value means that the analytic mass function overpredicts the simulation result.

The results for two feeder simulations are presented in the bottom panels of Figure 1. As discussed in Appendix B, the errors on our truncated feeder mass function are large compared to the hierarchical case with comparable . Also, the error becomes large at relatively small even for the feeder scaling case with . Therefore, based on our error evaluation, we do not expect that our feeder mass function describes the simulation results well for the more massive halos or at higher redshifts for which . We find that the feeder mass function formula Eq.(5.1) fits well our simulation results for F70 and F122 (see the left panel of Figure 2 for the F122 simulation results). Consistent with the error analysis of feeder mass function, Eq.(5.1) fits to F215 with are not equally good. With a larger , the F677 simulation is clearly not well fit by our truncated feeder mass function (see Figure 3). Finally, Figure 1 also shows the difference between the simulations and Eq.(5.1) assuming , , and that the moments scale exactly as in Eq.(2.4). That is, the bottom-most panels illustrate how calibrating on simulations shifts the purely analytic expectations for the non-Gaussian mass function. The dominant effect among these three factors is that of ; change in affects the non-Gaussian mass functions starting at . On the other hand, using different scalings of higher moments only change the expressions starting at and the factors modify the mass functions at a few percent level at most (see Figure 4).

In Figure 2, the right hand panel shows the fractional difference between the semi-analytic expression Eq.(5.1) (with the cumulants measured in the realizations and , fit) and the simulation results for the two feeder models: F122 and F215. We have plotted the fractional difference as a function of and the plotted points include simulation results from all three redshifts . The result can be interpreted as the error in the semi-analytic approach and qualitatively correlates with our analytic error analysis of the PDF: (i) the error at low is small, (ii) the error increases for higher but the error is smaller for smaller .

On a different note, by looking at the non-Gaussian mass function results from simulations only, we can verify that the F70 model is more non-Gaussian than the H99 model (compare top and bottom plots in Figure 1), even though the skewness of F70 is smaller than the skewness of H99. Similarly, we also see that the F215 model has comparable amount of non-Gaussianity as the H500 model. This verifies that the non-Gaussian mass function is sensitive to the total non-Gaussianity and the scaling of higher moments in the initial conditions.

11.522.53Mass ()F122
Figure 2: Left: the simulation results and the mass function prediction Eq.(5.1) for F122 model: () with . Right: the fractional error in the semi-analytic predictions for F215 and F122 models compared to simulation results. The plotted points include results from all three redshift values that we have looked at i.e. .
123456Mass ()F677
123456Mass ()M384
Figure 3: The simulation results and semi-analytic mass function prediction for: (i) left: a feeder simulation with (), and (ii) right: a mixed scaling simulation with ().

The right hand panel of Figure 3 is a mixed scaling simulation with , and approximately a third of the contribution to the coming from the feeder component. For the theory curve we used both Eq.(5.1) and Eq.(5.2) for the corresponding components but forced the same . As with the case with other simulations containing feeder scaling, the high mass function results are not described well by the Edgeworth mass function.

All in all, our analysis of the truncation error and the Edgeworth fits to the mass functions from simulations were qualitatively consistent with each other. Hence, we find that the error evaluation of the PDF is a good indicator of the accuracy of the Edgeworth (or Petrov) mass function. However, to fit the simulations well we needed to rescale the analytic ratio of non-Gaussian to Gaussian mass function by an extra dependent parameter, , for both scalings. This rescaling seems reasonable to enforce the same total matter density regardless of level of non-Gaussianity. In Figure 4, we plot our best fit values as a function of for both scalings. We find a simple linear relation between and :

Figure 4: dependence of defined in Eq.(5.1) and Eq.(5.2) for feeder and hierarchical mass functions respectively. The data points are best fit obtained from simulations and to obtain the best fit lines we require for .
Figure 5: We plot results from our simulations for both hierarchical and feeder scalings as density plots. The quantity plotted is . The points overlayed on the plots are the data points used to obtain the density plots. The red solid line shows the maximum trusted when allowing for 20% error in the PDF from to some . The maximum trusted line for the hierarchical case lies outside of the plot range. See Appendix B for details of how these curves were obtained.
Figure 6: We plot the relative difference between the semi-analytic predictions (Eq.(5.2 or Eq.(5.1)) and the simulation results as a density plot. See Eq.(5.5) for the precise definition of the quantity plotted. We have omitted simulation results for which the uncertainty is larger than . The important observation to note here is that the relative difference is quite small i.e. , below the solid red line.

Analytic suggestions for non-Gaussian mass functions have also been obtained using excursion set methods [66, 67, 31, 68, 69, 70, 71]. It would be interesting to compare the predictions for additional corrections to Eq.(3.7) from those methods to these simulations, particularly the results for shown in Figure 4.

To finish up our discussion of mass function results, in Figure 5 we show the non Gaussian mass function results (from simulations) in a two dimensional density plot on a plane, for both hierarchical and feeder scalings. We plot the quantity . In the same plots, we overlay the maximum that we can trust the PDF 3.5, , as a function of ; see error analysis in Appendix B for details. In Figure 6, we have plots similar to those in Figure 5 but now we are plotting the relative difference of our mass function predictions (from 5.2 and 5.1) from the simulation results. The quantity plotted is:


These plots summarize our results for the mass function (for simulations with ). Each isolated region in the density plots represents one simulation, which can be easily mapped to the exact simulation (in Table 1) by looking at the value. Figure 5 simply shows our simulation mass function results. But from Figure 6, we we can see that the magnitude of relative difference of our semi-analytic mass functions to the simulation results is generally less than ten percent () for the values of and that are calculated to be trustworthy by evaluating the series truncation error of the PDF 3.5. This is telling us that our semi-analytical mass function fits are consistent with the truncation error analysis.

Bias and Stochasticity:

Next we present simulation results and fit to the analytical predictions for large scale bias and large scale stochasticity. Before considering the non-Gaussian results, we present the results from our Gaussian simulations in Figure 7. Both the bias and the stochastic bias (after the shot noise correction to ) are constant for , consistent with lowest order predictions for Gaussian fluctuations. Further, the large scale stochasticity parameter is predicted to approach unity for Gaussian initial conditions. Figure 7 (right plot) is consistent with large scale stochasticity being constant for Gaussian initial conditions; however, the constant value is slightly greater than unity at when halos in the mass range are used at . For our other samples, the value of increased when using halos of larger mass. The values of the fitted parameters (bias and stochasticity) for our Gaussian simulations are summarized in Table 2. We also note that the shot noise corrections are large for some of our halo samples. For example, the ratio of shot noise correction to the uncorrected halo power spectrum is largest for the halo sample at (Gaussian simulations) in Table 2; at a reference wavenumber , this ratio is . Typically, this ratio is smaller () for the samples at smaller redshift in Table 2. Also, note that for many of our non-Gaussian models, the fractional shot noise contribution decreases to even for the sample as there are more halos in the non-Gaussian samples at higher redshifts. Some progress has been made recently towards a better understanding and modeling of large scale stochasticity [72], beyond the usual shot noise correction, for Gaussian initial conditions. In this paper, however, we are focusing on the non-Gaussian effect only.

Mass range of halos used