N-body simulations with generic non-Gaussian initial conditions I: power spectrum and halo mass function

We address the issue of setting up generic non-Gaussian initial conditions for N-body simulations. We consider inflationary-motivated primordial non-Gaussianity where the perturbations in the Bardeen potential are given by a dominant Gaussian part plus a non-Gaussian part specified by its bispectrum. The approach we explore here is suitable for any bispectrum, i.e. it does not have to be of the so-called separable or factorizable form. The procedure of generating a non-Gaussian field with a given bispectrum (and a given power spectrum for the Gaussian component) is not univocal, and care must be taken so that higher-order corrections do not leave a too large signature on the power spectrum. This is so far a limiting factor of our approach. We then run N-body simulations for the most popular inflationary-motivated non-Gaussian shapes. The halo mass function and the non-linear power spectrum agree with theoretical analytical approximations proposed in the literature, even if they were so far developed and tested only for a particular shape (the local one). We plan to make the simulations outputs available to the community via the non-Gaussian simulations comparison project web site http://icc.ub.edu/~liciaverde/NGSCP.html.


Introduction
The leading theory for the origin of primordial perturbations is inflation. Along with the shape and amplitude of the primordial power spectrum and the signature of a stochastic background of gravity waves, primordial non-Gaussianity offers a window to probe inflation. Even the simplest inflationary models predict deviations from Gaussian initial conditions, arising from the (self)-interaction of the field during inflation; for the simplest, slow roll, single field model these deviations are expected to be unmeasurably small [1,2]. For a thorough review of inflationary non-Gaussianity see e.g., [3] and references therein.
An important development of the past few years is the realization that a large, potentially detectable amount of non-Gaussianity can be produced by inflation if any of the conditions giving the standard, single-field, slow-roll is violated. These are: a) only one field is responsible for generating the primordial perturbations b) canonical kinetic energy of the field c) slow roll d) the quantum field was in the adiabatic (Bunch-Davies) vacuum. In particular the violation of each of these conditions leaves its "signature" on the statistical properties of the initial perturbations (e.g., [4][5][6][7][8][9][10][11] and references therein).
Deviations from Gaussianity can be characterized by the dependence of the bispectrum signal 1 on the configuration (or shape) of the three k-vectors in its argument. For example local type of non-Gaussianity yields a bispectrum that is dominated by the squeezed configurations (where k 1 ≪ k 2 ≃ k 3 ) and is generated by models of inflation involving multiple fields. 2 This shape is in general associated with models of inflation where non-Gaussianity is created by non-linearities developed outside the horizon [6]. When non-Gaussianity is created JCAP10(2010)022 2 Creating a 3D non-Gaussian field with a given bispectrum The argument of [35] (see also [36]) which applies to a two-dimensional field expanded in spherical harmonics, can be generalized to apply to a three-dimensional field transformed to Fourier space. Pioneering work can be found in [39,40], here we follow a different route. Let's start from (2.1) Here B is intended to be the bispectrum of the Φ field. In this convention Φ is the Bardeen potential (i.e −1× the standard gravitational potential, on sub-horizon scales). This equation refers to some time deep in the matter-dominated era i.e., we use the f NL CMB convention. We can decompose the Φ field in a Gaussian and non-Gaussian parts. In Fourier space this reads and therefore with Φ * k the complex conjugate of Φ k , and use it in eq. (2.2), we recover the identity of eq. (2.1).
Note that this ansatz yields a non-Gaussian field with a specified bispectrum, but this operation does not specify the higher orders. It is often assumed that non-Gaussianity is specified by a bispectrum and that higher-order correlations (and the power spectrum) are not affected by it. Only very recently some authors have began to consider the effect of non-Gaussianity on the trispectrum for some specific cases [41][42][43]. In this paper, however, we neglect contribtutions from higher-order correlations.
It is important to bear in mind (as already pointed out in [36]) that this procedure is not unique: in other words there will be more than one -non-equivalent-expression for Φ NG k yielding the same bispectrum. It is instructive to use the local non-Gaussian case as a worked example. The local case is usually defined in real space as [15,44,45]: where Φ G denotes a Gaussian random field. In Fourier space this becomes: where the constant has been absorbed in the k = 0 mode which is ignored. The bispectrum for the local non-Gaussianity is: rather than using eq. (2.7). This extra "degree of freedom" was used in ref. [36] to produce non-Gaussian fields with a given bispectrum and with a minimal non-Gaussian contribution to the power spectrum. We return to this point below.

