Higher moments of primordial nonGaussianity and Nbody simulations
Abstract
We perform cosmological Nbody simulations with nonGaussian 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 nonGaussianity that need not be weak. This scenario allows us to adjust the relative importance of nonGaussian 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 semianalytic prescriptions for the nonGaussian 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 IGC14/21
Higher moments of primordial nonGaussianity and Nbody 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
Contents
1 Introduction
Recently, the Planck satellite mission reported tight constraints on primordial nonGaussianity 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 nonGaussianity that is about two orders of magnitude below these constraints. NonGaussianity 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, nonGaussian 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 nonGaussianity tells us about the dynamics of those fields. Furthermore, if nonGaussianity 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 nonlinearity 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 nontrivial to extract when nonGaussianity is weak, some constraints with current data are already affected by assumptions about higher order moments. For example, [13] reports constraints on primordial nonGaussianity using a sample of 237 Xray clusters from the ROSAT AllSky 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 Xray clusters with SZ detected clusters that have been separately used to constrain nonGaussianity [15, 16, 17]. In addition, the eROSITA survey [6] and third generation surveys detecting clusters via the SunyaevZel’dovich effect are expected to provide much larger samples of clusters in the next few years. It will be interesting to revisit the nonGaussian 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 nonGaussianity 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 nonGaussian models predicted from inflation.
In this paper, we will use a slight variant of the popular local ansatz of primordial nonGaussianity 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 semianalytic prescription of [12] for the nonGaussian 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 nonGaussianity 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 nonGaussianity 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 nonGaussian statistics thanks to a contribution from a nonlinear, local function of a Gaussian field :
(2.1) 
Here parametrizes the size of the nonlinear term and the level of nonGaussianity. NonGaussianity 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 BunchDavies vacuum [41, 42].
The twopoint 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 wellconstrained 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 :
(2.2) 
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
(2.3) 
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:
(2.4) 
where . In the single source case, this is far too nonGaussian to be consistent with observations. Therefore, we consider the following two source model [46, 34, 11]:
(2.5) 
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
(2.6) 
the ratio of the contribution of the Gaussian part of the field to the total Gaussian power. is defined by
(2.7) 
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 scaleindependent constant.
The power spectrum, bispectrum, and trispectrum in our model are:
(2.8)  
(2.9)  
(2.10)  
where we have defined as
(2.11) 
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
(2.12) 
In Equations (2.8), (2.9) and (2.10) we have included terms that are usually subdominant in the case of single field, weakly nonGaussian local ansatz. In our model, these terms are important when the field is strongly nonGaussian. To discuss the observational constraints on , let’s consider , and . Observational constraint from small nonGaussianity 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 nonGaussian contribution to the power can shift the spectral index slightly, so that when the field is strongly nonGaussian the measured spectral index is close to, but not identical to, the spectral index of the Gaussian components^{1}^{1}1The 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 Nbody 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 nonGaussianity. 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, nonGaussian 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 type^{2}^{2}2One can also argue on statistical grounds that our observable universe is unlikely to have local nonGaussianity of the type written in Eq.(2.5) with large if inflation lasts much longer than 55 efolds [49, 50]).. Here we are using the two field, local model primarily as a phenomenological tool, easy to implement in Nbody simulations, to generate the feedertype 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 nonGaussian cosmology.
3 Abundance and Clustering Statistics
In this section we present the analytic predictions for the effect of locally nonGaussian, 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
(3.1) 
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
(3.2) 
is the smoothing function (here the Fourier transform of the real space tophat). Note that we will generally suppress the dependence in and in this paper, and usually write and only. We can now compute the connected npoint functions of the smoothed density contrast in real space:
(3.3) 
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 nonGaussian mass functions in that we calculate the ratio of the number density in the presence of nonGaussianity, , to the number density for Gaussian initial conditions, , for a particular halo mass at some redshift using the PressSchechter 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:
(3.4) 
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 nonGaussian initial conditions we need a nonGaussian PDF. The Petrov expansion [53] (which generalizes the Edgeworth expansion [54, 55]) expresses a nonGaussian PDF as a series in the cumulants of the distribution. In terms of ’s, this is:
(3.5) 
where , and ’s are the Hermite polynomials defined as . The second sum is over the nonnegative integer members of the set that satisfy
(3.6) 
For each set . This series can be integrated term by term to obtain , and with Eq.(3.4) gives ratio of nonGaussian mass function to Gaussian mass function [13]
(3.7) 
where
(3.8) 
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 )
(3.9)  
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 nonnegative 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 nonGaussian 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 nonGaussian term in the power spectrum. As a result, on scales other than we need to make a distinction between for the nonGaussian cosmology and for the Gaussian cosmology. In that case, there is an extra (mass or scale dependent) factor in the ratio , where
(3.11) 
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 nonGaussianity. The standard PressSchechter constant of proportionality is two, but for the nonGaussian case it is reasonable to fix the constant by requiring . Gaussian and nonGaussian cosmologies with identical will have slightly different normalization factors. This factor shifts further away from 2 as the level of nonGaussianity 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 nonGaussianity 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 nonGaussianity 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 nonGaussianity introduces a new, scaledependent term relating the power spectrum of halos to the power spectrum of the linear dark matter field. On large scales (small wave numbers) the nonGaussian 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 nonGaussianity was first demonstrated in [33]. An analytic derivation capturing the first order effect of nonGaussian 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 Nbody 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 nonGaussianity, 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 twofield model Eq.(2.5) can be calculated using the peakbackground split formalism [19, 62]. Appendix C has the peakbackgroundsplit calculation that results in Eq.(C.10) as the expression of bias. In case of small, local nonGaussianity ( for the case with a second, strongly nonGaussian field), the expression for large scale bias reduces to:
Here we have used subscripts on the bias coefficients to emphasize that the scaleindependent term depends on all sources of the fluctuations (), while the scaledependent term depends only on those sources with a primordial nonGaussian component (). Recall that in the more frequently quoted expression with a single, nonGaussian source, these two bias coefficients are both equal to the Gaussian bias plus scaleindependent corrections proportional to the level of nonGaussianity. 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, scaledependent term is to lowest order the Gaussian bias for the nonGaussian 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 nonGaussian 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 ()^{3}^{3}3First, 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 peakbackgroundsplit method, we can also calculate the halohalo power spectrum in terms of the underlying matter distribution . We can then define the large scale stochasticity, or crosscorrelation coefficient as
(3.13) 
The calculations and the corresponding expressions for and are given in Appendix C. Peakbackground 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 nonGaussianity (i.e. . For small nonGaussianity, 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 nonGaussian 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 GADGET2 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) nonGaussian 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 nonGaussian field was obtained by adding the three components as in Eq. (2.5). We ensured of the generated nonGaussian 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 .
Nbody simulations:
The public version of the GADGET2 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 interparticle 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 (halomatter cross power spectrum and halohalo 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:
(4.1) 
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 Nbody 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 
5 Results and Discussion
Mass Functions:
The hierarchical scaling has been considered a number of times already with the prediction from the Edgeworthseries 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 nonGaussian to Gaussian mass function is:
(5.1)  
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 Nbody simulations.
Similarly, for hierarchical scaling of moments, the expression for the nonGaussian mass function including terms upto is:
(5.2)  
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 nonGaussianity. To compute the cumulants necessary to calculate the dimensionless moments in the above formulae, we have used the MonteCarlo 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 nonGaussianity. However, the normalization requirement (which has not been analytically enforced) suggests that should depend on the level of nonGaussianity. The top panel of Figure 1 has our mass function results for two simulations with hierarchical scalings.
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 bottommost panels illustrate how calibrating on simulations shifts the purely analytic expectations for the nonGaussian mass function. The dominant effect among these three factors is that of ; change in affects the nonGaussian 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 semianalytic 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 semianalytic 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 nonGaussian mass function results from simulations only, we can verify that the F70 model is more nonGaussian 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 nonGaussianity as the H500 model. This verifies that the nonGaussian mass function is sensitive to the total nonGaussianity and the scaling of higher moments in the initial conditions.
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 nonGaussian 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 nonGaussianity. In Figure 4, we plot our best fit values as a function of for both scalings. We find a simple linear relation between and :
(5.3)  
(5.4) 
Analytic suggestions for nonGaussian 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:
(5.5) 
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 semianalytic 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 semianalytical 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 nonGaussian 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 nonGaussian models, the fractional shot noise contribution decreases to even for the sample as there are more halos in the nonGaussian 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 nonGaussian effect only.
Mass range of halos used 