Nbody simulations with generic nonGaussian initial conditions II: halo bias
Abstract
We present Nbody simulations for generic nonGaussian initial conditions with the aim of exploring and modelling the scaledependent halo bias. This effect is evident on very large scales requiring large simulation boxes. In addition, the previously available prescription to implement generic nonGaussian initial conditions has been improved to keep under control higherorder terms which were spoiling the power spectrum on large scales. We pay particular attention to the differences between physical, inflationmotivated primordial bispectra and their factorizable templates, and to the operational definition of the nonGaussian halo bias (which has both a scaledependent and an approximately scaleindependent contributions). We find that analytic predictions for both the nonGaussian halo mass function and halo bias work well once a fudge factor (which was introduced before but still lacks convincing physical explanation) is calibrated on simulations. The halo bias remains therefore an extremely promising tool to probe primordial nonGaussianity and thus to give insights into the physical mechanism that generated the primordial perturbations. The simulation outputs and tables of the analytic predictions will be made publicly available via the nonGaussian comparison project web site http://icc.ub.edu/~liciaverde/NGSCP.html.
a]Christian Wagner a,b]and Licia Verde
Prepared for submission to JCAP
Nbody simulations with generic nonGaussian initial conditions II: halo bias

ICCUBIEEC, University of Barcelona, Barcelona 08028, Spain.