Special cases and initial conditions implementation
Let's start from eq. (2.1). The bispectrum depends on the shape of the triangle made by the three k vectors in its argument. The detailed dependence of the bispectrum on the k vectors (also called in brief bispectrum shape) can help discriminate among different inflationary models. For example, local shape non-Gaussianities were the first type to be considered [15,44,45] and are a direct consequence of the non-linear relation between the inflaton fluctuations and the curvature perturbations that couple to matter and radiation. While single-field slow-roll inflation generates this shape of non-Gaussianity, it has been shown that the amplitude is proportional to the square of the slow-roll parameter, which is very small [1,2]. In contrast, large local non-Gaussianities can be generated in e.g., curvaton models [5], where the curvature perturbation can evolve outside the horizon or inflationary models with multiple scalar fields. Non-canonical kinetic terms and higher derivative contributions to the inflaton potential can produce significant levels of non-Gaussianity of the equilateral type if the speed of sound in these models is much smaller than the speed of light, which can be realized in certain brane inflation scenarios [46][47][48][49]. These two shapes arise due to the interactions of the field driving inflation. On the other hand, modified initial states (i.e. not starting from a Bunch-Davies adiabatic vacuum) also introduce deviations from Gaussianity. These are characterized by a different shape [7,14] which we call enfolded. Finally [13] introduced the so-called orthogonal shape. The motivation behind this shape is outlined in [13]. In summary, the space of non-Gaussianities produced by the most general single-field models, where the inflaton fluctuations have an approximate shift symmetry, is spanned by linear combinations of two independent shapes: the equilateral and the orthogonal. We shall show below that these four bispectrum types can be obtained from a linear combination of just three "kernels", which we will call F local , F A and F B , speeding up our calculations. The expressions for the bispectra for the four most popular shapes mentioned above are reported below.
In the local case: and we have used B to denote the bispectrum for f NL = 1, and the dependence on the three k vectors is understood. In the equilateral case:

JCAP10(2010)022
where F eq = −P (k 1 )P (k 2 ) + 2cyc. − 2[P (k 1 )P (k 2 )P (k 3 )] 2/3 (3.4) In the last equality we have introduced and, for simplifying the notation, we have avoided writing down the explicit dependence of F on the three k vectors. For the template for factorized enfolded [14]: where Note that this is not the true enfolded configuration, but it is a factorizable template that approximates it well.
Finally, for the orthogonal shape: It is important to note that these four shapes are not independent. For example the orthogonal one can be obtained from a combination of the other two: The factorization introduced here speeds up the implementation of the initial conditions for N-body simulations of these four shapes: initial conditions can be produced directly in Fourier space and split in a Gaussian realization Φ G k and three different non-Gaussian realizations Φ NG k , one for the local case, one corresponding to F A and one to F B . From this one can build initial conditions for different shapes and different f NL .
In the following we give some details on how we generate the initial particle distribution. We start off from the publicly available code by Sirko [52] for Gaussian initial conditions and modify it appropriately. The difference between setting up Gaussian and non-Gaussian initial conditions is the extra non-Gaussian contribution, Φ NG k , to the gravitational potential. For a generic, non-factorizable bispectrum, we can compute this field by disctretizing eq. (2.4), i.e., for each grid point k in Fourier space, we loop over all grid points k ′ :

JCAP10(2010)022
where Φ G k is a random realization of a Gaussian field with the power spectrum given by P (k) ∝ k ns−4 and n s denotes the spectral index. If the number of grid points is given by N 3 g , the computational costs of the generation of such a generic non-Gaussian field Φ NG k scales as N 6 g . For a modest grid size of 256 3 , this results already in 3 × 10 14 evaluations of the summand in eq. (3.12), which take about 1000 hours on a single core of a present-day CPU. For example, the computation for a 512 3 grid would take approximately 10 days on 256 cores.
For factorizable shapes the process can be greatly sped up. In fact note that in this case eq. (2.4) can be written as a convolution (or a sum of convolutions) of two auxiliary fields which evaluation can be swiftly done resorting to Fourier transforms. In fact if the bispectrum can be written in a factorizable form as The integral can then be quickly performed as a multiplication in real space (note the similarity with eqs. (2.5) and (2.6) where instead of the Φ field we have two auxiliary fields G and Q).
While for factorizable shapes the direct summation approach (eq. (3.12)) and the convolution approach are mathematically equivalent in a practical implementation they may be affected by different numerical effects. Here we explore both implementations for the local type. In the other cases we implement the direct summation approach to explore whether this approach is viable and uncover possible bottlenecks or limitations.
Next, the linear density field δ k at z = 0 is derived from the potential Φ k through the transfer function and the Poisson equation: where D(z) is the linear growth function normalized to (1 + z) in the matter-dominated era, and T (k) is the transfer function obtained with CAMB [51] and normalized to unity on large scales. Ω m is the present-day matter fraction and H 0 the Hubble constant. The particles are then displaced from a regular grid according to the displacement field at the initial redshift, z i = 49, using the Zel'dovich Approximation. 4

