N-body simulations with generic non-Gaussian initial conditions II: Halo bias

We present N-body simulations for generic non-Gaussian initial conditions with the aim of exploring and modelling the scale-dependent halo bias. This effect is evident at very large scales requiring large simulation boxes. In addition, the previously available prescription to implement generic non-Gaussian initial conditions has been improved to keep under control higher-order terms which were spoiling the power spectrum on large scales. We pay particular attention to the differences between physical, inflation-motivated primordial bispectra and their factorizable templates, and to the operational definition of the non-Gaussian halo bias (which has both a scale-dependent and an approximately scale-independent contributions). We find that analytic predictions for both the non-Gaussian halo mass function and halo bias work well once a calibration factor (which was introduced before) is calibrated on simulations. The halo bias remains therefore an extremely promising tool to probe primordial non-Gaussianity 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 non-Gaussian comparison project web site http://icc.ub.edu/~liciaverde/NGSCP.html


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 non-Gaussianity depend on the inflationary model. The simplest inflationary model, the standard single-field slow-roll inflation, produces non-Gaussianities which are unmeasurably small [2,3] (for a thorough review on inflationary non-Gaussianity see [4,5]).
Other inflationary models [6][7][8][9][10][11][12], which violate one or more assumptions of the standard single-field slow-roll inflation, may potentially cause much larger non-Gaussianities. 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, B Φ (k 1 , k 2 , k 3 ). The bispectrum is the Fourier transform of the 3-point-correlation function and is identical to zero for a Gaussian field. The functional form of the bispectrum specifies the form of the non-Gaussianity, while its amplitude, parametrized by the parameter f NL , gives the level of non-Gaussianity.
Since the Universe is statistically homogeneous and isotropic, the wave vectors k 1 , k 2 , and k 3 have to form a triangle. The shape of the triangle, for which the bispectrum peaks, is conventionally used to characterize the form of the non-Gaussianity. The local type of non-Gaussianity, Φ(x) = Φ G (x) + f NL Φ G (x) 2 − Φ G (x) 2 , where Φ G (x) is gravitational potential in the Gaussian case, yields a bispectrum that peaks for squeezed triangles (e.g., k 1 k 2 k 3 ). Potentially large non-Gaussianities of the local type are generally produced by inflationary models with multiple fields, where non-Gaussianity is created by non-linearities which develop outside the horizon [8].
Models in which the initial state of the inflaton is different from the Bunch-Davies vacuum [9,11,16,17], can produce non-Gaussianity of the enfolded shape, i.e. the bispectrum is dominated by flattened or enfolded triangle configurations (k 1 k 2 + k 3 ).
Finally, there is the so-called orthogonal shape of non-Gaussianity, 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 non-Gaussianity, since the temperature perturbations are still very small and well described by linear theory. If a shape of the primordial bispectrum is assumed 1 , constraints on the corresponding f NL value can be derived from the observed fluctuations in the CMB. The current constraints for the most commonly used shapes are −10 < f local NL < 74, −214 < f equil NL < 266, −410 < f orth NL < 6 (95% CL) [1], and f enfl NL = 114 ± 72 (68% CL) [20] . The Planck mission [21] is expected to measure f NL with an accuracy of ∆f NL 5.
Primordial non-Gaussianity leaves also observable signatures in the large-scale 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 non-linear gravitational evolution and biasing, this makes the measurement of the primordial contribution difficult [24][25][26] (see, however, [27]).
Another probe of non-Gaussianity from large-scale 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 non-Gaussainity [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 non-Gaussianity. This probe measures the primordial non-Gaussianity on scales of the order of 10 Mpc, which is smaller than the scales accessible by CMB measurements. Intriguingly, several groups [30][31][32] derived estimates of f NL from the observations of very massive galaxy clusters, which are as large as f NL ∼ 400 (see, however, [33]).
The clustering of halos on large scales offers another possibility to probe non-Gaussianity [34,35]. Non-Gaussianity 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 scale-dependent halo bias on large scales, which can be very different from the scale-independent halo bias in the Gaussian case. The amplitude of this effect scales approximately linearly in f NL , while the scale dependence is set by the shape of non-Gaussianity. For example, the local type non-Gaussianity generates a halo bias which scales as k −2 , whereas the halo bias of equilateral types of non-Gaussianity 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, −27 < f local NL < 70 (95% CL) [37], are competitive with the CMB measurements. Fisher matrix forecasts for future large-scale structure surveys predict that measurements of the non-Gaussian bias will constrain f local NL with an uncertainty of ∆f NL ∼ 1 [40][41][42][43].
All of the three aforementioned large-scale structure probes of non-Gaussianity are highly affected by non-linear gravitational evolution. The standard tool for taking nonlinearities into account are N-body simulations. However, until recently, only the local type of non-Gaussianity was simulated with N-body 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 N-body simulations with non-local non-Gaussianities. Using this method we were able to simulate the non-linear power spectrum and the cluster mass function for non-local shapes, and again found consistency with theoretical expectations. The technique of [48] was, however, not suitable to study the effect of non-Gaussianity 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 non-Gaussian halo bias.
In this paper, we modify our ansatz used for the generation of non-Gaussian 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 non-Gaussian simulations with large volumes, which we then use to study -for the first time by means of N-body simulationsthe non-Gaussian halo bias for different types of non-Gaussianity. 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 non-Gaussian cluster mass function.
The outline of the paper is the following. First, we review the non-Gaussian halo bias and its theoretical modelling in Sec. 2. In Sec. 3, we compare, in the context of the non-Gaussian 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 non-Gaussian 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 non-Gaussian mass function and halo bias in Sec. 6. Finally, we conclude in Sec. 7.

Non-Gaussian halo bias
The halo bias describes the relation between the halo density contrast, δ h , and the underlying dark matter density contrast, δ m . In Fourier space we can write the relation as δ h (k) = b(k)δ m (k) + n(k), where the random variable n(k) models the stochasticity coming from shot noise and possibly from physical halo formation processes. For Gaussian initial conditions the bias is scale-independent on large scales, b M,z (k) = b G 1;M,z for k 0.1 [49][50][51]. Analytic predictions and results of N-body simulations show that the linear 2 Gaussian bias, b G 1;M,z , increases with redshift z and halo mass M . In addition the halo bias may depend on environment and mass accretion history [52,53]. Primordial non-Gaussianity of the local type induces a scale dependence of the halo bias on large scales, which scales as ∼ k −2 . This particular scaling was first found in N-body 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 peak-background split, respectively. In Ref. [54], the same scale dependence of the non-Gaussian halo bias was found by computing the ellipsoidal collapse threshold in the presence of small non-Gaussianities. The analytic predictions were tested with further and larger N-body simulations [44][45][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 [27,[55][56][57][58][59]. In summary the non-Gaussian bias can be modelled to leading order as (2.1) where ∆b 1;M,z (f NL ) is the difference in the linear bias, which is non-zero, since a non-vanishing f NL changes the halo number density and hence the linear halo bias for a fixed halo mass M . Mainly, this effect leads to a scale-independent shift relative to the Gaussian (f NL = 0) linear bias [44,59,60]. Additionally, it makes the scale-dependent second term in above equation slightly non-linear in f NL [59,61]. In practice, when fitting observational data, the linear bias is estimated from the data itself.
The amplitude of the non-Gaussian contribution is given by f NL , the linear Lagrangian bias [b 1;M,z (f NL ) − 1], and the redshift-dependent collapse threshold qδ c /D(z), where δ z = 1.686 is the spherical collapse threshold in an Einstein-de Sitter universe, D(z) is the growth function normalized to (1 + z) −1 at high redshift, and q can either be regarded as a fudge factor introduced to improve the fitting to N-body simulations [46] or as a correction of the spherical collapse threshold [54,62]. The scale dependence of the non-Gaussian bias is hidden in M M (k) and F M (k). M M (k) connects the density with the primordial potential. The linear matter density contrast smoothed on a mass scale M is related to the primordial potential by the Poisson equation, Here, Ω m and H 0 are the present day matter fraction and the Hubble constant, respectively. The transfer function T (k) models the physics until the end of recombination and is normalised to unity on large scales. The function W M (k) is the Fourier transform of the top-hat filter with a smoothing scale given by the mass M . As W M (k → 0) = 1, we see that M M (k) scales as k 2 on large scales. The "form factor", F M (k), depends on the shape of the non-Gaussianity through the bispectrum of the primordial potential B Φ (k 1 , k 2 , k 3 ), [35]: where µ is defined by k · k = kk µ, P Φ (k) is the primordial power spectrum, and σ 2 M denotes the variance of linear density fluctuations on the mass scale M at redshiht z = 0. Considering the limit of k → 0, we see that on large scales the form factor depends primarily on the squeezed limit of the bispectrum. For the local type of non-Gaussianity, the "form factor" becomes scale and mass independent on large scales, F M (k 1) = 2. Hence, we obtain the aforementioned k −2 scale dependence of the halo bias in the case of local non-Gaussianity. As we discuss below in Sec. 3, other shapes of non-Gaussianity lead to different and in general weaker scale dependences.
One higher-order correction to ∆b, which was introduced in [44], is a relatively small contribution coming from the fractional difference in the Gaussian and non-Gaussian 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 f NL .
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 non-Gaussian halo bias was pointed out in [37], and further investigated and tested with N-body simulations in [63]. The formation redshift of the halos, z f , changes the amplitude of the non-Gaussian correction where µ(z f ) varies roughly between −2 and +2 from very young to very old halos [63]. For recent major mergers, µ takes the value −1 − 1/δ c ≈ −1.6 [37,63].
So far we discussed the scale-dependent halo bias caused by the local type of non-Gaussianity with a scale-independent f NL . The effect of a scale-dependent f NL were investigated analytically and by means of N-body simulations in [70,71]. In this paper we are mainly interested in the scale-dependent halo bias of non-local types of non-Gaussianity. The prescription presented above is sufficiently general to model the effects of other shapes of non-Gaussianity. As discussed above, the shape dependence is hidden in the form factor F M (k). 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 non-Gaussianity in [36]. Recently, [72] generalized the peak-background split approach to non-local types of non-Gaussianities and derived a very similar but not identical expression for the form factor. In the next section we discuss the analytic predictions of the non-Gaussian halo bias for non-local shapes using the high peak formulation of [35].

Physical shapes vs. templates
The functional form of the primordial bispectrum, B Φ (k 1 , k 2 , k 3 ), 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 k 1 , k 2 and k 3 . Separability is important for computationally efficient analysis and simulation of both CMB maps [13,19] and large-scale structure [48,73]. Hence, most physical shapes have been approximated by so-called 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 so-called "cosine" shape correlator is used [8], which absolute value is always ≤ 1 with 1 indicating perfect correlation. This measure of similarity is useful for constraints on the bispectrum coming from CMB or large-scale structure analysis, as it probes the entire range of possible triangle configurations. However, these templates can be misleading for predictions of the non-Gaussian 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 non-Gaussian 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 template, which was introduced in [74] as a separable approximation to the shape functions of DBI inflation [14], inflation with higher-derivative terms [13], and ghost inflation [15]: (3. 4) and the dependence of F on the three wavenumbers k is understood. The primordial power spectrum is specified by the amplitude A 0 at k 0 and the spectral index n s : 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 non-Gaussianity. 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, Bπ and Bπ 3 Φ , which correspond to the interaction operatorsπ(∂ i π) 2 andπ 3 (see [18] for details), are both close to the equilateral template. Especially, the cosine of Bπ with the equilateral template is very close to one. For this reason, we denote this shape also by "SSZ equilateral". 5 We now turn to the question how well the equilateral template captures the behaviour of these physical shapes in the squeezed limit. Choosing k 1 → 0 and using the following notation, k ≡ k 1 , k ≡ k 2 , and k 3 = |k + k | = k 2 + k 2 + 2k kµ, we obtain for the equilateral template where in the last line a scale-invariant power spectrum was assumed, i.e. n s = 1. 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 scale-invariant power spectrum, the scaling reduces for all models to k −1 with a model-dependent amplitude. For example for the bispectrum of DBI inflation (see [8] for an explicit expression) we get while the scaling of the bispectrum corresponding toπ(∂ i π) 2 [18] is (3.10) Here, we used the equilateral normalization convention for all shapes, i.e. Fig. 1, the scale-dependent part of the non-Gaussian halo bias, for a scale-invariant 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 1.5 and 1.256, respectively, in agreement with the above expressions for the squeezed limit. The solid lines show F M (k)/M M (k) computed with a power spectrum with n s = 0.95. The upturn and the wiggles on smaller scales stem from the transfer function included in M M (k). Note, in contrast to the local type of non-Gaussianity, F M (k)/M M (k) for the equilateral shapes is mass-dependent even on large scale.
In summary, predicting the non-Gaussian 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 50% and 26%, respectively.
For other templates, as we discuss now, the situation becomes more problematic. The 11) and the enfolded 7 template [16] scale both as ∼ k −2 in the squeezed limit, while the corresponding physical shapes scale as ∼ k −1 and ∼ k −3 , respectively. In Fig. 2 we show again the scale-dependent part of the non-Gaussian bias, F M (k)/M M (k), computed from these commonly used templates and the corresponding physical shapes. For comparison, the scale dependence obtained from the local type of non-Gaussianity, B loc Φ = 2f loc NL F loc , and from the equilateral template are also shown. Here, we used in all cases a scale-invariant power spectrum. The solid and dashed lines correspond to halos of mass 3 × 10 14 M /h and 3 × 10 11 M /h, respectively. The scale dependence of the halo bias caused by modified-initial-state type of non-Gaussianity [16] is almost degenerate with the one caused by the local type of non-Gaussianity, i.e. both predictions are very similar when f loc NL is rescaled by 1/8. However, on scales of ∼ 0.01h/Mpc 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 shapesdoes not lead to a scale-dependent 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 ∼ k −3 or as ∼ k −1 in the squeezed limit, while the scaling of the enfolded and orthogonal templates is ∼ k −2 .
In [19] a mode expansion of the bispectrum in separable functions was introduced. This expansion can be applied to arbitrary non-separable 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 large-scale 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 ∼ k −2 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 6 The orthogonal template was introduced as a separable approximation to the linear combination of the shapes Bπ and Bπ 3 Φ that is orthogonal (i.e. has a vanishing cosine) to Bπ This template approximates the shape function produced by inflationary models with a modifications of the initial-state of the inflaton field and has its maximum for flattened triangles, i.e. k1 = k2 + k3.  bispectrum does not speed up our current implementation of generating the initial conditions for N-body 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 non-Gaussian 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.

Initial conditions
In [48] we introduced a prescription how to generate non-Gaussian initial conditions for a given bispectrum, i.e. how we find a realization of the potential field Φ, which has the desired bispectrum We shortly review our method. First, we decompose the potential field in Fourier space, Φ k , in a Gaussian and a non-Gaussian part The Gaussian field Φ G k is statistically completely described by the power spectrum As the non-Gaussian contribution is small, the leading order of the bispectrum is given by Using the following ansatz for Φ N G 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 non-Gaussian contributions to the power spectrum can become very large on large scales. 8 Naively, one expects Φ N G Φ N G to be much smaller than the Gaussian contribution Φ G Φ G . However in the case of certain bispectra, Φ N G Φ N G diverges more strongly for small wavenumbers than the Gaussian power spectrum. With the above ansatz for Φ N G k the non-Gaussian contribution to the power spectrum reads .
In practice, there is a minimum wavenumber k min , the largest mode which fits in the simulation box, and a maximum wavenumber k max , which corresponds to the smallest wavelength still resolved by the used grid size. These wavenumbers provide large and small-scale cutoffs for the above integral. Nevertheless, this spurious contribution to the power spectrum can be large enough to affect substantially the results of the N-body 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, k 1, the integral is dominated by the bispectrum in the squeezed limit. Using the scalings presented in Sec. 3, we see that the non-Gaussian contribution scales as f 2 NL kP Φ (k), f 2 NL k −1 P Φ (k), and f 2 NL k −3 P Φ (k) for the equilateral, orthogonal/enfolded templates, and the local type of non-Gaussianity, 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 For the local type of non-Gaussianity, we found in [48] that our ansatz becomes identical to the real space formulation Here we generalize this modification by writing it in the following way .
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 k → 0, since they are suppressed by (4.10) Note that neither Eq. (4.5) nor Eq. (4.9) specify the trispectrum of the non-Gaussian 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 higher-order correlations.
Next, we can convert the potential Φ k into the linear density field δ k with the help of the transfer function, which we compute with CAMB [77], and the Poisson equation: 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 second-order Lagrangian perturbation theory [78][79][80].

Simulations
In this paper we are mainly interested in the scale-dependent halo bias of different types of non-Gaussianity. 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 trade-off between size and resolution is a box size of about 2 Gpc/h and a particle load of one billion particles. This enables us to probe large scales (k ∼ 0.003 h/Mpc) and resolve halos with masses above 10 13 M /h (∼ 20 particles).
The downside of our modified ansatz for the initial non-Gaussian part of the potential, Φ N G k given in Eq. (4.9), is the non-separability 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 k mode the following sum has to be computed where k runs over all grid points in Fourier space. If N g denotes the number of grid points in one dimension, then the computational cost of computing Φ N G k for all k modes scales as N 6 g . This scaling renders the usage of large grid sizes infeasible. Even for a moderately small grid size of 400, the generation of the non-Gaussian initial conditions takes two days on 256 cores of present-day 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, N p , which samples only the Gaussian contribution to the potential, and one smaller grid for the non-Gaussian part. 11 By doing this, we neglect non-Gaussianities on scales smaller than the grid spacing of the non-Gaussian 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 non-Gaussian grid will, presumably, hardly be affected. The choice of the non-Gaussian grid size is therefore a trade-off between computational resources and the lowest halo mass, which is still unaffected by the spatially unresolved non-Gaussianity.  To quantify this effect in more detail, we run two simulations of the local type of non-Gaussianity with f NL = 250. 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 non-Gaussian grid of size 400 (local-1024-400), in the other, both grids are of size 1024 (local-1024-1024). This large non-Gaussian grid size of 1024 is possible because in the case of local non-Gaussianity the computation can be performed efficiently in real space,  The simulations are carried out with Gadget-2 [81] taking only the gravitational interaction into account. The box size is 1875 Mpc/h, and 1024 3 particles with a softening length of 40 kpc/h are used to sample the density field. The simulations are started at an initial redshift of z initial = 49 using second-order Lagrangian perturbation theory to displace the particles from a regular grid.
Halos are found in the particle outputs of the N-body 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 redshift-dependent virial overdensity. With the halo catalogues and particle snapshots at hand, we can compute the Fourier modes δ h k and δ k of the halo and matter density field, respectively. This is achieved by first assigning the halos or particles to a regular grid (512 3 ) using the cloud-in-cell (CIC) scheme and then performing an FFT.
In Fig. 3, we show the halo bias obtained from the two simulations local-1024-1024 (solid lines) and local-1024-400 (dashed lines) for different halo masses at redshift z = 0 and z = 1. We compute the halo bias by the ratio of the halo-matter cross power spectrum Table 1. N-body simulations. N g denotes the size of the grid used for the non-Gaussian part of the potential. The size of the particle grid is identical to the size of the grid used for the Gaussian part of the potential and is given by N p . to the matter power spectrum with P m (k) = |δ k | 2 and P hm (k) = Re(δ h k δ * k ) , where the average is taken over all modes which lie in a shell of radius k and thickness ∆k = 0.003 h/Mpc. As expected, the clustering of halos with a mass below ∼ 3 × 10 13 M /h is affected by the smaller grid size used for the non-Gaussian Φ N G k . The grid spacing of the non-Gaussian field is 4.7 Mpc/h. This corresponds approximately to the Lagrangian radius of halos of mass 3 × 10 13 M /h. 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 non-Gaussian grid, is the halo number density. In Fig. 4, we show the ratio of the mass functions derived from the two simulations local-1024-400 and local-1024-1024. We find that the smaller non-Gaussian grid size decreases the halo number density at the low-mass end by a few percent. For halos more massive than ∼ 4 × 10 13 M /h, the halo number densities of both simulations agree within the errors.
In App. A, we compare the results of two N-body simulations of the equilateral type, for which we used the two different formulations for the non-Gaussian part of the potential, Eq. (4.5) and Eq. (4.9). This comparison gives us another estimate of the minimum halo mass, ∼ 5 × 10 13 M /h, above which the results are not affected by the low spatial resolution of the non-Gaussianity.
We conclude that, for halos more massive than 5×10 13 M /h, even the simulations which are set up with a small non-Gaussian grid (N g = 400) model correctly the characteristics of non-Gaussianity, in particular the halo bias and the halo number density. As a rule of thumb, the non-Gaussian grid spacing should be smaller than the Lagrangian radius of the lowest-mass halo of interest.
For the main study of the non-Gaussian halo bias, we perform a suite of simulations of different shapes of non-Gaussianity and different f NL 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 non-separable 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 non-Gaussian 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 non-Gaussian halo bias against N-body simulations. Lastly, other groups (e.g. [83]) are running N-body simulations of this type, too, which makes a future comparison easier.
The chosen values for f NL 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 N-body simulations. To do this accurately, we need high signal-to-noise and therefore a relatively large value of f NL . This is in particular the case for the equilateral type of non-Gaussianity, for which the non-Gaussian halo bias has only a very weak scale dependence. Note that if we find that the analytic modelling of the non-Gaussian bias works for relatively large non-Gaussianity, it will also be valid for small values of f NL .

Results
In this section we present the outcome of our N-body simulations of different types on non-Gaussianity. Here, we used the standard templates (namely the local, equilateral and orthogonal templates) to create the initial conditions for the suite of N-body 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 k-dependence in the squeezed limit, this lends support to their validity.
Although the simulations were targeted to study the non-Gaussian halo bias, we first consider the effect of non-Gaussianity on the halo mass function. Afterwards, we analyse in detail the scale-dependent halo bias induced by the different types of non-Gaussianity.

Halo mass function
Primordial non-Gaussianity 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 cluster-sized 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 non-Gaussian to the Gaussian halo mass function, R N G (M, z), can be modelled analytically in the framework of the Press-Schechter formalism [28,29] and the excursion set approach [84,85]. For the comparison with our N-body results, we use the formulas of [28] (MVJ) and [29] (LV), who used the saddle-point technique and the Edgeworth expansion, respectively, to derive the following expressions   [28] (MVJ) and [29] (LV) using the best-fit fudge factor of q = 0.75 and q = 0.9, respectively. and where σ 2 M is the linear mass variance on a mass scale M at z = 0. The redshift dependence is modelled by the collapse threshold, ∆ c (z) = √ qδ c D(z = 0)/D(z), which includes a fudge factor, √ q 12 . The non-Gaussianity is described by the normalized skewness of the linear density field, δ, on a mass scale M at z = 0, S 3,M = δ 3 M /σ 4 M . The normalized skewness scales linearly in f NL . The magnitude and sign of S 3,M depend on the shape of non-Gaussianity, while the mass dependence is only weakly dependent on the type of non-Gaussianity.
In  Figure 6. Effect on the halo mass function induced by the orthogonal shape of non-Gaussianity. Symbols and lines have the same meaning as in Fig. 5. Due to small signal to noise, the best-fit value for q at z = 0 is only determined to about 20%. and z = 0 (bottom right). The filled and open symbols represent the results for f NL = 250 and f NL = 60, respectively. The corresponding analytic predictions of MVJ and LV (with the best-fit value for q) are depicted by the solid and dashed lines, respectively. Overall, the agreement between simulations and predictions is very good. Only at the low-mass end at z = 0, 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 q = 0.75 for MVJ and q = 0.9 for LV is in agreement with our findings in [48]. Other groups [45,46], who worked with friends-of-friends (FOF) halos, used q ≈ 0.75 for both predictions, while [44], who applied a spherical-overdensity (SO) halo finder as we do, found q ≈ 1 for LV and q ≈ 0.75 for MVJ (V. Desjacques, private communication). For the standard definitions of FOF and SO halos, i.e. using a linking length of 0.2 times the mean particle distance (FOF) and the redshift-dependend virial overdensity (SO), [23] found that for a given halo mass the effect of non-Gaussianity 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 observationally-selected halos (e.g., SZ-selected or X-ray selected clusters). Nevertheless, the fact that q is quite close to unity in both cases, suggests that this will not introduce a dominant systematic effect.
The mass function results from our non-Gaussian N-body simulations of the orthogonal The triangles correspond to two simulations with the same realization of the initial Gaussian field but differing in the method how the initial non-Gaussian contribution was computed (Eq. 4.5 (black) and Eq. 4.9 (red)). The green squares depict the results of a second simulation with a different realization of the initial Gaussian field. Dashed and solid lines show the analytic predictions of [29] and [28].
type are presented in Fig. 6. Here, the f NL values are −1000 (filled symbols) and −250 (open symbols). Note that these simulations used a rather small grid size (N g = 400) for the non-Gaussian initial density field. Hence, halos with masses below 5 × 10 13 M /h are expected to show a smaller non-Gaussian effect than predicted by the theory. This is indeed the case, as can be seen in the figure. More interestingly, even at the high-mass end the non-Gaussian effect is smaller than one could expect. To bring the theory in agreement with the N-body data, the fudge factor has to be quite small, decreasing from q ≈ 0.75 to q ≈ 0.5 with decreasing redshift. Lastly, we consider the increase in the halo mass function due to non-Gaussianity of the equilateral type. In Fig. 7, the results of three simulations with f NL = 1000 are shown. Two of them, "eql. a)", started from initial conditions computed with our modified ansatz, Eq. (4.9), using a non-Gaussian grid size of 400. The third simulation, "eql. b)" (black triangles), used initial conditions generated with our original method, Eq. (4.5), with N g = 1024 (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 N-body results agree with each other for halo masses above 5 × 10 13 M /h. The fall-off of "eql. a)" at the low-mass end is caused by the small non-Gaussian grid. Interestingly, we find that, in the case of the equilateral shape of non-Gaussianity, q is very close to one and does not change significantly with redshift.
In conclusion, on the mass scales probed, the non-Gaussian effect on the halo mass function can reasonably well modelled by the analytic predictions, if a shape-dependent fudge factor is taken into account. In the case of the orthogonal type, there is a hint that q is in addition redshift dependent.

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 fall-off in the mass function reduces the number of halos per bin quickly. In the following, we only consider mass bins with more than 250 halos.
In order to keep the statistical error as small as possible in our analysis, we compute and fit the effect of the non-Gaussian halo bias in the following way. First, we compute the difference between the non-Gaussian and Gaussian bias from the corresponding N-body simulations: .

(6.3)
As the non-Gaussian and the Gaussian simulation share the same realization of the initial Gaussian field, ∆b(k) is almost free of sample variance and, in addition, has smaller shot noise than b N G (k) and b G (k) individually. We estimate the error on ∆b(k) directly from the distribution of ∆b(k) in each k bin (for details, see App. B).
Following Eq. 2.1, we model ∆b(k) by where b G 1 is the linear halo bias obtained from the Gaussian simulation on large scales (see App. B). The scale-independent shift, ∆b I , can be predicted by the difference in the linear bias derived from the non-Gaussian and Gaussian mass functions using the peak-background split approach [44] with R N G (M ) being the ratio of the non-Gaussian to the Gaussian mass function. Here, however, we treat ∆b I as a free parameter and compare it later on with the prediction derived from the mass functions. We choose q 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 non-Gaussianity on the halo bias for halos of mass 1.2 − 2.4 × 10 14 M /h at z = 0. Note that we plot ∆b(k) + 0.1. As ∆b I is negative, this addition of 0.1 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 k max = 0.1 Mpc/h) and include all terms given above. The short-dashed lines neglect ∆b I appearing inside the square brackets in Eq. (6.4). The inclusion of this term makes the non-Gaussian bias non-linear in f NL [59,61], since ∆b I depends on f NL . The dot-dashed lines neglect ∆b I completely. This scale-independent bias shift becomes important on smaller scales (k > 0.02), for which the scale-dependent part becomes small.
An example of the measured non-Gaussian bias from the simulations of the orthogonal type is given in Fig. 9. Here, the halos have again a mass of 1.2 − 2.4 × 10 14 M /h, but were found in simulation outputs at z = 1. 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 ∼ k −1 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 non-Gaussian halo bias induced by the equilateral  Figure 10. Same as in Fig. 8, but for the equilateral shape of non-Gaussianity. Note the linear scale of the y-axis.
shape. In Fig. 10, the halo bias corresponding to two different mass bins at z = 0 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.
-0.5 0 0  Figure 11. Best-fit values (fitting all modes with k < 0.1 Mpc/h) of the fudge factor q (left) and the scale-independent shift, ∆b I , normalized by f NL (right). Different colours correspond to different redshifts (see key in right panels). The symbol coding is the same we used in the mass function plots (Fig. 5, 6, and 7), see text for details. On the left-hand side, all points in between two vertical dashed lines correspond to the same halo mass bin, but show different redshifts (coulor), realizations (symbol shape), and f NL values (open/filled symbol). The horizontal dashed lines show the q values estimated from the LV non-Gaussian mass function. The coloured solid lines on the right-hand side show the predictions for the scale-independent bias shift using Eq. 6.5.
After having discussed for each type of non-Gaussianity typical examples, we show the complete set of best-fit values of the fitting parameters in Fig. 11. On the left-hand side, the best-fit values of the fudge factor q are presented. Different colours correspond to different redshifts: z = 1.5 (red), z = 1 (green), z = 0.67 (blue), and z = 0 (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 f NL = 60 and f NL = −250, respectively. Filled symbols show the results for f NL = 250 (local), f NL = −1000 (orth.), and f NL = 1000 (eql.). For clarity, the points of each mass bin are spread over the range of each mass bin.
For the local type of non-Gaussianity, we recover within the error bars the same q 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 non-Gaussian mass function, q ∼ 0.6. This suppression of the non-Gaussian halo bias effect is clearly redshift and halo mass dependent. In the case of the equilateral type of non-Gaussianity, 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 non-Gaussian halo bias effect for the equilateral shape was not found for the effect on the mass function, for which the measured fudge factor was q ∼ 1. 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 right-hand side of Fig. 11, the best-fit values of ∆b I normalized by f NL are shown. The colour and symbol coding is the same as before. As open and filled symbols (corresponding to different f NL values) are consistent with each other, we can infer that the scale-independent shift is linear in f NL for the f NL values probed. The solid lines represent the predictions from Eq. (6.5) using the LV mass function ratio, R N G LV (M ), and taking the measured fudge factor q 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 best-fit values and the predicted bias shift is very reasonable.

Conclusions
N-body simulations are an indispensable tool to test, develop and calibrate any statistical analysis of large-scale structure. In this paper, we have further explored the issue of setting up generic non-Gaussian initial conditions for N-body simulations. In a previous paper [48] we began addressing this issue, focussing -as we do here-on non-Gaussianities specified by a non-zero bispectrum. In [48] we focussed on the non-Gaussian corrections to the halo mass function and to the non-linear matter power spectrum, which, for the scales and redshfits of interest, could be reliably computed from simulation boxes of 600 Mpc/h on a side. Here we concentrate on the large-scale non-Gaussian halo bias. The scale-dependent bias is evident on very large scales thus needs multiple simulations of much larger boxes to be measured reliably for f NL values not too large as to push beyond the regime of validity of the analytical description of the effect.
The scale-dependent non-Gaussian halo bias has been recognized to be a very competitive approach to constrain primordial non-Gaussianity (e.g., [37,40,41,59]) yielding forecasted error bars on the local non-Gaussian parameter of ∆f loc NL ∼ 1 which makes this approach formally better than an ideal CMB experiment [86]. So far this effect has been explored with N-body simulations only for the so-called local type of non-Gaussianity. However, the shape of non-Gaussianity 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, higher-order contributions would leave an artificial signature on the large-scale 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 short-cuts that could be implemented for separable bispectra before, now cannot be applied in the same way.
Nevertheless we find a speed-up solution which involves sampling only the Gaussian contribution to the potential with a full resolution grid and using a smaller grid for the non-Gaussian 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 non-Gaussian 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 not-factorizable. 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 f NL 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 non-Guassian bias that scales as k −1 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 modified-initial-state type of non-Gaussianity is almost degenerate with the one caused by the local type of non-Gaussianity. The halo bias effect of the orthogonal physical shape is very similar to the effect of the equilateral shape.
We performed several simulations with non-Gaussian 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 k-dependence in the squeezed limit, this lends support to their validity.
Revisiting the mass function predictions, we find that the non-Gaussian effect on the halo mass function can be well modelled by the analytic predictions, if a shape-dependent fudge factor, q, 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 q is included: q effectively re-scales f NL . Interestingly, for the non-local types of non-Gaussianity, this fudge factor varies with redshift and halo mass.
Throughout we pay particular attention to the operational definition of non-Gaussian bias, taking into account that a non-zero f NL 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 (∼ k −1 ) or as the local shape (∼ k −3 ) in the squeezed limit (see, however, the quasi-single field models of [87,88], which have intermediate scalings). The equilateral shape leads to a small and almost scale-independent 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 f NL can be reinterpreted as limits on e.g., modified-initial-state type of non-Gaussianity, by an appropriate rescaling of f NL : f modin. NL 8f loc NL .
A Simulations of the equilateral type using two different approaches In this section we compare the results of two N-body simulations with f NL = 1000 non-Gaussian initial conditions of the equilateral type. Both simulations have a box size of 1875 Mpc/h and a particle load of 1024 3 . The simulations were started at z initial = 49 and carried out with Gadget-2 using the same numerical settings. In both cases, halos were identified with AHF [82]. The only difference between the simulations is the different prescription and implementation how the non-Gaussian part of the initial gravitational potential was generated. In the first case, "equilateral a)", we used our modified ansatz for Φ N G k , Eq. (4.9), and a grid size of 400 and 1024 for the non-Gaussian and Gaussian part, respectively. The initial conditions of the second simulation, "equilateral b)", were generated using our original ansatz, Eq. (4.5). The integral appearing in this ansatz is a sum of convolutions and can be computed efficiently with the help of FFTs. Hence, in this case, we were able to use an equal large grid size of 1024 for the non-Gaussian part, too. In both cases, we use the same realization of the Gaussian random field. Note although the two approaches are different at the level of Φ k , they both produce a non-Gaussian field, which has, at the leading order, the same power spectrum and bispectrum. However, high-order spectra, e.g. the trispectrum, are different.
In the remainder of this section, we consider observables which are affected by primordial non-Gaussianity, and analyse if the different ways of setting up the initial conditions lead to differences in these quantities. We start with the non-linear matter power spectrum. In Fig. 12, we show the ratio of the power spectrum measured from the two non-Gaussian simulations to the power spectrum obtained from a Gaussian simulation with the same realization. The cyan/blue and magenta/red lines correspond to redshift z = 0 and z = 1, respectively. The vertical black line shows the Nyquist frequency of the 400 3 grid used for the non-Gaussian part in "equilateral a)", which corresponds to the smallest scale still resolved by the grid. For both redshifts, significant differences between the non-Gaussian simulations only appear on scales smaller than the smallest scale sampled by the non-Gaussian grid used in "equilateral a)", i.e. on the right hand side of the vertical line. On the largest scales, the small number of modes leads to differences in the power spectra, since the term Φ G Φ N G ∝ Φ G Φ G Φ G does not completely vanish. The very good agreement on the intermediate scales confirms the theoretical expectation, that the non-Gaussian amplification of the power spectrum is at leading order caused by the primordial bispectrum.
Next, we show in Fig. 13 the halo bias for different halo masses at redshift z = 0. The halo bias is computed from the ratio of the halo-matter cross power spectrum to the matter power spectrum, see Eq. (5.2). The dashed and solid lines correspond to the simulations "equilateral a)" and "equilateral b)", respectively. As the effect of non-Gaussianity of the equilateral type on the halo bias is very weak (see Fig. 2), no effect of spatially unresolved non-Gaussianity, like we observed in the local case (see Fig. 3), is visible here. The small differences between the two simulations can be ascribed to the differences in the shot noise.
Lastly, we consider the halo mass function derived from the two non-Gaussian simulations. The ratio of the halo number densities measured from the simulations "equilateral a)" and "equilateral b)" as a function of halo mass is shown in Fig. 14. Similar to our test with the local type of non-Gaussianity (see Fig. 4), the number of halos with mass below 3×10 13 M /h is reduced by a few percent in the simulation for which a smaller non-Gaussian grid was used to set up the initial conditions. Halos with masses larger than ∼ 5 × 10 13 M /h are not affected by the smaller grid size. Both, the magnitude of the effect and the affected mass range are a bit larger than in the local case, which is probably due to the larger effect  on the halo mass function caused by the primordial non-Gaussianity with f eql NL = 1000 instead of f loc NL = 250 (see Fig. 5 and Fig. 7). As in the case of the non-linear matter power spectrum, no additional effects due to the differences in the initial non-Gaussian fields of "equilateral a)" and "equilateral b)" are noticeable. This strengthens the analytic modelling that the effects on the mass function due to non-Gaussianity are primarily caused by the primordial bispectrum.

B Details on the fitting of the non-Gaussian bias
We assume that the modes of the halo density field are related to the matter density by δ h = b(k)δ m + n(k), where n(k) is a noise variable which models the stochasticity in the halo formation process and the discrete nature of halos. The quantity which we want to measure from the simulations is the bias b(k). Assuming that the noise is uncorrelated with the density the bias is given by In practice, the average is taken over the finite number of modes per k bin. The term Re(nδ * m ) is therefore not exactly zero. Assuming Re(nδ * m ) has a Gaussian distribution with variance σ 2 nδ , the error of b(k) is then given by σ nδ / √ N modes / |δ m | 2 , where N modes is the number of independent modes.
In Fig. 15, we show the square root of the variance of the noise distribution for halos at z = 1 with masses of 1. respectively. The dashed lines show the Poisson noise prediction σ Poisson = 1/(2n h ), where n h is the halo number density. For the halos shown, the measured noise is in good agreement with the Poisson noise prediction. However, in general, depending on the halo mass and redshift, the actual noise can be smaller or larger than the Poisson noise and can also be scale-dependent (e.g., [89]). As the Gaussian and non-Gaussian simulations share the same realization of the Gaussian field, we expect that in the quantity the corresponding noise terms Re(n N G δ N G m * )/|δ N G m | and Re(n G δ G m * )/|δ G m | cancel each other partly. This is indeed what we find in the simulation data. The blue symbols in Fig. 15 represent σ n derived from the distribution of ∆ in each k-bin. The solid line is a linear fit to the data. When we fit the non-Gaussian halo bias models to the simulation data, we use this fit as our noise estimate. The resulting χ 2 /d.o.f. values are close to one.
Another ingredient in our fitting procedure of the non-Gaussian halo bias is the Gaussian linear bias, b G 1 . We estimate the linear bias by fitting the bias derived from the Gaussian simulations on large scales. In Fig. 16, the bias measurements for halos of different masses at redshift z = 1 are shown. The dashed lines depict the corresponding best fits of the linear bias. Only data points on large scales, k ≤ 0.07, were included in the fitting.
The remaining quantities needed for the fitting of the non-Gaussian bias (see Eq. (6.4)), like the growth function, the form factor, etc., are computed numerically for the given redshift,   type of non-Gaussianity, and halo mass 13 .
Having all quantities at hand, we fit the non-Gaussian bias derived from the simulation data with the model given in Eq. (6.4) out to wavenumbers k < k max . We use two different fitting ranges, a conservative k max = 0.03 Mpc/h and a more ambitious k max = 0.1 Mpc/h. Both fitting ranges yield results consistent with each other. However, in the case of the equilateral and orthogonal shape, the errors on the fitting parameters are significantly smaller (approximately by a factor of ten and two, respectively) when the larger fitting range is used. For the local type of non-Gaussianity, the larger fitting range reduces the error by approximately 30%.