ICREA, Institució Catalana de Recerca i Estudis Avançats.
Keywords: cosmological simulations, large scale structure, initial conditions
Contents
1 Introduction
Inflation, the standard theory for the origin of the primordial density fluctuations, predicts that the primordial density field can be described, to a good approximation, as a Gaussian random field. This nearly Gaussian nature of the primordial fluctuations is well established observationally (for current constraints see, e.g., [1]). However, small deviations from Gaussianity are predicted by inflation. The level and form of the nonGaussianity depend on the inflationary model. The simplest inflationary model, the standard singlefield slowroll inflation, produces nonGaussianities which are unmeasurably small [2, 3] (for a thorough review on inflationary nonGaussianity see [4, 5]).
Other inflationary models [6, 7, 8, 9, 10, 11, 12], which violate one or more assumptions of the standard singlefield slowroll inflation, may potentially cause much larger nonGaussianities. Hence, measurements of deviations from Gaussianity offer a unique window to distinguish between different inflationary models or to rule out classes of models.
Departures from Gaussianity are described at leading order by the bispectrum of the primordial fluctuations in the gravitational potential, . The bispectrum is the Fourier transform of the 3pointcorrelation function and is identical to zero for a Gaussian field. The functional form of the bispectrum specifies the form of the nonGaussianity, while its amplitude, parametrized by the parameter , gives the level of nonGaussianity.
Since the Universe is statistically homogeneous and isotropic, the wave vectors , , and have to form a triangle. The shape of the triangle, for which the bispectrum peaks, is conventionally used to characterize the form of the nonGaussianity. The local type of nonGaussianity, , where is gravitational potential in the Gaussian case, yields a bispectrum that peaks for squeezed triangles (e.g., ). Potentially large nonGaussianities of the local type are generally produced by inflationary models with multiple fields, where nonGaussianity is created by nonlinearities which develop outside the horizon [8].
Other inflationary models which involve higherderivative operators (see e.g. [13, 14, 15]) give rise to a bispectrum that has its maximum for equilateral triangles ().
Models in which the initial state of the inflaton is different from the BunchDavies vacuum [9, 11, 16, 17], can produce nonGaussianity of the enfolded shape, i.e. the bispectrum is dominated by flattened or enfolded triangle configurations ().
Finally, there is the socalled orthogonal shape of nonGaussianity, which cannot be easily described in terms of a triangle configuration [18]. This shape is —with respect to a specific scalar product (see [8] for an explicit definition)— orthogonal to the equilateral and local shape.
The cosmic microwave background (CMB) offers the cleanest probe of nonGaussianity, since the temperature perturbations are still very small and well described by linear theory. If a shape of the primordial bispectrum is assumed^{1}^{1}1For an approach without assuming a shape for the bispectrum, see [19]., constraints on the corresponding value can be derived from the observed fluctuations in the CMB. The current constraints for the most commonly used shapes are , , (95% CL) [1], and (68% CL) [20] . The Planck mission [21] is expected to measure with an accuracy of .
Primordial nonGaussianity leaves also observable signatures in the largescale structure of the Universe (for recent reviews, see [22, 23]). The bispectrum of the density field measured through observable tracers like galaxies, however, is dominated by nonlinear gravitational evolution and biasing, this makes the measurement of the primordial contribution difficult [24, 25, 26] (see, however, [27]).
Another probe of nonGaussianity from largescale structure observations, is the abundance of galaxy clusters or large voids in the Universe. Both of these objects originate from the tails of the nearly Gaussian initial density distribution and are therefore sensitive to primordial nonGaussainity [28, 29]. The signal depends on the skewness of the initial density distribution, which is given by an integral over the bispectrum, and is therefore less sensitive to the shape of the nonGaussianity. This probe measures the primordial nonGaussianity on scales of the order of , which is smaller than the scales accessible by CMB measurements. Intriguingly, several groups [30, 31, 32] derived estimates of from the observations of very massive galaxy clusters, which are as large as (see, however, [33]).
The clustering of halos on large scales offers another possibility to probe nonGaussianity [34, 35]. NonGaussianity couples small and large scales. The density peaks on small scales, which are the seeds for halo formation, are modulated by the large scale modes of the gravitational potential. This induces a scaledependent halo bias on large scales, which can be very different from the scaleindependent halo bias in the Gaussian case. The amplitude of this effect scales approximately linearly in , while the scale dependence is set by the shape of nonGaussianity. For example, the local type nonGaussianity generates a halo bias which scales as , whereas the halo bias of equilateral types of nonGaussianity is expected to have only a very mild scale dependence [36]. This technique was already used on observational data of galaxy surveys and quasar catalogues [37, 38, 39]. The derived constraints, (95% CL) [37], are competitive with the CMB measurements. Fisher matrix forecasts for future largescale structure surveys predict that measurements of the nonGaussian bias will constrain with an uncertainty of [40, 41, 42, 43].
All of the three aforementioned largescale structure probes of nonGaussianity are highly affected by nonlinear gravitational evolution. The standard tool for taking nonlinearities into account are Nbody simulations. However, until recently, only the local type of nonGaussianity was simulated with Nbody simulations and reasonable agreement with analytic predictions was found [34, 44, 45, 46, 47].
In [48], we presented a prescription for setting up initial conditions for Nbody simulations with nonlocal nonGaussianities. Using this method we were able to simulate the nonlinear power spectrum and the cluster mass function for nonlocal shapes, and again found consistency with theoretical expectations. The technique of [48] was, however, not suitable to study the effect of nonGaussianity on the halo bias, since —depending on the shape of the bispectrum— it introduced artificial power on large scales, which interferes with the expected signal of the nonGaussian halo bias.
In this paper, we modify our ansatz used for the generation of nonGaussian initial conditions in a way that these spurious contributions to the initial power spectrum are always kept under control. This enables us to set up and run nonGaussian simulations with large volumes, which we then use to study —for the first time by means of Nbody simulations— the nonGaussian halo bias for different types of nonGaussianity. As the simulated volume is larger than the one we were able to simulate in [48] and hence the statistical error on the number of massive halos is smaller, we use these simulations also to revisit the nonGaussian cluster mass function.
The outline of the paper is the following. First, we review the nonGaussian halo bias and its theoretical modelling in Sec. 2. In Sec. 3, we compare, in the context of the nonGaussian halo bias, the bispectra of inflationary models with the commonly used separable templates which approximate these bispectra. In Sec. 4, we review the method of setting up nonGaussian initial conditions for a given bispectrum and then present the modification of our ansatz. The practical implementation and numerical settings of our simulation suite are given in Sec. 5. We present and discuss our results on the nonGaussian mass function and halo bias in Sec. 6. Finally, we conclude in Sec. 7.
2 NonGaussian halo bias
The halo bias describes the relation between the halo density contrast, , and the underlying dark matter density contrast, . In Fourier space we can write the relation as , where the random variable models the stochasticity coming from shot noise and possibly from physical halo formation processes. For Gaussian initial conditions the bias is scaleindependent on large scales, for [49, 50, 51]. Analytic predictions and results of Nbody simulations show that the linear^{2}^{2}2The linear bias, , is called linear as it corresponds to the linear coefficient in the Taylor expansion of the halo density field: , where both density fields are smoothed on a sufficiently large scale. Gaussian bias, , increases with redshift and halo mass . In addition the halo bias may depend on environment and mass accretion history [52, 53].
Primordial nonGaussianity of the local type induces a scale dependence of the halo bias on large scales, which scales as . This particular scaling was first found in Nbody simulations by [34], who also gave a simple analytic understanding how this effect arises. Later [35, 37] rederived the analytic prediction in a more rigorous way using two different approaches, the high peak limit and the peakbackground split, respectively. In Ref. [54], the same scale dependence of the nonGaussian halo bias was found by computing the ellipsoidal collapse threshold in the presence of small nonGaussianities. The analytic predictions were tested with further and larger Nbody simulations [45, 44, 46, 47]. Reasonable agreement between theory and simulations was found if small corrections, which we discuss in detail below, are taken into account. In a number of subsequent papers further theoretical modelling was developed by applying (univariate and multivariate) local biasing in the framework of perturbation theory [55, 56, 57, 58, 59, 27]. In summary the nonGaussian bias can be modelled to leading order as
(2.1) 
where is the difference in the linear bias, which is nonzero, since a nonvanishing changes the halo number density and hence the linear halo bias for a fixed halo mass . Mainly, this effect leads to a scaleindependent shift relative to the Gaussian () linear bias [44, 59, 60]. Additionally, it makes the scaledependent second term in above equation slightly nonlinear in [59, 61]. In practice, when fitting observational data, the linear bias is estimated from the data itself.
The amplitude of the nonGaussian contribution is given by , the linear Lagrangian bias , and the redshiftdependent collapse threshold , where is the spherical collapse threshold in an Einsteinde Sitter universe, is the growth function normalized to at high redshift, and can either be regarded as a fudge factor introduced to improve the fitting to Nbody simulations [46] or as a correction of the spherical collapse threshold [54, 62]. The scale dependence of the nonGaussian bias is hidden in and . connects the density with the primordial potential. The linear matter density contrast smoothed on a mass scale is related to the primordial potential by the Poisson equation,
(2.2) 
Here, and are the present day matter fraction and the Hubble constant, respectively. The transfer function models the physics until the end of recombination and is normalised to unity on large scales. The function is the Fourier transform of the tophat filter with a smoothing scale given by the mass . As , we see that scales as on large scales.
The “form factor”, , depends on the shape of the nonGaussianity through the bispectrum of the primordial potential , [35]:
(2.3) 
where is defined by , is the primordial power spectrum, and denotes the variance of linear density fluctuations on the mass scale at redshiht . Considering the limit of , we see that on large scales the form factor depends primarily on the squeezed limit of the bispectrum. For the local type of nonGaussianity, the “form factor” becomes scale and mass independent on large scales, . Hence, we obtain the aforementioned scale dependence of the halo bias in the case of local nonGaussianity. As we discuss below in Sec. 3, other shapes of nonGaussianity lead to different and in general weaker scale dependences.
One higherorder correction to , which was introduced in [44], is a relatively small contribution coming from the fractional difference in the Gaussian and nonGaussian nonlinear matter power spectrum. This effect becomes larger at smaller scales but remains below a few percent on all scales for observationally allowed values of .
Some additional effects are also missed in Eq. (2.1): The impact of the mass accretion history of the halos —in particular recent major mergers— on the nonGaussian halo bias was pointed out in [37], and further investigated and tested with Nbody simulations in [63]. The formation redshift of the halos, , changes the amplitude of the nonGaussian correction by , where varies roughly between and from very young to very old halos [63]. For recent major mergers, takes the value [37, 63].
In addition, higherorder effects in the biasing description (e.g., terms of the order [55, 64, 59]) and in the primordial nonGaussianity (e.g., nonGaussianity at the cubic order quantified by [64, 59]) are neglected in Eq. (2.1).
Finally, general relativistic effects are not taken into account, which might become important on very large scales [65, 66, 67, 36, 68, 69].
So far we discussed the scaledependent halo bias caused by the local type of nonGaussianity with a scaleindependent . The effect of a scaledependent were investigated analytically and by means of Nbody simulations in [70, 71]. In this paper we are mainly interested in the scaledependent halo bias of nonlocal types of nonGaussianity. The prescription presented above is sufficiently general to model the effects of other shapes of nonGaussianity. As discussed above, the shape dependence is hidden in the form factor . The expression for the form factor, see Eq. (2.3), was derived by [35] in the framework of local biasing of density peaks and was evaluated for different types on nonGaussianity in [36]. Recently, [72] generalized the peakbackground split approach to nonlocal types of nonGaussianities and derived a very similar but not identical expression for the form factor. In the next section we discuss the analytic predictions of the nonGaussian halo bias for nonlocal shapes using the high peak formulation of [35].
3 Physical shapes vs. templates
The functional form of the primordial bispectrum, , predicted by inflationary models can be quite complex. In particular, the bispectrum is usually not separable, i.e. it cannot be written as a finite sum of terms which are factorizable as a product of functions of , and . Separability is important for computationally efficient analysis and simulation of both CMB maps [13, 19] and largescale structure [48, 73]. Hence, most physical shapes have been approximated by socalled templates, which are separable and constructed to effectively maximize the correlation between the physical shape and the template across all triangle configurations. In practice a socalled “cosine” shape correlator is used [8], which absolute value is always with 1 indicating perfect correlation. This measure of similarity is useful for constraints on the bispectrum coming from CMB or largescale structure analysis, as it probes the entire range of possible triangle configurations. However, these templates can be misleading for predictions of the nonGaussian halo bias, as in this case the shape is probed mostly in the squeezed limit. In this limit the templates do not necessarily need to be good approximations of the physical shapes to yield a cosine close to unity [20]. For the nonGaussian halo bias, however, the correct scaling of the templates in this regime is crucial for sensible predictions. Let us consider several templates in this regard. We start with the equilateral^{3}^{3}3This template peaks when the wavenumbers form a equilateral triangle. template, which was introduced in [74] as a separable approximation to the shape functions of DBI inflation [14], inflation with higherderivative terms [13], and ghost inflation [15]:
(3.1) 
with^{4}^{4}4The original expression of the equilateral template was proposed for a scaleinvariant power spectrum and given in terms of the wavenumbers . Here we generalize it to non scaleinvariant power spectra by substituting with with proper normalization. We use the same substitution ansatz in the case of other bispectra.
(3.2)  
(3.3)  
(3.4) 
and the dependence of on the three wavenumbers is understood. The primordial power spectrum is specified by the amplitude at and the spectral index :
(3.5) 
The equilateral template was also adopted in other models. More generally, the form of the bispectrum is interesting because of its direct connection to the (self) interactions of the inflaton. In the effective field theory of inflation [75], there is a direct connection between the coefficients of operators of the Lagrangian for the inflationary fluctuations and the shapes of primordial nonGaussianity. In [18], it was shown that the bispectrum of singlefield inflation with an approximate shift symmetry protecting the Goldstone boson is a linear combination of two different shapes. The two shapes, and , which correspond to the interaction operators and (see [18] for details), are both close to the equilateral template. Especially, the cosine of with the equilateral template is very close to one. For this reason, we denote this shape also by “SSZ equilateral”.^{5}^{5}5Previously, in the framework of general single field inflation, [9] obtained two independent shapes which span the same space as the two shapes of [18].
We now turn to the question how well the equilateral template captures the behaviour of these physical shapes in the squeezed limit. Choosing and using the following notation, , , and , we obtain for the equilateral template
(3.6)  
(3.7)  
(3.8) 
where in the last line a scaleinvariant power spectrum was assumed, i.e. . We find a similar scaling of the bispectrum in the squeezed limit also for the physical shapes generated by the aforementioned inflationary models. However, the exact scaling and its amplitude differ slightly from case to case. In the case of a scaleinvariant power spectrum, the scaling reduces for all models to with a modeldependent amplitude. For example for the bispectrum of DBI inflation (see [8] for an explicit expression) we get
(3.9) 
while the scaling of the bispectrum corresponding to [18] is
(3.10) 
Here, we used the equilateral normalization convention for all shapes, i.e. .
In Fig. 1, the scaledependent part of the nonGaussian halo bias, , is shown for different bispectra of the equilateral type. The dashed lines depict for a scaleinvariant power spectrum. Note that on very large scales the amplitude of the effect for the equilateral template and the physical models, namely DBI and SSZ eql., is different by a factor of and , respectively, in agreement with the above expressions for the squeezed limit. The solid lines show computed with a power spectrum with . The upturn and the wiggles on smaller scales stem from the transfer function included in . Note, in contrast to the local type of nonGaussianity, for the equilateral shapes is massdependent even on large scale.
In summary, predicting the nonGaussian halo bias with the equilateral template instead of the physical shapes results in approximately the correct scaling of the effect but a slightly smaller amplitude. The amplitude derived from the physical shapes of DBI and SSZ eql. is larger by and , respectively.
For other templates, as we discuss now, the situation becomes more problematic. The orthogonal^{6}^{6}6The orthogonal template was introduced as a separable approximation to the linear combination of the shapes and that is orthogonal (i.e. has a vanishing cosine) to template [18]
(3.11) 
and the enfolded^{7}^{7}7This template approximates the shape function produced by inflationary models with a modifications of the initialstate of the inflaton field and has its maximum for flattened triangles, i.e. . template [16]
(3.12) 
scale both as in the squeezed limit, while the corresponding physical shapes scale as and , respectively. In Fig. 2 we show again the scaledependent part of the nonGaussian bias, , computed from these commonly used templates and the corresponding physical shapes. For comparison, the scale dependence obtained from the local type of nonGaussianity, , and from the equilateral template are also shown. Here, we used in all cases a scaleinvariant power spectrum. The solid and dashed lines correspond to halos of mass and , respectively. The scale dependence of the halo bias caused by modifiedinitialstate type of nonGaussianity [16] is almost degenerate with the one caused by the local type of nonGaussianity, i.e. both predictions are very similar when is rescaled by . However, on scales of the different halo mass dependence lifts this degeneracy to some extent. The orthogonal shape of [18] (SSZ orth.) —having the same scaling in the squeezed limit as the equilateral shapes— does not lead to a scaledependent halo bias on large scales and is almost indistinguishable from the equilateral type. Only on smaller scales, the “orthogonal” halo bias has a slightly different scale dependence than the equilateral one.
In conclusion, in respect of the halo bias, the commonly used enfolded and orthogonal templates are not adequate to capture the behaviour of the underlying physical models. The physical shapes scale as or as in the squeezed limit, while the scaling of the enfolded and orthogonal templates is .
In [19] a mode expansion of the bispectrum in separable functions was introduced. This expansion can be applied to arbitrary nonseparable bispectra. In many cases already a small number of modes is sufficient to yield a large cosine with the original bispectrum. While this method of constructing separable templates for arbitrary bispectra is useful for efficient measurements of the bispectrum from the CMB or the largescale structure, it fails to reproduce the correct scaling of the bispectrum in the squeezed limit. By construction of the mode expansion, all templates derived in this way have a scaling of in the squeezed limit, irrespectively of the scaling of the original bispectrum.
New separable templates which are not only overall good approximations to the original bispectrum but also faithfully capture the scaling in the squeezed limit can most probably be constructed (e.g [20]). However, we will see in the next section that separability of the bispectrum does not speed up our current implementation of generating the initial conditions for Nbody simulation.
As a concluding remark, we note that similar considerations apply for the skewness specified by the bispectrum. The skewness is the relevant parameter used for the predictions of the nonGaussian halo mass function. Although having a large mean correlation over all triangle configurations, templates and corresponding physical shapes may lead to different values of the skewness. However, the differences in the skewness will be, in general, much smaller than the aforementioned differences in the scaling behaviour of the squeezed limit.
4 Initial conditions
In [48] we introduced a prescription how to generate nonGaussian initial conditions for a given bispectrum, i.e. how we find a realization of the potential field , which has the desired bispectrum
(4.1) 
We shortly review our method. First, we decompose the potential field in Fourier space, , in a Gaussian and a nonGaussian part
(4.2) 
The Gaussian field is statistically completely described by the power spectrum
(4.3) 
As the nonGaussian contribution is small, the leading order of the bispectrum is given by
(4.4) 
Using the following ansatz for
(4.5)  
with the complex conjugate of , one can easily show that
(4.6) 
Finally, we recover at leading order Eq. (4.1) by virtue of the symmetry property of the bispectrum.
One problem of the above ansatz is that the nonGaussian contributions to the power spectrum can become very large on large scales. ^{8}^{8}8On large scales, we require the nonGaussian contributions to the power spectrum to be negligible (or renormalizable), in order to be in agreement with CMB measurements. Naively, one expects to be much smaller than the Gaussian contribution . However in the case of certain bispectra, diverges more strongly for small wavenumbers than the Gaussian power spectrum. With the above ansatz for the nonGaussian contribution to the power spectrum reads
(4.7) 
In practice, there is a minimum wavenumber , the largest mode which fits in the simulation box, and a maximum wavenumber , which corresponds to the smallest wavelength still resolved by the used grid size. These wavenumbers provide large and smallscale cutoffs for the above integral. Nevertheless, this spurious contribution to the power spectrum can be large enough to affect substantially the results of the Nbody simulations. In particular, simulations with large box sizes, which are needed to predict the halo bias on large scales, suffer from these divergences. On large scales, , the integral is dominated by the bispectrum in the squeezed limit. Using the scalings presented in Sec. 3, we see that the nonGaussian contribution scales as , , and for the equilateral, orthogonal/enfolded templates, and the local type of nonGaussianity, respectively. Hence, in the case of the equilateral type the contribution is small, whereas in the other cases it starts to dominate the power spectrum on large scales.^{9}^{9}9In [73] the mode expansion technique of [19] was used to develop a computationally efficient method to generate initial conditions for Nbody simulations for a given bispectrum. However, due to the scaling in the squeezed limit of the expanded bispectrum, also in this case, nonGaussian contributions to the power spectrum of the form arise on large scales.
For the local type of nonGaussianity, we found in [48] that our ansatz becomes identical to the real space formulation , which does not lead to the aforementioned divergences, when we use in Eq. (4.5) instead of . Here we generalize this modification by writing it in the following way
(4.8) 
Inserting this expression in Eq. (4.5) yields the modified ansatz for (see also [72])
(4.9) 
Note that the integrand in the modified ansatz is in general not separable. This means that the integration cannot be written in a simple way as a sum of convolutions, which can be efficiently computed by Fast Fourier Transformations (FFTs). However, a possible avenue to still factorize our modified ansatz is the use of Schwinger parameters (see [76] for details). This is in contrast to our original ansatz, which is separable as long as the bispectrum is separable. However, the modified ansatz does not lead to contributions to the power spectrum, which diverge for , since they are suppressed by :^{10}^{10}10Nevertheless, without cutoffs the integral can still be divergent, e.g. in the local case this would lead to a renormalization of the amplitude. In practice, the cutoffs given by the box size and grid size render this contribution negligible.
(4.10) 
Note that neither Eq. (4.5) nor Eq. (4.9) specify the trispectrum of the nonGaussian field at the leading order. A possible inclusion of the trispectrum in the generation of the initial conditions was proposed in [73]. In this paper, however, we neglect contributions from higherorder correlations.
Next, we can convert the potential into the linear density field with the help of the transfer function, which we compute with CAMB [77], and the Poisson equation:
(4.11) 
Finally, the density is sampled with particles. The positions and velocities of the particles are derived from the displacement field at the initial redshift using the Zel’dovich Approximation or secondorder Lagrangian perturbation theory [78, 79, 80].
5 Simulations
In this paper we are mainly interested in the scaledependent halo bias of different types of nonGaussianity. As this effect, in general, is larger on large scales, we need to simulate large volumes but still have sufficient mass and force resolution to resolve halos. Since we are not interested in halos of a particular mass, a reasonable tradeoff between size and resolution is a box size of about and a particle load of one billion particles. This enables us to probe large scales () and resolve halos with masses above ( particles).
The downside of our modified ansatz for the initial nonGaussian part of the potential, given in Eq. (4.9), is the nonseparability of the integrand. This precludes us from using FFTs to compute the integral. Instead, we perform the integration in Fourier space by direct summation, i.e. for each mode the following sum has to be computed
(5.1) 
where runs over all grid points in Fourier space. If denotes the number of grid points in one dimension, then the computational cost of computing for all modes scales as . This scaling renders the usage of large grid sizes infeasible. Even for a moderately small grid size of , the generation of the nonGaussian initial conditions takes two days on 256 cores of presentday CPUs. In order to still use a large particle load, we use two different grid sizes: One grid being as large as the particle grid, , which samples only the Gaussian contribution to the potential, and one smaller grid for the nonGaussian part.^{11}^{11}11To compute the nonGaussian part, we use Eq. (5.1) where we only sum over modes which lie within the smaller grid. By doing this, we neglect nonGaussianities on scales smaller than the grid spacing of the nonGaussian grid. This not only changes the clustering on these scales, but may also affect the halo clustering on large scales, as the formation of the halos occurs on small scales. This effect will depend on the halo mass; halos which form on scales larger than the resolution of the nonGaussian grid will, presumably, hardly be affected. The choice of the nonGaussian grid size is therefore a tradeoff between computational resources and the lowest halo mass, which is still unaffected by the spatially unresolved nonGaussianity.
To quantify this effect in more detail, we run two simulations of the local type of nonGaussianity with . The settings of the simulations are identical except for the way we set up the initial conditions. In one case, we use a Gaussian grid of size 1024 and a nonGaussian grid of size 400 (local1024400), in the other, both grids are of size 1024 (local10241024). This large nonGaussian grid size of 1024 is possible because in the case of local nonGaussianity the computation can be performed efficiently in real space, .
We adopt a flat CDM cosmology with cosmological parameters consistent with current observational constraints [1], namely (the presentday matter fraction), (presentday baryon fraction), (the Hubble constant), (the spectral index of the primordial power spectrum), and (the rms of the linear density fluctuations at in spheres of ).
The simulations are carried out with Gadget2 [81] taking only the gravitational interaction into account. The box size is , and particles with a softening length of are used to sample the density field. The simulations are started at an initial redshift of using secondorder Lagrangian perturbation theory to displace the particles from a regular grid.
Halos are found in the particle outputs of the Nbody simulations by the publicly available halo finder AHF [82]. This halo finder defines halos as gravitationally bound objects with a spherical overdensity equal to the redshiftdependent virial overdensity. With the halo catalogues and particle snapshots at hand, we can compute the Fourier modes and of the halo and matter density field, respectively. This is achieved by first assigning the halos or particles to a regular grid () using the cloudincell (CIC) scheme and then performing an FFT.
In Fig. 3, we show the halo bias obtained from the two simulations local10241024 (solid lines) and local1024400 (dashed lines) for different halo masses at redshift and . We compute the halo bias by the ratio of the halomatter cross power spectrum to the matter power spectrum
(5.2) 
with and , where the average is taken over all modes which lie in a shell of radius and thickness . As expected, the clustering of halos with a mass below is affected by the smaller grid size used for the nonGaussian . The grid spacing of the nonGaussian field is . This corresponds approximately to the Lagrangian radius of halos of mass . The small difference in the halo clustering for more massive halos is consistent with shot noise. At higher redshift, the noise increases due the smaller number of halos per mass bin.
The other quantity of interest, which is also affected by the smaller nonGaussian grid, is the halo number density. In Fig. 4, we show the ratio of the mass functions derived from the two simulations local1024400 and local10241024. We find that the smaller nonGaussian grid size decreases the halo number density at the lowmass end by a few percent. For halos more massive than , the halo number densities of both simulations agree within the errors.
In App. A, we compare the results of two Nbody simulations of the equilateral type, for which we used the two different formulations for the nonGaussian part of the potential, Eq. (4.5) and Eq. (4.9). This comparison gives us another estimate of the minimum halo mass, , above which the results are not affected by the low spatial resolution of the nonGaussianity.
We conclude that, for halos more massive than , even the simulations which are set up with a small nonGaussian grid () model correctly the characteristics of nonGaussianity, in particular the halo bias and the halo number density. As a rule of thumb, the nonGaussian grid spacing should be smaller than the Lagrangian radius of the lowestmass halo of interest.
type of nonGaussianity # realizations local 250 400 1024 1 local 250 1024 1024 2 local 60 1024 1024 3 equilateral template 1000 400 1024 2 equilateral template^{a} 1000 1024 1024 1 orthogonal template 1000 400 1024 2 orthogonal template 250 400 1024 3 Gaussian   1024 3
\@textsuperscripta  In this case, we use ansatz Eq. (4.5) for the nonGaussian . The convolution is computed with an FFT in real space. 
For the main study of the nonGaussian halo bias, we perform a suite of simulations of different shapes of nonGaussianity and different values, which are summarized in Tab. 1. The cosmology and simulation settings, like box size, particle load, etc., are the same as given above.
Although the orthogonal template is not a good approximation of the physical shape (SSZ orth.) in the context of halo bias (see Sec. 3) and although our ansatz is applicable to nonseparable bispectra, we still use this template for some of our simulations. This has the following reasons. First, the corresponding physical shape (SSZ orth.) causes a halo bias which is very similar to the equilateral type and, hence, hard to measure. Second, the scale dependence of the nonGaussian halo bias of the orthogonal template lies in between the local and the equilateral one, which makes it an interesting toy model to test the analytic predictions of the nonGaussian halo bias against Nbody simulations. Lastly, other groups (e.g. [83]) are running Nbody simulations of this type, too, which makes a future comparison easier.
The chosen values for are partly larger than allowed by current constraints coming from CMB data [1]. Here, however, we want primarily to compare theoretical predictions with the results of Nbody simulations. To do this accurately, we need high signaltonoise and therefore a relatively large value of . This is in particular the case for the equilateral type of nonGaussianity, for which the nonGaussian halo bias has only a very weak scale dependence. Note that if we find that the analytic modelling of the nonGaussian bias works for relatively large nonGaussianity, it will also be valid for small values of .
6 Results
In this section we present the outcome of our Nbody simulations of different types on nonGaussianity. Here, we used the standard templates (namely the local, equilateral and orthogonal templates) to create the initial conditions for the suite of Nbody simulations. The templates are useful as toy models to test and calibrate the analytical predictions. If the analytical expressions can correctly reproduce the bias behaviour for shapes that have such different dependence in the squeezed limit, this lends support to their validity.
Although the simulations were targeted to study the nonGaussian halo bias, we first consider the effect of nonGaussianity on the halo mass function. Afterwards, we analyse in detail the scaledependent halo bias induced by the different types of nonGaussianity.
6.1 Halo mass function
Primordial nonGaussianity affects significantly the number density of objects which correspond to the tails of the probability distribution function (PDF) of the initial density field [28] (for recent reviews, see [22, 23]). Galaxy clustersized halos originate from rare high peaks in the density field smoothed on the scale of the cluster mass. If the PDF has a positive skewness on this mass scale, more of these high peaks are initially present and give rise to more galaxy clusters.
The ratio of the nonGaussian to the Gaussian halo mass function, , can be modelled analytically in the framework of the PressSchechter formalism [28, 29] and the excursion set approach [84, 85]. For the comparison with our Nbody results, we use the formulas of [28] (MVJ) and [29] (LV), who used the saddlepoint technique and the Edgeworth expansion, respectively, to derive the following expressions
(6.1) 
and
(6.2) 
where is the linear mass variance on a mass scale at . The redshift dependence is modelled by the collapse threshold, , which includes a fudge factor, ^{12}^{12}12This fudge factor was introduced in [46] and interpreted as a modification to the collapse threshold which also applies to the Gaussian mass function. Here, however, we do not adopt this interpretation and consider only as a fudge factor in the ratio of the nonGaussian to the Gaussian mass function, which we calibrate using Nbody simulations.. The nonGaussianity is described by the normalized skewness of the linear density field, , on a mass scale at , . The normalized skewness scales linearly in . The magnitude and sign of depend on the shape of nonGaussianity, while the mass dependence is only weakly dependent on the type of nonGaussianity.
In Fig. 5, we show the fractional difference in the nonGaussian and Gaussian mass functions derived from the simulations of the local type of nonGaussianity. The four panels correspond to different redshifts, (top left), (top right), (bottom left), and (bottom right). The filled and open symbols represent the results for and , respectively. The corresponding analytic predictions of MVJ and LV (with the bestfit value for ) are depicted by the solid and dashed lines, respectively. Overall, the agreement between simulations and predictions is very good. Only at the lowmass end at , the data points do not follow the theory lines very closely. This is not of practical importance, as the difference is much smaller than observationally uncertainties in halo mass measurements. The applied fudge factor of for MVJ and for LV is in agreement with our findings in [48]. Other groups [46, 45], who worked with friendsoffriends (FOF) halos, used for both predictions, while [44], who applied a sphericaloverdensity (SO) halo finder as we do, found for LV and for MVJ (V. Desjacques, private communication). For the standard definitions of FOF and SO halos, i.e. using a linking length of times the mean particle distance (FOF) and the redshiftdependend virial overdensity (SO), [23] found that for a given halo mass the effect of nonGaussianity on the number density of FOF halos is smaller than on the number density of SO halos. Note however that it remains to be seen whether FOF or SO halos most closely match observationallyselected halos (e.g., SZselected or Xray selected clusters). Nevertheless, the fact that is quite close to unity in both cases, suggests that this will not introduce a dominant systematic effect.
The mass function results from our nonGaussian Nbody simulations of the orthogonal type are presented in Fig. 6. Here, the values are (filled symbols) and (open symbols). Note that these simulations used a rather small grid size () for the nonGaussian initial density field. Hence, halos with masses below are expected to show a smaller nonGaussian effect than predicted by the theory. This is indeed the case, as can be seen in the figure. More interestingly, even at the highmass end the nonGaussian effect is smaller than one could expect. To bring the theory in agreement with the Nbody data, the fudge factor has to be quite small, decreasing from to with decreasing redshift.
Lastly, we consider the increase in the halo mass function due to nonGaussianity of the equilateral type. In Fig. 7, the results of three simulations with are shown. Two of them, “eql. a)”, started from initial conditions computed with our modified ansatz, Eq. (4.9), using a nonGaussian grid size of . The third simulation, “eql. b)” (black triangles), used initial conditions generated with our original method, Eq. (4.5), with (see also App. A). Remember that in the case of the equilateral type, artificial contributions to the power spectrum are kept under control also when our original ansatz is used. As expected, the Nbody results agree with each other for halo masses above . The falloff of “eql. a)” at the lowmass end is caused by the small nonGaussian grid. Interestingly, we find that, in the case of the equilateral shape of nonGaussianity, is very close to one and does not change significantly with redshift.
In conclusion, on the mass scales probed, the nonGaussian effect on the halo mass function can reasonably well modelled by the analytic predictions, if a shapedependent fudge factor is taken into account. In the case of the orthogonal type, there is a hint that is in addition redshift dependent.
6.2 Halo bias
In order to explicitly explore the dependence on the halo mass, we divide the halos into five mass bins:
The number of halos per mass bin decreases with mass. In particular, at high redshift and for large masses the exponential falloff in the mass function reduces the number of halos per bin quickly. In the following, we only consider mass bins with more than halos.
In order to keep the statistical error as small as possible in our analysis, we compute and fit the effect of the nonGaussian halo bias in the following way. First, we compute the difference between the nonGaussian and Gaussian bias from the corresponding Nbody simulations:
(6.3) 
As the nonGaussian and the Gaussian simulation share the same realization of the initial Gaussian field, is almost free of sample variance and, in addition, has smaller shot noise than and individually. We estimate the error on directly from the distribution of in each bin (for details, see App. B).
Following Eq. 2.1, we model by
(6.4) 
where is the linear halo bias obtained from the Gaussian simulation on large scales (see App. B). The scaleindependent shift, , can be predicted by the difference in the linear bias derived from the nonGaussian and Gaussian mass functions using the peakbackground split approach [44]
(6.5) 
with being the ratio of the nonGaussian to the Gaussian mass function. Here, however, we treat as a free parameter and compare it later on with the prediction derived from the mass functions. We choose as the second free parameter. All other quantities in Eq. (6.4) are computed from the theory and are kept fixed.
In Fig. 8, we show as an example the effect of local nonGaussianity on the halo bias for halos of mass at . Note that we plot . As is negative, this addition of is needed to still make use of the logarithmic scale.
The different line types visualize the effect of the different terms in Eq. (6.4). The solid lines show the best fit to the data (using all modes up to ) and include all terms given above. The shortdashed lines neglect appearing inside the square brackets in Eq. (6.4). The inclusion of this term makes the nonGaussian bias nonlinear in [59, 61], since depends on . The dotdashed lines neglect completely. This scaleindependent bias shift becomes important on smaller scales (), for which the scaledependent part becomes small.
An example of the measured nonGaussian bias from the simulations of the orthogonal type is given in Fig. 9. Here, the halos have again a mass of , but were found in simulation outputs at . Consequently, the number of halos is smaller than in the previous figure and the residual shot noise is larger. The line types have the same meaning as before. On large scales, the halo bias scales as as predicted by the theory. Hence, with increasing wavenumber, the effect does not drop as rapidly as in the local case and extends to smaller scales.
Next, we present an example of the nonGaussian halo bias induced by the equilateral shape. In Fig. 10, the halo bias corresponding to two different mass bins at is shown. As expected, the scale dependence is very weak and in agreement with the theoretical predictions. In particular, the observed mass dependence of the effect is consistent with the model predictions.
After having discussed for each type of nonGaussianity typical examples, we show the complete set of bestfit values of the fitting parameters in Fig. 11. On the lefthand side, the bestfit values of the fudge factor are presented. Different colours correspond to different redshifts: (red), (green), (blue), and (magenta). Triangles, boxes, and circles depict the three different realizations of the initial Gaussian random field used for the generation of the initial conditions. In the case of the local and orthogonal type, the open symbols correspond to and , respectively. Filled symbols show the results for (local), (orth.), and (eql.). For clarity, the points of each mass bin are spread over the range of each mass bin.
For the local type of nonGaussianity, we recover within the error bars the same value, which was needed to bring the mass functions in agreement with analytic predictions. For the orthogonal shape, the effect is suppressed by a factor which is similar to the one we found for the nonGaussian mass function, . This suppression of the nonGaussian halo bias effect is clearly redshift and halo mass dependent. In the case of the equilateral type of nonGaussianity, the error bars are as expected very large. Nevertheless, also in this case there is a hint of a redshift and halo mass dependence of the effect. Note, however, that this suppression visible for the nonGaussian halo bias effect for the equilateral shape was not found for the effect on the mass function, for which the measured fudge factor was . These findings are very interesting and —if solidified by larger simulations— may help to lead to a better theoretical understanding of the halo biasing (see discussion in [72]).
On the righthand side of Fig. 11, the bestfit values of normalized by are shown. The colour and symbol coding is the same as before. As open and filled symbols (corresponding to different values) are consistent with each other, we can infer that the scaleindependent shift is linear in for the values probed. The solid lines represent the predictions from Eq. (6.5) using the LV mass function ratio, , and taking the measured fudge factor from the mass function into account. Keeping in mind that the LV mass function is not in all cases a good fit to the mass functions derived from our simulations (see for example bottom right panel of Fig. 6), the agreement between the bestfit values and the predicted bias shift is very reasonable.
7 Conclusions
Nbody simulations are an indispensable tool to test, develop and calibrate any statistical analysis of largescale structure. In this paper, we have further explored the issue of setting up generic nonGaussian initial conditions for Nbody simulations. In a previous paper [48] we began addressing this issue, focussing —as we do here— on nonGaussianities specified by a nonzero bispectrum. In [48] we focussed on the nonGaussian corrections to the halo mass function and to the nonlinear matter power spectrum, which, for the scales and redshfits of interest, could be reliably computed from simulation boxes of on a side. Here we concentrate on the largescale nonGaussian halo bias. The scaledependent bias is evident on very large scales thus needs multiple simulations of much larger boxes to be measured reliably for values not too large as to push beyond the regime of validity of the analytical description of the effect.
The scaledependent nonGaussian halo bias has been recognized to be a very competitive approach to constrain primordial nonGaussianity (e.g., [37, 40, 41, 59]) yielding forecasted error bars on the local nonGaussian parameter of which makes this approach formally better than an ideal CMB experiment [86]. So far this effect has been explored with Nbody simulations only for the socalled local type of nonGaussianity. However, the shape of nonGaussianity is a window into the physical mechanism of inflation and the generation of primordial perturbations, and it is of value to be able to simulate initial conditions for general bispectrum shapes.
The technique presented in [48] was not suitable to study the halo bias because, in general, higherorder contributions would leave an artificial signature on the largescale power spectrum (which are the scales where the scale dependence of the halo bias is evident). These components could be kept under control only in specific cases. Here the issue is solved and these components are always under control. This however comes at the expense of computational cost: the new expressions are not separable (even if the bispectrum itself is) meaning that computational shortcuts that could be implemented for separable bispectra before, now cannot be applied in the same way.
Nevertheless we find a speedup solution which involves sampling only the Gaussian contribution to the potential with a full resolution grid and using a smaller grid for the nonGaussian part. We assess the performance of this approach for the local case, where the computation can be performed in real space and is not computationally intensive. By using a smaller grid size for the nonGaussian part of the initial conditions, we can keep the computational time taken by the calculation of the initial condition comparable to the time needed to run the simulation.
The form of the bispectrum predicted by inflationary models can be complicated and in general is notfactorizable. As factorization is very important for an efficient analysis of CMB data, physical bispectrum shapes have been approximated by factorizable templates. We highlight here that, since the halo bias is dominated by squeezed configurations, what is a good template for CMB analysis may not be correct for the halo bias. We show that the equilateral template results in roughly the correct scale dependence of the halo bias for the corresponding physical models, but with a shift in amplitude; if unaccounted for, it could introduce a bias on the estimated of up to 50%. Fortunately, this shift can be quantified well. Enfolded and orthogonal templates on the other hand do not work as well. They both lead to a nonGuassian bias that scales as on large scales. This is different from the halo bias induced by the corresponding physical models: The scale dependence of the halo bias caused by modifiedinitialstate type of nonGaussianity is almost degenerate with the one caused by the local type of nonGaussianity. The halo bias effect of the orthogonal physical shape is very similar to the effect of the equilateral shape.
We performed several simulations with nonGaussian initial conditions of different shapes to further test and calibrate the analytic predictions. For this purpose we used the standard templates: these are useful toy models as they span a range of scalings in the squeezed regime. If the analytical expressions can correctly reproduce the bias behaviour for shapes that have such different dependence in the squeezed limit, this lends support to their validity.
Revisiting the mass function predictions, we find that the nonGaussian effect on the halo mass function can be well modelled by the analytic predictions, if a shapedependent fudge factor, , is taken into account. In the case of the orthogonal type, there is a hint that this factor depends on redshift and halo mass.
For the halo bias we also find that the analytic predictions work well once a a shapedependent normalization factor is included: effectively rescales . Interestingly, for the nonlocal types of nonGaussianity, this fudge factor varies with redshift and halo mass.
Throughout we pay particular attention to the operational definition of nonGaussian bias, taking into account that a nonzero modifies the mass function and thus modifies also the linear bias. However, a Gaussian distribution with the same mass function would have a linear halo bias modified in the same way.
In most inflationary models considered so far, the bispectrum scales either as the equilateral shape () or as the local shape () in the squeezed limit (see, however, the quasisingle field models of [87, 88], which have intermediate scalings). The equilateral shape leads to a small and almost scaleindependent halo bias and is thus not easily accessible with this technique. Other shapes that, on large scales, have the same scale dependence as the local one, may in principle be distinguished through the halo bias dependence on halo mass on intermediate scales. On large scales, however, the limits for the local can be reinterpreted as limits on e.g., modifiedinitialstate type of nonGaussianity, by an appropriate rescaling of : .
Acknowledgements
We would like to thank Lotfi Boubekeur for useful discussions. The simulations and computations were performed on the cluster Hipatia (UB computing facilities and ERC grant FP7 IDEAS Phys.LSS 240117). We are grateful to the Centro de Ciencias de Benasque Pedro Pascual where the very last stages of the work were carried out. CW is supported by MICINN grant AYA200803531. LV acknowledges support from FP7PEOPLE200743IRG n. 202182 and FP7 IDEAS Phys.LSS 240117 and MICINN grant AYA200803531.
Appendix A Simulations of the equilateral type using two different approaches
In this section we compare the results of two Nbody simulations with nonGaussian initial conditions of the equilateral type. Both simulations have a box size of and a particle load of