Testing the initial conditions
In this section we analyze the quality of the non-Gaussian initial conditions. First we consider the 1-point function, i.e. the probability distribution function of the density contrast δ at the position x. Especially, we compute the variance, δ 2 , and the skewness, δ 3 , as a function of the smoothing scale R and compare the results with the analytic predictions. Furthermore, we calculate the power spectrum of the non-Gaussian initial conditions and demonstrate that the deviations from the Gaussian case are in most cases small. The bispectrum is a more complex quantity than the power spectrum as it depends on two variables and all modes are highly correlated. For this reason, it will not be shown here but is left for a forthcoming  Figure 1. Variance of the linear density field at z = 0 as a function of the smoothing radius R. The non-linearity parameter is f NL = 100. The symbols show the variance derived from the 8 different realizations of each type of non-Gaussianity. The lines depict the theoretical predictions. For clarity, the case "local a)" and "local b)" are slightly shifted horizontally. All other symbols and lines fall on top of each other as expected, as they are derived from the same input power spectrum and therefore share the same Gaussian part.
paper. In addition, the skewness is a higher signal-to-noise quantity than the bispectrum: we decide therefore to use the skewness as a "metric" to analyze the quality of the non-Gaussian initial conditions. First we generate eight Gaussian realizations of the density field on a grid of size 256 3 in a box with a side of 600 Mpc/h. From these Gaussian realizations we produce for each of the previous mentioned types of non-Gaussianity, i.e. local, equilateral, enfolded, and orthogonal, a non-Gaussian density field with an f NL of −500, −250, −100, 100, 250, and 500.
For the local type, we compute the non-Gaussian contribution in three different ways: a) the traditional way by squaring the Gaussian density field in real space b) in Fourier space using our ansatz eq. (3.12) with the bispectrum given by eq. (2.7) c) using again eq. (3.12) but with B(k 1 , k 2 , k 3 ) −→ 6f NL P (k 2 )P (k 3 ), this recovers the traditional method in real space except for aliasing effects introduced by the finite grid size.
In figure 1 and figure 2 we show the variance and skewness of each realization (data points) for all types of non-Gaussianity with f NL = 100 and compare them to analytic predictions (solid lines). The analytic prediction for the skewness is obtained by integrating the bispectrum (see e.g., eq. (4.13) of ref. [18]). We have truncated the integrals at the maximum and minimum k sampled by the simulation box. The magnitude of the effect depends on box size, scale and type of non-Gaussianity. For our simulation settings, we find it always to be below 15%. Note that the skewness of the Gaussian field due to the finite volume and grid size is subtracted from the measured skewness of the non-Gaussian field. Lines show the theoretical predictions. Only the orthogonal case has a negative skewness for positive f NL . For clarity, the case "local a)" and "local b)" are slightly shifted horizontally.
The moments of the density are computed from the linear density field smoothed with a top-hat filter of radius R. In case of the skewness, the δ 3 of the Gaussian realization is subtracted from the total skewness in order to reduce the noise introduced by the finite volume and grid size. The variances agree well with the theory except for the type "local b)". The increased variance for this type is caused by the term Φ NG k Φ NG −k , which, in this case, is not negligible. This term gives rise to a P 2 (k) term multiplied with a divergent integral, which in our case of a finite box and grid is truncated. This is reminiscent of the discussion in section V of [36]; for more details see appendix A. Apart from the "local b)" case, the values of the skewness obtained from the initial conditions also agree well with the predictions. Only for the orthogonal case small deviations at larger radii are visible.
Next we consider the power spectrum of the initial particle distributions. We compute the power spectrum by first assigning the particle to a 512 3 grid using the cloud-in-cell scheme. Then we perform a Fast Fourier Transform and average the k-modes over spherical shells with a thickness of ∆k = 0.01 h/Mpc. The ratios of the power spectrum of the non-Gaussian and the Gaussian initial conditions are shown in figure 3. The eight different realizations are depicted by the different symbols. For clarity, we show the individual realization only for k ≤ 0.1 h/Mpc. The solid line represents the mean of the realizations. On large scales, the power spectrum of the "local b)" case deviates strongly from the Gaussian power spectrum. At the lowest wave number k = 0.01 h/Mpc the power spectrum is enhanced by almost an order of magnitude! This behavior explains the offset of the variance and can be traced back to the second-order term Φ NG k Φ NG −k , which is further explored in the appendix A. We ob-  Figure 3. Ratio of the power spectra derived from the particle distributions of the non-Gaussian (f NL = 100) and Gaussian initial conditions. For k < 0.1 h/Mpc, individual realizations are shown by symbols. The lines depict the arithmetic mean of the realizations. For clarity, the case "local a)" and "local b)" are slightly shifted horizontally. serve a similar effect for the orthogonal shape, although the deviations are much smaller and a possible systematic shift, i.e. enhancement or suppression of the power spectrum, is -if at all -at the few percent level. Nevertheless we discard the orthogonal shape from now on and refer to the appendix A for more details on this issue. The case "local a)" and "local c)" agree perfectly with each other as it is expected from their mathematical equivalence. Both of these implementations of the local type and the equilateral case (almost not visible) do not show deviations from the Gaussian power spectrum. The power spectrum of the enfolded shape has very small deviations from the Gaussian power spectrum, the mean only deviates at the sub-percent level. We also checked variance, skewness, and power spectrum for the other f NL values and found qualitatively similar results. The deviations from the Gaussian power spectrum scale roughly linear with f NL , except for the "local b" case for which they scale as f 2 NL (see appendix A for details).

Testing the simulation settings
In this section we perform several tests to assure that the settings -like force and mass resolution, box size, and initial redshift-of our simulations are good enough to assure that systematic effects are smaller than the statistical errors. Note that we only consider ratios of non-Gaussian to Gaussian quantities in this paper. In particular, we are interested in ratios of the non-linear power spectrum and the halos mass function. Most of the systematic effects due to numerical limitations are expected to cancel out when considering ratios (see e.g. [53] addition, since the non-Gaussian and Gaussian simulations are based on the same realization of the random density field, also the statistical error on the ratios is reduced vastly. The standard settings of our main set of simulations can be found in the first half of table 1. The second half of the table describes the simulations which we use to investigate the numerical uncertainties of our main set of simulations. The generation of initial conditions done with the method "local a)" -the traditional real-space implementation of local non-Gaussian initial conditions-is by far not as computationally intensive as the method in Fourier space described in section 2 and 3. We exploit this fact and use the case "local a)" to perform larger N-body simulations to explore the numerical limitations of the smaller simulations of the main set. In addition, we use the case "local a)" as a benchmark for the implementation of local non-Gaussianity in Fourier space.
All our N-body simulations are performed with the publicly available code Gadget-2 [54]. Our choice of cosmology is a flat ΛCDM cosmology, which is consistent with the seven-year WMAP results [55]. In particular, we choose the following values: Ω m = 0.27, Ω b h 2 = 0.023, h = 0.7, n s = 0.95, and σ 8 = 0.8, where Ω m is the matter density, Ω b the baryon density, h the Hubble parameter, n s the spectral index of the primordial power spectrum, and σ 8 the rms of linear density fluctuations at z = 0 in a sphere of 8 Mpc/h. First, we consider the ratio of the non-linear power spectra derived from the local non-Gaussian simulations with f NL = 100 and Gaussian simulations at z = 1 and z = 0. Our force and mass resolution is such that we expect numerical uncertainties in the power spectrum to be at the sub-percent level for k 0.3 Mpc/h (see [56]). As we look at the ratio of the power spectrum the actual uncertainty will be even smaller. The results are presented in figure 4. The black dashed and solid lines depict the average of the eight realizations of the case "local c)" at redshift z = 1 and z = 0, respectively. The magenta and cyan lines show the two different implementations of the local type of non-Gaussianity for a single realization. The agreement is excellent and demonstrates the functionality of our method in Fourier space. The ratios derived from simulations of a larger box (1200 Mpc/h, blue lines) and from simulations with a higher starting redshift (z i = 99, red lines) are consistent with the ones derived from our main set of simulations. Only in the case of higher force and mass resolution, we find statistically significant deviations on smaller scales, i.e. the ratio falls off for k 0.   table 1). The dashed and solid black lines show the mean ratio derived from the main set of our simulations of type "local c)". and k 0.3 Mpc/h at z = 1 and z = 0, respectively. The small reduction of the ratio is probably caused by the enhancement of non-linear power due to the better resolution.
The second quantity we are interested in is the ratio of the non-Gaussian and Gaussian halo mass function. To identify halos in the simulation data, we use the publicly available halo finder AHF [57,58]. AHF defines halos to be gravitationally bound objects which have a spherical overdensity (SO) given by the redshift dependent virial overdensity. The minimal number of particles in a halo is 20, thus in the main set of our simulations we find halos of mass 2 × 10 13 M ⊙ /h or higher. Although 20 particles are too few to reliably resolve all halos of the corresponding mass (about fifty to several hundred particles are needed to resolve SO halos reliably [58,59]), we find that the ratio of the mass functions is not affected by the low mass resolution. This is demonstrated in figure 5. The black dashed and solid lines show the averaged fractional difference in the non-Gaussian and Gaussian mass functions at redshift z = 1.5 and z = 0.67, respectively, derived from the eight simulations of the type "local c)". The results obtained from the simulations with an eight times higher mass resolution are depicted by the green symbols and are consistent with the lower resolution results. In addition, figure 5 shows that the ratio is not biased by the finite-volume (blue symbols) nor the initial redshift (red symbols). The very good agreement between the two different ways of setting up the local case (cyan and magenta symbols) reassures us of the correct implementation of the method described in section 2 and 3.
Overall the results of this section give us confidence that for our main set of simulations -including the other types of non-Gaussianity -systematic effects due to numerical limitations are within the statistical errors.

Results
Here, we present our findings derived from the non-Gaussian simulations of the local, equilateral, and enfolded type for f NL = −500, −250, −100, 100, 250, and 500. First we present the results for the non-linear matter power spectrum. Afterwards, we turn to the halos mass function. The investigation of the non-Gaussian halo bias effect e.g., [28,29] and the bispectrum is left to forthcoming work [50].

Non-linear power spectrum
High precision in the theoretical prediction of the non-linear matter power spectrum at k ∼ 1 h/Mpc and above is needed to derive unbiased results from the data of upcoming weaklensing surveys e.g., [60]. Future Lyman-alpha forest surveys will also access these scales and test them with small statistical error-bars. Non-Gaussianities in the initial conditions alter the non-linear power spectrum at the few percent level. In figure    Note that the maximum of the ratio is larger at higher redshifts and its positions is shifted to smaller scales. The shape, especially the peak location and peak height, is consistent with predictions of the halo model presented in [62].

Halo mass function
At the high-mass end, the halo mass function is very sensitive to a primordial skewness in the probability distribution function of the density field [63]. Hence, galaxy cluster surveys offer in principle the possibility to probe primordial non-Gaussianities. In figure 9 we present the ratios of the cumulative mass functions derived from the non-Gaussian and Gaussian simulations. We show the results for the local, equilateral, and enfolded shape with different f NL at redshift z = 1.5 and z = 0.67. In addition to the data points we plot the analytic predictions of [20] (MVJ, blue lines) and [18] (LV, red lines). In order to convert the analytic ratio of the non-Gaussian and Gaussian mass functions, r NG (M ), provided in [20] and [18] into the corresponding ratio of the cumulative mass functions, R NG (> M ), we use the fitting formula of [64] for the SO halo mass function, n Tinker (M ), i.e. (6.1) We checked that the integrated Tinker fit is a good fit to the cumulative mass function of the SO halos found in our Gaussian simulations.  Figure 7. Same as in figure 6 but for the equilateral type of non-Gaussianity.  Figure 8. Same as in figure 6 but for the enfolded type of non-Gaussianity.
Note that, compared to [21], we do not find that substituting the linear spherical collapse overdensity by δ c −→ √ qδ c (with q = 0.75) improves the agreement between the N-body data and the predictions. This is in accordance with the findings of [25,26], who argue that the differences are due to the different halo identification algorithms. While in this paper we use halos defined by spherical overdensities, [21] applied a Friends-of-Friends halo finder to their simulations.
There is an indication that using SO halos a value q ∼ 0.9 yield a slightly better fit to our simulations data. Given our simulations size, this effect is however not highly significant and only visible for high values of f NL . While the exact value of q is an important theoretical JCAP10(2010)022 issue, in practice it may not be crucial. Recall that a 10% uncertainty in q will propagate into a 10% error on f NL . For values of f NL 10, as expected from most models, this uncertainty is at most comparable to the expected statistical errors of future measurements.
At the high-mass end, we find that LV fits the data better for positive f NL , whereas MVJ gives better predictions for negative f NL . In general, for large positive f NL the theoretical predictions overestimate the increase of halo number density at the low-mass end. However, for all shapes, the qualitative behavior is captured well by the analytic models.

Conclusions and discussion
We have addressed the issue of setting up generic non-Gaussian initial conditions for N-body simulations. In the era of precision cosmology, N-body simulations have become a crucial tool to test, develop and calibrate any statistical analysis of large-scale structure surveys. While the approach of constraining primordial non-Gaussianity with large-scale structure and galaxy surveys has recently received renewed attention, as far as we know, N-body simulations have been run only for the so-called local-type of non-Gaussianity. The shape of non-Gaussianity, however, is crucial if one wants to use such a signal as a window into the generation mechanism for primordial perturbations.
Building on the expertise developed in the context of Cosmic Microwave Background non-Gaussianity, we have shown how to set up non-Gaussian initial conditions for any type of non-Gaussianity specified by a primordial bispectrum. Given the current cosmological constraints all non-Gaussian fields we consider are given by the sum of a (dominant) Gaussian component and a (subdominant) non-Gaussian one.
The implementation is based on direct summation in Fourier space and is significantly more computationally intensive than the workhorse local case defined in real space, as it scales as N 6 g where N 3 g is the number of grid points in the simulation box. We have investigated the numerical effects of such an approach by comparing the results of the local-type case obtained with the two approaches. For factorizable templates, the implementation can be sped up significantly: in fact the operation can be rewritten as a convolution of two suitably defined auxiliary fields. This is useful, however not the main goal of our investigation, as our aim was to develop an approach suitable for non-factorizable templates. This is relevant because for some physically motivated inflationary models the existing factorizable templates may not be a good approximation of the effect e.g., [67][68][69] and references therein.
The procedure of generating a non-Gaussian field with a given bispectrum (and a given power spectrum for the Gaussian component) is not univocal and care must be taken so that higher-order corrections do not leave a too large signature on the power spectrum especially on large scales. This is so far a limiting factor of our approach, as these components can be kept under control only in specific cases. More investigation is clearly needed, possibly along the lines proposed by [36].
We have then explored the effects of several popular forms of primordial non-Gaussianity on the halo-mass function and the non-linear power spectrum deferring the analysis of other statistics such as the non-Gaussian halo bias to a forthcoming paper. We confirm that the non-Gaussian correction to the halo mass function is determined by the primordial skewness and that a suitable combination of the different analytical approximations proposed in the literature (depending on the regime of applicability) offer a good fit to the simulations data. We also confirm that independently of the type of non-Gaussianity spherical-overdensity JCAP10(2010)022 halos need a q correction factor much closer to unity (if any correction at all) than Friendsof-Friends halos.
We believe that the methodology illustrated and developed here will be relevant for the on-going and planned efforts of constraining primordial non-Gaussianity from large-scale structure including surveys of galaxies, high-redshift clusters, Lyman-alpha forest and weak lensing. A Non-Gaussian contributions to the initial power spectrum As before we write the gravitational potential as the sum of a Gaussian and non-Gaussian field, Φ k = Φ G k + Φ NG k , where the non-Gaussian part is quadratic in Φ G . The power spectrum, defined by Φ k Φ q = (2π) 3 δ D (k+q)P (k), consists then of the following terms Φ The term Φ G Φ NG will vanish in theory as it involves an odd number of Gaussian fields. In practice, we have only a finite number of modes to perform the average. Hence, especially on larges scales where there are only a few modes in the box, this term is not exactly zero. We investigate this term further below.
For certain combinations of the F -functions many of the above terms cancel. For example, with B ∝ (F local − F B ) the P 2 (k)-terms vanish in the above limit. Moreover for the JCAP10(2010)022 equilateral case (see eq. (3.4)) all terms with r ≥ 1 cancel, whereas for the enfolded type the term with P 4/3 (k) does not vanish. Now let us return to the term Φ G Φ NG linear in f NL . In the discretized version this becomes: For small k, there are only a few modes in the box and the skewness of the realized Gaussian field is not completely vanishing. Depending on the choice for B this noise can be amplified significantly. In particular, for our ansatz of the orthogonal type (eq. (3.10)) this contributes substantially to the power spectrum and causes the deviations visible in figure 3.